Improve layers

docs/examples/01_basic_model.ipynb

@@ -137,7 +137,7 @@
     ")\n",
     "\n",
     "if add_northsea:\n",
-    "    ds = nlmod.mdims.add_northsea(ds, cachedir=cachedir)\n"
+    "    ds = nlmod.read.rws.add_northsea(ds, cachedir=cachedir)\n"
    ]
   },
   {

docs/examples/03_local_grid_refinement.ipynb

@@ -140,7 +140,7 @@
     ")\n",
     "\n",
     "if add_northsea:\n",
-    "    ds = nlmod.mdims.add_northsea(ds, cachedir=cachedir)\n"
+    "    ds = nlmod.read.rws.add_northsea(ds, cachedir=cachedir)\n"
    ]
   },
   {

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

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

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__)
@@ -233,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
@@ -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.
 
@@ -576,12 +634,9 @@ 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.attrs["_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))
@@ -770,7 +825,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)]