From 0edb8f0f6c2ecf557880a5df11c3e71fc2e20f3a Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Tue, 18 Jun 2024 09:55:05 +0200 Subject: [PATCH 1/4] move inflow points to confluence --- hydromt_sfincs/sfincs.py | 13 +++++++------ hydromt_sfincs/workflows/flwdir.py | 17 +++++++++++++---- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py index dafde8bf..65327bcb 100644 --- a/hydromt_sfincs/sfincs.py +++ b/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, diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py index e13e33fe..7387ae29 100644 --- a/hydromt_sfincs/workflows/flwdir.py +++ b/hydromt_sfincs/workflows/flwdir.py @@ -82,7 +82,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 +109,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 +144,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,6 +160,10 @@ 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) @@ -164,7 +174,7 @@ def river_source_points( 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 +184,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 From 6620019fc64b454cfeddd10aeed10d3592a0792d Mon Sep 17 00:00:00 2001 From: roeldegoede Date: Mon, 26 Aug 2024 15:21:35 +0200 Subject: [PATCH 2/4] fix of test_river_source_points --- tests/test_flwdir.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_flwdir.py b/tests/test_flwdir.py index b6d55b1b..5e9f5fd3 100644 --- a/tests/test_flwdir.py +++ b/tests/test_flwdir.py @@ -46,13 +46,13 @@ def test_river_source_points(hydrography, data_catalog): 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 + assert gdf_src["riv_idx"].values[0] == 2 gdf_src = river_source_points(src_type="outflow", **kwargs) assert gdf_src.index.size == 1 # this data only one river - assert gdf_src["riv_idx"].values[0] == 34 + assert gdf_src["riv_idx"].values[0] == 0 gdf_src = river_source_points(src_type="headwater", **kwargs) assert gdf_src.index.size == 2 - assert np.isin(38, gdf_src["riv_idx"].values) + assert np.isin(2, gdf_src["riv_idx"].values) # test errors with pytest.raises(ValueError, match="src_type must be either"): From cce6073c1252caa49cd8250ccedf82f2479b05e0 Mon Sep 17 00:00:00 2001 From: roeldegoede Date: Mon, 26 Aug 2024 16:25:40 +0200 Subject: [PATCH 3/4] force river_inflow buffer to be 100 to match old default value (and fix the test) --- tests/data/sfincs_test.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/data/sfincs_test.yml b/tests/data/sfincs_test.yml index e0f37222..a5089958 100644 --- a/tests/data/sfincs_test.yml +++ b/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 From 6a4f45c0556a10083e2dcb55e0a9cfbf1b959a9b Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Thu, 29 Aug 2024 10:05:08 +0200 Subject: [PATCH 4/4] update test and changelog --- docs/changelog.rst | 2 ++ hydromt_sfincs/workflows/flwdir.py | 3 +-- tests/test_flwdir.py | 8 ++++---- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index 29b356c4..aed8bcb5 100644 --- a/docs/changelog.rst +++ b/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 ----- diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py index 7387ae29..f59ff334 100644 --- a/hydromt_sfincs/workflows/flwdir.py +++ b/hydromt_sfincs/workflows/flwdir.py @@ -1,4 +1,5 @@ """Flow direction and river network workflows for SFINCS models.""" + import logging from typing import Tuple @@ -169,9 +170,7 @@ def river_source_points( # 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 imprecise river geometries diff --git a/tests/test_flwdir.py b/tests/test_flwdir.py index 5e9f5fd3..9e416661 100644 --- a/tests/test_flwdir.py +++ b/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] == 2 + 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] == 0 gdf_src = river_source_points(src_type="headwater", **kwargs) assert gdf_src.index.size == 2 - assert np.isin(2, gdf_src["riv_idx"].values) # test errors with pytest.raises(ValueError, match="src_type must be either"):