Multiple improvements (#289)

* Multiple improvements

Update modpath executable
Add nlmod.plot.geotop_lithok_on_map()
Add "concentration" to perioddata in nlmod.gwf.buy()
Add optional "layer" argument to extrapolate_ds, so you can use method for geotop-data with layer="z" as well

* Simplify fill_top_and_bottom

See comment from Onno in issue #275

* Remove the layer dimension in top returned from get_regis

* Also remove the layer dimension in top returned from get_geotop

* Fix codacy issues

* Make sure add_geotop_to_regis_layers works without a layer dimension in top of regis and geotop

And fix codacy issue

* Make sure insert_layer works without a layer dimension in ds

And fix some tests
Fix extrapolate_ds for 2d-data-arrays

* Fix some geotop-errors

* Change mapserver layer for level areas of Wetterkip Fryslan

* Change extent of test, as Wetterskip Frysland does not return any geometries anymore

Change to HHNK

* Update url for level areas of noorderzijlvest

And improve test

* Process comments from @dbrakenhoff and other fixes

fix split_layer_ds and combine_layer_ds, and add tests
small fixes to notebooks

* Process more comments and improve tests

* Revert change which causes failed tests everywhere

* Add tests from plot.flopy, gwt and hfb

* Update 11_grid_rotation.ipynb

* Make a copy of the dataset in most methods in nlmod.layers

* Add docstrings and minor change in `set_model_top`

* Process comments from @dbrakenhoff

docs/examples/04_modifying_layermodels.ipynb

@@ -268,9 +268,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "layer_names = []\n",
-    "colors_new = {}\n",
-    "\n",
     "for j, i in split_reindexer.items():\n",
     "    if j not in colors:\n",
     "        colors[j] = colors[i]"

docs/examples/11_grid_rotation.ipynb

@@ -239,10 +239,7 @@
     "oc = nlmod.gwf.oc(ds, gwf)\n",
     "\n",
     "# create recharge package\n",
-    "rch = nlmod.gwf.rch(ds, gwf)\n",
-    "\n",
-    "# create storage package\n",
-    "sto = nlmod.gwf.sto(ds, gwf)"
+    "rch = nlmod.gwf.rch(ds, gwf)"
    ]
   },
   {

docs/examples/12_layer_generation.ipynb

@@ -71,7 +71,7 @@
    "outputs": [],
    "source": [
     "f, ax = nlmod.plot.get_map(extent, figsize=5)\n",
-    "nlmod.plot.data_array(regis[\"top\"].max(\"layer\"), edgecolor=\"k\")"
+    "nlmod.plot.data_array(regis[\"top\"], edgecolor=\"k\")"
    ]
   },
   {

docs/examples/17_unsaturated_zone_flow.ipynb

@@ -95,7 +95,7 @@
     "y = np.floor(y / 100) * 100\n",
     "dx = dy = 100\n",
     "extent = [x, x + dx, y, y + dy]\n",
-    "regis = nlmod.read.regis.get_regis(extent, cachename=\"regis.nc\", cachedir=cachedir)"
+    "regis = nlmod.read.regis.get_regis(extent, drop_layer_dim_from_top=False, cachename=\"regis.nc\", cachedir=cachedir)"
    ]
   },
   {

nlmod/dims/base.py

@@ -179,7 +179,7 @@ def to_model_ds(
     return ds
 
 
-def extrapolate_ds(ds, mask=None):
+def extrapolate_ds(ds, mask=None, layer="layer"):
     """Fill missing data in layermodel based on nearest interpolation.
 
     Used for ensuring layer model contains data everywhere. Useful for
@@ -190,18 +190,20 @@ def extrapolate_ds(ds, mask=None):
     ----------
     ds : xarray.DataSet
         Model layer DataSet
-    mask: np.ndarray, optional
+    mask : np.ndarray, optional
         Boolean mask for each cell, with a value of True if its value needs to
         be determined. When mask is None, it is determined from the botm-
         variable. The default is None.
+    layer : str, optional
+        The name of the layer dimension. The default is 'layer'.
 
     Returns
     -------
     ds : xarray.DataSet
         filled layermodel
     """
     if mask is None:
-        mask = np.isnan(ds["botm"]).all("layer").data
+        mask = np.isnan(ds["botm"]).all(layer).data
     if not mask.any():
         # all of the model cells are is inside the known area
         return ds
@@ -225,9 +227,9 @@ def extrapolate_ds(ds, mask=None):
         data = ds[key].data
         if ds[key].dims == dims:
             if np.isnan(data[mask]).sum() > 0:  # do not update if no NaNs
-                data[mask] = data[~mask, i]
-        elif ds[key].dims == ("layer",) + dims:
-            for lay in range(len(ds["layer"])):
+                data[mask] = data[~mask][i]
+        elif ds[key].dims == (layer,) + dims:
+            for lay in range(len(ds[layer])):
                 if np.isnan(data[lay, mask]).sum() > 0:  # do not update if no NaNs
                     data[lay, mask] = data[lay, ~mask][i]
         else: