From 90d8774b2fd786b0c838bb96406c66ee783c2f90 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 12:03:12 +0200
Subject: [PATCH 1/9] Improvements in layer definition

Add update_ds_from_layer_ds
Remove first_active_layer from ds, calculate when needed
Add some deprecation warnings
Improve structured_da_to_ds so it works for 3d grids to vertex grid as well
Remove Nan layers from regis just after getting the data again
Improve documentation
---
 nlmod/mdims/mbase.py    |   2 +
 nlmod/mdims/mgrid.py    | 114 +++++++++++++++++++++++++++++++++++-----
 nlmod/mdims/mlayers.py  |  67 +++++++++++++++++------
 nlmod/mdims/resample.py |  12 ++++-
 nlmod/read/regis.py     |  14 ++++-
 5 files changed, 178 insertions(+), 31 deletions(-)

diff --git a/nlmod/mdims/mbase.py b/nlmod/mdims/mbase.py
index b7482644..0c1afcb9 100644
--- a/nlmod/mdims/mbase.py
+++ b/nlmod/mdims/mbase.py
@@ -182,6 +182,8 @@ def extrapolate_ds(ds, mask=None):
     if not mask.any():
         # all of the model cells are is inside the known area
         return ds
+    if mask.all():
+        raise (Exception("The model only contains NaNs"))
     if ds.gridtype == "vertex":
         x = ds.x.data
         y = ds.y.data
diff --git a/nlmod/mdims/mgrid.py b/nlmod/mdims/mgrid.py
index 97338f8d..4bd6b751 100644
--- a/nlmod/mdims/mgrid.py
+++ b/nlmod/mdims/mgrid.py
@@ -25,12 +25,18 @@
 from packaging import version
 
 from .. import util
-from .mlayers import set_idomain, get_first_active_layer_from_idomain
+from .mlayers import (
+    set_idomain,
+    get_first_active_layer,
+    fill_nan_top_botm_kh_kv,
+)
 from .resample import (
     get_resampled_ml_layer_ds_vertex,
     affine_transform_gdf,
     get_affine_world_to_mod,
+    structured_da_to_ds,
 )
+from .mbase import extrapolate_ds
 from .rdp import rdp
 
 logger = logging.getLogger(__name__)
@@ -283,6 +289,58 @@ def refine(
     return ds
 
 
+def update_ds_from_layer_ds(ds, layer_ds, method="nearest", **kwargs):
+    """
+    Add variables from a layer Dataset to a model Dataset.
+    Keep de grid-information from the model Dataset (x and y or icell2d), but update
+    the layer dimension when neccesary.
+
+    Parameters
+    ----------
+    ds : TYPE
+        DESCRIPTION.
+    layer_ds : TYPE
+        DESCRIPTION.
+    method : str
+        THe method used for resampling layer_ds to the grid of ds
+    **kwargs : TYPE
+        DESCRIPTION.
+
+    Returns
+    -------
+    ds : TYPE
+        DESCRIPTION.
+
+    """
+    if not layer_ds.layer.equals(ds.layer):
+        # do not change the original Dataset
+        layer_ds = layer_ds.copy()
+        # update layers in ds
+        drop_vars = []
+        for var in ds.data_vars:
+            if "layer" in ds[var].dims:
+                if var not in layer_ds.data_vars:
+                    logger.info(
+                        f"Variable {var} is dropped, as it has dimension layer, but is not defined in layer_ds"
+                    )
+                drop_vars.append(var)
+        if len(drop_vars) > 0:
+            ds = ds.drop_vars(drop_vars)
+        ds = ds.assign_coords({"layer": layer_ds.layer})
+    if method in ["nearest", "linear"]:
+        layer_ds = layer_ds.interp(
+            x=ds.x, y=ds.y, method="nearest", kwargs={"fill_value": None}
+        )
+        for var in layer_ds.data_vars:
+            ds[var] = layer_ds[var]
+    else:
+        for var in layer_ds.data_vars:
+            ds[var] = structured_da_to_ds(layer_ds[var], ds, method=method)
+    ds = extrapolate_ds(ds)
+    ds = fill_nan_top_botm_kh_kv(ds, **kwargs)
+    return ds
+
+
 def col_to_list(col_in, ds, cellids):
     """Convert array data in ds to a list of values for specific cells.
 
@@ -342,9 +400,7 @@ def col_to_list(col_in, ds, cellids):
     return col_lst
 
 
-def lrc_to_reclist(
-    layers, rows, columns, cellids, ds, col1=None, col2=None, col3=None
-):
+def lrc_to_reclist(layers, rows, columns, cellids, ds, col1=None, col2=None, col3=None):
     """Create a reclist for stress period data from a set of cellids.
 
     Used for structured grids.
@@ -561,7 +617,9 @@ def da_to_reclist(
         if only_active_cells:
             cellids = np.where((mask) & (ds["idomain"] == 1))
             ignore_cells = np.sum((mask) & (ds["idomain"] != 1))
-            logger.info(f'ignore {ignore_cells} out of {np.sum(mask)} cells because idomain is inactive')
+            logger.info(
+                f"ignore {ignore_cells} out of {np.sum(mask)} cells because idomain is inactive"
+            )
         else:
             cellids = np.where(mask)
 
@@ -575,16 +633,15 @@ def da_to_reclist(
             return lrc_to_reclist(layers, rows, columns, cellids, ds, col1, col2, col3)
     else:
         if first_active_layer:
-            if "first_active_layer" not in ds:
-                ds["first_active_layer"] = get_first_active_layer_from_idomain(
-                    ds["idomain"]
-                )
-            cellids = np.where((mask) & (ds["first_active_layer"] != ds.nodata))
-            layers = col_to_list("first_active_layer", ds, cellids)
+            fal = get_first_active_layer(ds)
+            cellids = np.where((mask) & (fal != fal["_FillValue"]))
+            layers = col_to_list(fal, ds, cellids)
         elif only_active_cells:
             cellids = np.where((mask) & (ds["idomain"][layer] == 1))
             ignore_cells = np.sum((mask) & (ds["idomain"][layer] != 1))
-            logger.info(f'ignore {ignore_cells} out of {np.sum(mask)} cells because idomain is inactive')
+            logger.info(
+                f"ignore {ignore_cells} out of {np.sum(mask)} cells because idomain is inactive"
+            )
             layers = col_to_list(layer, ds, cellids)
         else:
             cellids = np.where(mask)
@@ -763,7 +820,38 @@ def add_info_to_gdf(
     min_total_overlap=0.5,
     geom_type="Polygon",
 ):
-    """ "Add information from gdf_from to gdf_to"""
+    """
+    Add information from gdf_from to gdf_to
+
+    Parameters
+    ----------
+    gdf_to : TYPE
+        DESCRIPTION.
+    gdf_from : TYPE
+        DESCRIPTION.
+    columns : TYPE, optional
+        DESCRIPTION. The default is None.
+    desc : TYPE, optional
+        DESCRIPTION. The default is "".
+    silent : TYPE, optional
+        DESCRIPTION. The default is False.
+    min_total_overlap : TYPE, optional
+        DESCRIPTION. The default is 0.5.
+    geom_type : TYPE, optional
+        DESCRIPTION. The default is "Polygon".
+
+    Raises
+    ------
+
+        DESCRIPTION.
+
+    Returns
+    -------
+    gdf_to : TYPE
+        DESCRIPTION.
+
+    """
+
     gdf_to = gdf_to.copy()
     if columns is None:
         columns = gdf_from.columns[~gdf_from.columns.isin(gdf_to.columns)]
diff --git a/nlmod/mdims/mlayers.py b/nlmod/mdims/mlayers.py
index 6f4b9de4..5a0ed12b 100644
--- a/nlmod/mdims/mlayers.py
+++ b/nlmod/mdims/mlayers.py
@@ -656,6 +656,9 @@ def add_kh_kv_from_ml_layer_to_dataset(
     some model dataset, such as regis, also have 'c' and 'kd' values. These
     are ignored at the moment
     """
+    DeprecationWarning(
+        "add_kh_kv_from_ml_layer_to_dataset is deprecated. Please use update_ds_from_layer_ds instead."
+    )
     ds.attrs["anisotropy"] = anisotropy
     ds.attrs["fill_value_kh"] = fill_value_kh
     ds.attrs["fill_value_kv"] = fill_value_kv
@@ -786,7 +789,7 @@ def fill_top_bot_kh_kv_at_mask(ds, fill_mask, gridtype="structured"):
     Parameters
     ----------
     ds : xr.DataSet
-        model dataset, should contain 'first_active_layer'
+        model dataset
     fill_mask : xr.DataArray
         1 where a cell should be replaced by masked value.
     gridtype : str, optional
@@ -889,7 +892,24 @@ def fill_top_and_bottom(ds):
     return ds
 
 
-def set_idomain(ds, nodata=-999, remove_nan_layers=True):
+def set_idomain(ds, remove_nan_layers=True):
+    """
+    Set idmomain in a model Dataset
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        The model Dataset.
+    remove_nan_layers : bool, optional
+        Removes layers which only contain inactive cells. The default is True.
+
+    Returns
+    -------
+    ds : TYPE
+        DESCRIPTION.
+
+    """
+    """Set idomain in a model ds"""
     # set idomain with a default of -1 (pass-through)
     ds["idomain"] = xr.full_like(ds["botm"], -1, int)
     # set idomain of cells  with a positive thickness to 1
@@ -902,13 +922,31 @@ def set_idomain(ds, nodata=-999, remove_nan_layers=True):
         # only keep layers with at least one active cell
         ds = ds.sel(layer=(ds["idomain"] > 0).any(ds["idomain"].dims[1:]))
     # TODO: set idomain above/below the first/last active layer to 0
-    ds["first_active_layer"] = get_first_active_layer_from_idomain(
-        ds["idomain"], nodata=nodata
-    )
-    ds.attrs["nodata"] = nodata
+    # TODO: remove 'active' and replace by logic of keeping inactive cells in idomain
     return ds
 
 
+def get_first_active_layer(ds, **kwargs):
+    """
+    Get the first active layer in each cell from a model ds
+
+    Parameters
+    ----------
+    ds : xr.DataSet
+        DESCRIPTION.
+    **kwargs : dict
+        Kwargs are passed on to get_first_active_layer_from_idomain.
+
+    Returns
+    -------
+    first_active_layer : xr.DataArray
+        raster in which each cell has the zero based number of the first
+        active layer. Shape can be (y, x) or (icell2d)
+
+    """
+    return get_first_active_layer_from_idomain(ds["idomain"], **kwargs)
+
+
 def get_first_active_layer_from_idomain(idomain, nodata=-999):
     """get the first active layer in each cell from the idomain.
 
@@ -926,7 +964,7 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999):
         raster in which each cell has the zero based number of the first
         active layer. Shape can be (y, x) or (icell2d)
     """
-    logger.info("get first active modellayer for each cell in idomain")
+    logger.debug("get first active modellayer for each cell in idomain")
 
     first_active_layer = xr.where(idomain[0] == 1, 0, nodata)
     for i in range(1, idomain.shape[0]):
@@ -935,7 +973,7 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999):
             i,
             first_active_layer,
         )
-
+    first_active_layer.attrs["_FillValue"] = nodata
     return first_active_layer
 
 
@@ -954,7 +992,8 @@ def add_northsea(ds, cachedir=None):
     ds.update(rws.get_northsea(ds, cachedir=cachedir, cachename="sea_ds.nc"))
 
     # fill top, bot, kh, kv at sea cells
-    fill_mask = (ds["first_active_layer"] == ds.nodata) * ds["northsea"]
+    fal = get_first_active_layer(ds)
+    fill_mask = (fal == fal["_FillValue"]) * ds["northsea"]
     ds = fill_top_bot_kh_kv_at_mask(ds, fill_mask, gridtype=ds.attrs["gridtype"])
 
     # add bathymetry noordzee
@@ -970,11 +1009,7 @@ def add_northsea(ds, cachedir=None):
     ds = jarkus.add_bathymetry_to_top_bot_kh_kv(ds, ds["bathymetry"], fill_mask)
 
     # update idomain on adjusted tops and bots
-    ds["thickness"] = calculate_thickness(ds)
-    ds["idomain"] = update_idomain_from_thickness(
-        ds["idomain"], ds["thickness"], ds["northsea"]
-    )
-    ds["first_active_layer"] = get_first_active_layer_from_idomain(ds["idomain"])
+    ds = set_idomain(ds)
     return ds
 
 
@@ -1004,7 +1039,9 @@ def update_idomain_from_thickness(idomain, thickness, mask):
         raster with adjusted idomain of each cell. dimensions should be
         (layer, y, x) or (layer, icell2d).
     """
-
+    DeprecationWarning(
+        "update_idomain_from_thickness is deprecated. Please use set_idomain instead."
+    )
     for ilay, thick in enumerate(thickness):
         if ilay == 0:
             mask1 = (thick == 0) * mask
diff --git a/nlmod/mdims/resample.py b/nlmod/mdims/resample.py
index 35aa3d43..bf8922b1 100644
--- a/nlmod/mdims/resample.py
+++ b/nlmod/mdims/resample.py
@@ -14,6 +14,8 @@
 from shapely.affinity import affine_transform
 from affine import Affine
 
+from ..util import get_da_from_da_ds
+
 logger = logging.getLogger(__name__)
 
 
@@ -533,7 +535,11 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
     elif ds.gridtype == "vertex":
         # assume the grid is a quadtree grid, where cells are refined by splitting them
         # in 4
-        da_out = xr.full_like(ds["area"], nodata)
+        dims = list(da.dims)
+        dims.remove("y")
+        dims.remove("x")
+        dims.append("icell2d")
+        da_out = get_da_from_da_ds(ds, dims=tuple(dims), data=nodata)
         for area in np.unique(ds["area"]):
             dx = dy = np.sqrt(area)
             x, y = get_xy_mid_structured(ds.extent, dx, dy)
@@ -545,7 +551,9 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
             da_temp = da_temp.rio.write_crs(da.rio.crs)
             da_temp = da.rio.reproject_match(da_temp, resampling)
             mask = ds["area"] == area
-            da_out[mask] = da_temp.sel(y=ds["y"][mask], x=ds["x"][mask])
+            da_out.loc[dict(icell2d=mask)] = da_temp.sel(
+                y=ds["y"][mask], x=ds["x"][mask]
+            )
     else:
         raise (Exception(f"Gridtype {ds.gridtype} not supported"))
 
diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py
index 10b231c1..6e8b5e31 100644
--- a/nlmod/read/regis.py
+++ b/nlmod/read/regis.py
@@ -75,7 +75,12 @@ def get_combined_layer_models(
 
 
 @cache.cache_netcdf
-def get_regis(extent, botm_layer="AKc", variables=("top", "botm", "kh", "kv")):
+def get_regis(
+    extent,
+    botm_layer="AKc",
+    variables=("top", "botm", "kh", "kv"),
+    remove_nan_layers=True,
+):
     """get a regis dataset projected on the modelgrid.
 
     Parameters
@@ -91,6 +96,9 @@ def get_regis(extent, botm_layer="AKc", variables=("top", "botm", "kh", "kv")):
         a tuple of the variables to keep from the regis Dataset. Possible
         entries in the list are 'top', 'botm', 'kD', 'c', 'kh', 'kv', 'sdh' and
         'sdv'. The default is ("top", "botm", "kh", "kv").
+    remove_nan_layers : bool, optional
+        When True, layers which contain only NaNs for the botm array are removed.
+        The default is True.
 
     Returns
     -------
@@ -117,6 +125,10 @@ def get_regis(extent, botm_layer="AKc", variables=("top", "botm", "kh", "kv")):
     # rename bottom to botm, as it is called in FloPy
     ds = ds.rename_vars({"bottom": "botm"})
 
+    if remove_nan_layers:
+        # only keep layers with at least one active cell
+        ds = ds.sel(layer=~(np.isnan(ds["botm"])).all(ds["botm"].dims[1:]))
+
     # slice data vars
     ds = ds[list(variables)]
 

From 7f2092beab3c92421b0ee54862215f7bcff5deae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 12:37:28 +0200
Subject: [PATCH 2/9] Create 12_layer_generation.ipynb

---
 docs/examples/12_layer_generation.ipynb | 193 ++++++++++++++++++++++++
 1 file changed, 193 insertions(+)
 create mode 100644 docs/examples/12_layer_generation.ipynb

diff --git a/docs/examples/12_layer_generation.ipynb b/docs/examples/12_layer_generation.ipynb
new file mode 100644
index 00000000..77f773bb
--- /dev/null
+++ b/docs/examples/12_layer_generation.ipynb
@@ -0,0 +1,193 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "163f60be",
+   "metadata": {},
+   "source": [
+    "# Layer generation\n",
+    "This notebook contains two workflows to generate a model dataset and add layer infomration to it:\n",
+    "\n",
+    "- Get data\n",
+    "- Regrid\n",
+    "\n",
+    "or\n",
+    "\n",
+    "- Create grid\n",
+    "- Add data to ds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4ef9418d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import nlmod\n",
+    "from shapely.geometry import LineString\n",
+    "import logging\n",
+    "logging.basicConfig(level=logging.INFO)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "13fb7319",
+   "metadata": {},
+   "source": [
+    "## Get REGIS data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "67d3b65b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "extent = [100000, 101000, 400000, 401000]\n",
+    "\n",
+    "# downlaod regis\n",
+    "regis = nlmod.read.regis.get_regis(extent)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e96dbc81",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "f, ax = nlmod.plot.get_map(extent, figsize=5)\n",
+    "nlmod.plot.da(regis['top'].max('layer'), edgecolor=\"k\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "49a18b1b",
+   "metadata": {},
+   "source": [
+    "## Define some properties of the grid\n",
+    "We choose a resolution of the calculation grid (delr and delc) larger than the resolution of regis (100 meter)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "73677bcd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delr = 200.0\n",
+    "delc = 200.0\n",
+    "\n",
+    "line = LineString([(extent[0], extent[2]), (extent[1], extent[3])])\n",
+    "refinement_features = [([line], \"line\", 3)]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d22c22ff",
+   "metadata": {},
+   "source": [
+    "## Regrid data\n",
+    "create a model dataset from regis and refine."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb9f87d",
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "ds = nlmod.mbase.to_model_ds(regis, delr=delr, delc=delc)\n",
+    "ds = nlmod.mgrid.refine(ds, model_ws=\"model_12\", refinement_features=refinement_features)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cf9fbfd6",
+   "metadata": {},
+   "source": [
+    "When we plot the model top, we see that all information has a resolution of 200 meter, also in the smaller cells. The original resolution of the regis data is 100 meter however. So information is lost in the process."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6f46ad4d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f, ax = nlmod.plot.get_map(extent, figsize=5)\n",
+    "nlmod.plot.da(ds[\"top\"], ds=ds, edgecolor=\"k\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e976498d",
+   "metadata": {},
+   "source": [
+    "## Add data to existing grid\n",
+    "We can also first create a grid first, and then warp the information from regis to this grid. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "af0c984f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = nlmod.mbase.get_default_ds(extent, delr=delr, delc=delc)\n",
+    "ds = nlmod.mgrid.refine(ds, model_ws=\"model_12\", refinement_features=refinement_features)\n",
+    "ds = nlmod.mgrid.update_ds_from_layer_ds(ds, regis, method=\"average\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "98b6d820",
+   "metadata": {},
+   "source": [
+    "When we plot the model top, we now see that in the cells with an equal resolution to that of Regis (or smaller) no infomration is lost is lost."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f65ea923",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f, ax = nlmod.plot.get_map(extent, figsize=5)\n",
+    "nlmod.plot.da(ds[\"top\"], ds=ds, edgecolor=\"k\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From cf52c3acc415dcb40e328b6419d7e2b14ccd9d05 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 13:13:57 +0200
Subject: [PATCH 3/9] Remove nodata from knmi.py

Move add_northsea
Fix tests
Fix noteboooks
---
 docs/examples/01_basic_model.ipynb           |  2 +-
 docs/examples/03_local_grid_refinement.ipynb |  2 +-
 nlmod/mdims/mgrid.py                         |  2 +-
 nlmod/mdims/mlayers.py                       | 37 ---------------
 nlmod/read/knmi.py                           | 49 ++++++--------------
 nlmod/read/rws.py                            | 38 +++++++++++++++
 tests/test_004_northsea.py                   |  9 ++--
 tests/test_007_run_notebooks.py              |  5 ++
 8 files changed, 66 insertions(+), 78 deletions(-)

diff --git a/docs/examples/01_basic_model.ipynb b/docs/examples/01_basic_model.ipynb
index 446a743f..a187d3af 100644
--- a/docs/examples/01_basic_model.ipynb
+++ b/docs/examples/01_basic_model.ipynb
@@ -137,7 +137,7 @@
     ")\n",
     "\n",
     "if add_northsea:\n",
-    "    ds = nlmod.mdims.add_northsea(ds)\n"
+    "    ds = nlmod.read.rws.add_northsea(ds)\n"
    ]
   },
   {
diff --git a/docs/examples/03_local_grid_refinement.ipynb b/docs/examples/03_local_grid_refinement.ipynb
index b06283cb..0c531d05 100644
--- a/docs/examples/03_local_grid_refinement.ipynb
+++ b/docs/examples/03_local_grid_refinement.ipynb
@@ -140,7 +140,7 @@
     ")\n",
     "\n",
     "if add_northsea:\n",
-    "    ds = nlmod.mdims.add_northsea(ds)\n"
+    "    ds = nlmod.read.rws.add_northsea(ds)\n"
    ]
   },
   {
diff --git a/nlmod/mdims/mgrid.py b/nlmod/mdims/mgrid.py
index 4bd6b751..0bf8003b 100644
--- a/nlmod/mdims/mgrid.py
+++ b/nlmod/mdims/mgrid.py
@@ -634,7 +634,7 @@ def da_to_reclist(
     else:
         if first_active_layer:
             fal = get_first_active_layer(ds)
-            cellids = np.where((mask) & (fal != fal["_FillValue"]))
+            cellids = np.where((mask) & (fal != fal._FillValue))
             layers = col_to_list(fal, ds, cellids)
         elif only_active_cells:
             cellids = np.where((mask) & (ds["idomain"][layer] == 1))
diff --git a/nlmod/mdims/mlayers.py b/nlmod/mdims/mlayers.py
index 5a0ed12b..798505e5 100644
--- a/nlmod/mdims/mlayers.py
+++ b/nlmod/mdims/mlayers.py
@@ -5,7 +5,6 @@
 import xarray as xr
 
 from . import resample
-from ..read import jarkus, rws
 
 logger = logging.getLogger(__name__)
 
@@ -977,42 +976,6 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999):
     return first_active_layer
 
 
-def add_northsea(ds, cachedir=None):
-    """a) get cells from modelgrid that are within the northsea, add data
-    variable 'northsea' to ds
-    b) fill top, bot, kh and kv add northsea cell by extrapolation
-    c) get bathymetry (northsea depth) from jarkus. Add datavariable
-    bathymetry to model dataset"""
-
-    logger.info(
-        "nan values at the northsea are filled using the bathymetry from jarkus"
-    )
-
-    # find grid cells with northsea
-    ds.update(rws.get_northsea(ds, cachedir=cachedir, cachename="sea_ds.nc"))
-
-    # fill top, bot, kh, kv at sea cells
-    fal = get_first_active_layer(ds)
-    fill_mask = (fal == fal["_FillValue"]) * ds["northsea"]
-    ds = fill_top_bot_kh_kv_at_mask(ds, fill_mask, gridtype=ds.attrs["gridtype"])
-
-    # add bathymetry noordzee
-    ds.update(
-        jarkus.get_bathymetry(
-            ds,
-            ds["northsea"],
-            cachedir=cachedir,
-            cachename="bathymetry_ds.nc",
-        )
-    )
-
-    ds = jarkus.add_bathymetry_to_top_bot_kh_kv(ds, ds["bathymetry"], fill_mask)
-
-    # update idomain on adjusted tops and bots
-    ds = set_idomain(ds)
-    return ds
-
-
 def update_idomain_from_thickness(idomain, thickness, mask):
     """get new idomain from thickness in the cells where mask is 1 (or True).
 
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index 7644906c..03e08942 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -7,12 +7,13 @@
 
 from .. import cache, util
 from ..mdims.resample import get_affine_mod_to_world
+from ..mdims.mlayers import get_first_active_layer
 
 logger = logging.getLogger(__name__)
 
 
 @cache.cache_netcdf
-def get_recharge(ds, method="linear", nodata=None):
+def get_recharge(ds, method="linear"):
     """add multiple recharge packages to the groundwater flow model with knmi
     data by following these steps:
 
@@ -40,19 +41,12 @@ def get_recharge(ds, method="linear", nodata=None):
         If 'linear', calculate recharge by subtracting evaporation from precipitation.
         If 'separate', add precipitation as 'recharge' and evaporation as 'evaporation'.
         The defaults is 'linear'.
-    nodata : int, optional
-        if the first_active_layer data array in ds has this value,
-        it means this cell is inactive in all layers. If nodata is None the
-        nodata value in ds is used.
-        the default is None.
 
     Returns
     -------
     ds : xr.DataSet
         dataset with spatial model data including the rch raster
     """
-    if nodata is None:
-        nodata = ds.nodata
 
     start = pd.Timestamp(ds.time.attrs["start"])
     # include the end day in the time series.
@@ -71,7 +65,7 @@ def get_recharge(ds, method="linear", nodata=None):
     shape = [len(ds_out[dim]) for dim in dims]
 
     locations, oc_knmi_prec, oc_knmi_evap = get_knmi_at_locations(
-        ds, start=start, end=end, nodata=nodata
+        ds, start=start, end=end
     )
 
     # add closest precipitation and evaporation measurement station to each cell
@@ -162,18 +156,13 @@ def _add_ts_to_ds(timeseries, loc_sel, variable, ds):
             ds[variable].loc[loc_sel.index, :] = model_recharge.values
 
 
-def get_locations_vertex(ds, nodata=-999):
+def get_locations_vertex(ds):
     """get dataframe with the locations of the grid cells of a vertex grid.
 
     Parameters
     ----------
     ds : xr.DataSet
         dataset containing relevant model grid information
-    nodata : int, optional
-        if the first_active_layer data array in ds has this value,
-        it means this cell is inactive in all layers. If nodata is None the
-        nodata value in ds is used.
-        the default is None
 
     Returns
     -------
@@ -182,7 +171,8 @@ def get_locations_vertex(ds, nodata=-999):
         includes the columns: x, y and layer
     """
     # get active locations
-    icell2d_active = np.where(ds["first_active_layer"] != nodata)[0]
+    fal = get_first_active_layer(ds)
+    icell2d_active = np.where(fal=fal._FillValue)[0]
 
     # create dataframe from active locations
     x = ds["x"].sel(icell2d=icell2d_active)
@@ -191,7 +181,7 @@ def get_locations_vertex(ds, nodata=-999):
         # transform coordinates into real-world coordinates
         affine = get_affine_mod_to_world(ds)
         x, y = affine * (x, y)
-    layer = ds["first_active_layer"].sel(icell2d=icell2d_active)
+    layer = fal.sel(icell2d=icell2d_active)
     locations = pd.DataFrame(
         index=icell2d_active, data={"x": x, "y": y, "layer": layer}
     )
@@ -200,18 +190,13 @@ def get_locations_vertex(ds, nodata=-999):
     return locations
 
 
-def get_locations_structured(ds, nodata=-999):
+def get_locations_structured(ds):
     """get dataframe with the locations of the grid cells of a structured grid.
 
     Parameters
     ----------
     ds : xr.DataSet
         dataset containing relevant model grid information
-    nodata : int, optional
-        if the first_active_layer data array in ds has this value,
-        it means this cell is inactive in all layers. If nodata is None the
-        nodata value in ds is used.
-        the default is None
 
     Returns
     -------
@@ -221,16 +206,15 @@ def get_locations_structured(ds, nodata=-999):
     """
 
     # store x and y mids in locations of active cells
-    rows, columns = np.where(ds["first_active_layer"] != nodata)
+    fal = get_first_active_layer(ds)
+    rows, columns = np.where(fal=fal._FillValue)
     x = np.array([ds["x"].data[col] for col in columns])
     y = np.array([ds["y"].data[row] for row in rows])
     if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:
         # transform coordinates into real-world coordinates
         affine = get_affine_mod_to_world(ds)
         x, y = affine * (x, y)
-    layers = [
-        ds["first_active_layer"].data[row, col] for row, col in zip(rows, columns)
-    ]
+    layers = [fal.data[row, col] for row, col in zip(rows, columns)]
 
     locations = hpd.ObsCollection(
         pd.DataFrame(
@@ -241,7 +225,7 @@ def get_locations_structured(ds, nodata=-999):
     return locations
 
 
-def get_knmi_at_locations(ds, start="2010", end=None, nodata=-999):
+def get_knmi_at_locations(ds, start="2010", end=None):
     """get knmi data at the locations of the active grid cells in ds.
 
     Parameters
@@ -252,11 +236,6 @@ def get_knmi_at_locations(ds, start="2010", end=None, nodata=-999):
         start date of measurements that you want, The default is '2010'.
     end :  str or datetime, optional
         end date of measurements that you want, The default is None.
-    nodata : int, optional
-        if the first_active_layer data array in ds has this value,
-        it means this cell is inactive in all layers. If nodata is None the
-        nodata value in ds is used.
-        the default is None
 
     Raises
     ------
@@ -274,9 +253,9 @@ def get_knmi_at_locations(ds, start="2010", end=None, nodata=-999):
     """
     # get locations
     if ds.gridtype == "structured":
-        locations = get_locations_structured(ds, nodata=nodata)
+        locations = get_locations_structured(ds)
     elif ds.gridtype == "vertex":
-        locations = get_locations_vertex(ds, nodata=nodata)
+        locations = get_locations_vertex(ds)
     else:
         raise ValueError("gridtype should be structured or vertex")
 
diff --git a/nlmod/read/rws.py b/nlmod/read/rws.py
index 5f5511b7..f32923fa 100644
--- a/nlmod/read/rws.py
+++ b/nlmod/read/rws.py
@@ -11,6 +11,8 @@
 
 from .. import cache, mdims, util
 
+from . import jarkus
+
 logger = logging.getLogger(__name__)
 
 
@@ -129,3 +131,39 @@ def get_northsea(ds, da_name="northsea"):
     ds_out = mdims.gdf_to_bool_dataset(ds, swater_zee, modelgrid, da_name)
 
     return ds_out
+
+
+def add_northsea(ds, cachedir=None):
+    """a) get cells from modelgrid that are within the northsea, add data
+    variable 'northsea' to ds
+    b) fill top, bot, kh and kv add northsea cell by extrapolation
+    c) get bathymetry (northsea depth) from jarkus. Add datavariable
+    bathymetry to model dataset"""
+
+    logger.info(
+        "nan values at the northsea are filled using the bathymetry from jarkus"
+    )
+
+    # find grid cells with northsea
+    ds.update(get_northsea(ds, cachedir=cachedir, cachename="sea_ds.nc"))
+
+    # fill top, bot, kh, kv at sea cells
+    fal = mdims.get_first_active_layer(ds)
+    fill_mask = (fal == fal._FillValue) * ds["northsea"]
+    ds = mdims.fill_top_bot_kh_kv_at_mask(ds, fill_mask, gridtype=ds.attrs["gridtype"])
+
+    # add bathymetry noordzee
+    ds.update(
+        jarkus.get_bathymetry(
+            ds,
+            ds["northsea"],
+            cachedir=cachedir,
+            cachename="bathymetry_ds.nc",
+        )
+    )
+
+    ds = jarkus.add_bathymetry_to_top_bot_kh_kv(ds, ds["bathymetry"], fill_mask)
+
+    # update idomain on adjusted tops and bots
+    ds = mdims.set_idomain(ds)
+    return ds
diff --git a/tests/test_004_northsea.py b/tests/test_004_northsea.py
index 43e7b3ea..cc3d628d 100644
--- a/tests/test_004_northsea.py
+++ b/tests/test_004_northsea.py
@@ -65,7 +65,8 @@ def test_fill_top_bot_kh_kv_seamodel():
     ds = test_001_model.test_get_ds_from_cache("basic_sea_model")
     ds.update(nlmod.read.rws.get_northsea(ds))
 
-    fill_mask = (ds["first_active_layer"] == ds.nodata) * ds["northsea"]
+    fal = nlmod.mdims.get_first_active_layer(ds)
+    fill_mask = (fal == fal._FillValue) * ds["northsea"]
     ds = nlmod.mdims.fill_top_bot_kh_kv_at_mask(ds, fill_mask)
 
     return ds
@@ -77,7 +78,8 @@ def test_fill_top_bot_kh_kv_nosea():
     ds = test_001_model.test_get_ds_from_cache("small_model")
     ds.update(nlmod.read.rws.get_northsea(ds))
 
-    fill_mask = (ds["first_active_layer"] == ds.nodata) * ds["northsea"]
+    fal = nlmod.mdims.get_first_active_layer(ds)
+    fill_mask = (fal == fal._FillValue) * ds["northsea"]
     ds = nlmod.mdims.fill_top_bot_kh_kv_at_mask(ds, fill_mask)
 
     return ds
@@ -114,7 +116,8 @@ def test_add_bathymetrie_to_top_bot_kh_kv_seamodel():
     ds.update(nlmod.read.rws.get_northsea(ds))
     ds.update(nlmod.read.jarkus.get_bathymetry(ds, ds["northsea"]))
 
-    fill_mask = (ds["first_active_layer"] == ds.nodata) * ds["northsea"]
+    fal = nlmod.mdims.get_first_active_layer(ds)
+    fill_mask = (fal == fal._FillValue) * ds["northsea"]
 
     ds = nlmod.read.jarkus.add_bathymetry_to_top_bot_kh_kv(
         ds, ds["bathymetry"], fill_mask
diff --git a/tests/test_007_run_notebooks.py b/tests/test_007_run_notebooks.py
index 57a5eabc..e1b54b91 100644
--- a/tests/test_007_run_notebooks.py
+++ b/tests/test_007_run_notebooks.py
@@ -74,3 +74,8 @@ def test_run_notebook_10_modpath():
 @pytest.mark.notebooks
 def test_run_notebook_11_grid_rotation():
     _run_notebook(nbdir, "11_grid_rotation.ipynb")
+
+
+@pytest.mark.notebooks
+def test_run_notebook_12_layer_generation():
+    _run_notebook(nbdir, "12_layer_generation.ipynb")

From 3e3407c93440aaccf4d2fb95243ba9b3cad2e35e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 13:17:55 +0200
Subject: [PATCH 4/9] Fix bug

---
 nlmod/read/knmi.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index 03e08942..d5cd6444 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -172,7 +172,7 @@ def get_locations_vertex(ds):
     """
     # get active locations
     fal = get_first_active_layer(ds)
-    icell2d_active = np.where(fal=fal._FillValue)[0]
+    icell2d_active = np.where(fal == fal._FillValue)[0]
 
     # create dataframe from active locations
     x = ds["x"].sel(icell2d=icell2d_active)
@@ -207,7 +207,7 @@ def get_locations_structured(ds):
 
     # store x and y mids in locations of active cells
     fal = get_first_active_layer(ds)
-    rows, columns = np.where(fal=fal._FillValue)
+    rows, columns = np.where(fal == fal._FillValue)
     x = np.array([ds["x"].data[col] for col in columns])
     y = np.array([ds["y"].data[row] for row in rows])
     if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:

From 5bb6495cbca69ab76a69f9619ce5747dc7305d4e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 14:12:10 +0200
Subject: [PATCH 5/9] Change == in != in knmi

---
 nlmod/read/knmi.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index d5cd6444..e226c882 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -172,7 +172,7 @@ def get_locations_vertex(ds):
     """
     # get active locations
     fal = get_first_active_layer(ds)
-    icell2d_active = np.where(fal == fal._FillValue)[0]
+    icell2d_active = np.where(fal != fal._FillValue)[0]
 
     # create dataframe from active locations
     x = ds["x"].sel(icell2d=icell2d_active)
@@ -207,7 +207,7 @@ def get_locations_structured(ds):
 
     # store x and y mids in locations of active cells
     fal = get_first_active_layer(ds)
-    rows, columns = np.where(fal == fal._FillValue)
+    rows, columns = np.where(fal != fal._FillValue)
     x = np.array([ds["x"].data[col] for col in columns])
     y = np.array([ds["y"].data[row] for row in rows])
     if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:

From 59299f495b98039914226f6d9b1c07803fc1b1e6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 14:20:32 +0200
Subject: [PATCH 6/9] Add rotated again in get_patches in plot.py

---
 nlmod/visualise/plots.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nlmod/visualise/plots.py b/nlmod/visualise/plots.py
index e7cd2249..3ce1d420 100644
--- a/nlmod/visualise/plots.py
+++ b/nlmod/visualise/plots.py
@@ -387,7 +387,7 @@ def da(da, ds=None, ax=None, rotated=False, **kwargs):
         if isinstance(ds, list):
             patches = ds
         else:
-            patches = get_patches(ds)
+            patches = get_patches(ds, rotated=rotated)
         pc = PatchCollection(patches, **kwargs)
         pc.set_array(da)
         ax.add_collection(pc)

From e22cd4d054b468cfd0f4af944b48b491f068b08b Mon Sep 17 00:00:00 2001
From: dbrakenhoff <d.brakenhoff@artesia-water.nl>
Date: Fri, 28 Oct 2022 14:20:42 +0200
Subject: [PATCH 7/9] fix floating docstring

---
 nlmod/mdims/mlayers.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/nlmod/mdims/mlayers.py b/nlmod/mdims/mlayers.py
index 798505e5..83a2efc4 100644
--- a/nlmod/mdims/mlayers.py
+++ b/nlmod/mdims/mlayers.py
@@ -908,7 +908,6 @@ def set_idomain(ds, remove_nan_layers=True):
         DESCRIPTION.
 
     """
-    """Set idomain in a model ds"""
     # set idomain with a default of -1 (pass-through)
     ds["idomain"] = xr.full_like(ds["botm"], -1, int)
     # set idomain of cells  with a positive thickness to 1

From 64b8955b79a6c12457cf2f86adcdcc6154eb5717 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 14:38:43 +0200
Subject: [PATCH 8/9] change fal._FillValue in fal.attrs["_FillValue"]

---
 nlmod/mdims/mgrid.py | 4 ++--
 nlmod/read/knmi.py   | 4 ++--
 nlmod/read/rws.py    | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/nlmod/mdims/mgrid.py b/nlmod/mdims/mgrid.py
index 0bf8003b..79d8ce43 100644
--- a/nlmod/mdims/mgrid.py
+++ b/nlmod/mdims/mgrid.py
@@ -239,7 +239,7 @@ def refine(
         # create a modelgrid with only one layer, to speed up Gridgen
         top = ds["top"].values
         botm = ds["botm"].values[[0]]
-        modelgrid = modelgrid_from_ds(ds, nlay=1, top=top, botm=botm)
+        modelgrid = modelgrid_from_ds(ds, rotated=False, nlay=1, top=top, botm=botm)
         g = Gridgen(modelgrid, model_ws=model_ws, exe_name=exe_name)
 
     ds_has_rotation = "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0
@@ -634,7 +634,7 @@ def da_to_reclist(
     else:
         if first_active_layer:
             fal = get_first_active_layer(ds)
-            cellids = np.where((mask) & (fal != fal._FillValue))
+            cellids = np.where((mask) & (fal != fal.attrs["_FillValue"]))
             layers = col_to_list(fal, ds, cellids)
         elif only_active_cells:
             cellids = np.where((mask) & (ds["idomain"][layer] == 1))
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index e226c882..666b3860 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -172,7 +172,7 @@ def get_locations_vertex(ds):
     """
     # get active locations
     fal = get_first_active_layer(ds)
-    icell2d_active = np.where(fal != fal._FillValue)[0]
+    icell2d_active = np.where(fal != fal.attrs["_FillValue"])[0]
 
     # create dataframe from active locations
     x = ds["x"].sel(icell2d=icell2d_active)
@@ -207,7 +207,7 @@ def get_locations_structured(ds):
 
     # store x and y mids in locations of active cells
     fal = get_first_active_layer(ds)
-    rows, columns = np.where(fal != fal._FillValue)
+    rows, columns = np.where(fal != fal.attrs["_FillValue"])
     x = np.array([ds["x"].data[col] for col in columns])
     y = np.array([ds["y"].data[row] for row in rows])
     if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:
diff --git a/nlmod/read/rws.py b/nlmod/read/rws.py
index f32923fa..4bf8f9a1 100644
--- a/nlmod/read/rws.py
+++ b/nlmod/read/rws.py
@@ -149,7 +149,7 @@ def add_northsea(ds, cachedir=None):
 
     # fill top, bot, kh, kv at sea cells
     fal = mdims.get_first_active_layer(ds)
-    fill_mask = (fal == fal._FillValue) * ds["northsea"]
+    fill_mask = (fal == fal.attrs["_FillValue"]) * ds["northsea"]
     ds = mdims.fill_top_bot_kh_kv_at_mask(ds, fill_mask, gridtype=ds.attrs["gridtype"])
 
     # add bathymetry noordzee

From 5410712e11df3a80ddc5b57e500a9e9a87730ddb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Fri, 28 Oct 2022 14:49:25 +0200
Subject: [PATCH 9/9] Fix regis tests

---
 nlmod/read/regis.py            | 13 +++++++++++--
 tests/test_002_regis_geotop.py | 13 +++++--------
 2 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py
index 6e8b5e31..21a73318 100644
--- a/nlmod/read/regis.py
+++ b/nlmod/read/regis.py
@@ -20,7 +20,11 @@
 
 @cache.cache_netcdf
 def get_combined_layer_models(
-    extent, regis_botm_layer="AKc", use_regis=True, use_geotop=True
+    extent,
+    regis_botm_layer="AKc",
+    use_regis=True,
+    use_geotop=True,
+    remove_nan_layers=True,
 ):
     """combine layer models into a single layer model.
 
@@ -42,6 +46,9 @@ def get_combined_layer_models(
         True if part of the layer model should be REGIS. The default is True.
     use_geotop : bool, optional
         True if part of the layer model should be geotop. The default is True.
+    remove_nan_layers : bool, optional
+        When True, layers which contain only NaNs for the botm array are removed.
+        The default is True.
 
     Returns
     -------
@@ -55,7 +62,9 @@ def get_combined_layer_models(
     """
 
     if use_regis:
-        regis_ds = get_regis(extent, regis_botm_layer)
+        regis_ds = get_regis(
+            extent, regis_botm_layer, remove_nan_layers=remove_nan_layers
+        )
     else:
         raise ValueError("layer models without REGIS not supported")
 
diff --git a/tests/test_002_regis_geotop.py b/tests/test_002_regis_geotop.py
index a4c9b03d..91538920 100644
--- a/tests/test_002_regis_geotop.py
+++ b/tests/test_002_regis_geotop.py
@@ -15,7 +15,7 @@ def test_get_regis(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
 
     regis_ds = regis.get_regis(extent)
 
-    assert regis_ds.dims["layer"] == 132
+    assert regis_ds.dims["layer"] == 20
 
     return regis_ds
 
@@ -23,12 +23,12 @@ def test_get_regis(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
 # @pytest.mark.skip(reason="too slow")
 def test_get_regis_botm_layer_BEk1(
     extent=[98700.0, 99000.0, 489500.0, 489700.0],
-    botm_layer="BEk1",
+    botm_layer="MSc",
 ):
     # extent, nrow, ncol = regis.fit_extent_to_regis(extent, delr, delc)
     regis_ds = regis.get_regis(extent, botm_layer)
-    assert regis_ds.dims["layer"] == 18
-    assert regis_ds.layer.values[-1] == "BEk1"
+    assert regis_ds.dims["layer"] == 15
+    assert regis_ds.layer.values[-1] == botm_layer
     return regis_ds
 
 
@@ -54,10 +54,7 @@ def test_get_regis_geotop_keep_all_layers(
     extent=[98600.0, 99000.0, 489400.0, 489700.0], delr=100.0, delc=100.0
 ):
     regis_geotop_ds = regis.get_combined_layer_models(
-        extent, use_regis=True, use_geotop=True
-    )
-    nlmod.mdims.to_model_ds(
-        regis_geotop_ds, delr=delr, delc=delc, remove_nan_layers=False
+        extent, use_regis=True, use_geotop=True, remove_nan_layers=False
     )
     assert regis_geotop_ds.dims["layer"] == 135
     return regis_geotop_ds