Bugfix subbasin (#263)

* remove deprecated calls to np.bool #261

* earlier error message if subbasin is too small #236

* add test

* update changelog

* bugfix test

* address review comments

docs/changelog.rst

@@ -30,6 +30,7 @@ Fixed
 - Wrong dtype for wflow_subcatch map. PR #247
 - **setup_gauges**: Allow snapping to river/mask for snap_uparea method. PR #248
 - Removed building a wflow model without a config file in the build notebook.
+- Deprecated np.bool and earlier error message for subbasin delination. PR #263
 
 v0.5.0 (February 2024)
 ======================

hydromt_wflow/wflow.py

@@ -417,6 +417,14 @@ def setup_rivers(
         """
         self.logger.info("Preparing river maps.")
 
+        # Check that river_upa threshold is bigger than the maximum uparea in the grid
+        if river_upa > float(self.grid[self._MAPS["uparea"]].max()):
+            raise ValueError(
+                f"river_upa threshold {river_upa} should be larger than the maximum \
+uparea in the grid {float(self.grid[self._MAPS['uparea']].max())} in order to create \
+river cells."
+            )
+
         rivdph_methods = ["gvf", "manning", "powlaw"]
         if rivdph_method not in rivdph_methods:
             raise ValueError(f'"{rivdph_method}" unknown. Select from {rivdph_methods}')

hydromt_wflow/workflows/basemaps.py

@@ -85,7 +85,7 @@ def hydrography(
     outidx = None
     if "mask" not in ds.coords and xy is None:
         ds.coords["mask"] = xr.Variable(
-            dims=ds.raster.dims, data=np.ones(ds.raster.shape, dtype=np.bool)
+            dims=ds.raster.dims, data=np.ones(ds.raster.shape, dtype=bool)
         )
     elif "mask" not in ds.coords:
         # NOTE if no subbasin mask is provided calculate it here
@@ -99,9 +99,27 @@ def hydrography(
         # NOTE: this mask is passed on from get_basin_geometry method
         logger.debug("Mask in dataset assumed to represent subbasins.")
     ncells = np.sum(ds["mask"].values)
-    logger.debug(f"(Sub)basin at original resolution has {ncells} cells.")
-
     scale_ratio = int(np.round(res / ds.raster.res[0]))
+
+    if ncells < 4:
+        raise ValueError(
+            "(Sub)basin at original resolution should at least consist of two cells on "
+            f"each axis and the total number of cells is {ncells}. "
+            "Consider using a larger domain or higher spatial resolution. "
+            "For subbasin models, consider a (higher) threshold on for example "
+            "upstream area or stream order to snap the outlet."
+        )
+    elif ncells < 100 and scale_ratio > 1:
+        logger.warning(
+            f"(Sub)basin at original resolution is small and has {ncells} cells. "
+            "This may results in errors later when upscaling flow directions. "
+            "If so, consider using a larger domain or higher spatial resolution. "
+            "For subbasin models, consider a (higher) threshold on for example "
+            "upstream area or stream order to snap the outlet."
+        )
+    else:
+        logger.debug(f"(Sub)basin at original resolution has {ncells} cells.")
+
     if scale_ratio > 1:  # upscale flwdir
         if flwdir is None:
             # NOTE initialize with mask is FALSE
@@ -154,7 +172,7 @@ def hydrography(
             mask_int.raster.set_nodata(-1)  # change nodata value
             ds_out.coords["mask"] = mask_int.raster.reproject_like(
                 da_flw, method="nearest"
-            ).astype(np.bool)
+            ).astype(bool)
             basins = ds_out["mask"].values.astype(np.int32)
             logger.warning(
                 "The basin delineation might be wrong as no original resolution outlets"
@@ -251,9 +269,10 @@ def hydrography(
     logger.debug(f"Outlet coordinates ({len(xy_pit[0])}/{npits}): {xy_pit_str}.")
     if np.any(np.asarray(ds_out.raster.shape) == 1):
         raise ValueError(
-            "The output extent should at consist of two cells on each axis. "
-            "Consider using a larger domain or higher spatial resolution. "
-            "For subbasin models, consider a (higher) threshold to snap the outlet."
+            "The output extent at model resolution should at least consist of two "
+            "cells on each axis. Consider using a larger domain or higher spatial "
+            "resolution. For subbasin models, consider a (higher) threshold to snap "
+            "the outlet."
         )
     return ds_out, flwdir_out
 

hydromt_wflow/workflows/river.py

@@ -521,7 +521,7 @@ def _width_fit(
     p0=[0.15, 0.65],
     logger=logger,  # rhine uparea based
 ):
-    outliers = np.full(np.sum(mask), False, dtype=np.bool)
+    outliers = np.full(np.sum(mask), False, dtype=bool)
     a, b = None, None
     # check if sufficient data
     if np.sum(mask) > 10:

tests/test_model_methods.py

@@ -45,6 +45,28 @@ def test_setup_basemaps(tmpdir):
 
     assert mod.grid["wflow_subcatch"].dtype == "int32"
 
+    # Test for too small basins
+    region = {"subbasin": [12.572061, 46.601984]}
+
+    with pytest.raises(ValueError) as error:  # noqa PT011
+        mod.setup_basemaps(
+            region=region,
+            hydrography_fn=hydrography,
+        )
+    assert str(error.value).startswith(
+        "(Sub)basin at original resolution should at least consist of two cells"
+    )
+
+    region = {"subbasin": [12.572061, 46.600359]}
+    with pytest.raises(ValueError) as error:  # noqa PT011
+        mod.setup_basemaps(
+            region=region,
+            hydrography_fn=hydrography,
+        )
+    assert str(error.value).startswith(
+        "The output extent at model resolution should at least consist of two cells on"
+    )
+
 
 def test_setup_grid(example_wflow_model):
     # Tests on setup_grid_from_raster