move inflow points to confluence if within buffer from model boundary (…

…#202)

* move inflow points to confluence

* fix of test_river_source_points

* force river_inflow buffer to be 100 to match old default value (and fix the test)

* update test and changelog

---------

Co-authored-by: roeldegoede <[email protected]>

docs/changelog.rst

@@ -17,6 +17,8 @@ Added
 Changed
 -------
 - improved subgrid tables are saved as NetCDF (#160)
+- In `SfincsModel.setup_river_inflow`, in case of a confluence within a user-defined buffer of the
+  model boundary the confluence rather than both tributaries is selected as inflow point. (#202)
 
 Fixed
 -----

hydromt_sfincs/sfincs.py

@@ -794,6 +794,7 @@ def setup_river_inflow(
         self,
         rivers: Union[str, Path, gpd.GeoDataFrame] = None,
         hydrography: Union[str, Path, xr.Dataset] = None,
+        buffer: float = 200,
         river_upa: float = 10.0,
         river_len: float = 1e3,
         river_width: float = 500,
@@ -814,7 +815,8 @@ def setup_river_inflow(
         Discharge is set to zero at these points, but can be updated
         using the `setup_discharge_forcing` or `setup_discharge_forcing_from_grid` methods.
 
-        Note: this method assumes the rivers are directed from up- to downstream.
+        Note: this method assumes the rivers are directed from up- to downstream. Use
+        `reverse_river_geom=True` if the rivers are directed from downstream to upstream.
 
         Adds model layers:
 
@@ -831,6 +833,10 @@ def setup_river_inflow(
             Path, data source name, or a xarray raster object for hydrography data.
 
             * Required layers: ['uparea', 'flwdir'].
+        buffer: float, optional
+            Buffer around the model region boundary to define in/outflow points [m],
+            by default 200 m. We suggest to use a buffer of at least twice the hydrography
+            resolution. Inflow points are moved to a downstreawm confluence if within the buffer.
         river_upa : float, optional
             Minimum upstream area threshold for rivers [km2], by default 10.0
         river_len: float, optional
@@ -892,11 +898,6 @@ def setup_river_inflow(
         elif hydrography is None:
             raise ValueError("Either hydrography or rivers must be provided.")
 
-        # estimate buffer based on model resolution
-        buffer = self.reggrid.dx
-        if self.crs.is_geographic:
-            buffer = buffer * 111111.0
-
         # get river inflow / headwater source points
         gdf_src = workflows.river_source_points(
             gdf_riv=gdf_riv,

hydromt_sfincs/workflows/flwdir.py

@@ -1,4 +1,5 @@
 """Flow direction and river network workflows for SFINCS models."""
+
 import logging
 from typing import Tuple
 
@@ -82,7 +83,7 @@ def river_source_points(
     gdf_riv: gpd.GeoDataFrame,
     gdf_mask: gpd.GeoDataFrame,
     src_type: str = "inflow",
-    buffer: float = 100,
+    buffer: float = 200,
     river_upa: float = 10,
     river_len: float = 1e3,
     da_uparea: xr.DataArray = None,
@@ -109,7 +110,8 @@ def river_source_points(
         If 'outflow', return points where the river flows out of the model domain.
         If 'headwater', return all headwater (including inflow) points within the model domain.
     buffer: float, optional
-        Buffer around gdf_mask to select river source points, by default 100 m.
+        Buffer around gdf_mask to select river source points, by default 200 m.
+        Inflow points are moved to a downstream confluence if within the buffer.
     river_upa : float, optional
         Minimum upstream area threshold for rivers [km2], by default 10.0
     river_len: float, optional
@@ -143,6 +145,11 @@ def river_source_points(
 
     # clip river to model gdf_mask
     gdf_riv = gdf_riv.to_crs(gdf_mask.crs).clip(gdf_mask.unary_union)
+    # keep only lines
+    gdf_riv = gdf_riv[
+        [t.endswith("LineString") for t in gdf_riv.geom_type]
+    ].reset_index(drop=True)
+
     # filter river network based on uparea and length
     if "uparea" in gdf_riv.columns:
         gdf_riv = gdf_riv[gdf_riv["uparea"] >= river_upa]
@@ -154,17 +161,19 @@ def river_source_points(
         )
         return gpd.GeoDataFrame()
 
+    # remove lines that fully are within the buffer of the mask boundary
+    bnd = gdf_mask.boundary.buffer(buffer).unary_union
+    gdf_riv = gdf_riv[~gdf_riv.within(bnd)]
+
     # get source points 1m before the start/end of the river
     # 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 = -1 if reverse_river_geom else 1
     gdf_up = gdf_riv.interpolate(dx).to_frame("geometry")
-    gdf_up["riv_idx"] = gdf_riv.index
     gdf_ds = gdf_riv.interpolate(-dx).to_frame("geometry")
-    gdf_ds["riv_idx"] = gdf_riv.index
 
     # get points that do not intersect with up/downstream end of other river segments
-    # use a small buffer of 5m around these points to account for dx and avoid issues with inprecise river geometries
+    # use a small buffer of 5m around these points to account for dx and avoid issues with imprecise river geometries
     if src_type in ["inflow", "headwater"]:
         pnts_ds = gdf_ds.buffer(5).unary_union
         gdf_pnt = gdf_up[~gdf_up.intersects(pnts_ds)].reset_index(drop=True)
@@ -174,7 +183,6 @@ def river_source_points(
 
     # get buffer around gdf_mask, in- and outflow points should be within this buffer
     if src_type in ["inflow", "outflow"]:
-        bnd = gdf_mask.boundary.buffer(buffer).unary_union
         gdf_pnt = gdf_pnt[gdf_pnt.intersects(bnd)].reset_index(drop=True)
 
     # log numer of source points

tests/data/sfincs_test.yml

@@ -56,6 +56,7 @@ setup_drainage_structures:
 
 setup_river_inflow:
   hydrography: merit_hydro
+  buffer: 100
   river_upa: 10
   river_len: 1000
   keep_rivers_geom: True

tests/test_flwdir.py

@@ -39,20 +39,20 @@ def test_river_source_points(hydrography, data_catalog):
     assert gdf_src.index.size == 6
     gdf_src = river_source_points(src_type="outflow", da_uparea=da_uparea, **kwargs)
     assert gdf_src.index.size == 1
-    assert np.isin(["geometry", "riv_idx", "uparea"], gdf_src.columns).all()
+    assert np.isin(["geometry", "uparea"], gdf_src.columns).all()
+    np.allclose(gdf_src.geometry[0].coords[:], [(322650.3, 5044385.7)])
 
     # test reverse oriented line
     gdf_riv = data_catalog.get_geodataframe("rivers_lin2019_v1")
     kwargs = dict(gdf_riv=gdf_riv, gdf_mask=gdf_mask, reverse_river_geom=True)
     gdf_src = river_source_points(src_type="inflow", **kwargs)
     assert gdf_src.index.size == 1  # this data only one river
-    assert gdf_src["riv_idx"].values[0] == 38
+    np.allclose(gdf_src.geometry[0].coords[:], [(322650.3, 5044385.7)])
     gdf_src = river_source_points(src_type="outflow", **kwargs)
+    np.allclose(gdf_src.geometry[0].coords[:], [(322554.0, 5044434.7)])
     assert gdf_src.index.size == 1  # this data only one river
-    assert gdf_src["riv_idx"].values[0] == 34
     gdf_src = river_source_points(src_type="headwater", **kwargs)
     assert gdf_src.index.size == 2
-    assert np.isin(38, gdf_src["riv_idx"].values)
 
     # test errors
     with pytest.raises(ValueError, match="src_type must be either"):