several fixes for river in- and outflow points. Ref #135 (#136)

* fix value error in setup_river_inflow. Ref #135 * Fix river_boundary_points, add reverse option * Update changelog.rst * Update hydromt_sfincs/workflows/flwdir.py Co-authored-by: DirkEilander <[email protected]> * Edit setup_river)outflow, fix type error * add warning for CRS and comments to make it easier to read river_boundary_points --------- Co-authored-by: DirkEilander <[email protected]>
Deltares · Nov 9, 2023 · 4915572 · 4915572
1 parent 2b2bc18
commit 4915572
Show file tree

Hide file tree

Showing 3 changed files with 44 additions and 9 deletions.
diff --git a/docs/changelog.rst b/docs/changelog.rst
@@ -11,16 +11,19 @@ v1.0.2 (unreleased)
 
 Added
 -----
+- 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
 
 Changed
 -------
 - Changed `setup_cn_infiltration_with_kr` into `setup_cn_infiltration_with_ks` since saturated hydraulic conductivity (ks) is used instead of recovery rate (kr) PR #126
 
+
 Fixed
 -----
 - writing COG files in `SfincsModel.setup_subgrid` (the COG driver settings were wrong) PR #117
 - a constant offset in the `datasets_dep` argument to `SfincsModel.setup_subgrid` and `SfincsModel.setup_dep` was ignored PR #119
+- Bugfixes in workflows.river_boundary_points to make sure function also works with geoDataFrame #PR 136
 - mismatch between gis data and the model grid causing issues while reading the model PR #128
 - `utils.downscale_floodmap` now also works for large (rotated) grids PR #145
 

diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py
@@ -19,7 +19,7 @@
 from hydromt.vector import GeoDataArray, GeoDataset
 from hydromt.workflows.forcing import da_to_timedelta
 from pyproj import CRS
-from shapely.geometry import box, LineString
+from shapely.geometry import LineString, box
 
 from . import DATADIR, plots, utils, workflows
 from .regulargrid import RegularGrid
@@ -756,6 +756,7 @@ def setup_river_inflow(
         merge: bool = False,
         first_index: int = 1,
         keep_rivers_geom: bool = False,
+        reverse_river_geom: bool = False,
     ):
         """Setup discharge (src) points where a river enters the model domain.
 
@@ -799,8 +800,9 @@ def setup_river_inflow(
             First index for the river source points, by default 1.
         keep_rivers_geom: bool, optional
             If True, keep a geometry of the rivers "rivers_inflow" in geoms. By default False.
-        buffer: int, optional
-            Buffer [no. of cells] around model domain, by default 10.
+        reverse_river_geom: bool, optional
+            If True, assume that segments in 'rivers' are drawn from downstream to upstream.
+            Only used if 'rivers' is not None, By default False
 
         See Also
         --------
@@ -817,7 +819,11 @@ def setup_river_inflow(
             )
             da_flwdir = ds["flwdir"]
             da_uparea = ds["uparea"]
-        elif rivers == "rivers_outflow" and rivers in self.geoms:
+        elif (
+            isinstance(rivers, str)
+            and rivers == "rivers_outflow"
+            and rivers in self.geoms
+        ):
             # reuse rivers from setup_river_in/outflow
             gdf_riv = self.geoms[rivers]
         elif rivers is not None:
@@ -836,6 +842,8 @@ def setup_river_inflow(
             river_len=river_len,
             river_upa=river_upa,
             inflow=True,
+            reverse_river_geom=reverse_river_geom,
+            logger=self.logger,
         )
         n = len(gdf_src.index)
         self.logger.info(f"Found {n} river inflow points.")
@@ -878,6 +886,7 @@ def setup_river_outflow(
         keep_rivers_geom: bool = False,
         reset_bounds: bool = False,
         btype: str = "outflow",
+        reverse_river_geom: bool = False,
     ):
         """Setup open boundary cells (mask=3) where a river flows
         out of the model domain.
@@ -921,6 +930,9 @@ def setup_river_outflow(
             by default False.
         btype: {'waterlevel', 'outflow'}
             Boundary type
+        reverse_river_geom: bool, optional
+            If True, assume that segments in 'rivers' are drawn from downstream to upstream.
+            Only used if rivers is not None, By default False
 
         See Also
         --------
@@ -936,7 +948,11 @@ def setup_river_outflow(
             )
             da_flwdir = ds["flwdir"]
             da_uparea = ds["uparea"]
-        elif rivers == "rivers_inflow" and rivers in self.geoms:
+        elif (
+            isinstance(rivers, str)
+            and rivers == "rivers_inflow"
+            and rivers in self.geoms
+        ):
             # reuse rivers from setup_river_in/outflow
             gdf_riv = self.geoms[rivers]
         elif rivers is not None:
@@ -956,6 +972,8 @@ def setup_river_outflow(
             river_len=river_len,
             river_upa=river_upa,
             inflow=False,
+            reverse_river_geom=reverse_river_geom,
+            logger=self.logger,
         )
 
         if len(gdf_out) > 0:

diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py
@@ -70,6 +70,8 @@ def river_boundary_points(
     river_upa: float = 10,
     river_len: float = 1e3,
     inflow: bool = True,
+    reverse_river_geom: bool = False,
+    logger: logging.Logger = logger,
 ) -> Tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
     """Returns the locations where a river flows in (`inflow=True`)
     or out (`inflow=False`) of the model region.
@@ -95,6 +97,9 @@ def river_boundary_points(
         Mimimum river length [m] within the model domain to define river cells, by default 1000 m.
     inflow: bool, optional
         If True, return inflow otherwise outflow boundary points, by default True.
+    reverse_river_geom: bool, optional
+        If True, assume that segments in 'rivers' are drawn from downstream to upstream.
+        Only used if 'rivers' is not None, By default False
 
     Returns
     -------
@@ -105,6 +110,11 @@ def river_boundary_points(
         region.geometry.type == "Polygon"
     ):
         raise ValueError("Boundary must be a GeoDataFrame of LineStrings.")
+    if res > 0.01 and region.crs.is_geographic:
+        # provide warning
+        logger.warning(
+            "The region crs is geographic, while the resolution seems to be in meters."
+        )
 
     if gdf_riv is None and (da_flwdir is None or da_uparea is None):
         raise ValueError("Either gdf_riv or da_flwdir and da_uparea must be provided.")
@@ -117,11 +127,15 @@ def river_boundary_points(
         gdf_riv = gdf_riv.to_crs(region.crs).clip(region.unary_union)
         # filter river network based on uparea and length
         if "uparea" in gdf_riv.columns:
-            gdf_riv = [gdf_riv["uparea"] >= river_upa]
+            gdf_riv = gdf_riv[gdf_riv["uparea"] >= river_upa]
         if "rivlen" in gdf_riv.columns:
             gdf_riv = gdf_riv[gdf_riv["rivlen"] > river_len]
-
+    # a positive dx results in a point near the start of the line (inflow)
+    # a negative dx results in a point near the end of the line (outflow)
     dx = res / 5 if inflow else -res / 5
+    if reverse_river_geom:
+        dx = -dx
+
     # move point a bit into the model domain
     gdf_pnt = gdf_riv.interpolate(dx).to_frame("geometry")
     # keep points on boundary cells
@@ -133,8 +147,8 @@ def river_boundary_points(
         gdf_pnt["uparea"] = da_uparea.raster.sample(gdf_pnt).values
         gdf_pnt = gdf_pnt.sort_values("uparea", ascending=False).reset_index(drop=True)
     if "rivwth" in gdf_riv.columns:
-        gdf_out = hydromt.gis_utils.nearest_merge(
-            gdf_out, gdf_riv, columns=["rivwth"], max_dist=10
+        gdf_pnt = hydromt.gis_utils.nearest_merge(
+            gdf_pnt, gdf_riv, columns=["rivwth"], max_dist=10
         )
 
     return gdf_pnt, gdf_riv