add setup_river_floodplain method

docs/api.rst

@@ -30,6 +30,7 @@ Setup components
    WflowModel.setup_config
    WflowModel.setup_basemaps
    WflowModel.setup_rivers
+   WflowModel.setup_river_floodplain
    WflowModel.setup_lakes
    WflowModel.setup_reservoirs
    WflowModel.setup_glaciers

docs/changelog.rst

@@ -12,12 +12,11 @@ Unreleased
 
 Added
 -----
-
+- New setup_river_floodplain method `PR #123 <https://github.com/Deltares/hydromt_wflow/pull/123>`_
 
 Changed
 -------
 
-
 Fixed
 -----
 - Bugfix with wrong nodata value in the hydrography method which caused errors for model which where not based on (sub)basins `PR #144 <https://github.com/Deltares/hydromt_wflow/pull/144>`_

docs/user_guide/wflow_model_setup.rst

@@ -38,6 +38,8 @@ a specific method see its documentation.
       - This component sets the region of interest and res (resolution in degrees) of the model.
     * - :py:func:`~WflowModel.setup_rivers`
       - This component sets the all river parameter maps.
+    * - :py:func:`~WflowModel.setup_river_floodplain`
+        This component sets a map with floodplain volume per flood depth.
     * - :py:func:`~WflowModel.setup_lakes`
       - This component generates maps of lake areas and outlets as well as parameters with average lake area, depth a discharge values.
     * - :py:func:`~WflowModel.setup_reservoirs`

examples/wflow_build.ini

@@ -25,6 +25,13 @@ smooth_len       = 5000             # length over which to smooth river depth an
 # river_routing   = local-inertial  # {'kinematic-wave', 'local-inertial'}
 # land_routing    = kinematic-wave  # {'kinematic-wave', 'local-inertial'}
 
+[setup_river_floodplain]
+hydrography_fn   = merit_hydro      # source hydrography data, should correspond to hydrography_fn in setup_basemaps
+river_upa        = 30               # minimum upstream area threshold for the river map [km2], default = 30 km2
+flood_depths     = [0.5, 1.0, 1.5, 2.0, 2.5] # flood depths at which a volume is derived, by default [0.5,1.0,1.5,2.0,2.5]
+floodplain_1d    = True             # Value for model.floodplain_1d setting in toml file, by default True.
+
+
 [setup_reservoirs]
 reservoirs_fn   = hydro_reservoirs  # source for reservoirs based on GRAND: {hydro_reservoirs}
 min_area        = 1.0           # minimum lake area to consider [km2]

examples/wflow_piave_clip/wflow_sbm.toml

@@ -34,6 +34,7 @@ kw_land_tstep = 3600
 thicknesslayers = [ 100, 300, 800,]
 river_routing = "local-inertial"
 land_routing = "local-inertial"
+floodplain_1d = true
 
 [output]
 path = "run_default/output.nc"
@@ -137,6 +138,9 @@ q_av = "q_river"
 [output.lateral.land]
 h = "h_land"
 
+[input.lateral.river.floodplain]
+volume = "floodplain_volume"
+
 [input.lateral.river.reservoir]
 area = "ResSimpleArea"
 areas = "wflow_reservoirareas"

examples/wflow_piave_subbasin/wflow_sbm.toml

@@ -34,6 +34,7 @@ kw_land_tstep = 3600
 thicknesslayers = [ 100, 300, 800,]
 river_routing = "local-inertial"
 land_routing = "local-inertial"
+floodplain_1d = true
 
 [output]
 path = "run_default/output.nc"
@@ -143,6 +144,9 @@ volume = "volume_reservoir"
 [state.lateral.river.lake]
 waterlevel = "waterlevel_lake"
 
+[input.lateral.river.floodplain]
+volume = "floodplain_volume"
+
 [input.lateral.river.reservoir]
 area = "ResSimpleArea"
 areas = "wflow_reservoirareas"

hydromt_wflow/wflow.py

@@ -13,6 +13,7 @@
 import pyproj
 import toml
 import codecs
+import pyflwdir
 from pyflwdir import core_d8, core_ldd, core_conversion
 from dask.diagnostics import ProgressBar
 import logging
@@ -392,6 +393,64 @@ def setup_rivers(
         self.staticgeoms.pop("rivers", None)  # remove old rivers if in staticgeoms
         self.rivers  # add new rivers to staticgeoms
 
+    def setup_river_floodplain(
+        self,
+        hydrography_fn,
+        river_upa: float = 30,
+        flood_depths: List = [0.5, 1.0, 1.5, 2.0, 2.5],
+        floodplain_1d: bool = True,
+    ):
+        """
+        This component adds a map with floodplain volume per flood depth,
+        which is used in the wflow 1D floodplain schematisation.
+
+        Requires `setup_rivers` to be executed beforehand.
+
+        Adds model layers:
+
+        * **floodplain_volume** map: river mask [-]
+
+        Parameters
+        ----------
+        hydrography_fn : str, Path
+            Name of data source for hydrography data.
+            Must be same as setup_basemaps for consistent results.
+
+            * Required variables: ['flwdir', 'uparea', 'elevtn']
+        river_upa : float
+            minimum upstream area threshold for drain in the HAND.
+        flood_depths : tuple of float, optional
+            flood depths at which a volume is derived, by default [0.5,1.0,1.5,2.0,2.5]
+        floodplain_1d : bool
+            Value for model.floodplain_1d setting in toml file, by default True.
+        """
+        if not hasattr(pyflwdir.FlwdirRaster, "ucat_volume"):
+            self.logger.warning("This method requires pyflwdir >= 0.5.6")
+            return
+
+        self.logger.info("Preparing 1D river floodplain_volume map.")
+
+        # read data
+        ds_hydro = self.data_catalog.get_rasterdataset(
+            hydrography_fn, geom=self.region, buffer=10
+        )
+
+        # get river floodplain volume
+        inv_rename = {v: k for k, v in self._MAPS.items() if v in self.staticmaps}
+        da_fldpln = workflows.river_floodplain_volume(
+            ds=ds_hydro,
+            ds_model=self.staticmaps.rename(inv_rename),
+            river_upa=river_upa,
+            flood_depths=flood_depths,
+            logger=self.logger,
+        )
+        self.set_staticmaps(da_fldpln, "floodplain_volume")
+
+        # update config
+        self.logger.debug(f'Update wflow config model.floodplain_1d="{floodplain_1d}"')
+        self.set_config("model.floodplain_1d", floodplain_1d)
+        self.set_config("input.lateral.river.floodplain.volume", "floodplain_volume")
+
     def setup_hydrodem(
         self,
         elevtn_map="wflow_dem",

hydromt_wflow/workflows/river.py

@@ -15,7 +15,7 @@
 
 logger = logging.getLogger(__name__)
 
-__all__ = ["river", "river_bathymetry", "river_width"]
+__all__ = ["river", "river_bathymetry", "river_width", "river_floodplain_volume"]
 
 RESAMPLING = {"climate": "nearest", "precip": "average"}
 NODATA = {"discharge": -9999}
@@ -61,7 +61,7 @@ def river(
         in window based smoothing of river length, by default 0.0. 
         The smoothing is skipped if min_riverlen_ratio = 0.
     channel_dir: {"up", "down"}
-        flow direcition in which to calculate (subgrid) river length and width
+        flow direction in which to calculate (subgrid) river length and width
 
     Returns:
     ds_out: xr.Dataset
@@ -74,7 +74,7 @@ def river(
         raise ValueError(f"One or more variables missing from ds: {dvars}.")
     if ds_model is not None and not np.all([v in ds_model for v in dvars_model]):
         raise ValueError(f"One or more variables missing from ds_model: {dvars_model}")
-    # sort sugrid
+    # sort subgrid
     subgrid = True
     if ds_model is None or not np.all([v in ds_model for v in ["x_out", "y_out"]]):
         subgrid = False
@@ -304,6 +304,101 @@ def river_bathymetry(
     return ds_model[["rivwth", "rivdph"]]
 
 
+def river_floodplain_volume(
+    ds,
+    ds_model,
+    river_upa=30,
+    flood_depths=[0.5, 1.0, 1.5, 2.0, 2.5],
+    dtype=np.float64,
+    logger=logger,
+):
+    """Calculate the floodplain volume at given flood depths based on a (subgrid) HAND map.
+
+    Parameters
+    ----------
+    ds: xr.Dataset
+        hydrography dataset containing "flwdir", "uparea", "elevtn" variables;
+    ds_model: xr.Dataset, optional
+        Model dataset with output grid, must contain "rivmsk", "rivwth", "rivlen"
+        for subgrid must contain "x_out", "y_out".
+    river_upa: float
+        minimum threshold to define the river when calculating HAND, by default 30 [km2]
+    flood_depths : list of float, optional
+        flood depths at which a volume is derived, by default [0.5,1.0,1.5,2.0,2.5]
+    dtype: numpy.dtype, optional
+        output dtype
+
+    Returns
+    -------
+    da_out : xr.DataArray with dims (flood_depth, y, x)
+        Floodplain volume [m3]
+    """
+    if not ds.raster.crs == ds_model.raster.crs:
+        raise ValueError("Hydrography dataset CRS does not match model CRS")
+
+    # check data variables.
+    dvars = ["flwdir", "elevtn", "uparea"]
+    dvars_model = ["rivmsk", "rivwth", "rivlen"]
+    if not np.all([v in ds for v in dvars]):
+        raise ValueError(f"One or more variables missing from ds: {dvars}.")
+    if not np.all([v in ds_model for v in dvars_model]):
+        raise ValueError(f"One or more variables missing from ds_model: {dvars_model}")
+
+    # initialize flood directions
+    flwdir = flw.flwdir_from_da(ds["flwdir"], mask=None)
+
+    # get river cell coordinates
+    riv_mask = ds_model["rivmsk"].values > 0
+    if "x_out" in ds_model and "y_out" in ds_model:  # use subgrid
+        idxs_out = ds.raster.xy_to_idx(
+            xs=ds_model["x_out"].values,
+            ys=ds_model["y_out"].values,
+            mask=riv_mask,
+            nodata=flwdir._mv,
+        )
+    else:
+        # get cell index of river cells
+        idxs_out = np.arange(ds_model.raster.size).reshape(ds_model.raster.shape)
+        idxs_out = np.where(riv_mask, idxs_out, flwdir._mv)
+
+    # calculate hand
+    hand = flwdir.hand(
+        drain=ds["uparea"].values >= river_upa, elevtn=ds["elevtn"].values
+    )
+
+    # calculate volume at user defined flood depths
+    # note that the ucat_map is not save currently but is required if
+    # we want to create a flood map by postprocessing
+    flood_depths = np.atleast_1d(flood_depths).ravel().astype(dtype)  # force 1D
+    _, ucat_vol = flwdir.ucat_volume(idxs_out=idxs_out, hand=hand, depths=flood_depths)
+    # force minimum volume based on river width * length
+    rivlen = np.maximum(ds_model["rivlen"].values, 1)  # avoid zero division errors
+    min_vol = ds_model["rivwth"].values * rivlen * flood_depths[0]
+    ucat_vol[0, :, :] = np.maximum(ucat_vol[0, :, :], min_vol).round(0)
+    fldwth = ucat_vol[0, :, :] / rivlen / flood_depths[0]
+    for i in range(1, ucat_vol.shape[0]):
+        dh = flood_depths[i] - flood_depths[i - 1]
+        min_vol = ucat_vol[i - 1, ...] + (fldwth * rivlen * dh)
+        ucat_vol[i, :, :] = np.maximum(ucat_vol[i, :, :], min_vol).round(0)
+        fldwth = (ucat_vol[i, :, :] - ucat_vol[i - 1, :, :]) / rivlen / dh
+
+    # return xarray DataArray
+    da_out = (
+        xr.DataArray(
+            coords={"flood_depth": flood_depths, **ds_model.raster.coords},
+            dims=("flood_depth", *ds_model.raster.dims),
+            data=ucat_vol,
+        )
+        .reset_coords(drop=True)
+        .where(ds_model["rivmsk"] > 0, -9999.0)
+    )
+    da_out.raster.set_nodata(-9999.0)
+    da_out.raster.set_crs(ds_model.raster.crs)
+
+    # TODO return a second dataset with hand and ucat_map variables for postprocessing
+    return da_out
+
+
 # TODO: methods below are redundant after version v0.1.4 and will be removed in
 # future versions
 
@@ -328,9 +423,7 @@ def river_width(
 ):
     nopars = a is None or b is None  # no manual a, b parameters
     fit = fit or (nopars and predictor not in ["discharge"])  # fit power-law on the fly
-    nowth = (
-        f"{rivwth_name}{obs_postfix}" not in ds_like
-    )  # no obseved with in staticmaps
+    nowth = f"{rivwth_name}{obs_postfix}" not in ds_like  # no observed width
     fill = fill and nowth == False  # fill datagaps and masked areas (lakes/res) in obs
 
     if nowth and fit:

tests/data/wflow_piave_build_subbasin.ini

@@ -21,6 +21,12 @@ elevtn_map      = wflow_dem     # {'dem_subgrid', 'wflow_dem'}
 river_routing   = local-inertial  # {'kinematic-wave', 'local-inertial'}
 land_routing    = local-inertial  # {'kinematic-wave', 'local-inertial'}
 
+[setup_river_floodplain]
+hydrography_fn   = merit_hydro      # source hydrography data, should correspond to hydrography_fn in setup_basemaps
+river_upa        = 30               # minimum upstream area threshold for the river map [km2], default = 30 km2
+flood_depths     = [0.5, 1.0, 1.5, 2.0, 2.5] # flood depths at which a volume is derived, by default [0.5,1.0,1.5,2.0,2.5]
+floodplain_1d    = True             # Value for model.floodplain_1d setting in toml file, by default True.
+
 [setup_reservoirs]
 reservoirs_fn = hydro_reservoirs
 min_area = 0.0