import logging
from pathlib import Path
from typing import TYPE_CHECKING, Union
import numpy as np
import pandas as pd
import xarray as xr
import xugrid as xu
from hydromt import hydromt_step
from hydromt.model.components import ModelComponent
from hydromt_sfincs import DATADIR, workflows
from hydromt_sfincs.utils import fill_nan_in_mask
from hydromt_sfincs.components.infiltration_common import (
ALL_VARS,
DEFAULT_INFILTRATIONFILE,
VARIABLES,
clear_data,
configure,
configured_flavor,
flavor_variables,
get_attrs,
reset_config,
_require_lulc_modifiers,
)
from hydromt_sfincs.components.quadtree import SfincsQuadtreeMixin
if TYPE_CHECKING:
from hydromt_sfincs import SfincsModel
logger = logging.getLogger(f"hydromt.{__name__}")
[docs]
class SfincsQuadtreeInfiltration(SfincsQuadtreeMixin, ModelComponent):
"""SFINCS infiltration component for quadtree grids.
Unsuffixed ``create_*`` methods estimate parameters from HSG, optional
Ksat, and optional land-use modifiers. Methods with a ``_from_maps`` suffix
use final SFINCS parameter maps directly.
"""
[docs]
def __init__(self, model: "SfincsModel"):
super().__init__(model=model)
@property
def data(self) -> xu.UgridDataset:
"""Return only infiltration layers, retaining the parent mesh."""
data = self.model.quadtree_grid.data
return data[[name for name in self.grid_variables if name in data]]
@property
def mask(self) -> xu.UgridDataArray:
return self.model.quadtree_grid.mask
grid_variables = ALL_VARS
def _set_layers(
self,
layers: dict[str, Union[np.ndarray, xr.DataArray, xu.UgridDataArray]],
flavor: str,
) -> None:
self.clear()
grid = self.mask.grid
for name, layer in layers.items():
values = (
layer.values
if isinstance(layer, (xr.DataArray, xu.UgridDataArray))
else layer
)
da = xr.DataArray(
fill_nan_in_mask(
values, self.mask.values, name, VARIABLES[name].fill_value
),
dims=[grid.face_dimension],
)
uda = xu.UgridDataArray(da, grid)
uda = uda.astype(np.float32)
uda.name = name
uda.attrs.update(get_attrs(name))
self.model.quadtree_grid.set(uda, name=name)
configure(self.model.config, flavor=flavor, grid_type="quadtree")
[docs]
def read(self) -> None:
"""Read the infiltration layers from their own file."""
# TODO discuss whether we want to clear infiltration data first
# the grid defines the mesh these layers live on, so load it first
if self.model.quadtree_grid._data is None:
self.model.quadtree_grid.read(read_components=False)
flavor = configured_flavor(self.model.config, grid_type="quadtree")
if flavor is None or flavor == "con":
return
filename = self.model.config.get("inffile", abs_path=True)
if filename is None or not filename.is_file():
return
self.model.quadtree_grid._read_layers(filename)
[docs]
def write(self) -> None:
"""Write the infiltration layers to their own file."""
flavor = configured_flavor(self.model.config, grid_type="quadtree")
if flavor is None or flavor == "con":
return
self.model.quadtree_grid._write_layers(
flavor_variables(flavor), "inffile", DEFAULT_INFILTRATIONFILE
)
[docs]
@hydromt_step
def create_constant(
self,
qinf: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
lulc: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
reclass_table: Union[str, Path, pd.DataFrame, None] = None,
reproj_method: str = "average",
nrmax: int = 2000,
) -> None:
"""Create spatially varying constant infiltration rate.
Adds model layers to SfincsModel.quadtree_grid.data:
* **qinf** map: constant infiltration rate [mm/hr]
Parameters
----------
qinf : str, Path, or RasterDataset
Spatially varying infiltration rates [mm/hr]. If this is an xr.Dataset,
it must contain a ``qinf`` variable.
lulc : str, Path, or RasterDataset
Land-use dataset. Must be combined with ``reclass_table``.
reclass_table : str, Path, or pd.DataFrame
Reclassification table with a ``qinf`` column to convert land-use classes
to infiltration rates [mm/hr].
reproj_method : str, optional
Resampling method for reprojecting source data to quadtree blocks.
nrmax : int, optional
Maximum number of cells per quadtree block.
"""
logger.info("Creating constant spatially varying infiltration rate.")
# get infiltration data
if qinf is not None:
da_qinf = self.data_catalog.get_rasterdataset(
qinf,
bbox=self.model.bbox,
variables=["qinf"],
buffer=10,
)
# TODO check if this one is really necessary?
# elif lulc is not None and ksat is not None:
# da_ksat = self.data_catalog.get_rasterdataset(
# ksat, bbox=self.model.bbox, buffer=10
# )
# da_lulc = self.data_catalog.get_rasterdataset(
# lulc, bbox=self.model.bbox, buffer=10
# )
# if lulc_modifier_table is None:
# lulc_modifier_table = (
# Path(DATADIR)
# / "infiltration"
# / "nlcd_infiltration_modifiers.csv"
# )
# if isinstance(lulc_modifier_table, pd.DataFrame):
# df_modifiers = lulc_modifier_table.copy()
# else:
# df_modifiers = self.data_catalog.get_dataframe(
# lulc_modifier_table,
# source_kwargs={
# "driver": {"name": "pandas", "options": {"index_col": 0}}
# },
# )
# da_qinf = workflows.constant_infiltration_from_ksat_lulc(
# da_ksat,
# da_lulc,
# df_modifiers,
# da_mask=self.mask,
# factor_ksat=factor_ksat,
# )
elif lulc is not None:
# landuse/landcover should always be combined with mapping
if reclass_table is None:
raise IOError(
f"Infiltration mapping file should be provided for {lulc}"
)
da_lulc = self.data_catalog.get_rasterdataset(
lulc,
bbox=self.model.bbox,
buffer=10,
variables=["lulc"],
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
variables=["qinf"],
source_kwargs={
"driver": {"name": "pandas", "options": {"index_col": 0}}
},
)
# reclassify
da_qinf = da_lulc.raster.reclassify(df_map)["qinf"]
else:
raise ValueError(
"Either qinf or lulc must be provided when setting up constant infiltration."
)
# set nodata to nan before reprojecting/interpolating
da_qinf = da_qinf.raster.mask_nodata()
n_cells = self.mask.grid.n_face
qinf = np.full(n_cells, np.nan)
def compute_constant_infiltration(da_like, ilev=None):
# reproject infiltration data to model grid
return da_qinf.raster.reproject_like(da_like, method=reproj_method)
self.compute_quadtree(
compute_constant_infiltration,
qinf,
nrmax=nrmax,
)
self._set_layers({"qinf": qinf}, flavor="c2d")
[docs]
@hydromt_step
def create_runoff_coefficient(
self,
runoff: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
lulc: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
reclass_table: Union[str, Path, pd.DataFrame, None] = None,
reproj_method: str = "average",
nrmax: int = 2000,
) -> None:
"""Create a spatially varying constant run-off coefficient.
A fixed percentage of the rainfall runs off; the remainder infiltrates
(SFINCS ``inftype = r2d``). A value of 0 means all rainfall infiltrates,
100 means all rainfall runs off.
Adds model layers to SfincsModel.quadtree_grid.data:
* **runoff** map: run-off coefficient [%]
Parameters
----------
runoff : str, Path, or RasterDataset
Spatially varying run-off coefficients [%], between 0 and 100. If this
is an xr.Dataset, it must contain a ``runoff`` variable.
lulc : str, Path, or RasterDataset
Land-use dataset. Must be combined with ``reclass_table``.
reclass_table : str, Path, or pd.DataFrame
Reclassification table with a ``runoff`` column to convert land-use
classes to run-off coefficients [%].
reproj_method : str, optional
Resampling method for reprojecting source data to quadtree blocks.
nrmax : int, optional
Maximum number of cells per quadtree block.
"""
logger.info("Creating spatially varying run-off coefficient.")
if runoff is not None:
da_runoff = self.data_catalog.get_rasterdataset(
runoff,
bbox=self.model.bbox,
variables=["runoff"],
buffer=10,
)
elif lulc is not None:
if reclass_table is None:
raise IOError(
f"Infiltration mapping file should be provided for {lulc}"
)
da_lulc = self.data_catalog.get_rasterdataset(
lulc,
bbox=self.model.bbox,
buffer=10,
variables=["lulc"],
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
variables=["runoff"],
source_kwargs={
"driver": {"name": "pandas", "options": {"index_col": 0}}
},
)
da_runoff = da_lulc.raster.reclassify(df_map)["runoff"]
else:
raise ValueError("Either runoff or lulc must be provided.")
da_runoff = da_runoff.raster.mask_nodata()
runoff = np.full(self.mask.grid.n_face, np.nan)
def compute_runoff(da_like, ilev=None):
return da_runoff.raster.reproject_like(da_like, method=reproj_method)
self.compute_quadtree(compute_runoff, runoff, nrmax=nrmax)
if np.nanmin(runoff) < 0.0 or np.nanmax(runoff) > 100.0:
# apply cap when run-off coefficients are out of the valid range
runoff = np.clip(runoff, 0.0, 100.0)
logger.warning(
"Run-off coefficients are out of the valid range (0-100). Clipping to valid range."
)
self._set_layers({"runoff": runoff}, flavor="r2d")
[docs]
@hydromt_step
def create_cn(
self,
cn: Union[str, Path, xr.DataArray, xr.Dataset],
antecedent_moisture: Union[str, None] = "avg",
reproj_method: str = "med",
nrmax: int = 2000,
) -> None:
"""Create Curve Number infiltration without recovery for quadtree grids.
Adds model layers:
* **scs** map: potential maximum soil moisture retention [inch]
Parameters
----------
cn : str, Path, or RasterDataset
Curve number data. Dataset inputs must contain a ``cn`` variable, or
``cn_<antecedent_moisture>`` when ``antecedent_moisture`` is set.
antecedent_moisture : {'dry', 'avg', 'wet'}, optional
Antecedent runoff condition used to select the source variable. Set to
None when the input already holds adjusted curve numbers.
By default 'avg'.
reproj_method : str, optional
Resampling method for reprojecting curve number data to quadtree blocks.
By default 'med'.
nrmax : int, optional
Maximum number of cells per quadtree block, by default 2000.
"""
# Add logger info
logger.info(
f"Creating curve number values for SFINCS quadtree grid with antecedent moisture condition: {antecedent_moisture}."
)
# get data
da_org = self.model.data_catalog.get_rasterdataset(
cn, bbox=self.model.bbox, buffer=10
)
v = "cn"
if antecedent_moisture:
v = f"cn_{antecedent_moisture}"
if isinstance(da_org, xr.Dataset) and v in da_org.data_vars:
da_org = da_org[v]
elif isinstance(da_org, xr.Dataset):
raise ValueError(f"Could not find variable {v} in {cn}")
n_cells = self.mask.grid.n_face
scs = np.full(n_cells, np.nan)
def compute_cn_infiltration(da_like, ilev=None):
da_cn = da_org.raster.reproject_like(da_like, method=reproj_method)
da_scs = workflows.cn_to_s(da_cn).round(3)
return da_scs
self.compute_quadtree(
compute_cn_infiltration,
scs,
nrmax=nrmax,
)
self._set_layers({"scs": scs}, flavor="cna")
[docs]
@hydromt_step
def create_cn_from_landuse_hsg(
self,
lulc: Union[str, Path, xr.DataArray, xr.Dataset],
hsg: Union[str, Path, xr.DataArray, xr.Dataset],
reclass_table: Union[str, Path, pd.DataFrame],
antecedent_moisture: str = "avg",
reproj_method: str = "average",
nrmax: int = 2000,
) -> None:
"""Create Curve Number infiltration from land use and HSG.
S is computed per land-use pixel and area-averaged onto the quadtree
faces, as in :py:meth:`create_cn_with_recovery`. Open water (CN 0 in the
table) counts as S = 0, so a face that is 60% water gets 40% of the land
S. Unmapped pixels are skipped.
Adds model layers:
* **scs** map: potential maximum soil moisture retention [inch]
Parameters
----------
lulc : str, Path, or RasterDataset
Name of gridded land use map.
hsg : str, Path, or RasterDataset
Name of gridded hydrologic soil group map.
reclass_table : str, Path, or DataFrame
Reclassification table mapping land use and hydrologic soil groups to curve numbers.
antecedent_moisture : {'dry', 'avg', 'wet'}, optional
Antecedent runoff conditions.
By default `avg`
reproj_method : str, optional
Resampling method for S to quadtree blocks, by default 'average'
(area-weighted).
nrmax : int, optional
Maximum number of cells per quadtree block, by default 2000.
"""
da_lulc = self.data_catalog.get_rasterdataset(
lulc, bbox=self.model.bbox, buffer=10, variables=["lulc"]
)
da_hsg = self.data_catalog.get_rasterdataset(
hsg, bbox=self.model.bbox, buffer=10, variables=["hsg"]
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
source_kwargs={"driver": {"name": "pandas", "options": {"index_col": 0}}},
)
da_cn = workflows.curve_number_from_landuse_hsg(da_lulc, da_hsg, df_map)
da_cn = workflows.adjust_curve_number(
da_cn,
antecedent_moisture=antecedent_moisture,
)
# S per pixel; missing pixels stay NaN nodata so the average skips them
da_s = workflows.cn_to_s(da_cn, fill_cn=None)
scs = np.full(self.mask.grid.n_face, np.nan)
def compute_cn_infiltration(da_like, ilev=None):
return da_s.raster.reproject_like(da_like, method=reproj_method).round(3)
self.compute_quadtree(compute_cn_infiltration, scs, nrmax=nrmax)
self._set_layers({"scs": scs}, flavor="cna")
[docs]
@hydromt_step
def create_cn_with_recovery(
self,
lulc: Union[str, Path, xr.DataArray, xr.Dataset],
hsg: Union[str, Path, xr.DataArray, xr.Dataset],
ksat: Union[str, Path, xr.DataArray, xr.Dataset],
reclass_table: Union[str, Path, pd.DataFrame],
effective: float,
block_size: int = 2000,
) -> None:
"""Create Curve Number infiltration with recovery for quadtree grids.
Adds model layers:
* **smax** map: maximum soil moisture retention [m]
* **seff** map: effective soil moisture retention [m]
* **ks** map: saturated hydraulic conductivity [mm/hr]
Input data are read by the component; each quadtree block is passed to
:py:func:`hydromt_sfincs.workflows.curve_number_with_recovery`.
Parameters
----------
lulc : str, Path, or RasterDataset
Landuse/landcover data set.
hsg : str, Path, or RasterDataset
Hydrologic soil group map in integers.
ksat : str, Path, or RasterDataset
Saturated hydraulic conductivity [mm/hr].
reclass_table : str, Path, or DataFrame
Reclassification table relating land cover and soil type to curve numbers.
effective : float
Fraction of ``smax`` that is effective soil retention, e.g. 0.50 for 50%.
block_size : int, optional
Maximum number of cells per quadtree block, by default 2000.
"""
da_lulc = self.data_catalog.get_rasterdataset(
lulc, bbox=self.model.bbox, buffer=10
)
da_hsg = self.data_catalog.get_rasterdataset(
hsg, bbox=self.model.bbox, buffer=10
)
da_ksat = self.data_catalog.get_rasterdataset(
ksat, bbox=self.model.bbox, buffer=10
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
source_kwargs={"driver": {"name": "pandas", "options": {"index_col": 0}}},
)
outputs = {
name: np.full(self.mask.grid.n_face, np.nan)
for name in ("smax", "seff", "ks")
}
def compute_cn_recovery_block(da_like, ilev=None):
ds = workflows.curve_number_with_recovery(
da_lulc,
da_hsg,
da_ksat,
df_map,
effective=effective,
da_mask=da_like,
)
return tuple(ds[name] for name in outputs)
self.compute_quadtree(compute_cn_recovery_block, outputs, nrmax=block_size)
self._set_layers(outputs, flavor="cnb")
[docs]
@hydromt_step
def create_green_ampt(
self,
hsg: Union[str, Path, xr.DataArray, xr.Dataset],
ksat: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
lulc: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
reclass_table: Union[str, Path, pd.DataFrame, None] = None,
lulc_modifier_table: Union[str, Path, pd.DataFrame, None] = None,
dual_hsg: str = "drained",
reproj_method: str = "average",
) -> None:
"""Estimate Green-Ampt infiltration from HSG and optional landuse.
Adds model layers:
* **psi** map: wetting front suction head [mm]
* **sigma** map: soil moisture deficit [-]
* **ks** map: saturated hydraulic conductivity [mm/hr]
Parameters
----------
hsg : str, Path, or RasterDataset
Hydrologic soil group map. By default, values are reclassified with
the bundled ``hsg_green_ampt.csv`` table.
ksat : str, Path, or RasterDataset, optional
Saturated hydraulic conductivity map, used to derive ``ks``
(``ks = ksat``, in mm/hr). Required unless the reclass table
already provides a ``ks`` column. The bundled ``hsg_green_ampt.csv``
has only ``psi``/``sigma``, so ksat is required when using it;
omitting it raises a ValueError.
lulc : str, Path, or RasterDataset, optional
Land-use map used to apply infiltration modifiers. Its classes must
match the index of ``lulc_modifier_table``.
reclass_table : str, Path, or DataFrame, optional
Table mapping HSG classes to Green-Ampt parameters.
lulc_modifier_table : str, Path, or DataFrame, optional
Required with ``lulc``. Table of modifier factors keyed by the
land-cover class codes in the supplied dataset.
dual_hsg : {'native', 'drained'}, optional
How to handle dual HSG classes, by default 'drained'. Native dual
classes must be present in the reclassification table.
reproj_method : str, optional
Resampling method for reprojecting final parameter maps to quadtree
blocks. By default 'average'.
"""
_require_lulc_modifiers(lulc, lulc_modifier_table)
if reclass_table is None:
reclass_table = Path(DATADIR) / "infiltration" / "hsg_green_ampt.csv"
da_hsg = self.data_catalog.get_rasterdataset(
hsg, bbox=self.model.bbox, buffer=10, variables=["hsg"]
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
source_kwargs={"driver": {"name": "pandas", "options": {"index_col": 0}}},
)
da_ksat = None
if ksat is not None:
da_ksat = self.data_catalog.get_rasterdataset(
ksat, bbox=self.model.bbox, buffer=10
)
da_ksat = workflows.mask_invalid_ksat(da_ksat).raster.reproject_like(
da_hsg, method="average"
)
if lulc is not None:
da_lulc = self.data_catalog.get_rasterdataset(
lulc, bbox=self.model.bbox, buffer=10, variables=["lulc"]
)
da_lulc = da_lulc.raster.reproject_like(da_hsg, method="nearest")
df_modifiers = self.data_catalog.get_dataframe(
lulc_modifier_table,
source_kwargs={
"driver": {"name": "pandas", "options": {"index_col": 0}}
},
)
ds = workflows.green_ampt_from_hsg_landuse(
da_hsg,
da_lulc,
df_map,
df_modifiers,
da_ksat=da_ksat,
dual_hsg=dual_hsg,
)
else:
ds = workflows.green_ampt_from_hsg(
da_hsg,
df_map,
da_ksat=da_ksat,
dual_hsg=dual_hsg,
)
outputs = {
name: np.full(self.mask.grid.n_face, np.nan)
for name in ("psi", "sigma", "ks")
}
def compute_block(da_like, ilev=None):
return tuple(
ds[name].raster.reproject_like(da_like, method=reproj_method)
for name in outputs
)
self.compute_quadtree(compute_block, outputs)
self._set_layers(outputs, flavor="gai")
[docs]
@hydromt_step
def create_green_ampt_from_maps(
self,
psi: Union[str, Path, xr.DataArray, xr.Dataset],
sigma: Union[str, Path, xr.DataArray, xr.Dataset],
ks: Union[str, Path, xr.DataArray, xr.Dataset],
reproj_method: str = "average",
) -> None:
"""Create Green-Ampt infiltration from final parameter maps.
Adds model layers:
* **psi** map: wetting front suction head [mm]
* **sigma** map: soil moisture deficit [-]
* **ks** map: saturated hydraulic conductivity [mm/hr]
Parameters
----------
psi, sigma, ks : str, Path, or RasterDataset
Data with final Green-Ampt parameters. Dataset inputs must contain
variables named ``psi``, ``sigma``, and ``ks`` respectively.
reproj_method : str, optional
Resampling method for reprojecting raster inputs to quadtree blocks.
"""
names = ("psi", "sigma", "ks")
layers = {}
raster_sources = {}
for name, source in zip(names, (psi, sigma, ks)):
da = self.data_catalog.get_rasterdataset(
source, bbox=self.model.bbox, buffer=10, variables=[name]
)
raster_sources[name] = da.raster.mask_nodata()
if raster_sources:
outputs = {
name: np.full(self.mask.grid.n_face, np.nan) for name in raster_sources
}
def compute_block(da_like, ilev=None):
return tuple(
raster_sources[name].raster.reproject_like(
da_like, method=reproj_method
)
for name in outputs
)
self.compute_quadtree(compute_block, outputs)
layers.update(outputs)
self._set_layers(layers, flavor="gai")
[docs]
@hydromt_step
def create_horton(
self,
hsg: Union[str, Path, xr.DataArray, xr.Dataset],
ksat: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
lulc: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
reclass_table: Union[str, Path, pd.DataFrame, None] = None,
lulc_modifier_table: Union[str, Path, pd.DataFrame, None] = None,
dual_hsg: str = "drained",
reproj_method: str = "average",
) -> None:
"""Estimate Horton infiltration from HSG and optional landuse.
Adds model layers:
* **f0** map: initial infiltration capacity [mm/hr]
* **fc** map: asymptotic infiltration capacity [mm/hr]
* **kd** map: Horton decay coefficient [hr-1]
Parameters
----------
hsg : str, Path, or RasterDataset
Hydrologic soil group map. By default, values are reclassified with
the bundled ``hsg_horton.csv`` table.
ksat : str, Path, or RasterDataset, optional
Saturated hydraulic conductivity map, used to derive ``fc``
(``fc = ksat * fc_scale``, ksat in mm/hr). Required unless the reclass
table already provides a literal ``fc`` column. The bundled
``hsg_horton.csv`` has only ``fc_scale``/``f0_scale``/``kd`` (and
``fc_scale`` is applied only when ksat is given), so ksat is required
when using it; omitting it raises a ValueError.
lulc : str, Path, or RasterDataset, optional
Land-use map used to apply infiltration modifiers. Its classes must
match the index of ``lulc_modifier_table``.
reclass_table : str, Path, or DataFrame, optional
Table mapping HSG classes to Horton parameters.
lulc_modifier_table : str, Path, or DataFrame, optional
Required with ``lulc``. Table of modifier factors keyed by the
land-cover class codes in the supplied dataset.
dual_hsg : {'native', 'drained'}, optional
How to handle dual HSG classes, by default 'drained'. Native dual
classes must be present in the reclassification table.
reproj_method : str, optional
Resampling method for reprojecting final parameter maps to quadtree
blocks. By default 'average'.
"""
_require_lulc_modifiers(lulc, lulc_modifier_table)
if reclass_table is None:
reclass_table = Path(DATADIR) / "infiltration" / "hsg_horton.csv"
da_hsg = self.data_catalog.get_rasterdataset(
hsg, bbox=self.model.bbox, buffer=10, variables=["hsg"]
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
source_kwargs={"driver": {"name": "pandas", "options": {"index_col": 0}}},
)
da_ksat = None
if ksat is not None:
da_ksat = self.data_catalog.get_rasterdataset(
ksat, bbox=self.model.bbox, buffer=10
)
da_ksat = workflows.mask_invalid_ksat(da_ksat).raster.reproject_like(
da_hsg, method="average"
)
if lulc is not None:
da_lulc = self.data_catalog.get_rasterdataset(
lulc, bbox=self.model.bbox, buffer=10, variables=["lulc"]
)
da_lulc = da_lulc.raster.reproject_like(da_hsg, method="nearest")
df_modifiers = self.data_catalog.get_dataframe(
lulc_modifier_table,
source_kwargs={
"driver": {"name": "pandas", "options": {"index_col": 0}}
},
)
ds = workflows.horton_from_hsg_landuse(
da_hsg,
da_lulc,
df_map,
df_modifiers,
da_ksat=da_ksat,
dual_hsg=dual_hsg,
)
else:
ds = workflows.horton_from_hsg(
da_hsg,
df_map,
da_ksat=da_ksat,
dual_hsg=dual_hsg,
)
outputs = {
name: np.full(self.mask.grid.n_face, np.nan) for name in ("f0", "fc", "kd")
}
def compute_block(da_like, ilev=None):
return tuple(
ds[name].raster.reproject_like(da_like, method=reproj_method)
for name in outputs
)
self.compute_quadtree(compute_block, outputs)
self._set_layers(outputs, flavor="hor")
[docs]
@hydromt_step
def create_horton_from_maps(
self,
f0: Union[str, Path, xr.DataArray, xr.Dataset],
fc: Union[str, Path, xr.DataArray, xr.Dataset],
kd: Union[str, Path, xr.DataArray, xr.Dataset],
reproj_method: str = "average",
) -> None:
"""Create Horton infiltration from final parameter maps.
Adds model layers:
* **f0** map: initial infiltration capacity [mm/hr]
* **fc** map: asymptotic infiltration capacity [mm/hr]
* **kd** map: Horton decay coefficient [hr-1]
Parameters
----------
f0, fc, kd : str, Path, RasterDataset, or UgridDataArray
Data with final Horton parameters. Dataset inputs must contain
variables named ``f0``, ``fc``, and ``kd`` respectively.
reproj_method : str, optional
Resampling method for reprojecting raster inputs to quadtree blocks.
"""
names = ("f0", "fc", "kd")
layers = {}
raster_sources = {}
for name, source in zip(names, (f0, fc, kd)):
da = self.data_catalog.get_rasterdataset(
source,
bbox=self.model.bbox,
buffer=10,
variables=[name],
)
raster_sources[name] = da.raster.mask_nodata()
if raster_sources:
outputs = {
name: np.full(self.mask.grid.n_face, np.nan) for name in raster_sources
}
def compute_block(da_like, ilev=None):
return tuple(
raster_sources[name].raster.reproject_like(
da_like, method=reproj_method
)
for name in outputs
)
self.compute_quadtree(compute_block, outputs)
layers.update(outputs)
self._set_layers(layers, flavor="hor")
[docs]
@hydromt_step
def create_bucket(
self,
hsg: Union[str, Path, xr.DataArray, xr.Dataset],
ksat: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
lulc: Union[str, Path, xr.DataArray, xr.Dataset, None] = None,
reclass_table: Union[str, Path, pd.DataFrame, None] = None,
lulc_modifier_table: Union[str, Path, pd.DataFrame, None] = None,
dual_hsg: str = "drained",
bucket_loss: Union[float, None] = None,
reproj_method: str = "average",
) -> None:
"""Estimate bucket infiltration from HSG and optional landuse.
Adds model layers:
* **bucket_smax** map: bucket maximum storage [mm]
* **bucket_k** map: bucket drainage coefficient [hr-1]
* **bucket_loss** map: bucket loss fraction [-]
Parameters
----------
hsg : str, Path, or RasterDataset
Hydrologic soil group map. By default, values are reclassified with
the bundled ``hsg_bucket.csv`` table.
ksat : str, Path, or RasterDataset, optional
Saturated hydraulic conductivity map, used to derive ``bucket_k``
(via a residence time ``bucket_smax / ksat``, ksat in mm/hr).
Required unless the reclass table provides a ``bucket_k`` or
``residence_time_hr`` column. The bundled ``hsg_bucket.csv`` provides
neither, so ksat is required when using it; omitting it raises a
ValueError. (``bucket_smax`` itself comes from the table and does not
need ksat.)
lulc : str, Path, or RasterDataset, optional
Land-use map used to apply infiltration modifiers. Its classes must
match the index of ``lulc_modifier_table``.
reclass_table : str, Path, or DataFrame, optional
Table mapping HSG classes to bucket parameters.
lulc_modifier_table : str, Path, or DataFrame, optional
Required with ``lulc``. Table of modifier factors keyed by the
land-cover class codes in the supplied dataset.
dual_hsg : {'native', 'drained'}, optional
How to handle dual HSG classes, by default 'drained'. Native dual
classes must be present in the reclassification table.
bucket_loss : float, optional
Uniform bucket loss fraction. Defaults to 0.0 without land use and
0.10 with land-use modifiers.
reproj_method : str, optional
Resampling method for reprojecting final parameter maps to quadtree
blocks. By default 'average'.
"""
_require_lulc_modifiers(lulc, lulc_modifier_table)
if reclass_table is None:
reclass_table = Path(DATADIR) / "infiltration" / "hsg_bucket.csv"
da_hsg = self.data_catalog.get_rasterdataset(
hsg, bbox=self.model.bbox, buffer=10, variables=["hsg"]
)
df_map = self.data_catalog.get_dataframe(
reclass_table,
source_kwargs={"driver": {"name": "pandas", "options": {"index_col": 0}}},
)
da_ksat = None
if ksat is not None:
da_ksat = self.data_catalog.get_rasterdataset(
ksat, bbox=self.model.bbox, buffer=10
)
da_ksat = workflows.mask_invalid_ksat(da_ksat).raster.reproject_like(
da_hsg, method="average"
)
if lulc is not None:
da_lulc = self.data_catalog.get_rasterdataset(
lulc, bbox=self.model.bbox, buffer=10, variables=["lulc"]
)
da_lulc = da_lulc.raster.reproject_like(da_hsg, method="nearest")
df_modifiers = self.data_catalog.get_dataframe(
lulc_modifier_table,
source_kwargs={
"driver": {"name": "pandas", "options": {"index_col": 0}}
},
)
ds = workflows.bucket_from_hsg_landuse(
da_hsg,
da_lulc,
df_map,
df_modifiers,
da_ksat=da_ksat,
dual_hsg=dual_hsg,
bucket_loss=0.10 if bucket_loss is None else bucket_loss,
)
else:
ds = workflows.bucket_from_hsg(
da_hsg,
df_map,
da_ksat=da_ksat,
bucket_loss=bucket_loss if np.isscalar(bucket_loss) else None,
dual_hsg=dual_hsg,
)
outputs = {
name: np.full(self.mask.grid.n_face, np.nan)
for name in flavor_variables("bkt")
}
def compute_block(da_like, ilev=None):
return tuple(
ds[name].raster.reproject_like(da_like, method=reproj_method)
for name in outputs
)
self.compute_quadtree(compute_block, outputs)
self._set_layers(outputs, flavor="bkt")
[docs]
@hydromt_step
def create_bucket_from_maps(
self,
bucket_smax: Union[str, Path, xr.DataArray, xr.Dataset],
bucket_k: Union[str, Path, xr.DataArray, xr.Dataset],
bucket_loss: Union[float, str, Path, xr.DataArray, xr.Dataset, None] = None,
reproj_method: str = "average",
) -> None:
"""Create bucket infiltration from final parameter maps.
Adds model layers:
* **bucket_smax** map: bucket maximum storage [mm]
* **bucket_k** map: bucket drainage coefficient [hr-1]
* **bucket_loss** map: bucket loss fraction [-]
Parameters
----------
bucket_smax, bucket_k : str, Path, RasterDataset
Data with final bucket parameters. Dataset inputs must contain
variables named ``bucket_smax`` and ``bucket_k`` respectively.
bucket_loss : float, str, Path, RasterDataset, optional
Uniform loss fraction or map with final bucket loss fractions.
Defaults to 0.0.
reproj_method : str, optional
Resampling method for reprojecting raster inputs to quadtree blocks.
"""
names = ("bucket_smax", "bucket_k")
layers = {}
raster_sources = {}
for name, source in zip(names, (bucket_smax, bucket_k)):
da = self.data_catalog.get_rasterdataset(
source,
bbox=self.model.bbox,
buffer=10,
variables=[name],
)
raster_sources[name] = da.raster.mask_nodata()
if bucket_loss is None or np.isscalar(bucket_loss):
layers["bucket_loss"] = np.full(
self.mask.grid.n_face,
np.float32(0.0 if bucket_loss is None else bucket_loss),
)
else:
da_loss = self.data_catalog.get_rasterdataset(
bucket_loss, bbox=self.model.bbox, buffer=10, variables=["bucket_loss"]
)
raster_sources["bucket_loss"] = da_loss.raster.mask_nodata()
if raster_sources:
outputs = {
name: np.full(self.mask.grid.n_face, np.nan) for name in raster_sources
}
def compute_block(da_like, ilev=None):
return tuple(
raster_sources[name].raster.reproject_like(
da_like, method=reproj_method
)
for name in outputs
)
self.compute_quadtree(compute_block, outputs)
layers.update(outputs)
self._set_layers(layers, flavor="bkt")
def clear(self) -> None:
"""Clear all infiltration layers from the model."""
self.model.quadtree_grid._data = clear_data(
self.model.quadtree_grid.data, keep=()
)
reset_config(self.model.config)
