Merge pull request #109 from Deltares/resolution_structures

adding argument to provide high resolut dem to determine structure height

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
 -------

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

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])

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):