Merge pull request #90 from ArtesiaWater/dev

Major update
gwmod · Sep 23, 2022 · 33e0d48 · 33e0d48
2 parents ce84793 + 7171e54
commit 33e0d48
Show file tree

Hide file tree

Showing 98 changed files with 7,966 additions and 4,684 deletions.
diff --git a/.gitignore b/.gitignore
@@ -144,5 +144,4 @@ cython_debug/
 nlmod/bin/
 flowchartnlmod.pptx
 tests/data/
-examples/*/
-data/nhflo/
+docs/examples/*/
diff --git a/.readthedocs.yml b/.readthedocs.yml
@@ -5,20 +5,22 @@
 # Required
 version: 2
 
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-20.04
+  tools:
+    python: "3.9"
+
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
-  configuration: docs/source/conf.py
+  configuration: docs/conf.py
 
 # Optionally build your docs in additional formats such as PDF and ePub
-formats: all
+# formats: all
 
-# Optionally set the version of Python and requirements required to build your docs
+# Optionally declare the Python requirements required to build your docs
 python:
-  version: "3.7"
   install:
     - requirements: docs/requirements.txt
-    - method: pip
+    - method: setuptools
       path: .
-
-build:
-  image: latest
diff --git a/README.md b/README.md
@@ -27,7 +27,7 @@ Install the module with pip:
 
 `pip install nlmod`
 
-`nlmod` has many dependencies `xarray`, `flopy`, `rasterio`, `owslib`, `hydropandas`, `netcdf4`, `pyshp`, `rtree`, `openpyxl` and `matplotlib`.
+`nlmod` has many dependencies `xarray`, `flopy`, `rasterio`, `rioxarray`, `owslib`, `hydropandas`, `netcdf4`, `pyshp`, `rtree`, `openpyxl` and `matplotlib`.
 
 When using pip the dependencies are automatically installed. Some dependencies are notoriously hard to install on certain platforms. 
 Please see the [dependencies](https://github.com/ArtesiaWater/hydropandas#dependencies) section of the `hydropandas` package for more information on how to install these packages manually. 

diff --git a/docs/Makefile b/docs/Makefile
@@ -5,7 +5,7 @@
 # from the environment for the first two.
 SPHINXOPTS    ?=
 SPHINXBUILD   ?= sphinx-build
-SOURCEDIR     = source
+SOURCEDIR     = .
 BUILDDIR      = build
 
 # Put it first so that "make" without argument is like "make help".

diff --git a/docs/source/conf.py → docs/conf.py b/docs/source/conf.py → docs/conf.py
@@ -47,7 +47,6 @@
     "IPython.sphinxext.ipython_console_highlighting",  # lowercase didn't work
     "sphinx.ext.autosectionlabel",
     "nbsphinx",
-    "nbsphinx_link",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -85,3 +84,6 @@
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
 html_static_path = ["_static"]
+
+# Allow errors in notebooks, so we can see the error online
+nbsphinx_allow_errors = True
diff --git a/docs/source/examples.rst → docs/examples.rst b/docs/source/examples.rst → docs/examples.rst
diff --git a/examples/01_basic_model.ipynb → docs/examples/01_basic_model.ipynb b/examples/01_basic_model.ipynb → docs/examples/01_basic_model.ipynb
@@ -17,6 +17,7 @@
    "metadata": {},
    "source": [
     "### Contents<a name=\"TOC\"></a>\n",
+    "0. [Download MODFLOW-binaries](#binaries)\n",
     "1. [Create model](#create)\n",
     "2. [Run model](#run)\n",
     "3. [Visualise](#visualise)"
@@ -50,6 +51,24 @@
     "logging.basicConfig(level=logging.INFO)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### [0. Download MODFLOW-binariesl](#TOC)<a name=\"binaries\"></a>\n",
+    "To run MODFLOW, we need to download the MODFLOW-excecutables. We do this with the following code:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if not nlmod.util.check_presence_mfbinaries():\n",
+    "    nlmod.util.download_mfbinaries()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -58,7 +77,7 @@
     "\n",
     "With the code below we create a modflow model with the name 'IJmuiden'. This model has the following properties :\n",
     "- an extent that covers part of the Northsea, Noordzeekanaal and the small port city IJmuiden.\n",
-    "- a structured grid based on the subsurface models [Regis](https://www.dinoloket.nl/regis-ii-het-hydrogeologische-model) and [Geotop](https://www.dinoloket.nl/detaillering-van-de-bovenste-lagen-met-geotop). The Regis layers that are not present within the extent are removed. In this case we use b'MSz1' as the bottom layer of the model. Use `nlmod.regis.get_layer_names()` to get all the layer names of Regis. All Regis layers below this layer are not used in the model. Geotop is used to replace the holoceen layer in Regis because there is no kh or kv defined for the holoceen in Regis. Part of the model is in the North sea. Regis and Geotop have no data there. Therefore the Regis and Geotop layers are extrapolated from the shore and the seabed is added using bathymetry data from [Jarkus](https://www.openearth.nl/rws-bathymetry/2018.html).\n",
+    "- a structured grid based on the subsurface models [Regis](https://www.dinoloket.nl/regis-ii-het-hydrogeologische-model) and [Geotop](https://www.dinoloket.nl/detaillering-van-de-bovenste-lagen-met-geotop). The Regis layers that are not present within the extent are removed. In this case we use 'MSz1' as the bottom layer of the model. Use `nlmod.read.regis.get_layer_names()` to get all the layer names of Regis. All Regis layers below this layer are not used in the model. Geotop is used to replace the holoceen layer in Regis because there is no kh or kv defined for the holoceen in Regis. Part of the model is in the North sea. Regis and Geotop have no data there. Therefore the Regis and Geotop layers are extrapolated from the shore and the seabed is added using bathymetry data from [Jarkus](https://www.openearth.nl/rws-bathymetry/2018.html).\n",
     "- starting heads of 1 in every cell.\n",
     "- the model is a steady state model of a single time step.\n",
     "- big surface water bodies (Northsea, IJsselmeer, Markermeer, Noordzeekanaal) within the extent are added as a general head boundary. The surface water bodies are obtained from a [shapefile](..\\data\\opp_water.shp).\n",
@@ -80,14 +99,14 @@
     "# model settings\n",
     "model_ws = 'model1'\n",
     "model_name = 'IJmuiden'\n",
+    "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n",
     "extent = [95000., 105000., 494000., 500000.]\n",
     "delr = 100.\n",
     "delc = 100.\n",
     "steady_state = True\n",
     "start_time = '2015-1-1'\n",
-    "gridtype = 'structured'\n",
     "use_regis = True\n",
-    "regis_botm_layer = b'MSz1'\n",
+    "regis_botm_layer = 'MSz1'\n",
     "use_geotop = True\n",
     "add_northsea = True\n",
     "starting_head = 1.0"
@@ -99,29 +118,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# create empty model dataset\n",
-    "model_ds = nlmod.mdims.get_empty_model_ds(model_name, model_ws)\n",
-    "\n",
-    "# add time discretisation\n",
-    "model_ds = nlmod.mdims.set_model_ds_time(model_ds,\n",
-    "                                         start_time=start_time,\n",
-    "                                         steady_state=steady_state, perlen=365*5)\n",
-    "\n",
-    "# add spatial discretisation\n",
-    "extent, nrow, ncol = nlmod.read.regis.fit_extent_to_regis(extent,\n",
-    "                                                          delr,\n",
-    "                                                          delc)\n",
-    "\n",
-    "# layer model\n",
-    "layer_model = nlmod.read.regis.get_combined_layer_models(extent, delr, delc,\n",
+    "layer_model = nlmod.read.regis.get_combined_layer_models(extent,\n",
     "                                                         use_regis=use_regis,\n",
     "                                                         regis_botm_layer=regis_botm_layer,\n",
     "                                                         use_geotop=use_geotop,\n",
-    "                                                         cachedir=model_ds.cachedir,\n",
+    "                                                         cachedir=cachedir,\n",
     "                                                         cachename='combined_layer_ds.nc')\n",
     "\n",
-    "# create modflow packages\n",
-    "sim, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds)"
+    "# create a model ds by changing grid of layer_model\n",
+    "ds = nlmod.mdims.to_model_ds(layer_model, model_name, model_ws, delr=delr, delc=delc)\n",
+    "\n",
+    "# add time discretisation\n",
+    "ds = nlmod.mdims.set_ds_time(ds, start_time=start_time, steady_state=steady_state, perlen=365*5)\n",
+    "\n",
+    "if add_northsea:\n",
+    "    ds = nlmod.mdims.add_northsea(ds)"
    ]
   },
   {
@@ -130,25 +141,29 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# update model_ds from layer model\n",
-    "model_ds = nlmod.mdims.update_model_ds_from_ml_layer_ds(model_ds,\n",
-    "                                                        layer_model,\n",
-    "                                                        keep_vars=['x', 'y'],\n",
-    "                                                        gridtype=gridtype,\n",
-    "                                                        add_northsea=add_northsea,\n",
-    "                                                        cachedir=model_ds.cachedir)\n",
+    "# create simulation \n",
+    "sim = nlmod.gwf.sim(ds)\n",
+    "\n",
+    "# create time discretisation\n",
+    "tdis = nlmod.gwf.tdis(ds, sim)\n",
+    "\n",
+    "# create groundwater flow model\n",
+    "gwf = nlmod.gwf.gwf(ds, sim)\n",
+    "\n",
+    "# create ims\n",
+    "ims = nlmod.gwf.ims(sim)\n",
     "\n",
     "# Create discretization\n",
-    "nlmod.mfpackages.dis_from_model_ds(model_ds, gwf)\n",
+    "nlmod.gwf.dis(ds, gwf)\n",
     "\n",
     "# create node property flow\n",
-    "nlmod.mfpackages.npf_from_model_ds(model_ds, gwf)\n",
+    "nlmod.gwf.npf(ds, gwf)\n",
     "\n",
     "# Create the initial conditions package\n",
-    "nlmod.mfpackages.ic_from_model_ds(model_ds, gwf, starting_head=starting_head)\n",
+    "nlmod.gwf.ic(ds, gwf, starting_head=starting_head)\n",
     "\n",
     "# Create the output control package\n",
-    "nlmod.mfpackages.oc_from_model_ds(model_ds, gwf);"
+    "nlmod.gwf.oc(ds, gwf);"
    ]
   },
   {
@@ -159,12 +174,12 @@
    "source": [
     "# voeg grote oppervlaktewaterlichamen toe o.b.v. rws shape\n",
     "da_name = 'rws_oppwater'\n",
-    "rws_ds = nlmod.read.rws.get_surface_water(model_ds,\n",
+    "rws_ds = nlmod.read.rws.get_surface_water(ds,\n",
     "                                           da_name,\n",
-    "                                           cachedir=model_ds.cachedir,\n",
+    "                                           cachedir=ds.cachedir,\n",
     "                                           cachename=da_name)\n",
-    "model_ds.update(rws_ds)\n",
-    "ghb = nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name)"
+    "ds.update(rws_ds)\n",
+    "ghb = nlmod.gwf.ghb(ds, gwf, da_name)"
    ]
   },
   {
@@ -174,15 +189,15 @@
    "outputs": [],
    "source": [
     "# surface level drain\n",
-    "ahn_ds = nlmod.read.ahn.get_ahn(model_ds, cachedir=model_ds.cachedir, cachename='ahn')\n",
-    "model_ds.update(ahn_ds)\n",
+    "ahn_ds = nlmod.read.ahn.get_ahn(ds, cachedir=ds.cachedir, cachename='ahn')\n",
+    "ds.update(ahn_ds)\n",
     "\n",
-    "drn = nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf)\n",
+    "drn = nlmod.gwf.surface_drain_from_ds(ds, gwf)\n",
     "\n",
     "\n",
     "# add constant head cells at model boundaries\n",
-    "model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain']))    \n",
-    "chd = nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head')"
+    "ds.update(nlmod.gwf.constant_head.chd_at_model_edge(ds, ds['idomain']))    \n",
+    "chd = nlmod.gwf.chd(ds, gwf, head='starting_head')"
    ]
   },
   {
@@ -192,18 +207,18 @@
    "outputs": [],
    "source": [
     "# add knmi recharge to the model datasets\n",
-    "knmi_ds = nlmod.read.knmi.get_recharge(model_ds, cachedir=model_ds.cachedir, cachename='recharge')\n",
-    "model_ds.update(knmi_ds)\n",
+    "knmi_ds = nlmod.read.knmi.get_recharge(ds, cachedir=ds.cachedir, cachename='recharge')\n",
+    "ds.update(knmi_ds)\n",
     "\n",
     "# create recharge package\n",
-    "rch = nlmod.mfpackages.rch_from_model_ds(model_ds, gwf)"
+    "rch = nlmod.gwf.rch(ds, gwf)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A big part of the model data is stored in the variable `model_ds` which is an `xarray.Dataset`. The data is shown below."
+    "A big part of the model data is stored in the variable `ds` which is an `xarray.Dataset`. The data is shown below."
    ]
   },
   {
@@ -212,7 +227,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model_ds"
+    "ds"
    ]
   },
   {
@@ -221,7 +236,7 @@
    "source": [
     "### [2. Write and Run](#TOC)<a name=\"run\"></a>\n",
     "Now that we've created all the modflow packages we need to write them to modflow files. You always have to write the modflow data to the model workspace before you can run the model. You can write the model files and run the model using the function `nlmod.util.write_and_run_model()` as shown below. This function has two additional options:\n",
-    "1. Write the model dataset to the disk if `write_model_ds` is `True`. This makes it easier and faster to load model data if you ever need it. \n",
+    "1. Write the model dataset to the disk if `write_ds` is `True`. This makes it easier and faster to load model data if you ever need it. \n",
     "2. Write a copy of this Jupyter Notebook to the same directory as the modflow files if `nb_path` is the name of this Jupyter Notebook. It can be useful to have a copy of the script that created the modflow files, together with the files. "
    ]
   },
@@ -231,7 +246,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nlmod.util.write_and_run_model(gwf, model_ds, write_model_ds=True, nb_path='01_basic_model.ipynb')"
+    "nlmod.gwf.write_and_run_model(gwf, ds, write_ds=True, nb_path='01_basic_model.ipynb')"
    ]
   },
   {
@@ -240,7 +255,7 @@
    "source": [
     "### [3. Visualise](#TOC)<a name=\"visualise\"></a>\n",
     "\n",
-    "Using the `model_ds` and `gwf` variables it is quite easy to visualise model data. Below the modelgrid together with the surface water is shown."
+    "Using the `ds` and `gwf` variables it is quite easy to visualise model data. Below the modelgrid together with the surface water is shown."
    ]
   },
   {
@@ -249,8 +264,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ax = nlmod.visualise.plots.plot_modelgrid(model_ds, gwf)\n",
-    "ax.figure.savefig(os.path.join(model_ds.figdir, 'mgrid_swater.png'), bbox_inches='tight')"
+    "ax = nlmod.visualise.plots.plot_modelgrid(ds, gwf)\n",
+    "ax.figure.savefig(os.path.join(ds.figdir, 'mgrid_swater.png'), bbox_inches='tight')"
    ]
   },
   {
@@ -267,27 +282,34 @@
    "outputs": [],
    "source": [
     "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 11))\n",
-    "model_ds['ahn'].plot(ax=axes[0][0])\n",
-    "model_ds['bot'][0].plot(ax=axes[0][1])\n",
-    "model_ds['idomain'][0].plot(ax=axes[1][0])\n",
-    "model_ds['chd'][0].plot(ax=axes[1][1])\n",
+    "ds['ahn'].plot(ax=axes[0][0])\n",
+    "ds['botm'][0].plot(ax=axes[0][1])\n",
+    "ds['idomain'][0].plot(ax=axes[1][0])\n",
+    "ds['chd'][0].plot(ax=axes[1][1])\n",
     "for axes1 in axes:\n",
     "    for ax in axes1:\n",
     "        ax.axis('scaled')\n",
     "\n",
-    "fig.savefig(os.path.join(model_ds.figdir, 'ahn_bot_idom_chd.png'), bbox_inches='tight')\n",
+    "fig.savefig(os.path.join(ds.figdir, 'ahn_bot_idom_chd.png'), bbox_inches='tight')\n",
     "\n",
     "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 11))\n",
-    "model_ds['bathymetry'].plot(ax=axes[0][0])\n",
-    "model_ds['northsea'].plot(ax=axes[0][1])\n",
-    "model_ds['kh'][1].plot(ax=axes[1][0])\n",
-    "model_ds['recharge'].plot(ax=axes[1][1])\n",
+    "ds['bathymetry'].plot(ax=axes[0][0])\n",
+    "ds['northsea'].plot(ax=axes[0][1])\n",
+    "ds['kh'][1].plot(ax=axes[1][0])\n",
+    "ds['recharge'].plot(ax=axes[1][1])\n",
     "\n",
     "for axes1 in axes:\n",
     "    for ax in axes1:\n",
     "        ax.axis('scaled')\n",
-    "fig.savefig(os.path.join(model_ds.figdir, 'bath_nsea_kh_top.png'), bbox_inches='tight')"
+    "fig.savefig(os.path.join(ds.figdir, 'bath_nsea_kh_top.png'), bbox_inches='tight')"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -303,7 +325,7 @@
   },
   "anaconda-cloud": {},
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },