Real wflow basins

docs/changelog.rst

@@ -11,6 +11,7 @@ Unreleased
 
 Added
 -----
+- If applicable, basins geometry based on the higher resolution DEM is stored seperately under **basins_highres**
 - New function **setup_1dmodel_connection** to connect wflow to 1D river model (eg Delft3D FM 1D, HEC-RAS, etc.) `PR #210 <https://github.com/Deltares/hydromt_wflow/pull/210>`_
 - New setup method for the **KsatHorFrac** parameter **setup_ksathorfarc** to up-downscale existing ksathorfrac maps. `PR #255 <https://github.com/Deltares/hydromt_wflow/pull/255>`_
 - new function **setup_pet_forcing** to reproject existing pet data rather than computing from other meteo data. PR #257
@@ -21,6 +22,7 @@ Added
 
 Changed
 -------
+- Basins geometry (**basins**) is now based on the actual wflow model resolution basins, instead of based on the supplied DEM
 - **setup_soilmaps**: the user can now supply variable sbm soil layer thicknesses. The Brooks Corey coefficient is then computed as weighted average over the sbm layers. PR #242
 - **setup_outlets**: the IDs of wflow_subcatch are used to define the outlets IDs rather than [1:n]. PR #247
 

hydromt_wflow/wflow.py

@@ -18,6 +18,7 @@
 import pandas as pd
 import pyflwdir
 import pyproj
+import shapely
 import toml
 import xarray as xr
 from dask.diagnostics import ProgressBar
@@ -214,6 +215,21 @@ def setup_basemaps(
         # retrieve global data (lazy!)
         ds_org = self.data_catalog.get_rasterdataset(hydrography_fn)
 
+        # Check on resolution (degree vs meter) depending on ds_org res/crs
+        scale_ratio = int(np.round(res / ds_org.raster.res[0]))
+        if scale_ratio < 1:
+            raise ValueError(
+                f"The model resolution {res} should be \
+larger than the {hydrography_fn} resolution {ds_org.raster.res[0]}"
+            )
+        if ds_org.raster.crs.is_geographic:
+            if res > 1:  # 111 km
+                raise ValueError(
+                    f"The model resolution {res} should be smaller than 1 degree \
+(111km) for geographic coordinate systems. "
+                    "Make sure you provided res in degree rather than in meters."
+                )
+
         # get basin geometry and clip data
         kind, region = hydromt.workflows.parse_region(region, logger=self.logger)
         xy = None
@@ -239,25 +255,17 @@ def setup_basemaps(
         if geom is not None and geom.crs is None:
             raise ValueError("wflow region geometry has no CRS")
 
+        # Set the basins geometry
         ds_org = ds_org.raster.clip_geom(geom, align=res, buffer=10)
         ds_org.coords["mask"] = ds_org.raster.geometry_mask(geom)
         self.logger.debug("Adding basins vector to geoms.")
-        self.set_geoms(geom, name="basins")
 
-        # Check on resolution (degree vs meter) depending on ds_org res/crs
-        scale_ratio = int(np.round(res / ds_org.raster.res[0]))
-        if scale_ratio < 1:
-            raise ValueError(
-                f"The model resolution {res} should be \
-larger than the {hydrography_fn} resolution {ds_org.raster.res[0]}"
-            )
-        if ds_org.raster.crs.is_geographic:
-            if res > 1:  # 111 km
-                raise ValueError(
-                    f"The model resolution {res} should be smaller than 1 degree \
-(111km) for geographic coordinate systems. "
-                    "Make sure you provided res in degree rather than in meters."
-                )
+        # Set name based on scale_factor
+        basins_suffix = ""
+        if scale_ratio != 1:
+            basins_suffix = "_highres"
+        self.set_geoms(geom, name=f"basins{basins_suffix}")
+
         # setup hydrography maps and set staticmap attribute with renamed maps
         ds_base, _ = workflows.hydrography(
             ds=ds_org,
@@ -285,6 +293,8 @@ def setup_basemaps(
         # Rename and add to grid
         rmdict = {k: v for k, v in self._MAPS.items() if k in ds_base.data_vars}
         self.set_grid(ds_base.rename(rmdict))
+        # Call basins once to set it
+        self.basins
 
         # setup topography maps
         ds_topo = workflows.topography(
@@ -1251,7 +1261,7 @@ def setup_gauges(
             kwargs.update(crs=code)
             gdf_gauges = self.data_catalog.get_geodataframe(
                 gauges_fn,
-                geom=self.basins,
+                geom=self.basins_highres,
                 assert_gtype="Point",
                 handle_nodata=NoDataStrategy.IGNORE,
                 **kwargs,
@@ -1260,15 +1270,15 @@ def setup_gauges(
             if self.data_catalog[gauges_fn].data_type == "GeoDataFrame":
                 gdf_gauges = self.data_catalog.get_geodataframe(
                     gauges_fn,
-                    geom=self.basins,
+                    geom=self.basins_highres,
                     assert_gtype="Point",
                     handle_nodata=NoDataStrategy.IGNORE,
                     **kwargs,
                 )
             elif self.data_catalog[gauges_fn].data_type == "GeoDataset":
                 da = self.data_catalog.get_geodataset(
                     gauges_fn,
-                    geom=self.basins,
+                    geom=self.basins_highres,
                     handle_nodata=NoDataStrategy.IGNORE,
                     **kwargs,
                 )
@@ -1422,7 +1432,7 @@ def setup_areamap(
         """
         self.logger.info(f"Preparing '{col2raster}' map from '{area_fn}'.")
         gdf_org = self.data_catalog.get_geodataframe(
-            area_fn, geom=self.basins, dst_crs=self.crs
+            area_fn, geom=self.basins_highres, dst_crs=self.crs
         )
         if gdf_org.empty:
             self.logger.warning(
@@ -1787,7 +1797,7 @@ def _setup_waterbodies(self, waterbodies_fn, wb_type, min_area=0.0, **kwargs):
             kwargs.update(predicate="contains")
         gdf_org = self.data_catalog.get_geodataframe(
             waterbodies_fn,
-            geom=self.basins,
+            geom=self.basins_highres,
             handle_nodata=NoDataStrategy.IGNORE,
             **kwargs,
         )
@@ -2042,7 +2052,7 @@ def setup_glaciers(self, glaciers_fn="rgi", min_area=1):
         self.logger.info("Preparing glacier maps.")
         gdf_org = self.data_catalog.get_geodataframe(
             glaciers_fn,
-            geom=self.basins,
+            geom=self.basins_highres,
             predicate="intersects",
             handle_nodata=NoDataStrategy.IGNORE,
         )
@@ -3163,14 +3173,22 @@ def read_staticgeoms(
     def write_geoms(
         self,
         geom_fn: str = "staticgeoms",
+        precision: int = 5,
     ):
         """Write geoms in <root/geom_fn> in GeoJSON format."""
         # to write use self.geoms[var].to_file()
         if not self._write:
             raise IOError("Model opened in read-only mode")
         if self.geoms:
             self.logger.info("Writing model staticgeom to file.")
+            grid_size = 10 ** (-precision)
             for name, gdf in self.geoms.items():
+                # TODO change to geopandas functionality once geopandas 1.0.0 comes
+                # See https://github.com/geopandas/geopandas/releases/tag/v1.0.0-alpha1
+                gdf.geometry = shapely.set_precision(
+                    gdf.geometry,
+                    grid_size=grid_size,
+                )
                 fn_out = join(self.root, geom_fn, f"{name}.geojson")
                 gdf.to_file(fn_out, driver="GeoJSON")
 
@@ -3684,13 +3702,21 @@ def basins(self):
                 .set_index("value")
                 .sort_index()
             )
-            gdf.index.name = self._MAPS["basins"]
             self.set_geoms(gdf, name="basins")
         else:
             self.logger.warning(f"Basin map {self._MAPS['basins']} not found in grid.")
             gdf = None
         return gdf
 
+    @property
+    def basins_highres(self):
+        """Returns a high resolution basin(s) geometry."""
+        if "basins_highres" in self.geoms:
+            gdf = self.geoms["basins_highres"]
+        else:
+            gdf = self.basins
+        return gdf
+
     @property
     def rivers(self):
         """Return a river geometry as a geopandas.GeoDataFrame.