From 6a116fc7ad33a1f7238b31bee06e878d80976b40 Mon Sep 17 00:00:00 2001 From: roeldegoede <83765910+roeldegoede@users.noreply.github.com> Date: Thu, 9 Nov 2023 13:25:09 +0100 Subject: [PATCH] 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 Co-authored-by: Dirk Eilander --- docs/api.rst | 1 + docs/changelog.rst | 1 + docs/user_guide/sfincs_model_setup.rst | 10 +- hydromt_sfincs/sfincs.py | 59 +++++++++++- hydromt_sfincs/workflows/__init__.py | 1 + hydromt_sfincs/workflows/storage_volume.py | 106 +++++++++++++++++++++ tests/test_1model_class.py | 73 ++++++++++++++ 7 files changed, 247 insertions(+), 4 deletions(-) create mode 100644 hydromt_sfincs/workflows/storage_volume.py diff --git a/docs/api.rst b/docs/api.rst index 1e28a77d..5e67ad26 100644 --- 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 index 31d60765..aab25c48 100644 --- 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 index e5067c5b..3b5a3179 100644 --- 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 index e74bba2c..ef67734d 100644 --- 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 index 4ea63132..f51f7eca 100644 --- 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 new file mode 100644 index 00000000..e9d48d3c --- /dev/null +++ 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 index fd29fa77..bc1741de 100644 --- 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+")