Profile water storage as calculated from SSURGO
Overview
This example illustrates a component-percent weighted mean water storage, based on total profile storage values. The comppct_r column is used as a weight to adjust the influence of the profile water storage values for each component: larger weights (larger components) result in more influence. When performing weighted means, use care in the selection of an appropriate weighting parameter. A selection of a weighting variable not direclty related to the property which is being averaged can result in an effect known as Simpson's paradox.
Some background on Soil Water Holding Capacity and Irrigation Management.
The general formula for computing profile water storage is defined as follows. Note that the awc_r column in the SSURGO database (chorizon table) is pre-adjusted to compensate for coarse fragments (see relevant SSURGO metadata below).
Calculate the water storage within a given component
component_whc = profile_sum_awc = sum( (hzdepb_r - hzdept_r) * awc_r )
Calculate the weighted mean value of the profile water storage, for a given map unit
weighted mean = sum(weights * x) / sum(weights)
mapunit_whc = sum( comppct_r * component_whc ) / sum( comppct_r )
Note: weights corresponding to non-soil components must be filtered out
From the SSURGO 2.1 table definitions:
The amount of water that an increment of soil depth, inclusive of fragments, can store that is available to plants. AWC is expressed as a volume fraction, and is commonly estimated as the difference between the water contents at 1/10 or 1/3 bar (field capacity) and 15 bars
(permanent wilting point) tension and adjusted for salinity, and fragments.
1. Reduce to 1:1 relationship with cokey by aggrigating AWC horizon data in the chorizon table, and linking to the component table.
SELECT mukey, compname, comppct_r, a.* FROM component
JOIN
(
SELECT cokey, sum( (hzdepb_r - hzdept_r) * awc_r) AS component_whc, sum((hzdepb_r - hzdept_r)) AS depth
FROM chorizon WHERE areasymbol = 'ca113'
GROUP BY cokey
) AS a
ON component.cokey = a.cokey
WHERE component.areasymbol = 'ca113'
ORDER BY mukey ;
Query Results
mukey | compname | comppct_r | cokey | component_whc | depth
--------+---------------------+-----------+----------------+---------------+-------
459225 | Millsholm | 30 | 459225:624008 | 6.82 | 58
459225 | Dibble | 55 | 459225:624007 | 12.36 | 86
459226 | Millsholm | 40 | 459226:624013 | 6.82 | 58
459226 | Dibble | 45 | 459226:624012 | 12.36 | 86
459227 | Millsholm | 30 | 459227:624019 | 6.82 | 58
459227 | Dibble | 60 | 459227:624018 | 12.36 | 99
...
2. Reduce to 1:1 relationship with mukey by aggreigating based on component percent:
SELECT mukey,
-- note that weights from non-soil components must be removed
-- otherwise, weighted mean values will be too low
SUM(comppct_r * component_whc) / SUM(comppct_r) AS comppct_weighted_whc_cm
FROM component
JOIN
(
SELECT cokey, sum( (hzdepb_r - hzdept_r) * awc_r) AS component_whc,
sum((hzdepb_r - hzdept_r)) AS depth
FROM chorizon
WHERE areasymbol = 'ca113'
GROUP BY cokey
) AS a
USING (cokey)
WHERE component.areasymbol = 'ca113'
-- filter out components that are missing soils data
AND a.component_whc IS NOT NULL
GROUP BY mukey ;
mukey | comppct_weighted_whc_cm
--------+-------------------------
459225 | 10
459226 | 10
459227 | 11
...
Link to polygons and setup access permissions:
-- create the new table with both geometry and attributes
CREATE TABLE yolo_whc AS
SELECT ogc_fid, wkb_geometry AS wkb_geometry, b.mukey, b.comppct_weighted_whc_cm
FROM mapunit_poly
-- use LEFT JOIN to include non-soil polygons in the result set
-- alternatively use JOIN to ignore non-soil polygons
LEFT JOIN
(
SELECT mukey,
-- note that weights from non-soil components must be removed
-- otherwise, weighted mean values will be too low
SUM(comppct_r * component_whc) / SUM(comppct_r) AS comppct_weighted_whc_cm
FROM component
JOIN
(
SELECT cokey, sum( (hzdepb_r - hzdept_r) * awc_r) AS component_whc,
sum((hzdepb_r - hzdept_r)) AS depth
FROM chorizon
WHERE areasymbol = 'ca113'
GROUP BY cokey
) AS a
USING (cokey)
WHERE component.areasymbol = 'ca113'
-- filter out components that are missing soils data
AND a.component_whc IS NOT NULL
GROUP BY mukey
) AS b
-- JOIN constraint
USING (mukey)
-- optional constraint to limit geometry search in mapunit_poly table
WHERE mapunit_poly.areasymbol = 'ca113' ;
Create indexes and register the new geometry:
-- create attribute and spatial index:
CREATE UNIQUE INDEX yolo_whc_idx ON yolo_whc (ogc_fid) ;
CREATE INDEX yolo_whc_spatial_idx ON yolo_whc USING gist (wkb_geometry gist_geometry_ops);
-- register in geometry_columns table:
INSERT INTO geometry_columns VALUES ('','public','yolo_whc','wkb_geometry',2,9001,'POLYGON');
-- optional:
-- fix permissions
-- GRANT SELECT on table yolo_whc to [user] ;
-- cleanup
vacuum analyze yolo_whc ;
Yolo County Water Holding Capacity Map: Profile total water holding capacity as computed by component percentage weighted average.
4. Dump to shapefile with ogr2ogr:
ogr2ogr ywhc.shp PG:"dbname=ssurgo_combined user=xxxx password=xxxx" yolo_whc
