remove top3d from calculate_thickness, fix for #65

docs/examples/03_local_grid_refinement.ipynb

@@ -71,9 +71,7 @@
     "model_ws = \"model3\"\n",
     "model_name = \"IJm_planeten\"\n",
     "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n",
-    "refine_shp_fname = os.path.join(\n",
-    "    nlmod.NLMOD_DATADIR, \"shapes\", \"planetenweg_ijmuiden\"\n",
-    ")\n",
+    "refine_shp_fname = os.path.join(\"data\", \"planetenweg_ijmuiden\")\n",
     "levels = 2\n",
     "extent = [95000.0, 105000.0, 494000.0, 500000.0]\n",
     "delr = 100.0\n",

docs/examples/07_resampling.ipynb

@@ -23,7 +23,8 @@
     "4. [Locally refined grid to structured grid](#4)\n",
     "5. [Fill nan values](#5)\n",
     "6. [Vector to grid](#6)\n",
-    "7. [Real world example](#7)"
+    "7. [Grid to reclist (stress period data)](#7)\n",
+    "8. [Real world example](#8)"
    ]
   },
   {
@@ -61,7 +62,12 @@
     "print(f'nlmod version: {nlmod.__version__}')\n",
     "\n",
     "# toon informatie bij het aanroepen van functies\n",
-    "logging.basicConfig(level=logging.INFO)"
+    "logging.basicConfig(level=logging.INFO)\n",
+    "\n",
+    "# verwijder ShapelyDeprecationWarning\n",
+    "import warnings\n",
+    "from shapely.errors import ShapelyDeprecationWarning\n",
+    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)"
    ]
   },
   {
@@ -272,7 +278,10 @@
    "outputs": [],
    "source": [
     "for method in ['nearest', 'linear', 'cubic', 'average', 'min']:\n",
-    "    struc_out = resample.structured_da_to_ds(struc2d, ds, method=method)\n",
+    "    if method == 'nearest':\n",
+    "        struc_out = resample.structured_da_to_ds(struc2d, ds, method=method)\n",
+    "    else:\n",
+    "        struc_out = resample.structured_da_to_ds(struc2d.astype(float), ds, method=method)        \n",
     "    compare_structured_data_arrays(struc2d, struc_out, method)"
    ]
   },
@@ -515,7 +524,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## [6. Vector to grid](#TOC)<a name=\"5\"></a>\n",
+    "## [6. Vector to grid](#TOC)<a name=\"6\"></a>\n",
     "\n",
     "Vector data can be points, lines or polygons often saved as shapefiles and visualised using GIS software. A common operation is to project vector data on a modelgrid. For example to add a surface water line to a grid. Here are some functions in `nlmod` to project vector data on a modelgrid."
    ]
@@ -824,7 +833,183 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## [7. Real world example](#TOC)<a name=\"6\"></a>\n",
+    "## [7. Grid to reclist](#TOC)<a name=\"7\"></a>\n",
+    "For some modflow packages (drn, riv, ghb, wel) you need to specify stress_period_data to create them using flopy. This stress_period_data consists of reclists (also called lrcd for a structured grid) for every time step. \n",
+    "\n",
+    "The function `da_to_reclist` can be used to convert grid data (both structured and vertex) to a reclist. This function has many arguments:\n",
+    "- `mask`, boolean DataArray to determine which cells should be added to the reclist. Can be 2d or 3d.\n",
+    "- `layer`, if `mask` is a 2d array the value of `layer` is used in the reclist. If `mask` is 3d or `first_active_layer` is True the `layer` argument is ignored.\n",
+    "- `only_active_cells`, if True only add cells with an idomain of 1 to the reclist\n",
+    "- `first_active_layer`, if True use the first active layer, obtained from the idomain, as the layer for each cell.\n",
+    "- `col1`,`col2` and `col3`, The column data of the reclist.\n",
+    "\n",
+    "The examples below show the result of each argument."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# add layer dimension\n",
+    "if 'layer' not in ds.dims:\n",
+    "    ds = ds.expand_dims({'layer':range(3)})\n",
+    "\n",
+    "# create some data arrays\n",
+    "ds['da1'] = ('layer', 'y','x'), np.random.randint(0,10,(ds.dims['layer'], ds.dims['y'],ds.dims['x']))\n",
+    "ds['da2'] = ('y','x'), np.random.randint(0,10,(ds.dims['y'],ds.dims['x']))\n",
+    "ds['da3'] = ('y','x'), np.random.randint(0,10,(ds.dims['y'],ds.dims['x']))\n",
+    "\n",
+    "# add a nodata value\n",
+    "ds.attrs['nodata'] = -999\n",
+    "\n",
+    "# create an idomain of ones except for the first cell which is zero\n",
+    "idomain = np.ones((ds.dims['layer'], ds.dims['y'],ds.dims['x']))\n",
+    "idomain[0,0,0] = 0\n",
+    "ds['idomain'] = ('layer','y','x'), idomain"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Mask and layer\n",
+    "If `mask` is a 2d array, the `layer` argument specifies the layer that is used in the reclist."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# structured 2d grid to reclist\n",
+    "mask2d = ds['da2'] == ds['da2'][0,0]\n",
+    "reclist1 = nlmod.mdims.da_to_reclist(ds, mask2d, col1=ds['da1'][0], col2='da2', layer=0, only_active_cells=False)\n",
+    "reclist1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If the `mask` is three dimensional the `layer` argument is ignored."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a 3dmask \n",
+    "mask3d = ds['da1'] == ds['da1'].values[0,0,0]\n",
+    "\n",
+    "# use this mask to create the reclist\n",
+    "reclist2 = nlmod.mdims.da_to_reclist(ds, mask3d, col1='da1',col2=100, layer=0, only_active_cells=False)\n",
+    "reclist2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Only active cells\n",
+    "With `only_active_cells=True` we make sure only active cells end up in the reclist. Which cells are active is based on the `idomain` in the model dataset."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Only return the cells with an active idomain \n",
+    "reclist3 = nlmod.mdims.da_to_reclist(ds, mask3d, col1='da1',col2=100, only_active_cells=True)\n",
+    "reclist3"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# also possible for a 2d grid\n",
+    "mask2d = ds['da2'] == ds['da2'][0,0]\n",
+    "reclist1 = nlmod.mdims.da_to_reclist(ds, mask2d, col1=ds['da1'][0], col2='da2', layer=0, only_active_cells=True)\n",
+    "reclist1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### First active_layer\n",
+    "Use `first_active_layer=True` to add the first active layer to the reclist. The first active layer is obtained from the idomain."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a reclist with col1 (str), col2 (DataArray), col3 (int)\n",
+    "reclist4 = nlmod.mdims.da_to_reclist(ds, mask2d, col1='da2',col2='da3', first_active_layer=True)\n",
+    "reclist4"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Reclist columns\n",
+    "The `col1`, `col2` and `col3` arguments specify what data should be listed in the reclist. The types can be `str`,`xarray.DataArray`,`None` or other. If the value is a `str` the corresponding DataArray from the Dataset is used to get data for the reclist. If the value is an `xarray.DataArray` the DataArray is used. If the value is `None` the column is not added to the reclist and if the value is from another type the value is used for every record in the reclist.\n",
+    "\n",
+    "Be aware that if `mask` is a 3d array, the DataArrays of the column should also be 3d."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create a reclist with col1 (str), col2 (DataArray), col3 (int)\n",
+    "reclist5 = nlmod.mdims.da_to_reclist(ds, mask3d, col1=ds['idomain'],col2='da1',col3=9, layer=0, only_active_cells=False)\n",
+    "reclist5"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Vertex model to reclist"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# add some random DataArray to the vertex dataset\n",
+    "da_vert = np.random.randint(0,10,(dsv['area'].shape))\n",
+    "dsv['da_vert'] = ('icell2d'), da_vert\n",
+    "\n",
+    "# create rec list from a vertex dataset\n",
+    "mask_vert = dsv['da_vert'] == dsv['da_vert'][0]\n",
+    "reclist6 = nlmod.mdims.da_to_reclist(dsv, mask_vert, col1='da_vert',col2=2330, only_active_cells=False)\n",
+    "reclist6"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## [8. Real world example](#TOC)<a name=\"8\"></a>\n",
     "In this example we will resample the values of the dutch Digital Terrain Model (DTM) from AHN4 to a structured grid and a vertex grid, using several methods. First we will download the AHN-information."
    ]
   },
@@ -930,9 +1115,9 @@
   },
   "anaconda-cloud": {},
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "nlmod",
    "language": "python",
-   "name": "python3"
+   "name": "nlmod"
   },
   "language_info": {
    "codemirror_mode": {

nlmod/gwf/gwf.py

@@ -27,7 +27,7 @@ def gwf(ds, sim, **kwargs):
     ----------
     ds : xarray.Dataset
         dataset with model data. Should have the dimension 'time' and the
-        attributes: model_name, mfversion, model_ws, time_units, start_time,
+        attributes: model_name, mfversion, model_ws, time_units, start,
         perlen, nstp, tsmult
     sim : flopy MFSimulation
         simulation object.
@@ -267,7 +267,7 @@ def ghb(ds, gwf, da_name, pname="ghb", **kwargs):
         ghb package
     """
 
-    ghb_rec = mdims.da_to_rec_list(
+    ghb_rec = mdims.da_to_reclist(
         ds,
         ds[f"{da_name}_cond"] != 0,
         col1=f"{da_name}_peil",
@@ -411,7 +411,7 @@ def chd(ds, gwf, chd="chd", head="starting_head", pname="chd", **kwargs):
         chd package
     """
     # get the stress_period_data
-    chd_rec = mdims.da_to_rec_list(ds, ds[chd] != 0, col1=head)
+    chd_rec = mdims.da_to_reclist(ds, ds[chd] != 0, col1=head)
 
     chd = flopy.mf6.ModflowGwfchd(
         gwf,
@@ -448,7 +448,7 @@ def surface_drain_from_ds(ds, gwf, surface_drn_cond=1000, pname="drn", **kwargs)
 
     ds.attrs["surface_drn_cond"] = surface_drn_cond
     mask = ds["ahn"].notnull()
-    drn_rec = mdims.da_to_rec_list(
+    drn_rec = mdims.da_to_reclist(
         ds,
         mask,
         col1="ahn",

nlmod/gwf/recharge.py

@@ -46,7 +46,7 @@ def model_datasets_to_rch(gwf, ds, pname="rch", **kwargs):
         mask = ds["rch_name"] != ""
         recharge = "rch_name"
 
-    spd = mdims.da_to_rec_list(
+    spd = mdims.da_to_reclist(
         ds,
         mask,
         col1=recharge,
@@ -137,7 +137,7 @@ def model_datasets_to_evt(
         mask = ds["evt_name"] != ""
         rate = "evt_name"
 
-    spd = mdims.da_to_rec_list(
+    spd = mdims.da_to_reclist(
         ds,
         mask,
         col1=surface,

nlmod/gwf/surface_water.py

@@ -777,7 +777,7 @@ def gdf_to_seasonal_pkg(
         **kwargs,
     )
     # add timeseries for the seasons 'winter' and 'summer'
-    tmin = pd.to_datetime(ds.time.start_time)
+    tmin = pd.to_datetime(ds.time.start)
     if tmin.month in summer_months:
         ts_data = [(0.0, 0.0, 1.0)]
     else: