ENH: added setup_storage_volume to account for green-infratructure (#101

) * added setup_storage_volume to account for green-infra * code optimization and "vol" added to _maps * write observations station names as strings again * fixing single_height in setup_greeninfra * add test for storage volume * extend test for location of the point-source in grid * move storage_volume intelligence to workflow * clean-up code * update docs * reformat with black * black reformat test --------- Co-authored-by: GundulaW <[email protected]> Co-authored-by: Dirk Eilander <[email protected]>
Deltares · Nov 9, 2023 · 6a116fc · 6a116fc
1 parent 94a5daa
commit 6a116fc
Show file tree

Hide file tree

Showing 7 changed files with 247 additions and 4 deletions.
diff --git a/docs/api.rst b/docs/api.rst
@@ -46,6 +46,7 @@ Setup components
    SfincsModel.setup_observation_lines
    SfincsModel.setup_structures
    SfincsModel.setup_drainage_structures
+   SfincsModel.setup_storage_volume
    SfincsModel.setup_waterlevel_forcing
    SfincsModel.setup_waterlevel_bnd_from_mask
    SfincsModel.setup_discharge_forcing

diff --git a/docs/changelog.rst b/docs/changelog.rst
@@ -11,6 +11,7 @@ v1.0.2 (unreleased)
 
 Added
 -----
+- `SfincsModel.setup_storage_volume` to account for green-infrastructure PR #101
 - Added reverse_river_geom keyword argument in workflows.river_boundary_points #PR 136
 - the COG files that are written automatically contain overviews for faster visualization PR #144
 

diff --git a/docs/user_guide/sfincs_model_setup.rst b/docs/user_guide/sfincs_model_setup.rst
@@ -103,10 +103,14 @@ Grid setup methods
      - This component adds a spatially varying constant infiltration rate map (qinffile) to the model grid.
    * - :py:func:`~hydromt_sfincs.SfincsModel.setup_cn_infiltration`
      - This component adds a potential maximum soil moisture retention map (scsfile) to the model grid based on a gridded curve number map.
-   * - :py:func:`~hydromt_sfincs.SfincsModel.setup_cn_infiltration_with_kr`
-     - This component adds a three layers related to the curve number (maximum and effective infiltration capacity; seff and smax) and recovery rate (kr) to the model grid based on landcover, Hydrological Similarity Group and saturated hydraulic conductivity (Ksat).
+   * - :py:func:`~hydromt_sfincs.SfincsModel.setup_cn_infiltration_with_ks`
+     - This component adds a three layers related to the curve number (maximum and effective infiltration capacity; seff and smax) and saturated hydraulic conductivity (ks, to account for recovery) to the model 
+     grid based on landcover, Hydrological Similarity Group and saturated hydraulic conductivity (Ksat).
+   * - :py:funct:`~hydromt_sfincs.SfincsModel.setup_storage_volume~
+     - This component adds a storage volume map (volfile) to the model grid to account for green-infrastructure.
    * - :py:func:`~hydromt_sfincs.SfincsModel.setup_subgrid`
-     - This component generates subgrid tables (sbgfile) for the model grid based on a list of elevation and Manning roughness datasets     
+     - This component generates subgrid tables (sbgfile) for the model grid based on a list of elevation and Manning roughness datasets
+
 Geoms setup methods
 -------------------
 

diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py
@@ -64,7 +64,7 @@ class SfincsModel(GridModel):
         ),
     }
     _FORCING_SPW = {"spiderweb": "spw"}  # TODO add read and write functions
-    _MAPS = ["msk", "dep", "scs", "manning", "qinf", "smax", "seff", "kr"]
+    _MAPS = ["msk", "dep", "scs", "manning", "qinf", "smax", "seff", "ks", "vol"]
     _STATES = ["rst", "ini"]
     _FOLDERS = []
     _CLI_ARGS = {"region": "setup_grid_from_region", "res": "setup_grid_from_region"}
@@ -79,6 +79,7 @@ class SfincsModel(GridModel):
         },
         "qinf": {"standard_name": "infiltration rate", "unit": "mm.hr-1"},
         "manning": {"standard_name": "manning roughness", "unit": "s.m-1/3"},
+        "vol": {"standard_name": "storage volume", "unit": "m3"},
         "bzs": {"standard_name": "waterlevel", "unit": "m+ref"},
         "bzi": {"standard_name": "wave height", "unit": "m"},
         "dis": {"standard_name": "discharge", "unit": "m3.s-1"},
@@ -1595,6 +1596,62 @@ def setup_drainage_structures(
         self.set_geoms(gdf_structures, "drn")
         self.set_config("drnfile", f"sfincs.drn")
 
+    def setup_storage_volume(
+        self,
+        storage_locs: Union[str, Path, gpd.GeoDataFrame],
+        volume: Union[float, List[float]] = None,
+        height: Union[float, List[float]] = None,
+        merge: bool = True,
+    ):
+        """Setup storage volume.
+
+        Adds model layer:
+        * **vol** map: storage volume for green infrastructure
+
+        Parameters
+        ----------
+        storage_locs : str, Path
+            Path, data source name, or geopandas object to storage location polygon or point geometry file.
+            Optional "volume" or "height" attributes can be provided to set the storage volume.
+        volume : float, optional
+            Storage volume [m3], by default None
+        height : float, optional
+            Storage height [m], by default None
+        merge : bool, optional
+            If True, merge with existing storage volumes, by default True.
+
+        """
+
+        # read, clip and reproject
+        gdf = self.data_catalog.get_geodataframe(
+            storage_locs,
+            geom=self.region,
+            buffer=10,
+        ).to_crs(self.crs)
+
+        if self.grid_type == "regular":
+            # if merge, add new storage volumes to existing ones
+            if merge and "vol" in self.grid:
+                da_vol = self.grid["vol"]
+            else:
+                da_vol = xr.full_like(self.mask, 0, dtype=np.float64)
+
+            # add storage volumes form gdf to da_vol
+            da_vol = workflows.add_storage_volume(
+                da_vol,
+                gdf,
+                volume=volume,
+                height=height,
+                logger=self.logger,
+            )
+
+            # set grid
+            mname = "vol"
+            da_vol.attrs.update(**self._ATTRS.get(mname, {}))
+            self.set_grid(da_vol, name=mname)
+            # update config
+            self.set_config(f"{mname}file", f"sfincs.{mname[:3]}")
+
     ### FORCING
     def set_forcing_1d(
         self,

diff --git a/hydromt_sfincs/workflows/__init__.py b/hydromt_sfincs/workflows/__init__.py
@@ -7,3 +7,4 @@
 from .merge import *
 from .tiling import *
 from .curvenumber import *
+from .storage_volume import *
diff --git a/hydromt_sfincs/workflows/storage_volume.py b/hydromt_sfincs/workflows/storage_volume.py
@@ -0,0 +1,106 @@
+import logging
+from typing import Union, List
+
+import geopandas as gpd
+import numpy as np
+import xarray as xr
+
+logger = logging.getLogger(__name__)
+
+
+def add_storage_volume(
+    da_vol: xr.DataArray,
+    gdf: gpd.GeoDataFrame,
+    volume: Union[float, List[float]] = None,
+    height: Union[float, List[float]] = None,
+    logger=logger,
+) -> xr.DataArray:
+    """Add storage volume to a grid based on a GeoDataFrame with storage locations.
+
+    Parameters
+    ----------
+    da_vol : xr.DataArray
+        DataArray with the grid to which the storage volume should be added.
+    gdf : gpd.GeoDataFrame
+        GeoDataFrame with the storage locations (polygon or point geometry file).
+        Optional "volume" or "height" attributes can be provided to set the storage volume.
+    volume : Union[float, List[float]], optional
+        Volume of the storage locations [m3], by default None.
+    height : Union[float, List[float]], optional
+        Height of the storage locations [m], by default None.
+
+    Returns
+    -------
+    xr.DataArray
+        DataArray with the grid including the storage volume.
+
+    """
+
+    # loop over the gdf rows and rasterize each geometry
+    for i, row in gdf.iterrows():
+        # create a gdf with only the current row
+        single_gdf = gpd.GeoDataFrame(gdf.loc[[i]]).reset_index(drop=True)
+
+        single_vol = single_gdf.get("volume", np.nan)
+        single_height = single_gdf.get("height", np.nan)
+
+        # check if volume or height is provided in the gdf or as input
+        if np.isnan(float(single_vol)):  # volume not provided or nan
+            if np.isnan(float(single_height)):  # height not provided or nan
+                if volume is not None:  # volume provided as input (list)
+                    single_vol = volume if not isinstance(volume, list) else volume[i]
+                elif height is not None:  # height provided as input (list)
+                    single_height = (
+                        height if not isinstance(height, list) else height[i]
+                    )
+                else:  # no volume or height provided
+                    logger.warning(
+                        f"No volume or height provided for storage location {i}"
+                    )
+                    continue
+            else:  # height provided in gdf
+                single_height = single_height[0]
+        else:  # volume provided in gdf
+            single_vol = single_vol[0]
+
+        # check if gdf has point or polyhon geometry
+        if single_gdf.geometry.type[0] == "Point":
+            # get x and y coordinate of the point in crs of the grid
+            x, y = single_gdf.geometry.iloc[0].coords[0]
+
+            # check if the grid is rotated
+            if da_vol.raster.rotation != 0:
+                # rotate the point
+                x, y = ~da_vol.raster.transform * (x, y)
+                # select the grid cell nearest to the point
+                closest_point = da_vol.reindex(x=da_vol.x, y=da_vol.y).sel(
+                    x=x, y=y, method="nearest"
+                )
+            else:
+                # select the grid cell nearest to the point
+                closest_point = da_vol.sel(x=x, y=y, method="nearest")
+
+            # add the volume to the grid cell
+            if not np.isnan(float(single_vol)):
+                da_vol.loc[
+                    dict(x=closest_point.x.item(), y=closest_point.y.item())
+                ] += single_vol
+            else:
+                logger.warning(f"No volume provided for storage location of type Point")
+
+        elif single_gdf.geometry.type[0] == "Polygon":
+            # rasterize the geometry
+            area = da_vol.raster.rasterize_geometry(single_gdf, method="area")
+            total_area = area.sum().values
+
+            # calculate the volume per cell and add it to the grid
+            if not np.isnan(float(single_vol)):
+                da_vol += area / total_area * single_vol
+            elif not np.isnan(float(single_height)):
+                da_vol += area * single_height
+            else:
+                logger.warning(
+                    f"No volume or height provided for storage location of type Polygon"
+                )
+
+    return da_vol
diff --git a/tests/test_1model_class.py b/tests/test_1model_class.py
@@ -3,8 +3,10 @@
 from os.path import isfile, join
 
 import numpy as np
+import geopandas as gpd
 import pandas as pd
 import pytest
+from shapely.geometry import Polygon, Point
 import xarray as xr
 from geopandas.testing import assert_geodataframe_equal
 from hydromt.cli.cli_utils import parse_config
@@ -192,6 +194,77 @@ def test_drainage_structures(tmpdir):
     assert len(mod.geoms["drn"].index) == nr_drainage_structures * 2
 
 
+def test_storage_volume(tmpdir):
+    tmp_root = str(tmpdir.join("storage_volume_test"))
+
+    # create two arbitrary but overlapping polygons
+    coords1 = [
+        (3.5, 3.5),
+        (6.5, 3.5),
+        (6.5, 6.5),
+        (4.5, 6.5),
+        (4.5, 7.25),
+        (6.5, 7.25),
+        (6.5, 8),
+        (3, 8),
+    ]
+    poly1 = Polygon(coords1)
+    # second polygon which overlaps aprtly with the first but is smaller
+    coords2 = [(6, 3), (7, 3), (7, 4), (6, 4)]
+    poly2 = Polygon(coords2)
+    # create a geodataframe with the two polygons
+    gdf = gpd.GeoDataFrame({"geometry": [poly1, poly2]}, crs=4326)
+    gdf["volume"] = [None, 1000]
+
+    # also create an arbitrary point
+    point = Point(5, 6)
+    point_gdf = gpd.GeoDataFrame({"geometry": [point]}, crs=4326)
+    point_gdf["volume"] = 20
+
+    # create a sfincs model
+    mod = SfincsModel(root=tmp_root, mode="w+")
+    mod.setup_grid_from_region(region={"bbox": [0, 0, 10, 10]}, res=20000, crs="utm")
+
+    # test setup_storage_volume with polygons
+    # one polygon has no volume specifed, the other has a volume of 1000
+    # the non-specified gets the volume of the input argument
+    mod.setup_storage_volume(storage_locs=gdf, volume=10000)
+
+    assert mod.grid["vol"].sum() == 11000
+
+    # test setup_storage_volume with points
+    mod.setup_storage_volume(storage_locs=point_gdf, merge=True)
+
+    assert mod.grid["vol"].sum() == 11020
+
+    # now redo the tests with a rotated grid
+    config = mod.config.copy()
+    mod = SfincsModel(root=tmp_root, mode="w+")
+
+    # get the config from the first model and add a rotation
+    config["rotation"] = 10
+    mod.config.update(config)
+    mod.update_grid_from_config()
+
+    # test setup_storage_volume with
+    # drop volume column from gdf
+    gdf = gdf.drop(columns=["volume"])
+    mod.setup_storage_volume(storage_locs=gdf, volume=[350, 800])
+
+    # check if the volumes are correct
+    assert np.isclose(mod.grid["vol"].sum(), 1150)
+
+    # drop volume column from gdf
+    point_gdf = point_gdf.drop(columns=["volume"])
+    mod.setup_storage_volume(storage_locs=point_gdf, volume=34.5, merge=False)
+
+    assert np.isclose(mod.grid["vol"].sum(), 34.5)
+
+    # check index of the point with maximum volume
+    index = mod.grid["vol"].argmax()
+    assert index == 1601
+
+
 def test_observations(tmpdir):
     root = TESTMODELDIR
     mod = SfincsModel(root=root, mode="r+")