diff --git a/docs/changelog.rst b/docs/changelog.rst index 57249db2..0403f85b 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -26,6 +26,7 @@ New - Added the option to use landuse/landcover data combined with a reclass table to `SfincsModel.setup_constant_infiltration`. PR #103 - Enabled to provide locations only (so no timeseries) for `SfincsModel.setup_waterlevel_forcing` and `SfincsModel.setup_discharge_forcing` PR #104 - New optional buffer argument in `SfincsModel.setup_discharge_forcing` to select gauges around boundary only. PR #104 +- New functionality within `SfincsModel.setup_structures` to use high resolution dem for weir elevation Changed ------- diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py index d38a0039..980d2a32 100644 --- a/hydromt_sfincs/sfincs.py +++ b/hydromt_sfincs/sfincs.py @@ -1320,6 +1320,8 @@ def setup_structures( self, structures: Union[str, Path, gpd.GeoDataFrame], stype: str, + dep: Union[str, Path, xr.DataArray] = None, + buffer: float = None, dz: float = None, merge: bool = True, **kwargs, @@ -1340,6 +1342,12 @@ def setup_structures( "par1" defaults to 0.6 if gdf has no "par1" column. stype : {'thd', 'weir'} Structure type. + dep : str, Path, xr.DataArray, optional + Path, data source name, or xarray raster object ('elevtn') describing the depth in an + alternative resolution which is used for sampling the weir. + buffer : float, optional + If provided, describes the distance from the centerline to the foot of the structure. + This distance is supplied to the raster.sample as the window (wdw). merge : bool, optional If True, merge with existing'stype' structures, by default True. dz: float, optional @@ -1356,20 +1364,44 @@ def setup_structures( "thd": ["name", "geometry"], "weir": ["name", "z", "par1", "geometry"], } - assert stype in cols + assert stype in cols, f"stype must be one of {list(cols.keys())}" gdf = gdf_structures[ [c for c in cols[stype] if c in gdf_structures.columns] ] # keep relevant cols structs = utils.gdf2linestring(gdf) # check if it parsed correct # sample zb values from dep file and set z = zb + dz - if stype == "weir" and dz is not None: - elv = self.grid["dep"] + if stype == "weir" and (dep is not None or dz is not None): + if dep is None or dep == "dep": + assert "dep" in self.grid, "dep layer not found" + elv = self.grid["dep"] + else: + elv = self.data_catalog.get_rasterdataset( + dep, geom=self.region, buffer=5, variables=["elevtn"] + ) + + # calculate window size from buffer + if buffer is not None: + res = abs(elv.raster.res[0]) + if elv.raster.crs.is_geographic: + res = res * 111111.0 + window_size = int(np.ceil(buffer / res)) + else: + window_size = 0 + self.logger.debug(f"Sampling elevation with window size {window_size}") + structs_out = [] for s in structs: pnts = gpd.points_from_xy(x=s["x"], y=s["y"]) - zb = elv.raster.sample(gpd.GeoDataFrame(geometry=pnts, crs=self.crs)) - s["z"] = zb.values + float(dz) + zb = elv.raster.sample( + gpd.GeoDataFrame(geometry=pnts, crs=self.crs), wdw=window_size + ) + if zb.ndim > 1: + zb = zb.max(axis=1) + + s["z"] = zb.values + if dz is not None: + s["z"] += float(dz) structs_out.append(s) gdf = utils.linestring2gdf(structs_out, crs=self.crs) # Else function if you define elevation of weir diff --git a/hydromt_sfincs/utils.py b/hydromt_sfincs/utils.py index 86a87a19..6431e14d 100644 --- a/hydromt_sfincs/utils.py +++ b/hydromt_sfincs/utils.py @@ -424,9 +424,9 @@ def gdf2linestring(gdf: gpd.GeoDataFrame) -> List[Dict]: feat = item.drop("geometry").dropna().to_dict() # check geom line = item.geometry - if line.type == "MultiLineString" and len(line.geoms) == 1: + if line.geom_type == "MultiLineString" and len(line.geoms) == 1: line = line.geoms[0] - if line.type != "LineString": + if line.geom_type != "LineString": raise ValueError("Invalid geometry type, only LineString is accepted.") xyz = tuple(zip(*line.coords[:])) feat["x"], feat["y"] = list(xyz[0]), list(xyz[1]) diff --git a/tests/test_1model_class.py b/tests/test_1model_class.py index 920e60d2..2384d101 100644 --- a/tests/test_1model_class.py +++ b/tests/test_1model_class.py @@ -136,6 +136,9 @@ def test_structs(tmpdir): assert "weirfile" in mod.config mod.write_geoms() assert isfile(join(mod.root, "sfincs.weir")) + # test with buffer + mod.setup_structures(fn_thd_gis, stype="weir", buffer=5, dep="dep", merge=False) + assert len(mod.geoms["weir"].index) == 2 def test_drainage_structures(tmpdir):