This page was generated from soil.ipynb. Interactive online version:
Soil Storage Capacity#
from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pygeohydro as gh import pygeoutils as geoutils from pynhd import NLDI
Some hydrological models have a parameter related the max soil storage capacity of watersheds. A simple way of estimating this parameter is to use porosity and thickness properties of a watershed soil. In this tutorial, we retrieve these data from two soil datasets that PyGeoHydro supports. PyGeoHydro has two functions for getting soil data:
soil_gnatsgo: You can access all soil properties that are available on gNATSGO dataset
soil_properties: Porosity, available water capacity, and field capacity properties are available through this function.
Let’s start by getting the basin geometry of the USGS station 11092450 using PyNHD:
basin = NLDI().get_basins("11092450") geom = basin.geometry
Next, we get the soil porosity. Note that
soil_properties return the requested property for the whole US. But, thanks to
xarray, the dataset is loaded lazily, so we don’t have to worry about memory usage. This function downloads the requested dataset(s) each of which is around 300 MB, so the first time that you run this function it can take a few minutes, depending on your internet speed. However, after it gets downloaded once, next calls to this function uses the downloaded
file, so it will be quick. By default, downloaded files are saved to
./cache folder, but you can specify a custom directory to save the
If you run this function without any arguments, it downloads all three soil properties and returns them as a single dataset. You can pass one or more of soil properties:
por for porosity,
awc for available water capacity, and
fc for field capacity.
soil = gh.soil_properties("por")
Now, we can mask this large dataset to our basin geometry using PyGeoUtils.
soil = geoutils.xarray_geomask(soil, geom, basin.crs) soil = soil.where(soil.porosity > soil.porosity.rio.nodata) soil["porosity"] = soil.porosity.rio.write_nodata(np.nan) _ = soil.porosity.plot()
Now that we have the porosity, we need to get soil thickness. We use Microsoft’s Planetary Computer service using
soil_gnatsgo function. This function has two additional dependencies that need to be installed:
planetary-computer. They can be installed using
mamba as follows:
mamba -c conda-forge install pystac-client planetary-computer
tk0_999a later to get the total soil thickness. You can check out the full description the layer and other available layers from the dataset’s webpage. Also, we convert the thickness unit from cm to mm.
thickness = gh.soil_gnatsgo("tk0_999a", geom, basin.crs).tk0_999a thickness = thickness.where(thickness < 2e6, drop=False) * 10 thickness = thickness.rio.write_nodata(np.nan) thickness.attrs["units"] = "mm" thickness.attrs["long_name"] = "soil thickness" thickness.name = "thickness"
soil.rio.resolution() we can see that the resolution of the thickness dataset is 10 m whereas that of the porosity is around 90 m. So, we should resample the dataset with finer resolution, thickness, to the coarser one. For this purpose we use
rioxarray.reproject_match function. For resampling, we use the weighted averaging function of
rasterio. Therefore, we need to 5 as the resampling method. The full list of supported resampling operations
can be found on
soil["thickness"] = thickness.rio.reproject_match(soil, resampling=5) _ = soil["thickness"].plot()
Considering that porosity is in mm/m and thickness is in mm, we can compute the storage capacity as follows:
soil["storage"] = soil.porosity * soil.thickness * 1e-3 soil["storage"] = soil["storage"].rio.write_nodata(np.nan) soil["storage"].attrs["units"] = "mm" soil["storage"].attrs["long_name"] = "soil storage capacity" soil["storage"].plot() plt.savefig(Path("_static", "soil_storage.png"), dpi=150)