From 9eab0611b8f4f470a7f8c288da5fe21c362ce386 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Wed, 3 Jan 2024 08:56:39 +0100 Subject: [PATCH 1/6] [WIP] start with including headwater points --- docs/api.rst | 2 +- hydromt_sfincs/sfincs.py | 113 ++++++++++++---- hydromt_sfincs/workflows/flwdir.py | 210 ++++++++++++++++++++--------- tests/conftest.py | 15 +++ tests/test_flwdir.py | 78 +++++++++++ 5 files changed, 323 insertions(+), 95 deletions(-) create mode 100644 tests/test_flwdir.py diff --git a/docs/api.rst b/docs/api.rst index 5e67ad26..1d85ab5c 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -154,7 +154,7 @@ SFINCS workflows workflows.merge_dataarrays workflows.burn_river_rect workflows.snap_discharge - workflows.river_boundary_points + workflows.river_source_points workflows.river_centerline_from_hydrography workflows.landuse workflows.cn_to_s diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py index 4d4e8d84..f71d5bfd 100644 --- a/hydromt_sfincs/sfincs.py +++ b/hydromt_sfincs/sfincs.py @@ -762,6 +762,8 @@ def setup_river_inflow( first_index: int = 1, keep_rivers_geom: bool = False, reverse_river_geom: bool = False, + src_type: str = "inflow", + mask_upstream_cells: bool = False, ): """Setup discharge (src) points where a river enters the model domain. @@ -808,13 +810,21 @@ def setup_river_inflow( 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 + src_type: {'inflow', 'headwater'}, optional + Source type, by default 'inflow' + If 'inflow', return points where the river flows into the model domain. + If 'headwater', return all headwater (including inflow) points within the model domain. + mask_upstream_cells: bool, optional + If True, mask upstream cells of the river source points, by default False. + Note that this requires hydrography data and is only used if `src_type='headwater'`. See Also -------- setup_discharge_forcing setup_discharge_forcing_from_grid """ - da_flwdir, da_uparea, gdf_riv = None, None, None + # get hydrography data + da_uparea = None if hydrography is not None: ds = self.data_catalog.get_rasterdataset( hydrography, @@ -822,9 +832,14 @@ def setup_river_inflow( variables=["uparea", "flwdir"], buffer=5, ) - da_flwdir = ds["flwdir"] - da_uparea = ds["uparea"] - elif ( + da_uparea = ds["uparea"] # reused in river_source_points + elif mask_upstream_cells: + raise ValueError( + "Masking of upstream cells requires hydrography data to be provided." + ) + + # get river centerlines + if ( isinstance(rivers, str) and rivers == "rivers_outflow" and rivers in self.geoms @@ -835,26 +850,54 @@ def setup_river_inflow( gdf_riv = self.data_catalog.get_geodataframe( rivers, geom=self.region ).to_crs(self.crs) - else: + elif hydrography is not None: + gdf_riv = workflows.river_centerline_from_hydrography( + da_flwdir = ds["flwdir"], + da_uparea = da_uparea, + river_upa=river_upa, + river_len=river_len, + gdf_mask=self.region + ) + elif hydrography is None: raise ValueError("Either hydrography or rivers must be provided.") - gdf_src, gdf_riv = workflows.river_boundary_points( - region=self.region, - res=self.reggrid.dx, + # 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, - da_flwdir=da_flwdir, - da_uparea=da_uparea, - river_len=river_len, + gdf_mask=self.region, + src_type=src_type, + buffer=buffer, river_upa=river_upa, - inflow=True, + river_len=river_len, + da_uparea=da_uparea, reverse_river_geom=reverse_river_geom, logger=self.logger, ) - n = len(gdf_src.index) - self.logger.info(f"Found {n} river inflow points.") - if n == 0: + if gdf_src.empty: return + # mask upstream cells + # FIXME basin mask should be based on shifted outlet cells + if mask_upstream_cells and src_type == "headwater": + gdf_bas, _ = workflows.basin_mask( + da_flwdir=ds["flwdir"], + gdf_outlet=gdf_src, + ) + self.setup_mask_active( + exclude_mask=gdf_bas, + reset_mask=False, + ) + # make sure that the river source points are not masked ! + self.setup_mask_active( + include_mask=gdf_src, + reset_mask=False, + ) + # set forcing src pnts gdf_src.index = gdf_src.index + first_index self.set_forcing_1d(gdf_locs=gdf_src.copy(), name="dis", merge=merge) @@ -943,7 +986,8 @@ def setup_river_outflow( -------- setup_mask_bounds """ - da_flwdir, da_uparea, gdf_riv = None, None, None + # get hydrography data + da_uparea = None if hydrography is not None: ds = self.data_catalog.get_rasterdataset( hydrography, @@ -951,9 +995,10 @@ def setup_river_outflow( variables=["uparea", "flwdir"], buffer=5, ) - da_flwdir = ds["flwdir"] - da_uparea = ds["uparea"] - elif ( + da_uparea = ds["uparea"] # reused in river_source_points + + # get river centerlines + if ( isinstance(rivers, str) and rivers == "rivers_inflow" and rivers in self.geoms @@ -964,22 +1009,36 @@ def setup_river_outflow( gdf_riv = self.data_catalog.get_geodataframe( rivers, geom=self.region ).to_crs(self.crs) + elif hydrography is not None: + gdf_riv = workflows.river_centerline_from_hydrography( + da_flwdir = ds["flwdir"], + da_uparea = da_uparea, + river_upa=river_upa, + river_len=river_len, + gdf_mask=self.region + ) else: raise ValueError("Either hydrography or rivers must be provided.") - # TODO reproject region and gdf_riv to utm zone if model crs is geographic - gdf_out, gdf_riv = workflows.river_boundary_points( - region=self.region, - res=self.reggrid.dx, + # 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_out = workflows.river_source_points( gdf_riv=gdf_riv, - da_flwdir=da_flwdir, - da_uparea=da_uparea, - river_len=river_len, + gdf_mask=self.region, + src_type="outflow", + buffer=buffer, river_upa=river_upa, - inflow=False, + river_len=river_len, + da_uparea=da_uparea, reverse_river_geom=reverse_river_geom, logger=self.logger, ) + if gdf_out.empty: + return if len(gdf_out) > 0: if "rivwth" in gdf_out.columns: diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py index 0756d8f2..3b3a7b61 100644 --- a/hydromt_sfincs/workflows/flwdir.py +++ b/hydromt_sfincs/workflows/flwdir.py @@ -11,8 +11,9 @@ logger = logging.getLogger(__name__) __all__ = [ - "river_boundary_points", + "river_source_points", "river_centerline_from_hydrography", + "basin_mask", ] @@ -21,7 +22,7 @@ def river_centerline_from_hydrography( da_uparea: xr.DataArray, river_upa: float = 10, river_len: float = 1e3, - mask: gpd.GeoDataFrame = None, + gdf_mask: gpd.GeoDataFrame = None, ) -> gpd.GeoDataFrame: """Returns the centerline of rivers based on a flow direction raster data (`da_flwdir`). @@ -37,7 +38,7 @@ def river_centerline_from_hydrography( river_len: float, optional Mimimum river length [m] within the model domain to define river cells, by default 1000 m. - mask: geopandas.GeoDataFrame, optional + gdf_mask: geopandas.GeoDataFrame, optional Polygon to clip river center lines before calculating the river length, by default None. @@ -47,108 +48,183 @@ def river_centerline_from_hydrography( River line vector data. """ # get river network from hydrography based on upstream area mask - flwdir = hydromt.flw.flwdir_from_da(da_flwdir, mask=da_uparea >= river_upa) - gdf_riv = gpd.GeoDataFrame.from_features(flwdir.streams(), crs=da_flwdir.raster.crs) + riv_mask=da_uparea >= river_upa + if not riv_mask.any(): + return gpd.GeoDataFrame() + flwdir = hydromt.flw.flwdir_from_da(da_flwdir, mask=riv_mask) + feats = flwdir.streams(uparea=da_uparea.values) + gdf_riv = gpd.GeoDataFrame.from_features(feats, crs=da_flwdir.raster.crs) # clip to mask and remove empty geometries - if mask is not None: - gdf_riv = gdf_riv.to_crs(mask.crs).clip(mask.unary_union) + if gdf_mask is not None and isinstance(gdf_mask, gpd.GeoDataFrame): + gdf_riv = gdf_riv.to_crs(gdf_mask.crs).clip(gdf_mask.unary_union) gdf_riv = gdf_riv[~gdf_riv.is_empty] # create river network from gdf to get distance from outlet 'rivlen' - gdf_riv["rivlen"] = gdf_riv.geometry.length + # length of river segments + if gdf_riv.crs.is_geographic: + gdf_riv["seglen"] = gdf_riv.to_crs("epsg:3857").geometry.length + else: + gdf_riv["seglen"] = gdf_riv.geometry.length + gdf_riv = gdf_riv[gdf_riv["seglen"] > 0] + if gdf_riv.empty or river_len == 0: + return gdf_riv + # accumulate to get river length from outlet flwdir = pyflwdir.from_dataframe(gdf_riv.set_index("idx"), ds_col="idx_ds") - gdf_riv["rivlen"] = flwdir.accuflux(gdf_riv["rivlen"].values, direction="down") + gdf_riv["rivdst"] = flwdir.accuflux(gdf_riv["seglen"].values, direction="down") + # get maximum river length from outlet (at headwater segments) for each river segment + gdf_riv["rivlen"] = flwdir.fillnodata(np.where(flwdir.n_upstream == 0, gdf_riv["rivdst"], 0), 0) + # filter river network based on total length gdf_riv = gdf_riv[gdf_riv["rivlen"] >= river_len] return gdf_riv -def river_boundary_points( - region: gpd.GeoDataFrame, - res: float, - gdf_riv: gpd.GeoDataFrame = None, - da_flwdir: xr.DataArray = None, - da_uparea: xr.DataArray = None, +def river_source_points( + gdf_riv: gpd.GeoDataFrame, + gdf_mask: gpd.GeoDataFrame, + src_type: str = 'inflow', + buffer: float = 100, river_upa: float = 10, river_len: float = 1e3, - inflow: bool = True, + da_uparea: xr.DataArray = None, 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. + or out (`inflow=False`) of the model gdf_mask. Rivers are based on either a river network vector data (`gdf_riv`) or a flow direction raster data (`da_flwdir`). Parameters ---------- - region: geopandas.GeoDataFrame - Polygon of model region of interest. - res: float - Model resolution [m]. - gdf_riv: geopandas.GeoDataFrame, optional + gdf_riv: geopandas.GeoDataFrame River network vector data, by default None. - da_flwdir: xarray.DataArray, optional - D8 flow direction raster data, by default None. - da_uparea: xarray.DataArray, optional - River upstream area raster data, by default None. + Requires 'uparea' and 'rivlen' attributes to + check for river length and upstream area thresholds. + gdf_mask: geopandas.GeoDataFrame + Polygon of model gdf_mask of interest. + src_type: ['inflow', 'outflow', 'headwater'], optional + Type of river source points to return, by default 'inflow'. + If 'inflow', return points where the river flows into the model domain. + 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. river_upa : float, optional Minimum upstream area threshold for rivers [km2], by default 10.0 river_len: float, optional 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. + da_uparea: xarray.DataArray, optional + River upstream area raster data, by default None. 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 ------- - gdf_src, gdf_riv: geopandas.GeoDataFrame - In-/outflow points and river line vector data. + gdf_pnt: geopandas.GeoDataFrame + Source points """ - if not isinstance(region, (gpd.GeoDataFrame, gpd.GeoSeries)) and np.all( - region.geometry.type == "Polygon" + # data checks + if not ( + isinstance(gdf_mask, (gpd.GeoDataFrame, gpd.GeoSeries)) + and np.all(np.isin(gdf_mask.geometry.type, ["Polygon", "MultiPolygon"])) ): - 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.") - elif gdf_riv is None: # get river network from hydrography - gdf_riv = river_centerline_from_hydrography( - da_flwdir, da_uparea, river_upa, river_len, mask=region - ) - else: - # clip river to model region - 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[gdf_riv["uparea"] >= river_upa] - if "rivlen" in gdf_riv.columns: - gdf_riv = gdf_riv[gdf_riv["rivlen"] > river_len] + raise TypeError("gdf_mask must be a GeoDataFrame of Polygons.") + if not ( + isinstance(gdf_riv, (gpd.GeoDataFrame, gpd.GeoSeries)) + and np.all(np.isin(gdf_riv.geometry.type, ["LineString", "MultiLineString"])) + ): + raise TypeError("gdf_riv must be a GeoDataFrame of LineStrings.") + if src_type not in ['inflow', 'outflow', 'headwater']: + raise ValueError("src_type must be either 'inflow', 'outflow', or 'headwater'.") + if gdf_mask.crs.is_geographic: # to pseudo mercator + gdf_mask = gdf_mask.to_crs("epsg:3857") + + # clip river to model gdf_mask + gdf_riv = gdf_riv.to_crs(gdf_mask.crs).clip(gdf_mask.unary_union) + # filter river network based on uparea and length + if "uparea" in gdf_riv.columns: + gdf_riv = gdf_riv[gdf_riv["uparea"] >= river_upa] + if "rivlen" in gdf_riv.columns: + gdf_riv = gdf_riv[gdf_riv["rivlen"] > river_len] + if gdf_riv.empty: + logger.warning("No rivers matching the uparea and rivlen thresholds found in gdf_riv.") + return gpd.GeoDataFrame() + + # 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 = 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 - bnd = region.boundary.buffer(res).unary_union # NOTE should be single geom - gdf_pnt = gdf_pnt[gdf_pnt.within(bnd)].reset_index(drop=True) + 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 + 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) + elif src_type == 'outflow': + pnts_up = gdf_up.buffer(5).unary_union + gdf_pnt = gdf_ds[~gdf_ds.intersects(pnts_up)].reset_index(drop=True) + + # 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 + logger.info(f"Found {gdf_pnt.index.size} {src_type} points.") # add uparea attribute if da_uparea is provided if da_uparea is not None: 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_pnt = hydromt.gis_utils.nearest_merge( - gdf_pnt, gdf_riv, columns=["rivwth"], max_dist=10 - ) - return gdf_pnt, gdf_riv + return gdf_pnt + +# NOTE: this function is should be available in hydromt +def basin_mask( + da_flwdir: xr.DataArray, + gdf_outlet: gpd.GeoDataFrame, + shift_nup: int = 0, +) -> tuple[gpd.GeoDataFrame, xr.DataArray]: + """Returns a basin mask of all cells upstream from the outlets. + + Parameters + ---------- + da_flwdir: xarray.DataArray + D8 flow direction raster data. + gdf_outlet: geopandas.GeoDataFrame + GeoDataFrame of outlet points. + + Returns + ------- + gdf_bas: geopandas.GeoDataFrame + GeoDataFrame of basin polygons. + da_bas: xarray.DataArray + Raster of basin polygons. + """ + if not ( + isinstance(gdf_outlet, (gpd.GeoDataFrame, gpd.GeoSeries)) + and np.all(gdf_outlet.geometry.type == "Point") + ): + raise TypeError("gdf_mask must be a GeoDataFrame of Points.") + flwdir = hydromt.flw.flwdir_from_da(da_flwdir) + if da_flwdir.raster.crs != gdf_outlet.crs: + gdf_outlet = gdf_outlet.to_crs(da_flwdir.raster.crs) + xy= gdf_outlet.geometry.x.values, gdf_outlet.geometry.y.values + # FIXME + # if shift_nup > 0: + # xy = flwdir.main_upstream(xy, n=shift_nup) + da_bas = xr.DataArray( + data=flwdir.basins(xy=xy).astype(np.int32), + dims=da_flwdir.raster.dims, + coords=da_flwdir.raster.coords, + ) + da_bas.raster.set_nodata(0) + da_bas.raster.set_crs(da_flwdir.raster.crs) + gdf_bas = da_bas.raster.vectorize() + gdf_bas['idx'] = gdf_bas.index + return gdf_bas, da_bas \ No newline at end of file diff --git a/tests/conftest.py b/tests/conftest.py index 247d693b..13d91a01 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -5,6 +5,7 @@ import pytest from hydromt_sfincs.sfincs import SfincsModel +from hydromt import DataCatalog TESTDATADIR = join(dirname(abspath(__file__)), "data") TESTMODELDIR = join(TESTDATADIR, "sfincs_test") @@ -37,3 +38,17 @@ def mod(tmpdir): mod.read() mod.set_root(str(tmpdir), mode="r+") return mod + +@pytest.fixture +def data_catalog(): + return DataCatalog('artifact_data') + +@pytest.fixture +def hydrography(data_catalog): + bbox = [12.64, 45.48, 12.82, 45.59] + ds_hydro = data_catalog.get_rasterdataset("merit_hydro", variables=['flwdir', 'uparea', 'basins'], bbox=bbox).load() + da_mask = (ds_hydro['basins'] == 210000039).astype(int) + da_mask.raster.set_crs(ds_hydro.raster.crs) + da_mask.raster.set_nodata(0) + gdf_mask = da_mask.raster.vectorize() + return ds_hydro['flwdir'], ds_hydro['uparea'], gdf_mask diff --git a/tests/test_flwdir.py b/tests/test_flwdir.py new file mode 100644 index 00000000..77255547 --- /dev/null +++ b/tests/test_flwdir.py @@ -0,0 +1,78 @@ +import pytest +from hydromt_sfincs.workflows.flwdir import basin_mask, river_centerline_from_hydrography, river_source_points +import geopandas as gpd +import numpy as np + +def test_river_centerline_from_hydrography(hydrography): + # get data + da_flwdir, da_uparea, gdf_mask = hydrography + # get river centerlines + gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, gdf_mask=gdf_mask) + # check + assert isinstance(gdf_riv, gpd.GeoDataFrame) + assert np.isin(['geometry', 'rivlen', 'uparea'], gdf_riv.columns).all() + assert np.isclose(gdf_riv['rivlen'].max(), 19665.50) + assert gdf_riv.index.size == 11 + # no rivers (uparea threshold too high) + gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, river_upa=1e6) + assert gdf_riv.empty + # no rivers (river length threshold too high) + gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, river_len=1e6) + assert gdf_riv.empty + +def test_river_source_points(hydrography, data_catalog): + # get data + da_flwdir, da_uparea, gdf_mask = hydrography + gdf_mask=gdf_mask.to_crs('EPSG:32633') + + # test with derived centerline + gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, gdf_mask=gdf_mask) + kwargs = dict(gdf_riv=gdf_riv, gdf_mask=gdf_mask) + gdf_src = river_source_points(src_type='inflow', **kwargs) + assert gdf_src.index.size == 3 + gdf_src = river_source_points(src_type='headwater', **kwargs) + 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() + + # 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 + 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 + 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"): + gdf_src = river_source_points(gdf_riv=gdf_riv, gdf_mask=gdf_mask, src_type='wrong') + with pytest.raises(TypeError, match="gdf_mask must be"): + gdf_src = river_source_points(gdf_riv=gdf_riv, gdf_mask=gdf_riv) + with pytest.raises(TypeError, match="gdf_riv must be"): + gdf_src = river_source_points(gdf_riv=gdf_mask, gdf_mask=gdf_mask) + + +def test_basin_mask(hydrography): + # get data and create outlet point + da_flwdir = hydrography[0] + points = gpd.points_from_xy([12.72375], [45.53625]) + gdf_outlet = gpd.GeoDataFrame(geometry=points, crs='EPSG:4326') + gdf_outlet = gdf_outlet.to_crs('EPSG:32633') + # test basin mask + gdf_bas, da_bas = basin_mask(da_flwdir, gdf_outlet) + assert gdf_bas.index.size == 1 + assert (gdf_bas['idx'].values == gdf_outlet.index.values).all() + assert gdf_bas.crs == da_flwdir.raster.crs + assert np.isin(np.unique(da_bas.values), [0, 1]).all() + # + gdf_outlet.to_file('outlet.geojson', driver='GeoJSON') + gdf_bas.to_file('basin.geojson', driver='GeoJSON') + # test errors + with pytest.raises(TypeError): + basin_mask(da_flwdir, gdf_outlet.buffer(1).to_frame("geometry")) \ No newline at end of file From 2ff08dd00ebc16ea7ade1599ddbcfd730453ac89 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Wed, 3 Jan 2024 09:19:12 +0100 Subject: [PATCH 2/6] update tests --- tests/conftest.py | 3 ++- .../data/sfincs_test/gis/rivers_inflow.geojson | 17 +++++++++-------- tests/data/sfincs_test/gis/src.geojson | 6 +++--- tests/data/sfincs_test/sfincs.src | 6 +++--- 4 files changed, 17 insertions(+), 15 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index 13d91a01..e25fc816 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,6 +1,7 @@ """add global fixtures""" from os.path import abspath, dirname, join +import numpy as np import pytest @@ -47,7 +48,7 @@ def data_catalog(): def hydrography(data_catalog): bbox = [12.64, 45.48, 12.82, 45.59] ds_hydro = data_catalog.get_rasterdataset("merit_hydro", variables=['flwdir', 'uparea', 'basins'], bbox=bbox).load() - da_mask = (ds_hydro['basins'] == 210000039).astype(int) + da_mask = (ds_hydro['basins'] == 210000039).astype(np.uint8) da_mask.raster.set_crs(ds_hydro.raster.crs) da_mask.raster.set_nodata(0) gdf_mask = da_mask.raster.vectorize() diff --git a/tests/data/sfincs_test/gis/rivers_inflow.geojson b/tests/data/sfincs_test/gis/rivers_inflow.geojson index 75169010..5221ab44 100644 --- a/tests/data/sfincs_test/gis/rivers_inflow.geojson +++ b/tests/data/sfincs_test/gis/rivers_inflow.geojson @@ -2,13 +2,14 @@ "type": "FeatureCollection", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32633" } }, "features": [ -{ "type": "Feature", "properties": { "idx": 16152, "idx_ds": 14884, "pit": false, "rivlen": 6648.7090264593971 }, "geometry": { "type": "LineString", "coordinates": [ [ 317317.852918574120849, 5044630.720567373558879 ], [ 317320.552235712530091, 5044723.300652908161283 ], [ 317323.251591718173586, 5044815.880749752745032 ], [ 317260.885173859074712, 5044910.358359371311963 ], [ 317263.585570152092259, 5045002.938477858901024 ], [ 317266.286005324334837, 5045095.518607657402754 ], [ 317268.986479375278577, 5045188.098748766817153 ], [ 317334.049404400866479, 5045186.201250261627138 ], [ 317399.112330119940452, 5045184.304428238421679 ], [ 317464.175256532384083, 5045182.408282697200775 ], [ 317529.238183637848124, 5045180.512813637033105 ], [ 317594.301111436099745, 5045178.618021056056023 ], [ 317662.058777295635082, 5045269.30406329780817 ], [ 317727.120743889186997, 5045267.410624660551548 ], [ 317792.182711175875738, 5045265.517862500622869 ], [ 317857.244679155293852, 5045263.625776817090809 ], [ 317922.306647827266715, 5045261.73436760995537 ], [ 317987.368617191677913, 5045259.843634877353907 ], [ 318052.430587247945368, 5045257.953578618355095 ], [ 318117.492557995952666, 5045256.064198832958937 ], [ 318182.554529435408767, 5045254.175495520234108 ], [ 318247.616501565964427, 5045252.287468680180609 ], [ 318312.678474387503229, 5045250.400118310004473 ], [ 318377.740447899850551, 5045248.513444410637021 ], [ 318442.802422102424316, 5045246.627446980215609 ], [ 318507.864396995282732, 5045244.742126018740237 ], [ 318572.92637257790193, 5045242.857481524348259 ], [ 318637.988348849990871, 5045240.973513497970998 ], [ 318703.050325811491348, 5045239.090221939608455 ], [ 318768.112303461995907, 5045237.207606843672693 ], [ 318833.174281801097095, 5045235.325668213889003 ], [ 318898.236260828853119, 5045233.44440604839474 ], [ 318963.298240544856526, 5045231.563820346258581 ], [ 319028.360220948699862, 5045229.683911105617881 ], [ 319093.422202040324919, 5045227.804678327403963 ], [ 319158.484183819324244, 5045225.926122009754181 ], [ 319223.546166285465006, 5045224.048242151737213 ], [ 319288.608149438456167, 5045222.171038753353059 ], [ 319353.670133278064895, 5045220.294511813670397 ], [ 319416.063369960756972, 5045125.83847651258111 ], [ 319483.794103015912697, 5045216.543487306684256 ], [ 319548.8560889136279, 5045214.668989739380777 ], [ 319613.918075496912934, 5045212.795168625190854 ], [ 319678.980062765418552, 5045210.922023965977132 ], [ 319744.042050718911923, 5045209.049555762670934 ], [ 319809.104039357160218, 5045207.177764010615647 ], [ 319874.166028679697774, 5045205.306648711673915 ], [ 319939.228018686524592, 5045203.436209863983095 ], [ 320004.290009377291426, 5045201.566447466611862 ], [ 320069.352000751474407, 5045199.69736152049154 ], [ 320131.755833059898578, 5045105.248756496235728 ], [ 320196.818788348115049, 5045103.381022477522492 ], [ 320261.88174431794323, 5045101.513964904472232 ], [ 320326.944700969266705, 5045099.647583779878914 ], [ 320392.007658301736228, 5045097.781879100948572 ], [ 320457.070616315118968, 5045095.916850868612528 ], [ 320522.13357500906568, 5045094.052499081008136 ], [ 320587.196534383459948, 5045092.188823738135397 ], [ 320652.259494437894318, 5045090.325824838131666 ], [ 320717.322455172252376, 5045088.463502380996943 ], [ 320782.385416586068459, 5045086.601856365799904 ], [ 320847.448378679226153, 5045084.740886791609228 ], [ 320912.511341451434419, 5045082.880593657493591 ], [ 320977.574304902402218, 5045081.020976963452995 ], [ 321042.637269031780306, 5045079.162036708556116 ], [ 321107.700233839452267, 5045077.303772890940309 ], [ 321172.763199325068854, 5045075.446185510605574 ], [ 321237.826165488222614, 5045073.589274567551911 ], [ 321302.889132328913547, 5045071.733040058985353 ], [ 321367.952099846675992, 5045069.877481986768544 ], [ 321433.015068041277118, 5045068.022600349038839 ], [ 321498.07803691260051, 5045066.168395143002272 ], [ 321563.141006460238714, 5045064.314866369590163 ], [ 321628.203976683900692, 5045062.462014029733837 ], [ 321693.26694758341182, 5045060.609838119708002 ], [ 321758.329919158364646, 5045058.758338641375303 ], [ 321823.392891408409923, 5045056.907515591010451 ], [ 321885.823730876029003, 5044962.477158785797656 ], [ 321950.887666980037466, 5044960.62768763396889 ], [ 322015.951603757217526, 5044958.77889291010797 ], [ 322081.015541207627393, 5044956.930774612352252 ], [ 322146.079479330859613, 5044955.08333274256438 ], [ 322213.770739075145684, 5045045.816782272420824 ], [ 322278.833716044609901, 5045043.970694209448993 ] ] } }, -{ "type": "Feature", "properties": { "idx": 10859, "idx_ds": 14884, "pit": false, "rivlen": 2839.9998192306093 }, "geometry": { "type": "LineString", "coordinates": [ [ 322781.434200728777796, 5046697.512981650419533 ], [ 322713.767283682827838, 5046606.773889691568911 ], [ 322646.098479583044536, 5046516.035487552173436 ], [ 322581.050901019130833, 5046517.878178010694683 ], [ 322578.427788468194194, 5046425.297775230370462 ], [ 322575.804713668418117, 5046332.717383888550103 ], [ 322573.181676620384678, 5046240.137003984302282 ], [ 322570.55867732455954, 5046147.556635519489646 ], [ 322632.988107778073754, 5046053.133583151735365 ], [ 322630.366146633517928, 5045960.553236592561007 ], [ 322692.798541183350608, 5045866.130880603566766 ], [ 322690.177618135814555, 5045773.55055595561862 ], [ 322752.612976667878684, 5045679.128896361216903 ], [ 322749.993091663403902, 5045586.548593632876873 ], [ 322747.373244373477064, 5045493.968302347697318 ], [ 322679.694303225143813, 5045403.229371787980199 ], [ 322612.013475374842528, 5045312.491131012327969 ], [ 322544.330760861688759, 5045221.753580015152693 ], [ 322479.269706231192686, 5045223.596960472874343 ], [ 322414.208652275206987, 5045225.441017352975905 ], [ 322346.52212744060671, 5045134.70551089476794 ], [ 322278.833716044609901, 5045043.970694209448993 ] ] } }, -{ "type": "Feature", "properties": { "idx": 896, "idx_ds": 10859, "pit": false, "rivlen": 4779.930632428006 }, "geometry": { "type": "LineString", "coordinates": [ [ 321482.641141249216162, 5047503.86443398706615 ], [ 321501.639048813842237, 5047475.111896035261452 ], [ 321566.676987534854561, 5047473.257716386578977 ], [ 321629.076994014380034, 5047378.823710976168513 ], [ 321626.439099021023139, 5047286.243220208212733 ], [ 321688.842070675455034, 5047191.809911049902439 ], [ 321753.882900070631877, 5047189.957757667638361 ], [ 321816.28879955841694, 5047095.525809776969254 ], [ 321881.330593178281561, 5047093.675008241087198 ], [ 321946.372387497220188, 5047091.824883125722408 ], [ 322008.782178246765397, 5046997.394971938803792 ], [ 322073.82493676955346, 5046995.546198667958379 ], [ 322138.86769598943647, 5046993.698101814836264 ], [ 322203.910455906065181, 5046991.850681381300092 ], [ 322268.953216519148555, 5046990.00393736269325 ], [ 322333.995977828453761, 5046988.157869760878384 ], [ 322399.038739833747968, 5046986.312478574924171 ], [ 322464.0815025344491, 5046984.467763802036643 ], [ 322529.124265930848196, 5046982.623725444078445 ], [ 322594.167030022246763, 5046980.780363499186933 ], [ 322659.209794808411971, 5046978.937677966430783 ], [ 322724.252560289227404, 5046977.095668844878674 ], [ 322789.295326464285608, 5046975.254336133599281 ], [ 322851.718643160886131, 5046880.833215908147395 ], [ 322849.09923068160424, 5046788.25276342779398 ], [ 322781.434200728777796, 5046697.512981650419533 ] ] } }, -{ "type": "Feature", "properties": { "idx": 8412, "idx_ds": 10859, "pit": false, "rivlen": 4801.4498644109563 }, "geometry": { "type": "LineString", "coordinates": [ [ 323915.856173059903085, 5047684.699420366436243 ], [ 323848.215789249981754, 5047593.948675882071257 ], [ 323783.179751736810431, 5047595.779179214499891 ], [ 323715.536518208216876, 5047505.029801967553794 ], [ 323650.49951922748005, 5047506.861659084446728 ], [ 323582.853436042438261, 5047416.113649071194232 ], [ 323517.815475604846142, 5047417.946859977208078 ], [ 323450.166542825230863, 5047327.200217196717858 ], [ 323382.515722643816844, 5047236.454264258034527 ], [ 323314.863015099545009, 5047145.709001156501472 ], [ 323247.208420231530908, 5047054.96442788746208 ], [ 323244.594709105964284, 5046962.383946596644819 ], [ 323176.937301789526828, 5046871.640075599774718 ], [ 323109.27800725447014, 5046780.896894428879023 ], [ 323106.662483467254788, 5046688.316449463367462 ], [ 323039.000376657466404, 5046597.573970547877252 ], [ 322971.336382734123617, 5046506.832181451842189 ], [ 322973.953756684320979, 5046599.412601556628942 ], [ 322911.525511753861792, 5046693.83233954757452 ], [ 322846.479855896090157, 5046695.672322393395007 ], [ 322781.434200728777796, 5046697.512981650419533 ] ] } }, -{ "type": "Feature", "properties": { "idx": 7081, "idx_ds": 8412, "pit": false, "rivlen": 5935.9603288002381 }, "geometry": { "type": "LineString", "coordinates": [ [ 324776.863968961988576, 5048216.456981705501676 ], [ 324709.241647010203451, 5048125.69735639449209 ], [ 324641.617437186010648, 5048034.938420946709812 ], [ 324576.586205770261586, 5048036.760802630335093 ], [ 324508.959145894623362, 5047946.003234423696995 ], [ 324443.926952959271148, 5047947.826969872228801 ], [ 324376.297043094295077, 5047857.070768904872239 ], [ 324311.263888650108129, 5047858.895858122967184 ], [ 324246.230734904180281, 5047860.721623728983104 ], [ 324181.19758185709361, 5047862.548065725713968 ], [ 324113.56289767078124, 5047771.794586765579879 ], [ 324048.528783124231268, 5047773.622382535599172 ], [ 323980.891249095089734, 5047682.870270816609263 ], [ 323915.856173059903085, 5047684.699420366436243 ] ] } }, -{ "type": "Feature", "properties": { "idx": 7059, "idx_ds": 8412, "pit": false, "rivlen": 6087.3256808398864 }, "geometry": { "type": "LineString", "coordinates": [ [ 323346.219378565321676, 5048256.675819529220462 ], [ 323408.636396655580029, 5048162.259956204332411 ], [ 323471.055378940887749, 5048067.844778727740049 ], [ 323533.476325376657769, 5047973.430287101306021 ], [ 323601.117818619997706, 5048064.177694470621645 ], [ 323666.149039515061304, 5048062.345166937448084 ], [ 323731.180261113564484, 5048060.513315798714757 ], [ 323728.572876859572716, 5047967.93270149640739 ], [ 323796.211483415099792, 5048058.682141054421663 ], [ 323858.637248025392182, 5047964.271026403643191 ], [ 323856.031827110797167, 5047871.69042157009244 ], [ 323788.392331925220788, 5047780.940329549834132 ], [ 323853.426443676406052, 5047779.109828205779195 ], [ 323915.856173059903085, 5047684.699420366436243 ] ] } }, -{ "type": "Feature", "properties": { "idx": 4226, "idx_ds": 7081, "pit": false, "rivlen": 10509.221184722774 }, "geometry": { "type": "LineString", "coordinates": [ [ 328209.030331759888213, 5048993.231382192112505 ], [ 328181.268359407375101, 5048955.922479469329119 ], [ 328113.704410359845497, 5048865.127519024536014 ], [ 328046.138572573894635, 5048774.333248464390635 ], [ 327978.570846088579856, 5048683.539667789824307 ], [ 327916.093371216848027, 5048777.908328552730381 ], [ 327851.070771580096334, 5048779.696883112192154 ], [ 327786.04817263816949, 5048781.486114014871418 ], [ 327721.025574391474947, 5048783.27602126263082 ], [ 327656.00297684018733, 5048785.066604854539037 ], [ 327590.980379984597676, 5048786.857864793390036 ], [ 327525.957783824880607, 5048788.649801079183817 ], [ 327463.488054458110128, 5048883.023199899122119 ], [ 327398.466422700206749, 5048884.816487887874246 ], [ 327335.999620497226715, 5048979.191247972659767 ], [ 327270.978953160112724, 5048980.985887664370239 ], [ 327205.958286522538401, 5048982.781203705817461 ], [ 327138.379902650369331, 5048891.996403333730996 ], [ 327073.358274384750985, 5048893.793073072098196 ], [ 327005.777039637905546, 5048803.009639935567975 ], [ 326940.7544497550698, 5048804.807663375511765 ], [ 326873.170364196063019, 5048714.0255974708125 ], [ 326870.608904660621192, 5048621.44484331086278 ], [ 326805.584390180418268, 5048623.244221447966993 ], [ 326737.996527747309301, 5048532.463535306043923 ], [ 326672.971051682892721, 5048534.264267148450017 ], [ 326605.380338601884432, 5048443.484948233701289 ], [ 326540.353900964139029, 5048445.287033787928522 ], [ 326475.327464022557251, 5048447.089795700274408 ], [ 326410.301027777604759, 5048448.89323397539556 ], [ 326345.274592229456175, 5048450.697348609566689 ], [ 326277.678141783049796, 5048359.921428979374468 ], [ 326212.650744670594577, 5048361.726897332817316 ], [ 326147.623348254826851, 5048363.533042050898075 ], [ 326082.59595253597945, 5048365.339863134548068 ], [ 326017.568557514343411, 5048367.147360583767295 ], [ 325955.11599359055981, 5048461.536240063607693 ], [ 325890.089562928362284, 5048463.345089253969491 ], [ 325825.063132965238765, 5048465.154614810831845 ], [ 325760.036703701538499, 5048466.964816735126078 ], [ 325697.588994490448385, 5048561.356408222578466 ], [ 325632.56352960865479, 5048563.167961890809238 ], [ 325567.538065428147092, 5048564.980191928334534 ], [ 325499.930993646383286, 5048474.212388136424124 ], [ 325434.904567883873824, 5048476.025971915572882 ], [ 325369.878142822592054, 5048477.840232067741454 ], [ 325304.851718462596182, 5048479.655168594792485 ], [ 325239.825294804293662, 5048481.470781493932009 ], [ 325174.798871847859118, 5048483.287070769816637 ], [ 325172.212485807249323, 5048390.706377064809203 ], [ 325107.185100576141849, 5048392.523343715816736 ], [ 325039.569441279105376, 5048301.760306535288692 ], [ 324974.541094485786743, 5048303.578626939095557 ], [ 324906.922584951040335, 5048212.816957004368305 ], [ 324841.893276606220752, 5048214.636631164699793 ], [ 324774.271916895755567, 5048123.876328473910689 ], [ 324776.863968961988576, 5048216.456981705501676 ] ] } }, -{ "type": "Feature", "properties": { "idx": 3269, "idx_ds": 7081, "pit": false, "rivlen": 6897.1489752897605 }, "geometry": { "type": "LineString", "coordinates": [ [ 324537.761280027509201, 5049060.525895957835019 ], [ 324540.111560188117437, 5049056.971507646143436 ], [ 324537.515320308331866, 5048964.390755049884319 ], [ 324534.919117764104158, 5048871.810013936832547 ], [ 324597.346479434694629, 5048777.406234257854521 ], [ 324529.726824684301391, 5048686.648566163145006 ], [ 324592.156187070882879, 5048592.244807437993586 ], [ 324654.587513500766363, 5048497.841734572313726 ], [ 324717.020803929539397, 5048403.439347565174103 ], [ 324714.427714325021952, 5048310.858672353439033 ], [ 324776.863968961988576, 5048216.456981705501676 ] ] } } +{ "type": "Feature", "properties": { "idx": 14884, "idx_ds": 16458, "pit": false, "uparea": 4464.71875, "seglen": 771.83868994732154, "rivdst": 771.83868994732154, "rivlen": 10509.221184746086 }, "geometry": { "type": "LineString", "coordinates": [ [ 322278.833716050314251, 5045043.970694203861058 ], [ 322341.271297740284353, 5044949.545065672136843 ], [ 322338.645939589070622, 5044856.964860216714442 ], [ 322401.08648500306299, 5044762.539928111247718 ], [ 322463.528993797488511, 5044668.115681963041425 ], [ 322525.973465927876532, 5044573.692121772095561 ], [ 322588.419901350163855, 5044479.269247545860708 ], [ 322650.868300019472372, 5044384.847059288993478 ] ] } }, +{ "type": "Feature", "properties": { "idx": 16152, "idx_ds": 14884, "pit": false, "uparea": 11.159213066101074, "seglen": 5876.8703365176552, "rivdst": 6648.7090264649769, "rivlen": 6648.7090264649769 }, "geometry": { "type": "LineString", "coordinates": [ [ 317317.852918575285003, 5044630.720567367970943 ], [ 317320.55223571381066, 5044723.300652902573347 ], [ 317323.251591719337739, 5044815.880749748088419 ], [ 317260.88517386035528, 5044910.358359365724027 ], [ 317263.58557015331462, 5045002.938477852381766 ], [ 317266.286005325498991, 5045095.518607651814818 ], [ 317268.986479376559146, 5045188.09874876216054 ], [ 317334.049404402030632, 5045186.201250256970525 ], [ 317399.112330121221021, 5045184.304428233765066 ], [ 317464.175256533664651, 5045182.408282692544162 ], [ 317529.2381836391869, 5045180.512813631445169 ], [ 317594.301111437613145, 5045178.61802105139941 ], [ 317662.058777297148481, 5045269.304063292220235 ], [ 317727.120743890875019, 5045267.41062465403229 ], [ 317792.182711177505553, 5045265.517862495034933 ], [ 317857.244679156923667, 5045263.625776811502874 ], [ 317922.306647829129361, 5045261.734367604367435 ], [ 317987.368617193424143, 5045259.843634870834649 ], [ 318052.430587249866221, 5045257.953578611835837 ], [ 318117.492557997698896, 5045256.064198826439679 ], [ 318182.554529437213205, 5045254.175495514646173 ], [ 318247.616501568001695, 5045252.287468674592674 ], [ 318312.678474389598705, 5045250.400118304416537 ], [ 318377.740447902062442, 5045248.513444404117763 ], [ 318442.802422104636207, 5045246.627446973696351 ], [ 318507.864396997552831, 5045244.742126012220979 ], [ 318572.926372580230236, 5045242.857481518760324 ], [ 318637.988348852319177, 5045240.973513492383063 ], [ 318703.05032581393607, 5045239.090221933089197 ], [ 318768.112303464440629, 5045237.207606838084757 ], [ 318833.174281803774647, 5045235.325668208301067 ], [ 318898.236260831588879, 5045233.444406042806804 ], [ 318963.298240547534078, 5045231.563820339739323 ], [ 319028.360220951493829, 5045229.683911100029945 ], [ 319093.422202043177094, 5045227.804678320884705 ], [ 319158.484183822292835, 5045225.926122003234923 ], [ 319223.546166288258974, 5045224.048242146149278 ], [ 319288.60814944136655, 5045222.171038747765124 ], [ 319353.670133281033486, 5045220.294511808082461 ], [ 319416.063369963841978, 5045125.838476506993175 ], [ 319483.794103018997703, 5045216.54348730109632 ], [ 319548.856088916771114, 5045214.668989732861519 ], [ 319613.918075500172563, 5045212.795168619602919 ], [ 319678.980062768678181, 5045210.922023960389197 ], [ 319744.04205072222976, 5045209.049555757082999 ], [ 319809.104039360536262, 5045207.177764004096389 ], [ 319874.166028683190234, 5045205.30664870608598 ], [ 319939.228018690133467, 5045203.436209858395159 ], [ 320004.290009380783886, 5045201.566447461023927 ], [ 320069.352000755257905, 5045199.697361513972282 ], [ 320131.755833063740283, 5045105.248756491579115 ], [ 320196.818788351723924, 5045103.381022471934557 ], [ 320261.881744321726728, 5045101.513964899815619 ], [ 320326.944700973108411, 5045099.647583774290979 ], [ 320392.007658305694349, 5045097.781879095360637 ], [ 320457.070616319077089, 5045095.916850863955915 ], [ 320522.133575013023801, 5045094.052499075420201 ], [ 320587.196534387534484, 5045092.188823733478785 ], [ 320652.259494442143477, 5045090.325824833475053 ], [ 320717.32245517661795, 5045088.46350237634033 ], [ 320782.385416590492241, 5045086.601856360211968 ], [ 320847.448378683708142, 5045084.740886786952615 ], [ 320912.511341455916408, 5045082.880593652836978 ], [ 320977.574304906884208, 5045081.020976958796382 ], [ 321042.637269036378711, 5045079.16203670296818 ], [ 321107.70023384410888, 5045077.303772886283696 ], [ 321172.763199329725467, 5045075.446185505948961 ], [ 321237.826165492879227, 5045073.589274562895298 ], [ 321302.88913233357016, 5045071.73304005432874 ], [ 321367.952099851565436, 5045069.877481982111931 ], [ 321433.015068046166562, 5045068.022600343450904 ], [ 321498.078036917606369, 5045066.168395137414336 ], [ 321563.141006465186365, 5045064.314866364002228 ], [ 321628.203976689022966, 5045062.462014024145901 ], [ 321693.266947588417679, 5045060.609838114120066 ], [ 321758.329919163370505, 5045058.758338635787368 ], [ 321823.392891413765028, 5045056.907515586353838 ], [ 321885.823730881325901, 5044962.47715878020972 ], [ 321950.88766698539257, 5044960.627687628380954 ], [ 322015.951603762689047, 5044958.778892903588712 ], [ 322081.015541212982498, 5044956.930774606764317 ], [ 322146.07947933638934, 5044955.083332736045122 ], [ 322213.77073908073362, 5045045.816782267764211 ], [ 322278.833716050314251, 5045043.970694203861058 ] ] } }, +{ "type": "Feature", "properties": { "idx": 10859, "idx_ds": 14884, "pit": false, "uparea": 4449.14501953125, "seglen": 2068.1611292829434, "rivdst": 2839.9998192302651, "rivlen": 10509.221184746086 }, "geometry": { "type": "LineString", "coordinates": [ [ 322781.434200734831393, 5046697.512981644831598 ], [ 322713.767283688823227, 5046606.773889687843621 ], [ 322646.098479588981718, 5046516.035487548448145 ], [ 322581.0509010249516, 5046517.878178006969392 ], [ 322578.42778847401496, 5046425.297775225713849 ], [ 322575.804713674238883, 5046332.717383882962167 ], [ 322573.181676626205444, 5046240.137003980576992 ], [ 322570.558677330496721, 5046147.556635513901711 ], [ 322632.988107783952728, 5046053.13358314614743 ], [ 322630.36614663945511, 5045960.553236587904394 ], [ 322692.798541189171374, 5045866.130880598910153 ], [ 322690.177618141751736, 5045773.55055595189333 ], [ 322752.612976673990488, 5045679.128896355628967 ], [ 322749.993091669399291, 5045586.548593627288938 ], [ 322747.373244379530661, 5045493.968302343040705 ], [ 322679.694303231080994, 5045403.229371782392263 ], [ 322612.013475380837917, 5045312.491131006740034 ], [ 322544.330760867451318, 5045221.753580009564757 ], [ 322479.269706236955244, 5045223.596960467286408 ], [ 322414.208652280969545, 5045225.441017346456647 ], [ 322346.522127446136437, 5045134.705510889180005 ], [ 322278.833716050314251, 5045043.970694203861058 ] ] } }, +{ "type": "Feature", "properties": { "idx": 896, "idx_ds": 10859, "pit": false, "uparea": 4036.349609375, "seglen": 1939.9308131965008, "rivdst": 4779.9306324267654, "rivlen": 4779.9306324267654 }, "geometry": { "type": "LineString", "coordinates": [ [ 321482.64114125672495, 5047503.864433980546892 ], [ 321501.639048818673473, 5047475.111896032467484 ], [ 321566.67698753986042, 5047473.257716383785009 ], [ 321629.076994019327685, 5047378.823710975237191 ], [ 321626.439099026028998, 5047286.243220204487443 ], [ 321688.842070680577308, 5047191.809911046177149 ], [ 321753.882900075754151, 5047189.957757663913071 ], [ 321816.288799563539214, 5047095.525809773243964 ], [ 321881.330593183636665, 5047093.675008237361908 ], [ 321946.372387502575293, 5047091.824883121065795 ], [ 322008.782178252236918, 5046997.394971935078502 ], [ 322073.82493677502498, 5046995.546198663301766 ], [ 322138.86769599490799, 5046993.698101812042296 ], [ 322203.910455911653116, 5046991.850681376643479 ], [ 322268.953216524794698, 5046990.00393735896796 ], [ 322333.995977834216319, 5046988.157869757153094 ], [ 322399.038739839335904, 5046986.312478571198881 ], [ 322464.081502540386282, 5046984.467763798311353 ], [ 322529.124265936668962, 5046982.623725440353155 ], [ 322594.167030028067529, 5046980.780363495461643 ], [ 322659.209794814349152, 5046978.937677962705493 ], [ 322724.252560295164585, 5046977.095668840222061 ], [ 322789.295326470397413, 5046975.254336130805314 ], [ 322851.718643166939728, 5046880.833215903490782 ], [ 322849.099230687657837, 5046788.252763424068689 ], [ 322781.434200734831393, 5046697.512981644831598 ] ] } }, +{ "type": "Feature", "properties": { "idx": 8412, "idx_ds": 10859, "pit": false, "uparea": 407.14508056640625, "seglen": 1961.4500451814827, "rivdst": 4801.4498644117475, "rivlen": 10509.221184746086 }, "geometry": { "type": "LineString", "coordinates": [ [ 323915.85617306700442, 5047684.699420362710953 ], [ 323848.215789256966673, 5047593.948675878345966 ], [ 323783.179751743678935, 5047595.779179209843278 ], [ 323715.53651821508538, 5047505.029801963828504 ], [ 323650.499519234348554, 5047506.861659081652761 ], [ 323582.853436049132142, 5047416.113649068400264 ], [ 323517.815475611481816, 5047417.946859973482788 ], [ 323450.166542831866536, 5047327.200217193923891 ], [ 323382.515722650277894, 5047236.454264254309237 ], [ 323314.863015106122475, 5047145.709001152776182 ], [ 323247.208420238108374, 5047054.964427883736789 ], [ 323244.59470911254175, 5046962.383946592919528 ], [ 323176.937301796046086, 5046871.640075594186783 ], [ 323109.278007260872982, 5046780.896894425153732 ], [ 323106.662483473774046, 5046688.316449458710849 ], [ 323039.000376663752832, 5046597.573970545083284 ], [ 322971.336382740293629, 5046506.832181447185576 ], [ 322973.953756690607406, 5046599.412601552903652 ], [ 322911.525511759973597, 5046693.832339542917907 ], [ 322846.479855902143754, 5046695.672322388738394 ], [ 322781.434200734831393, 5046697.512981644831598 ] ] } }, +{ "type": "Feature", "properties": { "idx": 7081, "idx_ds": 8412, "pit": false, "uparea": 395.60482788085938, "seglen": 1134.5104643920597, "rivdst": 5935.960328803807, "rivlen": 10509.221184746086 }, "geometry": { "type": "LineString", "coordinates": [ [ 324776.863968969788402, 5048216.45698170363903 ], [ 324709.241647018003277, 5048125.697356392629445 ], [ 324641.617437193635851, 5048034.938420943915844 ], [ 324576.586205777886789, 5048036.760802627541125 ], [ 324508.959145902073942, 5047946.003234419971704 ], [ 324443.926952966838144, 5047947.826969868503511 ], [ 324376.297043101862073, 5047857.070768902078271 ], [ 324311.263888657558709, 5047858.895858119241893 ], [ 324246.230734911630861, 5047860.721623726189137 ], [ 324181.19758186431136, 5047862.548065721988678 ], [ 324113.56289767799899, 5047771.794586761854589 ], [ 324048.528783131507225, 5047773.622382532805204 ], [ 323980.891249102249276, 5047682.870270811952651 ], [ 323915.85617306700442, 5047684.699420362710953 ] ] } }, +{ "type": "Feature", "properties": { "idx": 7059, "idx_ds": 8412, "pit": false, "uparea": 10.091493606567383, "seglen": 1285.8758164290573, "rivdst": 6087.325680840805, "rivlen": 6087.325680840805 }, "geometry": { "type": "LineString", "coordinates": [ [ 323346.21937857195735, 5048256.675819526426494 ], [ 323408.636396662099287, 5048162.259956203401089 ], [ 323471.055378947523423, 5048067.844778724946082 ], [ 323533.476325383235235, 5047973.430287098512053 ], [ 323601.117818626749795, 5048064.177694467827678 ], [ 323666.149039521929808, 5048062.345166934654117 ], [ 323731.180261120432988, 5048060.513315795920789 ], [ 323728.57287686644122, 5047967.932701493613422 ], [ 323796.211483422084711, 5048058.682141051627696 ], [ 323858.637248032377101, 5047964.271026399917901 ], [ 323856.031827117782086, 5047871.690421567298472 ], [ 323788.392331932263914, 5047780.940329547040164 ], [ 323853.426443683449179, 5047779.109828202053905 ], [ 323915.85617306700442, 5047684.699420362710953 ] ] } }, +{ "type": "Feature", "properties": { "idx": 4226, "idx_ds": 7081, "pit": false, "uparea": 10.004246711730957, "seglen": 4573.2608559422788, "rivdst": 10509.221184746086, "rivlen": 10509.221184746086 }, "geometry": { "type": "LineString", "coordinates": [ [ 328209.030331772693899, 5048993.231382194906473 ], [ 328181.268359418259934, 5048955.922479469329119 ], [ 328113.704410370555706, 5048865.127519022673368 ], [ 328046.138572584663052, 5048774.33324846252799 ], [ 327978.570846099290065, 5048683.539667786099017 ], [ 327916.093371227441821, 5048777.908328551799059 ], [ 327851.070771590690129, 5048779.696883110329509 ], [ 327786.048172648646869, 5048781.486114013008773 ], [ 327721.025574401952326, 5048783.276021261699498 ], [ 327656.002976850490086, 5048785.066604853607714 ], [ 327590.980379994900431, 5048786.857864792458713 ], [ 327525.957783835299779, 5048788.649801078252494 ], [ 327463.488054468354676, 5048883.023199897259474 ], [ 327398.466422710451297, 5048884.816487885080278 ], [ 327335.999620507354848, 5048979.191247972659767 ], [ 327270.978953170240857, 5048980.985887663438916 ], [ 327205.958286532608327, 5048982.781203703954816 ], [ 327138.379902660264634, 5048891.996403331868351 ], [ 327073.358274394646287, 5048893.793073071166873 ], [ 327005.777039647684433, 5048803.009639932774007 ], [ 326940.754449764848687, 5048804.80766337364912 ], [ 326873.170364205667283, 5048714.025597468018532 ], [ 326870.608904670225456, 5048621.444843308068812 ], [ 326805.58439019008074, 5048623.244221446104348 ], [ 326737.996527756971773, 5048532.463535303249955 ], [ 326672.97105169238057, 5048534.264267147518694 ], [ 326605.380338611488696, 5048443.484948230907321 ], [ 326540.353900973510463, 5048445.287033784203231 ], [ 326475.327464031986892, 5048447.089795696549118 ], [ 326410.301027786976192, 5048448.89323397167027 ], [ 326345.274592238711193, 5048450.697348606772721 ], [ 326277.678141792246606, 5048359.921428975649178 ], [ 326212.65074467973318, 5048361.726897329092026 ], [ 326147.623348263965454, 5048363.533042048104107 ], [ 326082.59595254494343, 5048365.3398631317541 ], [ 326017.568557523307391, 5048367.147360580042005 ], [ 325955.11599359952379, 5048461.536240060813725 ], [ 325890.089562937093433, 5048463.3450892502442 ], [ 325825.063132974028122, 5048465.154614808037877 ], [ 325760.036703710153233, 5048466.964816733263433 ], [ 325697.588994499179535, 5048561.356408221647143 ], [ 325632.563529617269523, 5048563.167961889877915 ], [ 325567.538065436645411, 5048564.980191926471889 ], [ 325499.930993654881604, 5048474.212388133630157 ], [ 325434.90456789243035, 5048476.025971912778914 ], [ 325369.878142830915749, 5048477.840232064947486 ], [ 325304.851718470919877, 5048479.655168591067195 ], [ 325239.82529481255915, 5048481.470781490206718 ], [ 325174.798871856066398, 5048483.287070767022669 ], [ 325172.21248581551481, 5048390.706377062015235 ], [ 325107.185100584290922, 5048392.523343712091446 ], [ 325039.569441287196241, 5048301.760306532494724 ], [ 324974.541094493819401, 5048303.578626936301589 ], [ 324906.922584958956577, 5048212.816957001574337 ], [ 324841.893276614253409, 5048214.636631162837148 ], [ 324774.271916903497186, 5048123.876328472048044 ], [ 324776.863968969788402, 5048216.45698170363903 ] ] } }, +{ "type": "Feature", "properties": { "idx": 3269, "idx_ds": 7081, "pit": false, "uparea": 372.01608276367188, "seglen": 961.18864648479268, "rivdst": 6897.1489752886, "rivlen": 6897.1489752886 }, "geometry": { "type": "LineString", "coordinates": [ [ 324537.761280037520919, 5049060.525895952247083 ], [ 324540.111560195684433, 5049056.971507644280791 ], [ 324537.515320315840654, 5048964.390755048952997 ], [ 324534.919117771671154, 5048871.810013936832547 ], [ 324597.346479442261625, 5048777.406234255060554 ], [ 324529.726824691868387, 5048686.648566162213683 ], [ 324592.156187078508083, 5048592.244807437062263 ], [ 324654.587513508391567, 5048497.841734569519758 ], [ 324717.020803937222809, 5048403.439347562380135 ], [ 324714.427714332705364, 5048310.858672351576388 ], [ 324776.863968969788402, 5048216.45698170363903 ] ] } } ] } diff --git a/tests/data/sfincs_test/gis/src.geojson b/tests/data/sfincs_test/gis/src.geojson index 027c2d24..544e09cf 100644 --- a/tests/data/sfincs_test/gis/src.geojson +++ b/tests/data/sfincs_test/gis/src.geojson @@ -2,8 +2,8 @@ "type": "FeatureCollection", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32633" } }, "features": [ -{ "type": "Feature", "properties": { "index": 1, "uparea": 4039.6552734375 }, "geometry": { "type": "Point", "coordinates": [ 321499.179280363314319, 5047478.834652632474899 ] } }, -{ "type": "Feature", "properties": { "index": 2, "uparea": 372.2264404296875 }, "geometry": { "type": "Point", "coordinates": [ 324539.39005017321324, 5049031.24278330616653 ] } }, -{ "type": "Feature", "properties": { "index": 3, "uparea": 10.424988746643066 }, "geometry": { "type": "Point", "coordinates": [ 328191.121168049343396, 5048969.163522810675204 ] } } +{ "type": "Feature", "properties": { "index": 1, "riv_idx": 0, "uparea": 4039.6552734375 }, "geometry": { "type": "Point", "coordinates": [ 321483.19241256051464, 5047503.030107934959233 ] } }, +{ "type": "Feature", "properties": { "index": 2, "riv_idx": 5, "uparea": 372.2264404296875 }, "geometry": { "type": "Point", "coordinates": [ 324538.312838675745297, 5049059.691759831272066 ] } }, +{ "type": "Feature", "properties": { "index": 3, "riv_idx": 2, "uparea": 10.424988746643066 }, "geometry": { "type": "Point", "coordinates": [ 328208.433359648974147, 5048992.429120215587318 ] } } ] } diff --git a/tests/data/sfincs_test/sfincs.src b/tests/data/sfincs_test/sfincs.src index bfe2e7ea..8c3ff568 100644 --- a/tests/data/sfincs_test/sfincs.src +++ b/tests/data/sfincs_test/sfincs.src @@ -1,3 +1,3 @@ -321499.2 5047478.8 -324539.4 5049031.2 -328191.1 5048969.2 +321483.2 5047503.0 +324538.3 5049059.7 +328208.4 5048992.4 From 7b272545bc63b6b7128d2a3d59daf0cb4d868a23 Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Tue, 12 Mar 2024 10:16:45 +0100 Subject: [PATCH 3/6] remove river based masking of active cells --- hydromt_sfincs/sfincs.py | 41 ++++++++-------------------------------- 1 file changed, 8 insertions(+), 33 deletions(-) diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py index fbbaffba..c3f9bac6 100644 --- a/hydromt_sfincs/sfincs.py +++ b/hydromt_sfincs/sfincs.py @@ -786,7 +786,6 @@ def setup_river_inflow( keep_rivers_geom: bool = False, reverse_river_geom: bool = False, src_type: str = "inflow", - mask_upstream_cells: bool = False, ): """Setup discharge (src) points where a river enters the model domain. @@ -837,9 +836,6 @@ def setup_river_inflow( Source type, by default 'inflow' If 'inflow', return points where the river flows into the model domain. If 'headwater', return all headwater (including inflow) points within the model domain. - mask_upstream_cells: bool, optional - If True, mask upstream cells of the river source points, by default False. - Note that this requires hydrography data and is only used if `src_type='headwater'`. See Also -------- @@ -856,10 +852,6 @@ def setup_river_inflow( buffer=5, ) da_uparea = ds["uparea"] # reused in river_source_points - elif mask_upstream_cells: - raise ValueError( - "Masking of upstream cells requires hydrography data to be provided." - ) # get river centerlines if ( @@ -875,11 +867,11 @@ def setup_river_inflow( ).to_crs(self.crs) elif hydrography is not None: gdf_riv = workflows.river_centerline_from_hydrography( - da_flwdir = ds["flwdir"], - da_uparea = da_uparea, + da_flwdir=ds["flwdir"], + da_uparea=da_uparea, river_upa=river_upa, river_len=river_len, - gdf_mask=self.region + gdf_mask=self.region, ) elif hydrography is None: raise ValueError("Either hydrography or rivers must be provided.") @@ -889,7 +881,7 @@ def setup_river_inflow( if self.crs.is_geographic: buffer = buffer * 111111.0 - # get river inflow / headwater source points + # get river inflow / headwater source points gdf_src = workflows.river_source_points( gdf_riv=gdf_riv, gdf_mask=self.region, @@ -904,23 +896,6 @@ def setup_river_inflow( if gdf_src.empty: return - # mask upstream cells - # FIXME basin mask should be based on shifted outlet cells - if mask_upstream_cells and src_type == "headwater": - gdf_bas, _ = workflows.basin_mask( - da_flwdir=ds["flwdir"], - gdf_outlet=gdf_src, - ) - self.setup_mask_active( - exclude_mask=gdf_bas, - reset_mask=False, - ) - # make sure that the river source points are not masked ! - self.setup_mask_active( - include_mask=gdf_src, - reset_mask=False, - ) - # set forcing src pnts gdf_src.index = gdf_src.index + first_index self.set_forcing_1d(gdf_locs=gdf_src.copy(), name="dis", merge=merge) @@ -1034,11 +1009,11 @@ def setup_river_outflow( ).to_crs(self.crs) elif hydrography is not None: gdf_riv = workflows.river_centerline_from_hydrography( - da_flwdir = ds["flwdir"], - da_uparea = da_uparea, + da_flwdir=ds["flwdir"], + da_uparea=da_uparea, river_upa=river_upa, river_len=river_len, - gdf_mask=self.region + gdf_mask=self.region, ) else: raise ValueError("Either hydrography or rivers must be provided.") @@ -1048,7 +1023,7 @@ def setup_river_outflow( if self.crs.is_geographic: buffer = buffer * 111111.0 - # get river inflow / headwater source points + # get river inflow / headwater source points gdf_out = workflows.river_source_points( gdf_riv=gdf_riv, gdf_mask=self.region, From cc3d46a73f7bc1adea867afd0fd5c599a097d67f Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Tue, 12 Mar 2024 16:09:55 +0100 Subject: [PATCH 4/6] run pre-commit --- hydromt_sfincs/workflows/flwdir.py | 49 +++++++++++++---------- tests/conftest.py | 12 ++++-- tests/test_flwdir.py | 64 +++++++++++++++++------------- 3 files changed, 71 insertions(+), 54 deletions(-) diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py index 3b3a7b61..33c7d409 100644 --- a/hydromt_sfincs/workflows/flwdir.py +++ b/hydromt_sfincs/workflows/flwdir.py @@ -48,7 +48,7 @@ def river_centerline_from_hydrography( River line vector data. """ # get river network from hydrography based on upstream area mask - riv_mask=da_uparea >= river_upa + riv_mask = da_uparea >= river_upa if not riv_mask.any(): return gpd.GeoDataFrame() flwdir = hydromt.flw.flwdir_from_da(da_flwdir, mask=riv_mask) @@ -71,7 +71,9 @@ def river_centerline_from_hydrography( flwdir = pyflwdir.from_dataframe(gdf_riv.set_index("idx"), ds_col="idx_ds") gdf_riv["rivdst"] = flwdir.accuflux(gdf_riv["seglen"].values, direction="down") # get maximum river length from outlet (at headwater segments) for each river segment - gdf_riv["rivlen"] = flwdir.fillnodata(np.where(flwdir.n_upstream == 0, gdf_riv["rivdst"], 0), 0) + gdf_riv["rivlen"] = flwdir.fillnodata( + np.where(flwdir.n_upstream == 0, gdf_riv["rivdst"], 0), 0 + ) # filter river network based on total length gdf_riv = gdf_riv[gdf_riv["rivlen"] >= river_len] return gdf_riv @@ -80,7 +82,7 @@ def river_centerline_from_hydrography( def river_source_points( gdf_riv: gpd.GeoDataFrame, gdf_mask: gpd.GeoDataFrame, - src_type: str = 'inflow', + src_type: str = "inflow", buffer: float = 100, river_upa: float = 10, river_len: float = 1e3, @@ -98,7 +100,7 @@ def river_source_points( ---------- gdf_riv: geopandas.GeoDataFrame River network vector data, by default None. - Requires 'uparea' and 'rivlen' attributes to + Requires 'uparea' and 'rivlen' attributes to check for river length and upstream area thresholds. gdf_mask: geopandas.GeoDataFrame Polygon of model gdf_mask of interest. @@ -126,18 +128,18 @@ def river_source_points( """ # data checks if not ( - isinstance(gdf_mask, (gpd.GeoDataFrame, gpd.GeoSeries)) + isinstance(gdf_mask, (gpd.GeoDataFrame, gpd.GeoSeries)) and np.all(np.isin(gdf_mask.geometry.type, ["Polygon", "MultiPolygon"])) ): raise TypeError("gdf_mask must be a GeoDataFrame of Polygons.") if not ( - isinstance(gdf_riv, (gpd.GeoDataFrame, gpd.GeoSeries)) + isinstance(gdf_riv, (gpd.GeoDataFrame, gpd.GeoSeries)) and np.all(np.isin(gdf_riv.geometry.type, ["LineString", "MultiLineString"])) ): raise TypeError("gdf_riv must be a GeoDataFrame of LineStrings.") - if src_type not in ['inflow', 'outflow', 'headwater']: + if src_type not in ["inflow", "outflow", "headwater"]: raise ValueError("src_type must be either 'inflow', 'outflow', or 'headwater'.") - if gdf_mask.crs.is_geographic: # to pseudo mercator + if gdf_mask.crs.is_geographic: # to pseudo mercator gdf_mask = gdf_mask.to_crs("epsg:3857") # clip river to model gdf_mask @@ -148,7 +150,9 @@ def river_source_points( if "rivlen" in gdf_riv.columns: gdf_riv = gdf_riv[gdf_riv["rivlen"] > river_len] if gdf_riv.empty: - logger.warning("No rivers matching the uparea and rivlen thresholds found in gdf_riv.") + logger.warning( + "No rivers matching the uparea and rivlen thresholds found in gdf_riv." + ) return gpd.GeoDataFrame() # get source points 1m before the start/end of the river @@ -156,21 +160,21 @@ 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_up["riv_idx"] = gdf_riv.index gdf_ds = gdf_riv.interpolate(-dx).to_frame("geometry") - gdf_ds['riv_idx'] = gdf_riv.index - + 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 - if src_type in ['inflow', 'headwater']: + # use a small buffer of 5m around these points to account for dx and avoid issues with inprecise 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) - elif src_type == 'outflow': + elif src_type == "outflow": pnts_up = gdf_up.buffer(5).unary_union gdf_pnt = gdf_ds[~gdf_ds.intersects(pnts_up)].reset_index(drop=True) - + # get buffer around gdf_mask, in- and outflow points should be within this buffer - if src_type in ['inflow', 'outflow']: + 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) @@ -184,11 +188,12 @@ def river_source_points( return gdf_pnt + # NOTE: this function is should be available in hydromt def basin_mask( da_flwdir: xr.DataArray, gdf_outlet: gpd.GeoDataFrame, - shift_nup: int = 0, + shift_nup: int = 0, ) -> tuple[gpd.GeoDataFrame, xr.DataArray]: """Returns a basin mask of all cells upstream from the outlets. @@ -207,14 +212,14 @@ def basin_mask( Raster of basin polygons. """ if not ( - isinstance(gdf_outlet, (gpd.GeoDataFrame, gpd.GeoSeries)) + isinstance(gdf_outlet, (gpd.GeoDataFrame, gpd.GeoSeries)) and np.all(gdf_outlet.geometry.type == "Point") ): raise TypeError("gdf_mask must be a GeoDataFrame of Points.") flwdir = hydromt.flw.flwdir_from_da(da_flwdir) if da_flwdir.raster.crs != gdf_outlet.crs: gdf_outlet = gdf_outlet.to_crs(da_flwdir.raster.crs) - xy= gdf_outlet.geometry.x.values, gdf_outlet.geometry.y.values + xy = gdf_outlet.geometry.x.values, gdf_outlet.geometry.y.values # FIXME # if shift_nup > 0: # xy = flwdir.main_upstream(xy, n=shift_nup) @@ -226,5 +231,5 @@ def basin_mask( da_bas.raster.set_nodata(0) da_bas.raster.set_crs(da_flwdir.raster.crs) gdf_bas = da_bas.raster.vectorize() - gdf_bas['idx'] = gdf_bas.index - return gdf_bas, da_bas \ No newline at end of file + gdf_bas["idx"] = gdf_bas.index + return gdf_bas, da_bas diff --git a/tests/conftest.py b/tests/conftest.py index e25fc816..c6220cbf 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -40,16 +40,20 @@ def mod(tmpdir): mod.set_root(str(tmpdir), mode="r+") return mod + @pytest.fixture def data_catalog(): - return DataCatalog('artifact_data') + return DataCatalog("artifact_data") + @pytest.fixture def hydrography(data_catalog): bbox = [12.64, 45.48, 12.82, 45.59] - ds_hydro = data_catalog.get_rasterdataset("merit_hydro", variables=['flwdir', 'uparea', 'basins'], bbox=bbox).load() - da_mask = (ds_hydro['basins'] == 210000039).astype(np.uint8) + ds_hydro = data_catalog.get_rasterdataset( + "merit_hydro", variables=["flwdir", "uparea", "basins"], bbox=bbox + ).load() + da_mask = (ds_hydro["basins"] == 210000039).astype(np.uint8) da_mask.raster.set_crs(ds_hydro.raster.crs) da_mask.raster.set_nodata(0) gdf_mask = da_mask.raster.vectorize() - return ds_hydro['flwdir'], ds_hydro['uparea'], gdf_mask + return ds_hydro["flwdir"], ds_hydro["uparea"], gdf_mask diff --git a/tests/test_flwdir.py b/tests/test_flwdir.py index 77255547..d3b5bc60 100644 --- a/tests/test_flwdir.py +++ b/tests/test_flwdir.py @@ -1,8 +1,13 @@ import pytest -from hydromt_sfincs.workflows.flwdir import basin_mask, river_centerline_from_hydrography, river_source_points +from hydromt_sfincs.workflows.flwdir import ( + basin_mask, + river_centerline_from_hydrography, + river_source_points, +) import geopandas as gpd import numpy as np + def test_river_centerline_from_hydrography(hydrography): # get data da_flwdir, da_uparea, gdf_mask = hydrography @@ -10,8 +15,8 @@ def test_river_centerline_from_hydrography(hydrography): gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, gdf_mask=gdf_mask) # check assert isinstance(gdf_riv, gpd.GeoDataFrame) - assert np.isin(['geometry', 'rivlen', 'uparea'], gdf_riv.columns).all() - assert np.isclose(gdf_riv['rivlen'].max(), 19665.50) + assert np.isin(["geometry", "rivlen", "uparea"], gdf_riv.columns).all() + assert np.isclose(gdf_riv["rivlen"].max(), 19665.50) assert gdf_riv.index.size == 11 # no rivers (uparea threshold too high) gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, river_upa=1e6) @@ -19,39 +24,42 @@ def test_river_centerline_from_hydrography(hydrography): # no rivers (river length threshold too high) gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, river_len=1e6) assert gdf_riv.empty - + + def test_river_source_points(hydrography, data_catalog): # get data da_flwdir, da_uparea, gdf_mask = hydrography - gdf_mask=gdf_mask.to_crs('EPSG:32633') - + gdf_mask = gdf_mask.to_crs("EPSG:32633") + # test with derived centerline gdf_riv = river_centerline_from_hydrography(da_flwdir, da_uparea, gdf_mask=gdf_mask) kwargs = dict(gdf_riv=gdf_riv, gdf_mask=gdf_mask) - gdf_src = river_source_points(src_type='inflow', **kwargs) + gdf_src = river_source_points(src_type="inflow", **kwargs) assert gdf_src.index.size == 3 - gdf_src = river_source_points(src_type='headwater', **kwargs) + gdf_src = river_source_points(src_type="headwater", **kwargs) assert gdf_src.index.size == 6 - gdf_src = river_source_points(src_type='outflow', da_uparea=da_uparea, **kwargs) + 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", "riv_idx", "uparea"], gdf_src.columns).all() - # test reverse oriented line - gdf_riv = data_catalog.get_geodataframe('rivers_lin2019_v1') + # 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 - 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 - gdf_src = river_source_points(src_type='headwater', **kwargs) + 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 + 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 + 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(38, gdf_src["riv_idx"].values) # test errors with pytest.raises(ValueError, match="src_type must be either"): - gdf_src = river_source_points(gdf_riv=gdf_riv, gdf_mask=gdf_mask, src_type='wrong') + gdf_src = river_source_points( + gdf_riv=gdf_riv, gdf_mask=gdf_mask, src_type="wrong" + ) with pytest.raises(TypeError, match="gdf_mask must be"): gdf_src = river_source_points(gdf_riv=gdf_riv, gdf_mask=gdf_riv) with pytest.raises(TypeError, match="gdf_riv must be"): @@ -62,17 +70,17 @@ def test_basin_mask(hydrography): # get data and create outlet point da_flwdir = hydrography[0] points = gpd.points_from_xy([12.72375], [45.53625]) - gdf_outlet = gpd.GeoDataFrame(geometry=points, crs='EPSG:4326') - gdf_outlet = gdf_outlet.to_crs('EPSG:32633') + gdf_outlet = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326") + gdf_outlet = gdf_outlet.to_crs("EPSG:32633") # test basin mask gdf_bas, da_bas = basin_mask(da_flwdir, gdf_outlet) assert gdf_bas.index.size == 1 - assert (gdf_bas['idx'].values == gdf_outlet.index.values).all() + assert (gdf_bas["idx"].values == gdf_outlet.index.values).all() assert gdf_bas.crs == da_flwdir.raster.crs assert np.isin(np.unique(da_bas.values), [0, 1]).all() - # - gdf_outlet.to_file('outlet.geojson', driver='GeoJSON') - gdf_bas.to_file('basin.geojson', driver='GeoJSON') + # + gdf_outlet.to_file("outlet.geojson", driver="GeoJSON") + gdf_bas.to_file("basin.geojson", driver="GeoJSON") # test errors with pytest.raises(TypeError): - basin_mask(da_flwdir, gdf_outlet.buffer(1).to_frame("geometry")) \ No newline at end of file + basin_mask(da_flwdir, gdf_outlet.buffer(1).to_frame("geometry")) From 756c7caf9f86c248308c4c81b185a264c2491840 Mon Sep 17 00:00:00 2001 From: roeldegoede Date: Wed, 10 Apr 2024 17:01:28 +0200 Subject: [PATCH 5/6] remove basin_mask functionality since not used --- hydromt_sfincs/workflows/flwdir.py | 46 ------------------------------ tests/test_flwdir.py | 21 -------------- 2 files changed, 67 deletions(-) diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py index 33c7d409..f720660f 100644 --- a/hydromt_sfincs/workflows/flwdir.py +++ b/hydromt_sfincs/workflows/flwdir.py @@ -187,49 +187,3 @@ def river_source_points( gdf_pnt = gdf_pnt.sort_values("uparea", ascending=False).reset_index(drop=True) return gdf_pnt - - -# NOTE: this function is should be available in hydromt -def basin_mask( - da_flwdir: xr.DataArray, - gdf_outlet: gpd.GeoDataFrame, - shift_nup: int = 0, -) -> tuple[gpd.GeoDataFrame, xr.DataArray]: - """Returns a basin mask of all cells upstream from the outlets. - - Parameters - ---------- - da_flwdir: xarray.DataArray - D8 flow direction raster data. - gdf_outlet: geopandas.GeoDataFrame - GeoDataFrame of outlet points. - - Returns - ------- - gdf_bas: geopandas.GeoDataFrame - GeoDataFrame of basin polygons. - da_bas: xarray.DataArray - Raster of basin polygons. - """ - if not ( - isinstance(gdf_outlet, (gpd.GeoDataFrame, gpd.GeoSeries)) - and np.all(gdf_outlet.geometry.type == "Point") - ): - raise TypeError("gdf_mask must be a GeoDataFrame of Points.") - flwdir = hydromt.flw.flwdir_from_da(da_flwdir) - if da_flwdir.raster.crs != gdf_outlet.crs: - gdf_outlet = gdf_outlet.to_crs(da_flwdir.raster.crs) - xy = gdf_outlet.geometry.x.values, gdf_outlet.geometry.y.values - # FIXME - # if shift_nup > 0: - # xy = flwdir.main_upstream(xy, n=shift_nup) - da_bas = xr.DataArray( - data=flwdir.basins(xy=xy).astype(np.int32), - dims=da_flwdir.raster.dims, - coords=da_flwdir.raster.coords, - ) - da_bas.raster.set_nodata(0) - da_bas.raster.set_crs(da_flwdir.raster.crs) - gdf_bas = da_bas.raster.vectorize() - gdf_bas["idx"] = gdf_bas.index - return gdf_bas, da_bas diff --git a/tests/test_flwdir.py b/tests/test_flwdir.py index d3b5bc60..b6d55b1b 100644 --- a/tests/test_flwdir.py +++ b/tests/test_flwdir.py @@ -1,6 +1,5 @@ import pytest from hydromt_sfincs.workflows.flwdir import ( - basin_mask, river_centerline_from_hydrography, river_source_points, ) @@ -64,23 +63,3 @@ def test_river_source_points(hydrography, data_catalog): gdf_src = river_source_points(gdf_riv=gdf_riv, gdf_mask=gdf_riv) with pytest.raises(TypeError, match="gdf_riv must be"): gdf_src = river_source_points(gdf_riv=gdf_mask, gdf_mask=gdf_mask) - - -def test_basin_mask(hydrography): - # get data and create outlet point - da_flwdir = hydrography[0] - points = gpd.points_from_xy([12.72375], [45.53625]) - gdf_outlet = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326") - gdf_outlet = gdf_outlet.to_crs("EPSG:32633") - # test basin mask - gdf_bas, da_bas = basin_mask(da_flwdir, gdf_outlet) - assert gdf_bas.index.size == 1 - assert (gdf_bas["idx"].values == gdf_outlet.index.values).all() - assert gdf_bas.crs == da_flwdir.raster.crs - assert np.isin(np.unique(da_bas.values), [0, 1]).all() - # - gdf_outlet.to_file("outlet.geojson", driver="GeoJSON") - gdf_bas.to_file("basin.geojson", driver="GeoJSON") - # test errors - with pytest.raises(TypeError): - basin_mask(da_flwdir, gdf_outlet.buffer(1).to_frame("geometry")) From 5f82ab644840e79e0bd30d2dc9dc30482eef7741 Mon Sep 17 00:00:00 2001 From: roeldegoede Date: Wed, 10 Apr 2024 17:10:05 +0200 Subject: [PATCH 6/6] fix tests --- hydromt_sfincs/workflows/flwdir.py | 1 - 1 file changed, 1 deletion(-) diff --git a/hydromt_sfincs/workflows/flwdir.py b/hydromt_sfincs/workflows/flwdir.py index f720660f..e13e33fe 100644 --- a/hydromt_sfincs/workflows/flwdir.py +++ b/hydromt_sfincs/workflows/flwdir.py @@ -13,7 +13,6 @@ __all__ = [ "river_source_points", "river_centerline_from_hydrography", - "basin_mask", ]