Soil Structural Data Comparison#
This notebook simply compares several different sets of soil structural data, clarifying how they can be used and noting their imperfections.
To run an integrated hydrology simulation, we need at least:
porosity [-]
permeability [m^2]
van Genuchten properties (n [-], alpha [Pa^-1], and residual saturation [-])
across the entire domain. These vary both spatially and with depth.
There are two ways to specify this data:
as a raster, with every cell having its own value
as an indexed field, which combines a “structural model” that identifies which media each cell is in, plus a lookup table mapping media to properties.
For a structural model, we typically define two or three layers: bedrock, a geologic layer, and a soil layer. For the horizontal component, it is most convenient to work with formations – horizontal shapes whose bodies are assumed to be homogeneous within the shape and layer.
We can work with fields of porosity and permeability. Currently ATS does not support fields of van Genuchten parameters, but this could be changed relatively simply.
There is no one right or best model! This notebook exists to let the user set their domain and then explore the various datasets to see which look the best on their domain.
# these can be turned on for development work
%load_ext autoreload
%autoreload 2
import watershed_workflow.ui
watershed_workflow.ui.setup_logging(1)
import numpy as np
from matplotlib import pyplot as plt
import watershed_workflow
import watershed_workflow.plot
import watershed_workflow.crs
import watershed_workflow.ui
import watershed_workflow.sources
# user input -- choose the HUC to work on
huc = '140200010204'
crs = watershed_workflow.crs.default_crs
# collect the HUCs from WBD
watershed = watershed_workflow.sources.ManagerWBD().getShapesByID([huc,], out_crs=crs)
# plot what we have so far -- an image of the HUC
watershed.exterior.plot(color='k')
2026-03-05 09:33:03,422 - root - INFO: fixing column: geometry
Available Products#
The following products are routinely available across the US (or globally).
GLHYMPS#
The Global Hydrogeology Maps product provides data by shapefile/formation, and provides:
porosity
permeability
depth-to-bedrock
This dataset is fairly complete, but there are some oddities – most notably much of the DTB field seems to be 0, and a significant subset of the porosity are 0 (mostly alluvial sediments, WW corrects this). But the data is pretty good and complete. However, it is missing WRM parameters (or equivalently, texture information).
Gleeson, T., Moosdorf, N., Hartmann, J., & van Beek, L. P. H. (2014). A glimpse beneath earth’s surface: GLobal HYdrogeology MaPS (GLHYMPS) of permeability and porosity. Geophysical Research Letters, 41(11), 3891–3898. https://doi.org/10.1029/2014GL059856
Huscroft, J., Gleeson, T., Hartmann, J., & Börker, J. (2018). Compiling and mapping global permeability of the unconsolidated and consolidated Earth: GLobal HYdrogeology MaPS 2.0 (GLHYMPS 2.0). Geophysical Research Letters, 45(4), 1897–1904. https://doi.org/10.1002/2017GL075860
# first let's look at GLYHMPS
# -- get the data
glhymps_mgr = watershed_workflow.sources.ManagerGLHYMPS()
gl = glhymps_mgr.getShapesByGeometry(watershed, out_crs=crs)
2026-03-05 09:33:04,108 - root - INFO: Shapes manager: buffered+snapped query box = (-11965520.459519299, 3933522.1022809995, -11864010.211718498, 4035032.3500817996)
2026-03-05 09:33:04,109 - root - INFO: from file: /Users/uec/code/watershed_workflow/data/soil_structure/GLHYMPS/GLHYMPS.shp
2026-03-05 09:33:11,950 - root - INFO: fixing column: geometry
# convert it to standard ATS datasets
gl = watershed_workflow.soil_properties.mangleGLHYMPSProperties(gl)
gl
| ID | name | source | permeability [m^2] | logk_stdev [-] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| 259 | 715639 | GLHYMPS-715639 | GLHYMPS | 6.309573e-16 | 2.50 | 0.19 | 0.000025 | 1.5 | 0.01 | MULTIPOLYGON (((-941668.093 1811291.15, -94199... |
| 278 | 715707 | GLHYMPS-715707 | GLHYMPS | 1.000000e-13 | 2.00 | 0.22 | 0.000294 | 1.5 | 0.01 | MULTIPOLYGON (((-942876.044 1814661.916, -9427... |
| 287 | 715766 | GLHYMPS-715766 | GLHYMPS | 3.162278e-13 | 1.80 | 0.09 | 0.000817 | 1.5 | 0.01 | POLYGON ((-947358.28 1813697.646, -948118.709 ... |
| 289 | 715779 | GLHYMPS-715779 | GLHYMPS | 1.000000e-13 | 2.00 | 0.22 | 0.000294 | 1.5 | 0.01 | POLYGON ((-944025.909 1813309.319, -944527.763... |
| 291 | 715796 | GLHYMPS-715796 | GLHYMPS | 6.309573e-16 | 2.50 | 0.19 | 0.000025 | 1.5 | 0.01 | MULTIPOLYGON (((-933660.561 1797293.054, -9335... |
| 894 | 726604 | GLHYMPS-726604 | GLHYMPS | 6.309573e-16 | 2.50 | 0.19 | 0.000025 | 1.5 | 0.01 | MULTIPOLYGON (((-941747.425 1812754.646, -9418... |
| 897 | 726608 | GLHYMPS-726608 | GLHYMPS | 6.309573e-16 | 2.50 | 0.19 | 0.000025 | 1.5 | 0.01 | MULTIPOLYGON (((-953864.501 1813999.086, -9542... |
| 917 | 726639 | GLHYMPS-726639 | GLHYMPS | 3.162278e-13 | 1.80 | 0.09 | 0.000817 | 1.5 | 0.01 | POLYGON ((-948778.786 1815596.525, -948024.977... |
| 918 | 726642 | GLHYMPS-726642 | GLHYMPS | 1.000000e-13 | 2.00 | 0.22 | 0.000294 | 1.5 | 0.01 | MULTIPOLYGON (((-944367.697 1813260.891, -9445... |
| 935 | 726664 | GLHYMPS-726664 | GLHYMPS | 1.000000e-13 | 2.00 | 0.22 | 0.000294 | 1.5 | 0.01 | POLYGON ((-952109.643 1813503.499, -951346.338... |
| 938 | 726667 | GLHYMPS-726667 | GLHYMPS | 3.162278e-13 | 1.80 | 0.09 | 0.000817 | 1.5 | 0.01 | POLYGON ((-940163.503 1810040.229, -940089.924... |
| 1042 | 730801 | GLHYMPS-730801 | GLHYMPS | 3.019952e-11 | 1.61 | NaN | NaN | 1.5 | 0.01 | POLYGON ((-945003.906 1820005.318, -944742.248... |
# -- print extent of properties
gl_logK = np.log10(gl['permeability [m^2]'][:])
print('logK', gl_logK.min(), gl_logK.max())
poro = gl['porosity [-]']
print('poro', poro.min(), poro.max())
logK -15.2 -10.52
poro 0.09 0.22
gl.crs
<Projected CRS: EPSG:5070>
Name: NAD83 / Conus Albers
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.
- bounds: (-124.79, 24.41, -66.91, 49.38)
Coordinate Operation:
- name: Conus Albers
- method: Albers Equal Area
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
# -- color a raster by the polygons to see the structure
gl_color = watershed_workflow.data.rasterizeGeoDataFrame(gl, 'ID', 10)
gl_color = gl_color.rio.clip(watershed.geometry, watershed.crs)
# plot the formations
plt.figure()
unique_inds, gl_cmap, gl_norm, gl_ticks, gl_labels = \
watershed_workflow.colors.createIndexedColormap(gl['ID'])
gl_color.plot.imshow(cmap=gl_cmap, norm=gl_norm)
plt.title("GLHMPS")
Text(0.5, 1.0, 'GLHMPS')
# generate glhymps porosity, permeability fields
# -- porosity
gl_poro = watershed_workflow.data.rasterizeGeoDataFrame(gl, 'porosity [-]', 10, nodata=np.nan)
gl_poro = gl_poro.rio.clip(watershed.geometry, watershed.crs)
# -- permeability
gl_perm = watershed_workflow.data.rasterizeGeoDataFrame(gl, 'permeability [m^2]', 10, nodata=np.nan)
gl_perm = gl_perm.rio.clip(watershed.geometry, watershed.crs)
# plot porosity and perm for glhymps
fig, axs = plt.subplots(1,2, figsize=(15,6))
gl_poro.plot.imshow(ax=axs[0])
gl_perm.plot.imshow(ax=axs[1])
axs[0].set_title('porosity [-]')
axs[1].set_title('log permeability [m^2]')
Text(0.5, 1.0, 'log permeability [m^2]')
Pelletier et al (2016)#
In Pelletier et al 2016, a global depth-to-bedrock raster was developed. This has provided a reliable source for data about deep bedrock.
Pelletier, Jon D., et al. “A gridded global data set of soil, intact regolith, and sedimentary deposit thicknesses for regional and global land surface modeling.” Journal of Advances in Modeling Earth Systems 8.1 (2016): 41-65.
pell_mgr = watershed_workflow.sources.ManagerPelletierDTB()
pell = pell_mgr.getDataset(watershed.geometry[0], watershed.crs)['band_1']
pell = pell.rio.clip(watershed.to_crs(pell.rio.crs).geometry)
print(pell)
2026-03-05 09:33:12,689 - root - INFO: from file: /Users/uec/code/watershed_workflow/data/soil_structure/PelletierDTB/Global_Soil_Regolith_Sediment_1304/data/average_soil_and_sedimentary-deposit_thickness.tif
2026-03-05 09:33:12,720 - root - INFO: Incoming shape area = 0.005521611592362775
2026-03-05 09:33:12,720 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.025000000000019895
2026-03-05 09:33:12,722 - root - INFO: ... buffered shape area = 0.015938643341995972
2026-03-05 09:33:12,722 - root - INFO: ... snapped bounding box = (-107.13333333341859, 38.80000000003088, -106.95000000008511, 38.92500000003098)
<xarray.DataArray 'band_1' (y: 8, x: 15)> Size: 480B
array([[nan, nan, nan, 0., 0., nan, nan, 0., 0., 1., 1., nan, nan,
nan, nan],
[nan, 2., 1., 0., 0., 0., 0., 0., 1., 1., 0., 0., nan,
nan, nan],
[nan, 0., 1., 1., 0., 0., 0., 1., 1., 1., 1., 1., 1.,
nan, 2.],
[ 1., 1., 0., 0., 1., 0., 1., 1., 1., 1., 1., 1., 1.,
0., nan],
[nan, 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., nan],
[nan, 1., 1., 1., 0., 0., 0., 0., 0., 1., 1., 1., 1.,
0., nan],
[nan, nan, 1., 1., 0., 0., 1., 1., nan, nan, 1., 1., 0.,
nan, nan],
[nan, nan, nan, 0., nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan]], dtype=float32)
Coordinates:
* y (y) float64 64B 38.89 38.88 38.87 38.86 38.85 38.85 38.84 38.83
* x (x) float64 120B -107.1 -107.1 -107.1 ... -107.0 -107.0 -107.0
spatial_ref int64 8B 0
Attributes:
DataType: Generic
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0
# plot Pelletier DTB
fig, ax = plt.subplots(1, 1)
pell.plot.imshow(ax=ax)
<matplotlib.image.AxesImage at 0x3202c88c0>
SSURGO#
The NRCS’s SSURGO/STATSGO products provide shape-based information about the top 2m of soil. This data is provided by formation, and includes multiple formations (map-units), in which each map unit has (potentially) multiple components and each component has (potentially) multiple horizons. Short of taking the most prevelant component, we must average, so we tend to average horizons and components and just apply this uniformly in the top 2m. Provided include a complete set of parameters:
porosity
permeability
bulk density
texture
and more…
Unfortunately, while map units seem to cover the area of interest, they don’t all have at least one horizon with valid data, so some formations end up with NaN values. So while SSURGO can be used as an “overlay” with good soil layer information, it cannot be used alone without a gap-filling plan.
Additionally, some locations have texture data but not porosity or permeability – we can use Rosetta to simulate these properties from the texture info.
Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture. Soil Survey Geographic (SSURGO) Database. https://sdmdataaccess.nrcs.usda.gov
Schaap, M. G., Leij, F. J., & van Genuchten, M. T. (2001). Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology, 251(3–4), 163–176. https://doi.org/10.1016/S0022-1694(01)00466-8
# ssurgo
nrcs_mgr = watershed_workflow.sources.ManagerNRCS()
nrcs = nrcs_mgr.getShapesByGeometry(watershed.geometry[0], watershed.crs)
nrcs
2026-03-05 09:33:12,874 - root - INFO: Shapes manager: buffered+snapped query box = (-107.11, 38.824, -106.974, 38.898)
2026-03-05 09:33:12,875 - root - INFO: Using superset cache: /Users/uec/code/watershed_workflow/data/soil_structure/SSURGO/SSURGO_-107.1100_38.8240_-106.9740_38.8980.gpkg
2026-03-05 09:33:12,890 - root - INFO: fixing column: geometry
| mukey | residual saturation [-] | Rosetta porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | Rosetta permeability [m^2] | thickness [m] | permeability [m^2] | porosity [-] | bulk density [g/cm^3] | total sand pct [%] | total silt pct [%] | total clay pct [%] | source | geometry | ID | name | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 498185 | 0.240269 | 0.502882 | 0.000092 | 1.319916 | 3.617805e-13 | 1.520000 | 2.279759e-13 | 0.265563 | 1.207632 | 32.946001 | 24.529703 | 42.524296 | NRCS | POLYGON ((-107.01831 38.88358, -107.01793 38.8... | 498185 | NRCS-498185 |
| 2 | 498205 | 0.181598 | 0.399993 | 0.000139 | 1.457559 | 4.938634e-13 | 1.520000 | 2.734080e-12 | 0.159737 | 1.413553 | 63.518421 | 22.988158 | 13.493421 | NRCS | MULTIPOLYGON (((-106.97469 38.86306, -106.9758... | 498205 | NRCS-498205 |
| 3 | 498206 | 0.181598 | 0.399993 | 0.000139 | 1.457559 | 4.938634e-13 | 1.520000 | 2.734080e-12 | 0.159737 | 1.413553 | 63.518421 | 22.988158 | 13.493421 | NRCS | MULTIPOLYGON (((-106.97863 38.87193, -106.9793... | 498206 | NRCS-498206 |
| 4 | 498208 | 0.208316 | 0.428990 | 0.000071 | 1.420938 | 2.830658e-13 | 1.520000 | 1.184189e-12 | 0.185986 | 1.302697 | 40.662797 | 38.201188 | 21.136015 | NRCS | MULTIPOLYGON (((-106.98546 38.87211, -106.9855... | 498208 | NRCS-498208 |
| 7 | 498231 | 0.213285 | 0.466663 | 0.000059 | 1.416784 | 3.782943e-13 | 1.070000 | 7.331979e-13 | 0.403218 | 1.197010 | 32.639846 | 41.078960 | 26.281194 | NRCS | MULTIPOLYGON (((-107.01992 38.86516, -107.0199... | 498231 | NRCS-498231 |
| 8 | 509477 | 0.175133 | 0.567740 | 0.000103 | 1.383180 | 2.438575e-12 | 1.520000 | 3.039095e-12 | 0.428779 | 0.807465 | 62.072368 | 17.973684 | 19.953947 | NRCS | MULTIPOLYGON (((-107.08069 38.82938, -107.0808... | 509477 | NRCS-509477 |
| 9 | 509479 | 0.230885 | 0.378365 | 0.000079 | 1.374540 | 1.031721e-13 | 1.250000 | 3.172415e-12 | NaN | 1.524444 | 39.279352 | 38.916498 | 21.804150 | NRCS | MULTIPOLYGON (((-107.07093 38.824, -107.07265 ... | 509479 | NRCS-509479 |
| 10 | 509481 | 0.230885 | 0.378365 | 0.000079 | 1.374540 | 1.031721e-13 | 1.250000 | 3.172415e-12 | NaN | 1.524444 | 39.279352 | 38.916498 | 21.804150 | NRCS | MULTIPOLYGON (((-107.10982 38.83356, -107.1094... | 509481 | NRCS-509481 |
| 11 | 509482 | 0.207877 | 0.369886 | 0.000092 | 1.397760 | 1.474479e-13 | 1.205000 | 4.654687e-12 | NaN | 1.530585 | 45.389069 | 37.890283 | 16.720648 | NRCS | MULTIPOLYGON (((-107.07674 38.82524, -107.0765... | 509482 | NRCS-509482 |
| 12 | 509513 | 0.219597 | 0.445168 | 0.000072 | 1.395746 | 2.727324e-13 | 1.520000 | 9.122162e-13 | 0.237423 | 1.280039 | 38.422252 | 35.742028 | 25.835720 | NRCS | MULTIPOLYGON (((-107.11 38.87229, -107.10976 3... | 509513 | NRCS-509513 |
| 13 | 509514 | 0.241801 | 0.474101 | 0.000079 | 1.340429 | 2.435161e-13 | 1.520000 | 2.375523e-13 | 0.271630 | 1.274496 | 30.592982 | 31.535307 | 37.871711 | NRCS | MULTIPOLYGON (((-107.00813 38.82402, -107.0080... | 509514 | NRCS-509514 |
| 14 | 509529 | 0.196811 | 0.356163 | 0.000157 | 1.406091 | 2.260935e-13 | 1.520000 | 5.850942e-12 | NaN | 1.607097 | 63.571711 | 22.369079 | 14.059211 | NRCS | MULTIPOLYGON (((-107.06266 38.87676, -107.0627... | 509529 | NRCS-509529 |
| 15 | 509532 | 0.199783 | 0.412306 | 0.000087 | 1.422578 | 2.966015e-13 | 1.520000 | 1.526682e-12 | 0.182238 | 1.356006 | 47.934946 | 34.134907 | 17.930147 | NRCS | MULTIPOLYGON (((-107.08744 38.82506, -107.0874... | 509532 | NRCS-509532 |
| 16 | 509544 | NaN | NaN | NaN | NaN | NaN | 1.520000 | 7.000000e-15 | NaN | NaN | NaN | NaN | NaN | NRCS | MULTIPOLYGON (((-107.04628 38.87972, -107.0465... | 509544 | NRCS-509544 |
| 17 | 509547 | NaN | NaN | NaN | NaN | NaN | 1.520000 | 4.230700e-11 | NaN | 2.030000 | NaN | NaN | NaN | NRCS | MULTIPOLYGON (((-107.06984 38.82436, -107.0695... | 509547 | NRCS-509547 |
| 18 | 509548 | 0.189314 | 0.404765 | 0.000104 | 1.433685 | 3.617371e-13 | 1.520000 | 2.006536e-12 | 0.193350 | 1.379123 | 53.676535 | 31.072368 | 15.251096 | NRCS | MULTIPOLYGON (((-107.04637 38.82679, -107.0460... | 509548 | NRCS-509548 |
| 20 | 509559 | 0.210870 | 0.420608 | 0.000078 | 1.410010 | 2.516507e-13 | 0.618125 | 7.504374e-13 | 0.151250 | 1.341513 | 43.240132 | 35.652961 | 21.106908 | NRCS | MULTIPOLYGON (((-107.06648 38.88147, -107.0667... | 509559 | NRCS-509559 |
| 21 | 509561 | 0.246665 | 0.467696 | 0.000081 | 1.332817 | 2.088819e-13 | 1.396471 | 1.590777e-13 | 0.218625 | 1.305591 | 30.228443 | 31.026630 | 38.744927 | NRCS | POLYGON ((-107.09952 38.8459, -107.09938 38.84... | 509561 | NRCS-509561 |
| 24 | 509577 | NaN | NaN | NaN | NaN | NaN | 1.520000 | NaN | NaN | NaN | NaN | NaN | NaN | NRCS | MULTIPOLYGON (((-107.0234 38.83592, -107.02335... | 509577 | NRCS-509577 |
| 27 | 509733 | 0.157366 | 0.400499 | 0.000221 | 1.762072 | 1.819618e-12 | 0.860000 | 5.262229e-12 | 0.449752 | 1.423158 | 82.200000 | 9.300000 | 8.500000 | NRCS | MULTIPOLYGON (((-106.97431 38.87336, -106.9743... | 509733 | NRCS-509733 |
| 30 | 509794 | 0.220962 | 0.383130 | 0.000140 | 1.367408 | 2.004353e-13 | 1.520000 | 1.554090e-12 | 0.415245 | 1.545395 | 60.571711 | 18.441447 | 20.986842 | NRCS | POLYGON ((-106.97543 38.8732, -106.97545 38.87... | 509794 | NRCS-509794 |
# plot ssurgo porosity, note white inside the domain is NaN values
fig, axs = plt.subplots(1,2, figsize=(12,5))
nrcs.plot('porosity [-]', ax=axs[0], legend=True)
axs[0].set_title('porosity [-]')
nrcs.plot('Rosetta porosity [-]', ax=axs[1], legend=True)
axs[1].set_title('Rosetta porosity [-]')
Text(0.5, 1.0, 'Rosetta porosity [-]')
import matplotlib.colors as mpc
# plot ssurgo porosity, note white inside the domain is NaN values
fig, axs = plt.subplots(1,2, figsize=(12,5))
nrcs.plot('permeability [m^2]', ax=axs[0], legend=True,
norm=mpc.LogNorm(vmin=1e-16, vmax=1e-10))
axs[0].set_title('permeability [m^2]')
nrcs.plot('Rosetta permeability [m^2]', ax=axs[1], legend=True,
norm=mpc.LogNorm(vmin=1e-16, vmax=1e-10))
axs[1].set_title('Rosetta permeability [m^2]')
Text(0.5, 1.0, 'Rosetta permeability [m^2]')
SoilGrids (2017)#
SoilGrids includes 250m, global rasters of:
depth to bedrock
bulk density
texture via percent sand/silt/clay
and more…
Effectively everything is available that is needed to run Rosetta to get a full suite of parameters. Vertically, bulk density and texture information are provided at 7 layers, the deepest starting at 2m. Data is complete globally, and doesn’t seem to have much missing data except under water bodies. The only downside is the raster nature of the data – to use this directly we would have to run every pixel through Rosetta and provide fields of Van Genuchten parameters. We could conceivably use the rasters of porosity and permeability, but it seems more convenient to work by formation instead of as a raster.
Note that here we use the 2017 product instead of the 2020 product as that has yet to include depth-to-bedrock. The depth-to-bedrock here seems fairly inaccurate – generally we have stopped using SoilGrids for this reason.
Therefore we keep this section here in case anyone wants to play with this product, but comment it out.
Hengl, T., et al. (2017). SoilGrids250m: Global gridded soil information based on machine learning. PLOS ONE, 12(2), e0169748. https://doi.org/10.1371/journal.pone.0169748
# SoilGrids
#
# NOTE: this may not work any more -- we may need to deprecate SoilGrids.
#sg_mgr = watershed_workflow.sources.ManagerSoilGrids2017()
#sg = sg_mgr.getDataset(watershed)
#print(sg)
SoilGrids has holes where there is standing water, but we still need soil in these places. Fill these using a nearest neighbor algorithm.
HydroFrame HF-Hydrodata CONUS2#
The HydroFrame project provides the CONUS2 subsurface parameterization via the hf_hydrodata package. The primary product is pf_indicator, a 3D integer grid (10 z-layers × y × x) on a 1-km Lambert Conformal Conic grid covering the continental US. Each integer maps to a hydrogeologic unit with known properties (porosity, permeability, van Genuchten parameters, specific storage).
Compared to SSURGO and GLHYMPS:
Coverage: CONUS-wide, no gaps, no missing data within the active domain
Resolution: 1 km (coarser than SSURGO but complete)
Depth: 10 layers to ~100 m depth
Properties: a full property set is available via
getIndicatorTable()— no Rosetta neededLimitation: requires a free HydroFrame account and PIN; PIN expires after a few days of non-use
This dataset is particularly well-suited as a bedrock/deep-subsurface layer to complement SSURGO’s shallow soil information.
Setup: add to ~/.watershed_workflowrc:
[HFHydrodata]
email = user@example.com
pin = 1234
Generate a PIN at https://hydrogen.princeton.edu/pin.
Maxwell, R. M., et al. (2015). The imprint of climate and geology on the residence times of groundwater. Geophysical Research Letters, 43(2), 701–708. https://doi.org/10.1002/2015GL066916
Yang, C., et al. (2023). CONUS2.0: A new integrated hydrologic simulation of the contiguous United States at high resolution. Geoscientific Model Development. https://doi.org/10.5194/gmd-2023-205
Leonarduzzi, E., et al. (2024). HydroData: A multi-scale hydrology dataset for the contiguous United States. Earth and Space Science, 11, e2023EA003021. https://doi.org/10.1029/2023EA003021
# Download the CONUS2 pf_indicator for the watershed
hfhd_mgr = watershed_workflow.sources.ManagerHFHydrodata()
hfhd_ds = hfhd_mgr.getDataset(watershed.geometry[0], watershed.crs)
hfhd_indicator = hfhd_ds['pf_indicator']
print(hfhd_indicator)
2026-03-05 09:33:13,937 - root - INFO: Incoming shape area = 0.005521611592362775
2026-03-05 09:33:13,937 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.026999999999999996
2026-03-05 09:33:13,939 - root - INFO: ... buffered shape area = 0.01692038260373745
2026-03-05 09:33:13,939 - root - INFO: ... snapped bounding box = (-107.136, 38.799, -106.94699999999999, 38.925)
<xarray.DataArray 'pf_indicator' (z: 10, y: 11, x: 17)> Size: 7kB
[1870 values with dtype=float32]
Coordinates:
* z (z) int32 40B 0 1 2 3 4 5 6 7 8 9
* y (y) float64 88B -7.55e+04 -7.45e+04 ... -6.65e+04 -6.55e+04
* x (x) float64 136B -8.515e+05 -8.505e+05 ... -8.355e+05
spatial_ref int64 8B 0
# Plot the indicator in the top-most layer (z=0 = shallowest)
# Reproject watershed boundary to CONUS2 LCC for overlay
watershed_lcc = watershed.to_crs(hfhd_mgr.native_crs_out)
# Build a discrete colormap and norm for the integer indicator values
present_ids = sorted(set(int(v) for v in hfhd_indicator.values.ravel() if not np.isnan(v)))
n = len(present_ids)
cmap = plt.cm.get_cmap('tab20', n)
import matplotlib.colors as mcolors
boundaries = [v - 0.5 for v in present_ids] + [present_ids[-1] + 0.5]
norm = mcolors.BoundaryNorm(boundaries, ncolors=n)
fig, axs = plt.subplots(1, 2, figsize=(14, 5))
for ax, z, title in [
(axs[0], 0, 'CONUS2 pf_indicator — top layer (z=0)'),
(axs[1], -1, 'CONUS2 pf_indicator — bottom layer (z=9)'),
]:
im = ax.imshow(
hfhd_indicator.isel(z=z).values,
cmap=cmap, norm=norm,
extent=[hfhd_indicator.x.values[0] - 500, hfhd_indicator.x.values[-1] + 500,
hfhd_indicator.y.values[0] - 500, hfhd_indicator.y.values[-1] + 500],
origin='upper',
)
watershed_lcc.boundary.plot(ax=ax, color='k', linewidth=1.5)
ax.set_title(title)
cbar = fig.colorbar(im, ax=ax, ticks=present_ids)
cbar.ax.set_yticklabels([str(v) for v in present_ids])
cbar.set_label('indicator')
plt.tight_layout()
/var/folders/fj/y8h3hc4x3fx9jtychppwtybr0000gp/T/ipykernel_54938/1298079568.py:8: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed in 3.11. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap()`` or ``pyplot.get_cmap()`` instead.
cmap = plt.cm.get_cmap('tab20', n)
# Show the property lookup table for the indicators present in this watershed
table = hfhd_mgr.getIndicatorTable()
present = sorted(hfhd_indicator.values.astype(int).ravel().tolist())
present_ids = sorted(set(present))
table.loc[table.index.isin(present_ids)]
| porosity | permeability_x | permeability_y | permeability_z | vg_alpha | vg_n | sres | specific_storage | |
|---|---|---|---|---|---|---|---|---|
| indicator | ||||||||
| 2 | 0.390 | 0.043630 | 0.043630 | 0.043630 | 3.467 | 2.738 | 0.00010 | 0.0001 |
| 3 | 0.387 | 0.015841 | 0.015841 | 0.015841 | 2.692 | 2.445 | 0.00010 | 0.0001 |
| 6 | 0.399 | 0.005009 | 0.005009 | 0.005009 | 1.122 | 2.479 | 0.00010 | 0.0001 |
| 7 | 0.384 | 0.005493 | 0.005493 | 0.005493 | 2.089 | 2.318 | 0.00010 | 0.0001 |
| 9 | 0.442 | 0.003387 | 0.003387 | 0.003387 | 1.585 | 2.413 | 0.00010 | 0.0001 |
| 19 | 0.050 | 0.005000 | 0.005000 | 0.000500 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 20 | 0.100 | 0.010000 | 0.010000 | 0.001000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 21 | 0.120 | 0.020000 | 0.020000 | 0.002000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 22 | 0.300 | 0.030000 | 0.030000 | 0.003000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 23 | 0.010 | 0.040000 | 0.040000 | 0.040000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 24 | 0.150 | 0.050000 | 0.050000 | 0.005000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 25 | 0.220 | 0.060000 | 0.060000 | 0.006000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
| 26 | 0.270 | 0.080000 | 0.080000 | 0.008000 | 0.500 | 2.500 | 0.00001 | 0.0001 |
# Derive porosity and permeability_x fields by mapping indicator -> table value
import numpy as np
import xarray as xr
ind_top = hfhd_indicator.isel(z=0)
poro_vals = table['porosity'].to_dict()
perm_vals = table['permeability_x'].to_dict()
poro_data = np.vectorize(poro_vals.get)(ind_top.values.astype(int))
perm_data = np.vectorize(perm_vals.get)(ind_top.values.astype(int))
poro_da = xr.DataArray(poro_data.astype(np.float32), coords=ind_top.coords, dims=ind_top.dims,
name='porosity')
perm_da = xr.DataArray(perm_data.astype(np.float32), coords=ind_top.coords, dims=ind_top.dims,
name='permeability_x [m/h]')
fig, axs = plt.subplots(1, 2, figsize=(14, 5))
poro_da.plot.imshow(ax=axs[0], cmap='YlOrBr', add_colorbar=True)
watershed_lcc.boundary.plot(ax=axs[0], color='k', linewidth=1.5)
axs[0].set_title('CONUS2 porosity — top layer (z=0)')
perm_da.plot.imshow(ax=axs[1], cmap='Blues', add_colorbar=True)
watershed_lcc.boundary.plot(ax=axs[1], color='k', linewidth=1.5)
axs[1].set_title('CONUS2 permeability_x [m/h] — top layer (z=0)')
plt.tight_layout()
SoilGrids 2.0#
SoilGrids 2.0 provides 250 m global gridded soil properties at six standard GlobalSoilMap depth intervals (0–5, 5–15, 15–30, 30–60, 60–100, 100–200 cm). Data are delivered via the ISRIC OGC WCS — no authentication required, CC-BY 4.0 license.
Available variables include bulk density (bdod), clay/sand/silt content, soil organic carbon (soc), pH in water (phh2o), coarse fragments (cfvo), nitrogen, and organic carbon density (ocd). The default request fetches the four texture variables needed for Rosetta (clay, sand, silt, bdod) and automatically computes van Genuchten parameters pixel-by-pixel.
Compared to SSURGO and GLHYMPS:
Coverage: global, no gaps within land areas
Resolution: 250 m raster (finer than CONUS2, coarser than SSURGO map units)
Depth: six layers to 2 m
Properties: full VG parameter set derived via Rosetta; no depth-to-bedrock in v2.0 (use Pelletier for that)
Limitation: purely statistical model — accuracy varies by region and depth
Poggio, L., et al. (2021). SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. SOIL, 7, 217–240. https://doi.org/10.5194/soil-7-217-2021
# Download SoilGrids 2.0 default variables (clay, sand, silt, bdod) for the watershed.
# Rosetta is run automatically because all four texture variables are present.
sg_mgr = watershed_workflow.sources.ManagerSoilGrids()
sg_ds = sg_mgr.getDataset(watershed.geometry[0], watershed.crs)
print(sg_ds)
2026-03-05 09:33:14,574 - root - INFO: Incoming shape area = 0.005521611592362775
2026-03-05 09:33:14,574 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.006
2026-03-05 09:33:14,575 - root - INFO: ... buffered shape area = 0.007736872189567648
2026-03-05 09:33:14,575 - root - INFO: ... snapped bounding box = (-107.114, 38.82, -106.97, 38.902)
2026-03-05 09:33:14,585 - root - INFO: Running Rosetta on SoilGrids texture rasters
2026-03-05 09:33:14,586 - root - INFO: Running Rosetta for van Genutchen parameters
2026-03-05 09:33:18,106 - root - INFO: ... done
<xarray.Dataset> Size: 639kB
Dimensions: (depth: 6, y: 41, x: 72)
Coordinates:
* depth (depth) float32 24B 0.025 0.1 0.225 ... 0.8 1.5
* y (y) float64 328B 38.9 38.9 38.9 ... 38.82 38.82
* x (x) float64 576B -107.1 -107.1 ... -107.0
spatial_ref int64 8B 0
Data variables:
clay (depth, y, x) float32 71kB 21.4 21.9 ... 10.7
sand (depth, y, x) float32 71kB 41.3 37.1 ... 61.2
silt (depth, y, x) float32 71kB 37.3 41.1 ... 28.1
bdod (depth, y, x) float32 71kB 1.07 1.05 ... 1.55
residual saturation [-] (depth, y, x) float32 71kB 0.1921 ... 0.1789
porosity [-] (depth, y, x) float32 71kB 0.488 ... 0.362
van Genuchten alpha [Pa^-1] (depth, y, x) float32 71kB 6.219e-05 ... 0.0...
van Genuchten n [-] (depth, y, x) float32 71kB 1.437 ... 1.445
permeability [m^2] (depth, y, x) float32 71kB 8.046e-13 ... 3.1...
Attributes:
source: https://maps.isric.org
variable: clay
name: SoilGrids 2.0
# Plot clay content and porosity at the shallowest and deepest layers
# Rows: clay (top), porosity (bottom). Columns: 0-5 cm (left), 100-200 cm (right).
watershed_ll = watershed.to_crs(sg_ds.rio.crs)
fig, axs = plt.subplots(2, 2, figsize=(14, 10))
plots = [
(0, 0, 'clay', 'Clay [%]', 'YlOrBr', 0),
(0, 1, 'clay', 'Clay [%]', 'YlOrBr', -1),
(1, 0, 'porosity [-]', 'Porosity [-]','Blues', 0),
(1, 1, 'porosity [-]', 'Porosity [-]','Blues', -1),
]
depth_titles = ['0–5 cm', '100–200 cm']
for row, col, var, label, cmap, z_idx in plots:
ax = axs[row, col]
sg_ds[var].isel(depth=z_idx).plot.imshow(ax=ax, cmap=cmap, add_colorbar=True)
watershed_ll.boundary.plot(ax=ax, color='k', linewidth=1.5)
ax.set_title(f'SoilGrids2 {label} — {depth_titles[col]}')
plt.tight_layout()
# Plot log10(permeability) at the shallowest and deepest layers
import matplotlib.colors as mpc
fig, axs = plt.subplots(1, 2, figsize=(14, 5))
for ax, z_idx, title in [
(axs[0], 0, 'SoilGrids2 permeability [m²] — 0–5 cm'),
(axs[1], -1, 'SoilGrids2 permeability [m²] — 100–200 cm'),
]:
perm = sg_ds['permeability [m^2]'].isel(depth=z_idx)
perm.plot.imshow(ax=ax, cmap='viridis',
norm=mpc.LogNorm(vmin=1e-14, vmax=1e-11),
add_colorbar=True)
watershed_ll.boundary.plot(ax=ax, color='k', linewidth=1.5)
ax.set_title(title)
plt.tight_layout()
POLARIS (30 m CONUS)#
POLARIS provides 30-m resolution, spatially continuous, probabilistic maps of soil hydraulic properties for the contiguous United States (CONUS). Unlike SSURGO polygon-based data or CONUS2 indicator classes, POLARIS delivers continuous float fields directly — van Genuchten parameters (θ_r, θ_s, α, n) and saturated hydraulic conductivity (K_sat) are pre-computed at each 30-m pixel without needing Rosetta.
Data are served as plain-HTTP GeoTIFFs tiled by 1-degree cells. No authentication is required.
Available variables include clay, sand, silt, organic matter (), bulk density (), pH, bubbling pressure (), and the full set of hydraulic parameters: (log₁₀ cm/hr), (log₁₀ cm⁻¹), , , .
Compared to SoilGrids 2.0:
Coverage: CONUS only (SoilGrids is global)
Resolution: 30 m raster (much finer than SoilGrids 250 m)
Depth: same six layers to 2 m
Properties: VG parameters provided directly — no Rosetta step needed
Units note: and are stored as log₁₀ values
License: CC BY-NC 4.0
Chaney, N.W., Minasny, B., Herman, J.D., Nauman, T.W., Brungard, C.W., Morgan, C.L.S., McBratney, A.B., Wood, E.F., and Yimam, Y. (2019). POLARIS soil properties: 30-m probabilistic maps of soil properties over the contiguous United States. Water Resources Research, 55, 2916–2938. https://doi.org/10.1029/2018WR022797
# Download POLARIS default variables (theta_s, theta_r, alpha, n, ksat) for the watershed.
# VG parameters are provided directly — no Rosetta step needed.
polaris_mgr = watershed_workflow.sources.ManagerPOLARIS()
polaris_ds = polaris_mgr.getDataset(watershed.geometry[0], watershed.crs)
print(polaris_ds)
2026-03-05 09:33:18,935 - root - INFO: Incoming shape area = 0.005521611592362775
2026-03-05 09:33:18,935 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.0008339999999999999
2026-03-05 09:33:18,936 - root - INFO: ... buffered shape area = 0.005821068049095576
2026-03-05 09:33:18,936 - root - INFO: ... snapped bounding box = (-107.10728399999999, 38.82687, -106.97634599999999, 38.895813999999994)
2026-03-05 09:33:18,937 - root - INFO: POLARIS: downloading theta_s/mean over 2 tile(s)
2026-03-05 09:33:18,937 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/0_5/lat3839_lon-108-107.tif
2026-03-05 09:33:25,318 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/0_5/lat3839_lon-107-106.tif
2026-03-05 09:33:30,080 - rasterio._env - WARNING: CPLE_IllegalArg in 5ed145d3-eb23-4bc7-be6d-bba6d7a8c7ee.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:30,416 - rasterio._env - WARNING: CPLE_IllegalArg in c322fb52-a548-4d7b-b48c-f64e8f30c4e5.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:31,070 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/5_15/lat3839_lon-108-107.tif
2026-03-05 09:33:33,071 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/5_15/lat3839_lon-107-106.tif
2026-03-05 09:33:35,474 - rasterio._env - WARNING: CPLE_IllegalArg in b3244f41-a2dc-48d0-bcaf-3fa6d1964238.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:35,816 - rasterio._env - WARNING: CPLE_IllegalArg in 20ae0b80-960d-4882-833e-a648eb91f545.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:36,476 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/15_30/lat3839_lon-108-107.tif
2026-03-05 09:33:38,097 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/15_30/lat3839_lon-107-106.tif
2026-03-05 09:33:41,597 - rasterio._env - WARNING: CPLE_IllegalArg in 783f3a0e-a132-4f37-9460-a0805175e019.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:41,943 - rasterio._env - WARNING: CPLE_IllegalArg in 9e242562-9155-49b6-ba68-49c1be391e1b.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:42,601 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/30_60/lat3839_lon-108-107.tif
2026-03-05 09:33:45,762 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/30_60/lat3839_lon-107-106.tif
2026-03-05 09:33:47,910 - rasterio._env - WARNING: CPLE_IllegalArg in 6514e9a4-3b22-4693-a378-6689983cd0f8.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:48,250 - rasterio._env - WARNING: CPLE_IllegalArg in 0111f121-6165-40e3-b4eb-5469d15a9125.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:48,898 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/60_100/lat3839_lon-108-107.tif
2026-03-05 09:33:51,381 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/60_100/lat3839_lon-107-106.tif
2026-03-05 09:33:57,328 - rasterio._env - WARNING: CPLE_IllegalArg in 9564eb3a-7fc2-4333-b8b0-75ed0f12e4b8.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:57,662 - rasterio._env - WARNING: CPLE_IllegalArg in 9fdf8cd9-9d25-40f0-9fbc-b2ec5fd010fb.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:33:58,310 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/100_200/lat3839_lon-108-107.tif
2026-03-05 09:34:04,626 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_s/mean/100_200/lat3839_lon-107-106.tif
2026-03-05 09:34:08,022 - rasterio._env - WARNING: CPLE_IllegalArg in 2df93964-f099-4a86-a33e-f7ad95411322.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:08,355 - rasterio._env - WARNING: CPLE_IllegalArg in e46beef1-6446-4b00-b062-2a34cdc3f164.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:09,128 - root - INFO: Written to: /Users/uec/code/watershed_workflow/data/soil_structure/POLARIS/POLARIS_theta_s_mean_-107.1073_38.8269_-106.9763_38.8958.nc
2026-03-05 09:34:09,138 - root - INFO: POLARIS: downloading theta_r/mean over 2 tile(s)
2026-03-05 09:34:09,138 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/0_5/lat3839_lon-108-107.tif
2026-03-05 09:34:10,964 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/0_5/lat3839_lon-107-106.tif
2026-03-05 09:34:12,708 - rasterio._env - WARNING: CPLE_IllegalArg in 38131a42-e40c-4b48-9975-7d32835cc806.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:13,019 - rasterio._env - WARNING: CPLE_IllegalArg in 1f73f0e9-45e8-46e1-914e-93a13857d5be.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:13,637 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/5_15/lat3839_lon-108-107.tif
2026-03-05 09:34:15,886 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/5_15/lat3839_lon-107-106.tif
2026-03-05 09:34:17,901 - rasterio._env - WARNING: CPLE_IllegalArg in c668f2cc-fcff-415e-8362-ae28f4ea86ce.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:18,216 - rasterio._env - WARNING: CPLE_IllegalArg in b79519c5-a438-42e4-9901-2b23defa314a.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:18,820 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/15_30/lat3839_lon-108-107.tif
2026-03-05 09:34:23,409 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/15_30/lat3839_lon-107-106.tif
2026-03-05 09:34:25,338 - rasterio._env - WARNING: CPLE_IllegalArg in 371ddb4c-69c8-4420-8e80-55f92d496f67.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:25,652 - rasterio._env - WARNING: CPLE_IllegalArg in 49dde3f5-83cb-4436-8445-91c5c68b767f.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:26,242 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/30_60/lat3839_lon-108-107.tif
2026-03-05 09:34:29,206 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/30_60/lat3839_lon-107-106.tif
2026-03-05 09:34:31,342 - rasterio._env - WARNING: CPLE_IllegalArg in ca868628-e48f-4a77-baaa-fd96e5728b90.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:31,653 - rasterio._env - WARNING: CPLE_IllegalArg in 524b5358-40e4-4cb2-965d-e6f775615f43.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:32,245 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/60_100/lat3839_lon-108-107.tif
2026-03-05 09:34:34,725 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/60_100/lat3839_lon-107-106.tif
2026-03-05 09:34:36,540 - rasterio._env - WARNING: CPLE_IllegalArg in f76469d4-bd72-4415-af01-54599634ad10.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:36,851 - rasterio._env - WARNING: CPLE_IllegalArg in afbdc086-202d-4745-99e6-adefa253c7db.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:37,449 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/100_200/lat3839_lon-108-107.tif
2026-03-05 09:34:39,827 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/theta_r/mean/100_200/lat3839_lon-107-106.tif
2026-03-05 09:34:42,558 - rasterio._env - WARNING: CPLE_IllegalArg in fe2bda71-ab7b-4c65-b5f1-e650bf6abc74.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:42,870 - rasterio._env - WARNING: CPLE_IllegalArg in 1794ee26-eee7-4c49-aa6e-de620977cb6a.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:43,470 - root - INFO: Written to: /Users/uec/code/watershed_workflow/data/soil_structure/POLARIS/POLARIS_theta_r_mean_-107.1073_38.8269_-106.9763_38.8958.nc
2026-03-05 09:34:43,479 - root - INFO: POLARIS: downloading alpha/mean over 2 tile(s)
2026-03-05 09:34:43,479 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/0_5/lat3839_lon-108-107.tif
2026-03-05 09:34:45,786 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/0_5/lat3839_lon-107-106.tif
2026-03-05 09:34:49,257 - rasterio._env - WARNING: CPLE_IllegalArg in f06faebc-a20b-40c6-b394-fdd4a2d15e6e.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:49,569 - rasterio._env - WARNING: CPLE_IllegalArg in aafcc737-eef1-4c3a-a049-711abc9d5f8a.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:50,189 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/5_15/lat3839_lon-108-107.tif
2026-03-05 09:34:52,715 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/5_15/lat3839_lon-107-106.tif
2026-03-05 09:34:54,483 - rasterio._env - WARNING: CPLE_IllegalArg in 84338732-3bcd-4100-b464-39cf77770ad9.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:54,799 - rasterio._env - WARNING: CPLE_IllegalArg in 910d016e-8097-4a84-a4a6-ba2eb2e1b58f.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:55,402 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/15_30/lat3839_lon-108-107.tif
2026-03-05 09:34:57,287 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/15_30/lat3839_lon-107-106.tif
2026-03-05 09:34:59,159 - rasterio._env - WARNING: CPLE_IllegalArg in 719a31d8-022a-48bd-8a60-9196c99631e0.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:34:59,474 - rasterio._env - WARNING: CPLE_IllegalArg in dc170bb9-1203-49f8-83d4-0caf10437129.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:00,074 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/30_60/lat3839_lon-108-107.tif
2026-03-05 09:35:02,879 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/30_60/lat3839_lon-107-106.tif
2026-03-05 09:35:05,126 - rasterio._env - WARNING: CPLE_IllegalArg in aef776c4-317f-456c-a355-144f0697f242.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:05,436 - rasterio._env - WARNING: CPLE_IllegalArg in 73a964dd-64e7-4973-9ebf-2c37c5f84c82.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:06,037 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/60_100/lat3839_lon-108-107.tif
2026-03-05 09:35:09,203 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/60_100/lat3839_lon-107-106.tif
2026-03-05 09:35:11,813 - rasterio._env - WARNING: CPLE_IllegalArg in b2f80c00-100e-451b-b8d7-ec71885f488a.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:12,125 - rasterio._env - WARNING: CPLE_IllegalArg in e3c3d646-be7d-408c-acff-f1aec6fc4827.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:12,725 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/100_200/lat3839_lon-108-107.tif
2026-03-05 09:35:14,499 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/alpha/mean/100_200/lat3839_lon-107-106.tif
2026-03-05 09:35:16,874 - rasterio._env - WARNING: CPLE_IllegalArg in 9dc3b0d4-ceba-4626-aff8-f4198c2d0501.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:17,185 - rasterio._env - WARNING: CPLE_IllegalArg in 8ec044b7-76cc-4ec6-9888-2e5a18a63929.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:17,781 - root - INFO: Written to: /Users/uec/code/watershed_workflow/data/soil_structure/POLARIS/POLARIS_alpha_mean_-107.1073_38.8269_-106.9763_38.8958.nc
2026-03-05 09:35:17,791 - root - INFO: POLARIS: downloading n/mean over 2 tile(s)
2026-03-05 09:35:17,791 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/0_5/lat3839_lon-108-107.tif
2026-03-05 09:35:20,596 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/0_5/lat3839_lon-107-106.tif
2026-03-05 09:35:22,869 - rasterio._env - WARNING: CPLE_IllegalArg in 176313d9-b6ea-4847-8ec9-055834c7260b.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:23,199 - rasterio._env - WARNING: CPLE_IllegalArg in 47772efe-af69-4d9b-bd2e-b000c56244f0.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:23,837 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/5_15/lat3839_lon-108-107.tif
2026-03-05 09:35:25,401 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/5_15/lat3839_lon-107-106.tif
2026-03-05 09:35:26,942 - rasterio._env - WARNING: CPLE_IllegalArg in 45a68171-6807-4869-923e-8f51c15eb9fb.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:27,278 - rasterio._env - WARNING: CPLE_IllegalArg in 783ea2e8-a9f4-4b40-86ec-2c00e5636b53.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:27,909 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/15_30/lat3839_lon-108-107.tif
2026-03-05 09:35:31,167 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/15_30/lat3839_lon-107-106.tif
2026-03-05 09:35:32,703 - rasterio._env - WARNING: CPLE_IllegalArg in 534ad0e4-eedb-4d34-9042-b3b1c9226312.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:33,041 - rasterio._env - WARNING: CPLE_IllegalArg in f4906d35-cef6-4d08-9d0d-f2246af6b7d2.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:33,672 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/30_60/lat3839_lon-108-107.tif
2026-03-05 09:35:43,592 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/30_60/lat3839_lon-107-106.tif
2026-03-05 09:35:45,625 - rasterio._env - WARNING: CPLE_IllegalArg in ae2050c6-1e8e-4f9d-b330-c288a15b2bea.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:45,957 - rasterio._env - WARNING: CPLE_IllegalArg in bab1d316-bb53-43f1-81e6-ddd103a2cc54.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:46,579 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/60_100/lat3839_lon-108-107.tif
2026-03-05 09:35:48,479 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/60_100/lat3839_lon-107-106.tif
2026-03-05 09:35:50,590 - rasterio._env - WARNING: CPLE_IllegalArg in a2863603-ccd4-404f-a244-84c4144be3d1.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:50,924 - rasterio._env - WARNING: CPLE_IllegalArg in 3c60330f-ae8c-459a-87fd-4fdf5d4e815c.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:51,550 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/100_200/lat3839_lon-108-107.tif
2026-03-05 09:35:53,342 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/n/mean/100_200/lat3839_lon-107-106.tif
2026-03-05 09:35:56,553 - rasterio._env - WARNING: CPLE_IllegalArg in 67524494-9109-4af7-9dbe-265db3f2b6c6.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:56,886 - rasterio._env - WARNING: CPLE_IllegalArg in 8b3b048c-d868-4ab7-9e4b-713237831f0e.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:35:57,513 - root - INFO: Written to: /Users/uec/code/watershed_workflow/data/soil_structure/POLARIS/POLARIS_n_mean_-107.1073_38.8269_-106.9763_38.8958.nc
2026-03-05 09:35:57,521 - root - INFO: POLARIS: downloading ksat/mean over 2 tile(s)
2026-03-05 09:35:57,521 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/0_5/lat3839_lon-108-107.tif
2026-03-05 09:36:00,799 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/0_5/lat3839_lon-107-106.tif
2026-03-05 09:36:04,073 - rasterio._env - WARNING: CPLE_IllegalArg in c11f2549-a9f2-401e-a9b0-e0c66a454225.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:04,404 - rasterio._env - WARNING: CPLE_IllegalArg in 816b5039-d7b1-43ac-b39c-e713918b2209.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:05,024 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/5_15/lat3839_lon-108-107.tif
2026-03-05 09:36:07,925 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/5_15/lat3839_lon-107-106.tif
2026-03-05 09:36:10,979 - rasterio._env - WARNING: CPLE_IllegalArg in ec8c2ff7-e5d4-4ba7-a05e-2819371980a6.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:11,274 - rasterio._env - WARNING: CPLE_IllegalArg in fc1d90fe-17cc-43c4-9291-d3977c62f1fb.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:11,845 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/15_30/lat3839_lon-108-107.tif
2026-03-05 09:36:15,095 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/15_30/lat3839_lon-107-106.tif
2026-03-05 09:36:18,527 - rasterio._env - WARNING: CPLE_IllegalArg in f41a4603-dc94-466d-b1ea-90eb6166cc48.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:18,824 - rasterio._env - WARNING: CPLE_IllegalArg in 444ae261-2bec-4df3-9abc-28113a8dda41.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:19,392 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/30_60/lat3839_lon-108-107.tif
2026-03-05 09:36:25,511 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/30_60/lat3839_lon-107-106.tif
2026-03-05 09:36:28,438 - rasterio._env - WARNING: CPLE_IllegalArg in 998a81e7-f248-4000-91ef-c27227cce58e.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:28,726 - rasterio._env - WARNING: CPLE_IllegalArg in a209c528-a3e7-4b1e-9d61-a077b297a817.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:29,302 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/60_100/lat3839_lon-108-107.tif
2026-03-05 09:36:31,129 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/60_100/lat3839_lon-107-106.tif
2026-03-05 09:36:33,843 - rasterio._env - WARNING: CPLE_IllegalArg in 88bf0498-e4cf-466d-914f-97c40a1da819.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:34,131 - rasterio._env - WARNING: CPLE_IllegalArg in 8f30d596-d4ea-43fc-9595-74a78dd3510d.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:34,693 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/100_200/lat3839_lon-108-107.tif
2026-03-05 09:36:36,801 - root - INFO: GET http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ksat/mean/100_200/lat3839_lon-107-106.tif
2026-03-05 09:36:38,573 - rasterio._env - WARNING: CPLE_IllegalArg in fc7ca4ab-1d2d-46c8-a907-dfe9d2eea86d.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:38,866 - rasterio._env - WARNING: CPLE_IllegalArg in b48a5da8-6158-403f-baa3-0c185fb9f586.tif: BLOCKXSIZE can only be used with TILED=YES
2026-03-05 09:36:39,442 - root - INFO: Written to: /Users/uec/code/watershed_workflow/data/soil_structure/POLARIS/POLARIS_ksat_mean_-107.1073_38.8269_-106.9763_38.8958.nc
<xarray.Dataset> Size: 14MB
Dimensions: (depth: 6, y: 249, x: 473)
Coordinates:
* depth (depth) float32 24B 0.025 0.1 0.225 0.45 0.8 1.5
* y (y) float64 2kB 38.9 38.9 38.9 38.89 ... 38.83 38.83 38.83
* x (x) float64 4kB -107.1 -107.1 -107.1 ... -107.0 -107.0 -107.0
spatial_ref int64 8B 0
Data variables:
theta_s (depth, y, x) float32 3MB ...
theta_r (depth, y, x) float32 3MB ...
alpha (depth, y, x) float32 3MB ...
n (depth, y, x) float32 3MB ...
ksat (depth, y, x) float32 3MB ...
Attributes:
source: http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0
variable: theta_s
stat: mean
name: POLARIS
# Plot theta_s (porosity) and van Genuchten n at shallowest and deepest layers.
# Rows: theta_s (top), n (bottom). Columns: 0-5 cm (left), 100-200 cm (right).
watershed_ll = watershed.to_crs(polaris_ds.rio.crs)
fig, axs = plt.subplots(2, 2, figsize=(14, 10))
plots = [
(0, 0, 'theta_s', 'Porosity θ_s [-]', 'Blues', 0),
(0, 1, 'theta_s', 'Porosity θ_s [-]', 'Blues', -1),
(1, 0, 'n', 'van Genuchten n [-]', 'RdYlGn', 0),
(1, 1, 'n', 'van Genuchten n [-]', 'RdYlGn', -1),
]
depth_titles = ['0–5 cm', '100–200 cm']
for row, col, var, label, cmap, z_idx in plots:
ax = axs[row, col]
polaris_ds[var].isel(depth=z_idx).plot.imshow(ax=ax, cmap=cmap, add_colorbar=True)
watershed_ll.boundary.plot(ax=ax, color='k', linewidth=1.5)
ax.set_title(f'POLARIS {label} — {depth_titles[col]}')
plt.tight_layout()
# Plot K_sat (stored as log10(cm/hr)) at shallowest and deepest layers.
import matplotlib.colors as mpc
fig, axs = plt.subplots(1, 2, figsize=(14, 5))
for ax, z_idx, title in [
(axs[0], 0, 'POLARIS K_sat [log₁₀ cm/hr] — 0–5 cm'),
(axs[1], -1, 'POLARIS K_sat [log₁₀ cm/hr] — 100–200 cm'),
]:
polaris_ds['ksat'].isel(depth=z_idx).plot.imshow(
ax=ax, cmap='viridis', add_colorbar=True,
cbar_kwargs={'label': 'log₁₀(K_sat) [log₁₀ cm/hr]'}
)
watershed_ll.boundary.plot(ax=ax, color='k', linewidth=1.5)
ax.set_title(title)
plt.tight_layout()
