Merge branch 'main' into biosphere

geb/HRUs.py

@@ -520,6 +520,19 @@ def __init__(self, data, model) -> None:
             self.scaling, axis=1
         )
         self.cell_size = self.data.grid.cell_size / self.scaling
+
+        # get lats and lons for subgrid
+        self.lon = np.linspace(
+            self.gt[0] + self.cell_size / 2,
+            self.gt[0] + self.cell_size * submask_width - self.cell_size / 2,
+            submask_width,
+        )
+        self.lat = np.linspace(
+            self.gt[3] + self.cell_size / 2,
+            self.gt[3] + self.cell_size * submask_height - self.cell_size / 2,
+            submask_height,
+        )
+
         if self.model.load_initial_data:
             self.land_use_type = np.load(
                 os.path.join(self.data.get_save_state_path(), "HRU.land_use_type.npz")

geb/cli.py

@@ -316,11 +316,7 @@ def decorator(func):
             "-d",
             type=str,
             multiple=True,
-            default=[
-                Path(os.environ.get("GEB_PACKAGE_DIR"))
-                / "examples"
-                / "data_catalog.yml"
-            ],
+            default=[Path(os.environ.get("GEB_PACKAGE_DIR")) / "data_catalog.yml"],
             help="""A list of paths to the data library YAML files. By default the data_catalog in the examples is used. If this is not set, defaults to data_catalog.yml""",
         )
         @click.option(

geb/setup/geb.py

@@ -1775,6 +1775,7 @@ def setup_groundwater(
         minimum_thickness_confined_layer=50,
         maximum_thickness_confined_layer=1000,
         intial_heads_source="GLOBGM",
+        force_one_layer=False,
     ):
         """
         Sets up the MODFLOW grid for GEB. This code is adopted from the GLOBGM
@@ -1829,7 +1830,7 @@ def setup_groundwater(
         confining_layer = self.snap_to_grid(confining_layer, self.grid)
         assert confining_layer.shape == self.grid.raster.shape
 
-        if not (confining_layer == 0).all():  # two-layer-model
+        if not (confining_layer == 0).all() and not force_one_layer:  # two-layer-model
             two_layers = True
         else:
             two_layers = False
@@ -1943,10 +1944,6 @@ def setup_groundwater(
             head_upper_layer = head_upper_layer - dem_globgm + dem
             assert head_upper_layer.shape == self.grid.raster.shape
 
-            # assert concistency of datasets. If one layer, this layer should be all nan
-            if not two_layers:
-                assert np.isnan(head_upper_layer).all()
-
             head_lower_layer = self.data_catalog.get_rasterdataset(
                 "head_lower_globgm", bbox=self.bounds, buffer=0
             ).compute()