diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9c2b3fdf..a2dd31f4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -17,13 +17,13 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.9] + python-version: [3.11] steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} @@ -32,12 +32,10 @@ jobs: python -m pip install --upgrade pip pip install -e .[ci] - - name: Lint with flake8 + - name: Download executables needed for tests + shell: bash -l {0} run: | - # stop the build if there are Python syntax errors or undefined names - flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics - # exit-zero treats all errors as warnings. - flake8 . --count --exit-zero --max-complexity=10 --max-line-length=80 --statistics + python -c "import nlmod; nlmod.util.download_mfbinaries()" - name: Run tests only env: diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index 72f83ba7..0890d7eb 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -9,22 +9,24 @@ on: jobs: deploy: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 + - name: Set up Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: '3.9' + - name: Install dependencies run: | python -m pip install --upgrade pip pip install build setuptools wheel + - name: build binary wheel and a source tarball run: | python -m build --sdist --wheel --outdir dist/ + - name: Publish a Python distribution to PyPI uses: pypa/gh-action-pypi-publish@release/v1 with: diff --git a/README.md b/README.md index 0286b4d9..5b3f23f5 100644 --- a/README.md +++ b/README.md @@ -18,9 +18,13 @@ groundwater models, makes models more reproducible and transparent. The functions in `nlmod` have four main objectives: -1. Create and adapt the temporal and spatial discretization of a MODFLOW model using an xarray Dataset (`nlmod.dims`). -2. Download and read data from external sources, project this data on the modelgrid and add this data to an xarray Dataset (`nlmod.read`). -3. Use data in an xarray Dataset to build modflow packages for both groundwater flow and transport models using FloPy (`nlmod.sim`, `nlmod.gwf` and `nlmod.gwt` for Modflow 6 and `nlmod.modpath` for Modpath). +1. Create and adapt the temporal and spatial discretization of a MODFLOW model using an + xarray Dataset (`nlmod.dims`). +2. Download and read data from external sources, project this data on the modelgrid and + add this data to an xarray Dataset (`nlmod.read`). +3. Use data in an xarray Dataset to build modflow packages for both groundwater flow + and transport models using FloPy (`nlmod.sim`, `nlmod.gwf` and `nlmod.gwt` for + Modflow 6 and `nlmod.modpath` for Modpath). 4. Visualise modeldata in Python (`nlmod.plot`) or GIS software (`nlmod.gis`). More information can be found on the documentation-website: diff --git a/docs/conf.py b/docs/conf.py index b2798667..27df1c88 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -10,10 +10,11 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # -from nlmod import __version__ import os import sys +from nlmod import __version__ + sys.path.insert(0, os.path.abspath(".")) diff --git a/docs/examples/00_model_from_scratch.ipynb b/docs/examples/00_model_from_scratch.ipynb index 3163a444..f452a9e5 100644 --- a/docs/examples/00_model_from_scratch.ipynb +++ b/docs/examples/00_model_from_scratch.ipynb @@ -20,10 +20,9 @@ "outputs": [], "source": [ "import flopy as fp\n", - "import matplotlib.pyplot as plt\n", - "import nlmod\n", - "import numpy as np\n", - "import pandas as pd" + "import pandas as pd\n", + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/01_basic_model.ipynb b/docs/examples/01_basic_model.ipynb index ca18a257..ed358fb3 100644 --- a/docs/examples/01_basic_model.ipynb +++ b/docs/examples/01_basic_model.ipynb @@ -18,12 +18,7 @@ "metadata": {}, "outputs": [], "source": [ - "import logging\n", - "import os\n", "\n", - "import flopy\n", - "import geopandas as gpd\n", - "import matplotlib.pyplot as plt\n", "import nlmod" ] }, diff --git a/docs/examples/02_surface_water.ipynb b/docs/examples/02_surface_water.ipynb index 1dcbba79..ab34565f 100644 --- a/docs/examples/02_surface_water.ipynb +++ b/docs/examples/02_surface_water.ipynb @@ -25,8 +25,9 @@ "import os\n", "\n", "import flopy\n", - "import rioxarray\n", "import matplotlib.pyplot as plt\n", + "import rioxarray\n", + "\n", "import nlmod" ] }, diff --git a/docs/examples/03_local_grid_refinement.ipynb b/docs/examples/03_local_grid_refinement.ipynb index 2ccf8c01..9fed4ef6 100644 --- a/docs/examples/03_local_grid_refinement.ipynb +++ b/docs/examples/03_local_grid_refinement.ipynb @@ -22,15 +22,12 @@ "source": [ "import os\n", "\n", - "import flopy\n", - "import numpy as np\n", - "import pandas as pd\n", "import geopandas as gpd\n", - "import hydropandas as hpd\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "from IPython.display import HTML\n", - "import nlmod\n", - "import warnings" + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/04_modifying_layermodels.ipynb b/docs/examples/04_modifying_layermodels.ipynb index 5a37f1e9..f7fe8687 100644 --- a/docs/examples/04_modifying_layermodels.ipynb +++ b/docs/examples/04_modifying_layermodels.ipynb @@ -20,11 +20,11 @@ "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", - "import nlmod\n", "import numpy as np\n", - "import pandas as pd\n", - "from nlmod.plot import DatasetCrossSection\n", - "from shapely.geometry import LineString" + "from shapely.geometry import LineString\n", + "\n", + "import nlmod\n", + "from nlmod.plot import DatasetCrossSection" ] }, { diff --git a/docs/examples/05_caching.ipynb b/docs/examples/05_caching.ipynb index 0c705248..f2e2de21 100644 --- a/docs/examples/05_caching.ipynb +++ b/docs/examples/05_caching.ipynb @@ -20,12 +20,10 @@ "outputs": [], "source": [ "import os\n", - "import flopy\n", - "import geopandas as gpd\n", - "import matplotlib.pyplot as plt\n", - "import nlmod\n", - "import numpy as np\n", - "import xarray as xr" + "\n", + "import xarray as xr\n", + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/06_gridding_vector_data.ipynb b/docs/examples/06_gridding_vector_data.ipynb index 11267282..56104eb7 100644 --- a/docs/examples/06_gridding_vector_data.ipynb +++ b/docs/examples/06_gridding_vector_data.ipynb @@ -18,17 +18,15 @@ "metadata": {}, "outputs": [], "source": [ - "import nlmod\n", + "import geopandas as gpd\n", + "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", - "\n", - "import matplotlib.pyplot as plt\n", - "\n", - "import geopandas as gpd\n", + "from IPython.display import display\n", "from shapely.geometry import LineString, Point\n", "from shapely.geometry import Polygon as shp_polygon\n", "\n", - "from IPython.display import display" + "import nlmod" ] }, { @@ -391,8 +389,9 @@ " ds = ds.expand_dims({\"layer\": range(3)})\n", "\n", "# create some data arrays\n", - "ds[\"da1\"] = (\"layer\", \"y\", \"x\"), np.random.randint(\n", - " 0, 10, (ds.sizes[\"layer\"], ds.sizes[\"y\"], ds.sizes[\"x\"])\n", + "ds[\"da1\"] = (\n", + " (\"layer\", \"y\", \"x\"),\n", + " np.random.randint(0, 10, (ds.sizes[\"layer\"], ds.sizes[\"y\"], ds.sizes[\"x\"])),\n", ")\n", "ds[\"da2\"] = (\"y\", \"x\"), np.random.randint(0, 10, (ds.sizes[\"y\"], ds.sizes[\"x\"]))\n", "ds[\"da3\"] = (\"y\", \"x\"), np.random.randint(0, 10, (ds.sizes[\"y\"], ds.sizes[\"x\"]))\n", diff --git a/docs/examples/07_resampling.ipynb b/docs/examples/07_resampling.ipynb index 093e046b..bb9a5ce9 100644 --- a/docs/examples/07_resampling.ipynb +++ b/docs/examples/07_resampling.ipynb @@ -18,25 +18,17 @@ "metadata": {}, "outputs": [], "source": [ - "import nlmod\n", - "from nlmod import resample\n", + "\n", + "import geopandas as gpd\n", + "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", - "import flopy\n", - "import warnings\n", - "\n", - "\n", "from matplotlib.colors import Normalize\n", - "from matplotlib.patches import Polygon\n", - "from matplotlib.collections import PatchCollection\n", - "import matplotlib.pyplot as plt\n", - "\n", - "import geopandas as gpd\n", - "from shapely.geometry import LineString, Point\n", - "from shapely.geometry import Polygon as shp_polygon\n", "from scipy.interpolate import RectBivariateSpline\n", + "from shapely.geometry import LineString, Point\n", "\n", - "from IPython.display import display" + "import nlmod\n", + "from nlmod import resample" ] }, { @@ -100,7 +92,7 @@ "outputs": [], "source": [ "ds[\"data_nan\"] = ds[\"data\"].copy()\n", - "ds[\"data_nan\"].data[0, 1] = np.NaN\n", + "ds[\"data_nan\"].data[0, 1] = np.nan\n", "\n", "fig, ax = plt.subplots()\n", "ax.set_aspect(\"equal\")\n", @@ -144,7 +136,7 @@ "outputs": [], "source": [ "dsv[\"data_nan\"] = dsv[\"data\"].copy()\n", - "dsv[\"data_nan\"][7] = np.NaN\n", + "dsv[\"data_nan\"][7] = np.nan\n", "\n", "fig, ax = plt.subplots()\n", "ax.set_aspect(\"equal\")\n", diff --git a/docs/examples/08_gis.ipynb b/docs/examples/08_gis.ipynb index cde523a3..40c56731 100644 --- a/docs/examples/08_gis.ipynb +++ b/docs/examples/08_gis.ipynb @@ -19,14 +19,12 @@ "outputs": [], "source": [ "import os\n", + "\n", "import flopy\n", - "import geopandas as gpd\n", - "import matplotlib.pyplot as plt\n", - "import nlmod\n", - "import numpy as np\n", "import xarray as xr\n", - "from IPython.display import FileLink, FileLinks\n", - "from shapely.geometry import Polygon" + "from IPython.display import FileLink\n", + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb index 1444756c..1a381db5 100644 --- a/docs/examples/09_schoonhoven.ipynb +++ b/docs/examples/09_schoonhoven.ipynb @@ -18,18 +18,17 @@ "outputs": [], "source": [ "import os\n", + "\n", "import flopy\n", + "import geopandas as gpd\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", - "import nlmod\n", "import numpy as np\n", - "import xarray as xr\n", "import pandas as pd\n", - "import hydropandas as hpd\n", - "import geopandas as gpd\n", - "from nlmod.plot import DatasetCrossSection\n", "from shapely.geometry import LineString, Point\n", - "import warnings" + "\n", + "import nlmod\n", + "from nlmod.plot import DatasetCrossSection" ] }, { diff --git a/docs/examples/10_modpath.ipynb b/docs/examples/10_modpath.ipynb index c70458cf..0a4e77c9 100644 --- a/docs/examples/10_modpath.ipynb +++ b/docs/examples/10_modpath.ipynb @@ -28,9 +28,10 @@ "\n", "import flopy\n", "import matplotlib.pyplot as plt\n", - "import nlmod\n", "import numpy as np\n", - "import xarray as xr" + "import xarray as xr\n", + "\n", + "import nlmod" ] }, { @@ -164,7 +165,7 @@ " marker=\"o\",\n", " color=\"red\",\n", ")\n", - "ax.set_title(f\"pathlines\")\n", + "ax.set_title(\"pathlines\")\n", "ax.legend(loc=\"upper right\")" ] }, @@ -265,7 +266,7 @@ " marker=\"o\",\n", " color=\"red\",\n", " )\n", - " ax.set_title(f\"pathlines\")\n", + " ax.set_title(\"pathlines\")\n", " ax.legend(loc=\"upper right\")\n", "\n", " if i == 1:\n", @@ -370,7 +371,7 @@ " marker=\"o\",\n", " color=\"red\",\n", " )\n", - " ax.set_title(f\"pathlines\")\n", + " ax.set_title(\"pathlines\")\n", " ax.legend(loc=\"upper right\")\n", "\n", " if i == 1:\n", @@ -477,7 +478,7 @@ " marker=\"o\",\n", " color=\"red\",\n", " )\n", - " ax.set_title(f\"pathlines\")\n", + " ax.set_title(\"pathlines\")\n", " ax.legend(loc=\"upper right\")\n", "\n", " if i == 1:\n", diff --git a/docs/examples/11_grid_rotation.ipynb b/docs/examples/11_grid_rotation.ipynb index 7d60bdc0..61f3d8d0 100644 --- a/docs/examples/11_grid_rotation.ipynb +++ b/docs/examples/11_grid_rotation.ipynb @@ -30,11 +30,10 @@ "outputs": [], "source": [ "import os\n", + "\n", "import matplotlib\n", - "import nlmod\n", - "import pandas as pd\n", - "import warnings\n", - "import flopy" + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/12_layer_generation.ipynb b/docs/examples/12_layer_generation.ipynb index a3185d8b..f61f1e02 100644 --- a/docs/examples/12_layer_generation.ipynb +++ b/docs/examples/12_layer_generation.ipynb @@ -26,8 +26,9 @@ "metadata": {}, "outputs": [], "source": [ - "import nlmod\n", - "from shapely.geometry import LineString" + "from shapely.geometry import LineString\n", + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/13_plot_methods.ipynb b/docs/examples/13_plot_methods.ipynb index ebc79dcc..7c027527 100644 --- a/docs/examples/13_plot_methods.ipynb +++ b/docs/examples/13_plot_methods.ipynb @@ -36,12 +36,13 @@ "outputs": [], "source": [ "import os\n", + "\n", + "import flopy\n", "import matplotlib.pyplot as plt\n", "import xarray as xr\n", - "import flopy\n", + "\n", "import nlmod\n", - "from nlmod.plot import DatasetCrossSection\n", - "import warnings" + "from nlmod.plot import DatasetCrossSection" ] }, { diff --git a/docs/examples/14_stromingen_example.ipynb b/docs/examples/14_stromingen_example.ipynb index 10da127f..562b4575 100644 --- a/docs/examples/14_stromingen_example.ipynb +++ b/docs/examples/14_stromingen_example.ipynb @@ -33,11 +33,13 @@ "outputs": [], "source": [ "import os\n", + "\n", "import flopy as fp\n", "import geopandas as gpd\n", - "import nlmod\n", "from pandas import date_range\n", "\n", + "import nlmod\n", + "\n", "nlmod.util.get_color_logger(\"INFO\")\n", "print(f\"nlmod version: {nlmod.__version__}\")" ] diff --git a/docs/examples/15_geotop.ipynb b/docs/examples/15_geotop.ipynb index 13500e10..28584803 100644 --- a/docs/examples/15_geotop.ipynb +++ b/docs/examples/15_geotop.ipynb @@ -16,10 +16,11 @@ "metadata": {}, "outputs": [], "source": [ - "from shapely.geometry import LineString\n", + "import matplotlib\n", "import matplotlib.pyplot as plt\n", - "import nlmod\n", - "import matplotlib" + "from shapely.geometry import LineString\n", + "\n", + "import nlmod" ] }, { diff --git a/docs/examples/16_groundwater_transport.ipynb b/docs/examples/16_groundwater_transport.ipynb index 34d70dc7..d72ea7f8 100644 --- a/docs/examples/16_groundwater_transport.ipynb +++ b/docs/examples/16_groundwater_transport.ipynb @@ -21,11 +21,12 @@ "outputs": [], "source": [ "# import packages\n", - "import nlmod\n", - "import pandas as pd\n", - "import xarray as xr\n", "import flopy as fp\n", "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import xarray as xr\n", + "\n", + "import nlmod\n", "\n", "# set up pretty logging and show package versions\n", "nlmod.util.get_color_logger(\"INFO\")\n", @@ -575,7 +576,7 @@ " ax.set_xlabel(\"x [m]\")\n", " ax.set_ylabel(\"elevation [m NAP]\")\n", " # convert to pandas timestamp for prettier printing\n", - " ax.set_title(f\"time = {pd.Timestamp(c.time.isel(time=time_idx).values)}\");" + " ax.set_title(f\"time = {pd.Timestamp(c.time.isel(time=time_idx).values)}\")" ] }, { diff --git a/docs/examples/17_unsaturated_zone_flow.ipynb b/docs/examples/17_unsaturated_zone_flow.ipynb index b9500d4a..76b21f03 100644 --- a/docs/examples/17_unsaturated_zone_flow.ipynb +++ b/docs/examples/17_unsaturated_zone_flow.ipynb @@ -28,9 +28,11 @@ "source": [ "# import packages\n", "import os\n", + "\n", + "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", + "\n", "import nlmod\n", "\n", "# set up pretty logging and show package versions\n", @@ -47,7 +49,11 @@ "source": [ "# Ignore warning about 1d layers, in numpy 1.26.1 and flopy 3.4.3\n", "import warnings\n", - "warnings.filterwarnings(\"ignore\", message=\"Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future.\")" + "\n", + "warnings.filterwarnings(\n", + " \"ignore\",\n", + " message=\"Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future.\",\n", + ")" ] }, { @@ -94,7 +100,9 @@ "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, drop_layer_dim_from_top=False, cachename=\"regis.nc\", cachedir=cachedir)" + "regis = nlmod.read.regis.get_regis(\n", + " extent, drop_layer_dim_from_top=False, cachename=\"regis.nc\", cachedir=cachedir\n", + ")" ] }, { diff --git a/docs/examples/cache_example.py b/docs/examples/cache_example.py index 5d33d7b0..7ca76792 100644 --- a/docs/examples/cache_example.py +++ b/docs/examples/cache_example.py @@ -1,11 +1,12 @@ -import nlmod import numpy as np import xarray as xr +import nlmod + @nlmod.cache.cache_netcdf() def func_to_create_a_dataset(number): - """create a dataarray as an example for the caching method. + """Create a dataarray as an example for the caching method. Parameters ---------- diff --git a/docs/examples/generate_logo.py b/docs/examples/generate_logo.py index c423df20..87bbfcb1 100644 --- a/docs/examples/generate_logo.py +++ b/docs/examples/generate_logo.py @@ -4,11 +4,14 @@ @author: ruben """ + import os -import nlmod + import art_tools import matplotlib.pyplot as plt +import nlmod + filled = False n = 2 dx = 10_000 diff --git a/docs/examples/run_all_notebooks.py b/docs/examples/run_all_notebooks.py index 86f5101e..a1602a31 100644 --- a/docs/examples/run_all_notebooks.py +++ b/docs/examples/run_all_notebooks.py @@ -1,11 +1,13 @@ # %% +import os +import time from glob import glob + import nbformat -from nbconvert.preprocessors import ExecutePreprocessor -import time import numpy as np +from nbconvert.preprocessors import ExecutePreprocessor + import nlmod -import os logger = nlmod.util.get_color_logger("INFO") @@ -31,7 +33,7 @@ logger.info(f"Running {notebook} succeeded in {seconds} seconds") except Exception as exception: logger.error(f"Running notebook failed: {notebook}") - elapsed_time[notebook] = np.NaN + elapsed_time[notebook] = np.nan exceptions[notebook] = exception # save results in notebook diff --git a/nlmod/__init__.py b/nlmod/__init__.py index d3905157..5582151b 100644 --- a/nlmod/__init__.py +++ b/nlmod/__init__.py @@ -1,9 +1,4 @@ -# -*- coding: utf-8 -*- -"""Created on Thu Jan 7 12:13:44 2021. - -@author: oebbe -""" - +# ruff: noqa: F401 E402 import os NLMOD_DATADIR = os.path.join(os.path.dirname(__file__), "data") diff --git a/nlmod/cache.py b/nlmod/cache.py index 1aa3672f..0ed0d441 100644 --- a/nlmod/cache.py +++ b/nlmod/cache.py @@ -720,8 +720,17 @@ def ds_contains( datavars.append("tsmult") if attrs_ds: - # set by `nlmod.base.to_model_ds()` and `nlmod.base.set_ds_attrs()`, excluding "created_on" - attrs_ds_required = ["model_name", "mfversion", "exe_name", "model_ws", "figdir", "cachedir", "transport"] + # set by `nlmod.base.to_model_ds()` and `nlmod.base.set_ds_attrs()`, + # excluding "created_on" + attrs_ds_required = [ + "model_name", + "mfversion", + "exe_name", + "model_ws", + "figdir", + "cachedir", + "transport", + ] attrs.extend(attrs_ds_required) # User-friendly error messages if missing from ds diff --git a/nlmod/dims/__init__.py b/nlmod/dims/__init__.py index d5892790..18c147c6 100644 --- a/nlmod/dims/__init__.py +++ b/nlmod/dims/__init__.py @@ -1,3 +1,4 @@ +# ruff: noqa: F401 F403 from . import base, grid, layers, resample, time from .attributes_encodings import * from .base import * diff --git a/nlmod/dims/attributes_encodings.py b/nlmod/dims/attributes_encodings.py index b93587fc..5b4e30b0 100644 --- a/nlmod/dims/attributes_encodings.py +++ b/nlmod/dims/attributes_encodings.py @@ -228,14 +228,13 @@ def get_encodings( def compute_scale_and_offset(min_value, max_value): - """ - Reduce precision of the dataset by storing it as int16 and thereby reducing the precision. - - Computes the scale_factor and offset for the dataset using a min_value and max_value to - transform the range of the dataset to the range of valid int16 values. The packed value - is computed as: + """Reduce precision of the dataset by storing it as int16. + + Computes the scale_factor and offset for the dataset using a min_value and max_value + to transform the range of the dataset to the range of valid int16 values. The packed + value is computed as: packed_value = (unpacked_value - add_offset) / scale_factor - + Parameters ---------- min_value : float @@ -259,8 +258,9 @@ def compute_scale_and_offset(min_value, max_value): def is_int16_allowed(vmin, vmax, dval_max): - """compute the loss of resolution by storing a float as int16 (`dval`). And - compare it with the maximum allowed loss of resolution (`dval_max`). + """Compute the loss of resolution by storing a float as int16 (`dval`). + + Compare it with the maximum allowed loss of resolution (`dval_max`). Parameters ---------- @@ -274,7 +274,8 @@ def is_int16_allowed(vmin, vmax, dval_max): Returns ------- bool - True if the loss of resolution is allowed, False otherwise""" + True if the loss of resolution is allowed, False otherwise + """ nsteps = 32766 + 32767 dval = (vmax - vmin) / nsteps return dval <= dval_max diff --git a/nlmod/dims/base.py b/nlmod/dims/base.py index c2b729b6..a56abe85 100644 --- a/nlmod/dims/base.py +++ b/nlmod/dims/base.py @@ -44,7 +44,6 @@ def set_ds_attrs( ds : xarray dataset model dataset. """ - if model_name is not None and len(model_name) > 16 and mfversion == "mf6": raise ValueError("model_name can not have more than 16 characters") ds.attrs["model_name"] = model_name @@ -338,7 +337,6 @@ def _get_structured_grid_ds( dictionary, and a coordinate reference system specified by `crs`, if provided. """ - if attrs is None: attrs = {} attrs.update({"gridtype": "structured"}) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 2aadff70..332c6cbc 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -39,19 +39,17 @@ affine_transform_gdf, get_affine_mod_to_world, get_affine_world_to_mod, - structured_da_to_ds, - get_delr, get_delc, + get_delr, + structured_da_to_ds, ) logger = logging.getLogger(__name__) def snap_extent(extent, delr, delc): - """ - snap the extent in such a way that an integer number of columns and rows fit - in the extent. The new extent is always equal to, or bigger than the - original extent. + """Snap the extent in such a way that an integer number of columns and rows fit in + the extent. The new extent is always equal to, or bigger than the original extent. Parameters ---------- @@ -96,7 +94,7 @@ def snap_extent(extent, delr, delc): def xy_to_icell2d(xy, ds): - """get the icell2d value of a point defined by its x and y coordinates. + """Get the icell2d value of a point defined by its x and y coordinates. Parameters ---------- @@ -118,7 +116,7 @@ def xy_to_icell2d(xy, ds): def xy_to_row_col(xy, ds): - """get the row and column values of a point defined by its x and y coordinates. + """Get the row and column values of a point defined by its x and y coordinates. Parameters ---------- @@ -142,7 +140,7 @@ def xy_to_row_col(xy, ds): def xyz_to_cid(xyz, ds=None, modelgrid=None): - """get the icell2d value of a point defined by its x and y coordinates. + """Get the icell2d value of a point defined by its x and y coordinates. Parameters ---------- @@ -231,7 +229,6 @@ def modelgrid_from_ds(ds, rotated=True, nlay=None, top=None, botm=None, **kwargs def modelgrid_to_vertex_ds(mg, ds, nodata=-1): """Add information about the calculation-grid to a model dataset.""" - warnings.warn( "'modelgrid_to_vertex_ds' is deprecated and will be removed in a" "future version, please use 'modelgrid_to_ds' instead", @@ -345,8 +342,9 @@ def get_dims_coords_from_modelgrid(mg): def gridprops_to_vertex_ds(gridprops, ds, nodata=-1): - """Gridprops is a dictionary containing keyword arguments needed to - generate a flopy modelgrid instance.""" + """Gridprops is a dictionary containing keyword arguments needed to generate a flopy + modelgrid instance. + """ _, xv, yv = zip(*gridprops["vertices"]) ds["xv"] = ("iv", np.array(xv)) ds["yv"] = ("iv", np.array(yv)) @@ -520,8 +518,8 @@ def refine( def ds_to_gridprops(ds_in, gridprops, method="nearest", icvert_nodata=-1): - """resample a dataset (xarray) on an structured grid to a new dataset with - a vertex grid. + """Resample a dataset (xarray) on an structured grid to a new dataset with a vertex + grid. Returns a dataset with resampled variables and the untouched variables. @@ -577,7 +575,7 @@ def ds_to_gridprops(ds_in, gridprops, method="nearest", icvert_nodata=-1): # add other variables for not_interp_var in not_interp_vars: ds_out[not_interp_var] = structured_da_to_ds( - da=ds_in[not_interp_var], ds=ds_out, method=method, nodata=np.NaN + da=ds_in[not_interp_var], ds=ds_out, method=method, nodata=np.nan ) has_rotation = "angrot" in ds_out.attrs and ds_out.attrs["angrot"] != 0.0 if has_rotation: @@ -601,8 +599,7 @@ def ds_to_gridprops(ds_in, gridprops, method="nearest", icvert_nodata=-1): def get_xyi_icell2d(gridprops=None, ds=None): - """Get x and y coordinates of the cell mids from the cellids in the grid - properties. + """Get x and y coordinates of the cell mids from the cellids in the grid properties. Parameters ---------- @@ -634,9 +631,9 @@ def get_xyi_icell2d(gridprops=None, ds=None): 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. + """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 ---------- @@ -723,7 +720,6 @@ def col_to_list(col_in, ds, cellids): col_lst : list raster values from ds presented in a list per cell. """ - if isinstance(col_in, str) and col_in in ds: col_in = ds[col_in] if isinstance(col_in, xr.DataArray): @@ -951,7 +947,6 @@ def cols_to_reclist(ds, cellids, *args, cellid_column=0): cellid_column : int, optional Adds the cellid ((layer, row, col) or (layer, icell2d)) to the reclist in this column number. Do not add cellid when cellid_column is None. The default is 0. - """ cols = [col_to_list(col, ds, cellids) for col in args] if cellid_column is not None: @@ -1079,8 +1074,7 @@ def da_to_reclist( def polygon_to_area(modelgrid, polygon, da, gridtype="structured"): - """create a grid with the surface area in each cell based on a polygon - value. + """Create a grid with the surface area in each cell based on a polygon value. Parameters ---------- @@ -1126,8 +1120,8 @@ def polygon_to_area(modelgrid, polygon, da, gridtype="structured"): def gdf_to_data_array_struc( gdf, gwf, field="VALUE", agg_method=None, interp_method=None ): - """Project vector data on a structured grid. Aggregate data if multiple - geometries are in a single cell. + """Project vector data on a structured grid. Aggregate data if multiple geometries + are in a single cell. Parameters ---------- @@ -1192,11 +1186,11 @@ def gdf_to_data_array_struc( def gdf_to_da( - gdf, ds, column, agg_method=None, fill_value=np.NaN, min_total_overlap=0.0, ix=None + gdf, ds, column, agg_method=None, fill_value=np.nan, min_total_overlap=0.0, ix=None ): - """Project vector data on a grid. Aggregate data if multiple - geometries are in a single cell. Supports structured and vertex grids. - This method replaces gdf_to_data_array_struc. + """Project vector data on a grid. Aggregate data if multiple geometries are in a + single cell. Supports structured and vertex grids. This method replaces + gdf_to_data_array_struc. Parameters ---------- @@ -1273,7 +1267,7 @@ def gdf_to_da( def interpolate_gdf_to_array(gdf, gwf, field="values", method="nearest"): - """interpolate data from a point gdf. + """Interpolate data from a point gdf. Parameters ---------- @@ -1462,9 +1456,9 @@ def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None): def gdf_to_bool_da(gdf, ds, ix=None, buffer=0.0, **kwargs): - """convert a GeoDataFrame with polygon geometries into a data array - corresponding to the modelgrid in which each cell is 1 (True) if one or - more geometries are (partly) in that cell. + """Convert a GeoDataFrame with polygon geometries into a data array corresponding to + the modelgrid in which each cell is 1 (True) if one or more geometries are (partly) + in that cell. Parameters ---------- @@ -1489,9 +1483,9 @@ def gdf_to_bool_da(gdf, ds, ix=None, buffer=0.0, **kwargs): def gdf_to_bool_ds(gdf, ds, da_name, keep_coords=None, ix=None, buffer=0.0, **kwargs): - """convert a GeoDataFrame with polygon geometries into a model dataset with - a data_array named 'da_name' in which each cell is 1 (True) if one or more - geometries are (partly) in that cell. + """Convert a GeoDataFrame with polygon geometries into a model dataset with a + data_array named 'da_name' in which each cell is 1 (True) if one or more geometries + are (partly) in that cell. Parameters ---------- @@ -1703,7 +1697,7 @@ def gdf_to_grid( def get_thickness_from_topbot(top, bot): - """get thickness from data arrays with top and bots. + """Get thickness from data arrays with top and bots. Parameters ---------- @@ -1745,9 +1739,9 @@ def get_thickness_from_topbot(top, bot): def get_vertices_arr(ds, modelgrid=None, vert_per_cid=4, epsilon=0, rotated=False): - """get vertices of a vertex modelgrid from a ds or the modelgrid. Only - return the 4 corners of each cell and not the corners of adjacent cells - thus limiting the vertices per cell to 4 points. + """Get vertices of a vertex modelgrid from a ds or the modelgrid. Only return the 4 + corners of each cell and not the corners of adjacent cells thus limiting the + vertices per cell to 4 points. This method uses the xvertices and yvertices attributes of the modelgrid. When no modelgrid is supplied, a modelgrid-object is created from ds. @@ -1809,9 +1803,9 @@ def get_vertices_arr(ds, modelgrid=None, vert_per_cid=4, epsilon=0, rotated=Fals def get_vertices(ds, vert_per_cid=4, epsilon=0, rotated=False): - """get vertices of a vertex modelgrid from a ds or the modelgrid. Only - return the 4 corners of each cell and not the corners of adjacent cells - thus limiting the vertices per cell to 4 points. + """Get vertices of a vertex modelgrid from a ds or the modelgrid. Only return the 4 + corners of each cell and not the corners of adjacent cells thus limiting the + vertices per cell to 4 points. This method uses the xvertices and yvertices attributes of the modelgrid. When no modelgrid is supplied, a modelgrid-object is created from ds. @@ -1842,7 +1836,6 @@ def get_vertices(ds, vert_per_cid=4, epsilon=0, rotated=False): vertices_da : xarray DataArray Vertex coördinates per cell with dimensions(cid, no_vert, 2). """ - # obtain vertices_arr = get_vertices_arr( @@ -1863,8 +1856,8 @@ def get_vertices(ds, vert_per_cid=4, epsilon=0, rotated=False): @cache.cache_netcdf(coords_2d=True) def mask_model_edge(ds, idomain=None): - """get data array which is 1 for every active cell (defined by idomain) at - the boundaries of the model (xmin, xmax, ymin, ymax). Other cells are 0. + """Get data array which is 1 for every active cell (defined by idomain) at the + boundaries of the model (xmin, xmax, ymin, ymax). Other cells are 0. Parameters ---------- @@ -1940,7 +1933,6 @@ def polygons_from_model_ds(model_ds): polygons : list of shapely Polygons list with polygon of each raster cell. """ - if model_ds.gridtype == "structured": delr = get_delr(model_ds) delc = get_delc(model_ds) diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py index a3c44d7b..c468582f 100644 --- a/nlmod/dims/layers.py +++ b/nlmod/dims/layers.py @@ -49,7 +49,7 @@ def calculate_thickness(ds, top="top", bot="botm"): raise ValueError("2d top should have same last dimension as bot") # subtracting floats can result in rounding errors. Mainly anoying for zero thickness layers. - thickness = thickness.where(~np.isclose(thickness, 0.), 0.) + thickness = thickness.where(~np.isclose(thickness, 0.0), 0.0) if isinstance(ds[bot], xr.DataArray): thickness.name = "thickness" @@ -69,8 +69,8 @@ def calculate_thickness(ds, top="top", bot="botm"): def calculate_transmissivity( ds, kh="kh", thickness="thickness", top="top", botm="botm" ): - """calculate the transmissivity (T) as the product of the horizontal - conductance (kh) and the thickness (D). + """Calculate the transmissivity (T) as the product of the horizontal conductance + (kh) and the thickness (D). Parameters ---------- @@ -94,7 +94,6 @@ def calculate_transmissivity( T : xarray.DataArray DataArray containing transmissivity (T). NaN where layer thickness is zero """ - if thickness in ds: thickness = ds[thickness] else: @@ -123,7 +122,7 @@ def calculate_transmissivity( def calculate_resistance(ds, kv="kv", thickness="thickness", top="top", botm="botm"): - """calculate vertical resistance (c) between model layers from the vertical + """Calculate vertical resistance (c) between model layers from the vertical conductivity (kv) and the thickness. The resistance between two layers is assigned to the top layer. The bottom model layer gets a resistance of infinity. @@ -149,7 +148,6 @@ def calculate_resistance(ds, kv="kv", thickness="thickness", top="top", botm="bo c : xarray.DataArray DataArray containing vertical resistance (c). NaN where layer thickness is zero """ - if thickness in ds: thickness = ds[thickness] else: @@ -224,7 +222,6 @@ def split_layers_ds( Dataset with new tops and bottoms taking into account split layers, and filled data for other variables. """ - layers = list(ds.layer.data) # Work on a shallow copy of split_dict @@ -569,7 +566,6 @@ def combine_layers_ds( Dataset with new tops and bottoms taking into account combined layers, and recalculated values for parameters (kh, kv, kD, c). """ - data_vars = [] for dv in [kh, kv, kD, c]: if dv is not None: @@ -644,7 +640,7 @@ def combine_layers_ds( def add_kh_kv_from_ml_layer_to_ds( ml_layer_ds, ds, anisotropy, fill_value_kh, fill_value_kv ): - """add kh and kv from a model layer dataset to the model dataset. + """Add kh and kv from a model layer dataset to the model dataset. Supports structured and vertex grids. @@ -811,8 +807,9 @@ def set_layer_thickness(ds, layer, thickness, change="botm", copy=True): def set_minimum_layer_thickness(ds, layer, min_thickness, change="botm", copy=True): - """Make sure layer has a minimum thickness by lowering the botm of the - layer where neccesary.""" + """Make sure layer has a minimum thickness by lowering the botm of the layer where + neccesary. + """ assert layer in ds.layer assert change == "botm", "Only change=botm allowed for now" if copy: @@ -834,8 +831,8 @@ def set_minimum_layer_thickness(ds, layer, min_thickness, change="botm", copy=Tr def remove_thin_layers( ds, min_thickness=0.1, update_thickness_every_layer=False, copy=True ): - """ - Remove cells with a thickness less than min_thickness (setting the thickness to 0) + """Remove cells with a thickness less than min_thickness (setting the thickness to + 0) The thickness of the removed cells is added to the first active layer below @@ -855,11 +852,11 @@ def remove_thin_layers( copy : bool, optional If copy=True, data in the return value is always copied, so the original Dataset is not altered. The default is True. + Returns ------- ds : xr.Dataset Dataset containing information about layers. - """ if "layer" in ds["top"].dims: msg = "remove_thin_layers does not support top with a layer dimension" @@ -904,8 +901,8 @@ def remove_thin_layers( def get_kh_kv(kh, kv, anisotropy, fill_value_kh=1.0, fill_value_kv=0.1, idomain=None): - """create kh en kv grid data for flopy from existing kh, kv and anistropy - grids with nan values (typically from REGIS). + """Create kh en kv grid data for flopy from existing kh, kv and anistropy grids with + nan values (typically from REGIS). fill nans in kh grid in these steps: 1. take kv and multiply by anisotropy, if this is nan: @@ -1018,7 +1015,6 @@ def fill_top_bot_kh_kv_at_mask(ds, fill_mask): ds : xr.DataSet model dataset with adjusted data variables: 'top', 'botm', 'kh', 'kv' """ - # zee cellen hebben altijd een top gelijk aan 0 ds["top"].values = np.where(fill_mask, 0, ds["top"]) @@ -1065,7 +1061,6 @@ def fill_nan_top_botm_kh_kv( 2. Remove inactive layers, with no positive thickness anywhere 3. Compute kh and kv, filling nans with anisotropy or fill_values """ - # 1 ds = remove_layer_dim_from_top(ds) @@ -1088,8 +1083,7 @@ def fill_nan_top_botm_kh_kv( def fill_nan_top_botm(ds): - """ - Remove Nans in non-existent layers in botm and top variables + """Remove Nans in non-existent layers in botm and top variables. The NaNs are removed by setting the value to the top and botm of higher/lower layers that do exist. @@ -1122,8 +1116,7 @@ def fill_nan_top_botm(ds): def set_nan_top_and_botm(ds, copy=True): - """ - Sets Nans for non-existent layers in botm and top variables + """Sets Nans for non-existent layers in botm and top variables. Nans are only added to top when it contains a layer dimension. @@ -1159,8 +1152,7 @@ def remove_layer_dim_from_top( return_inconsistencies=False, copy=True, ): - """ - Change top from 3d to 2d, removing NaNs in top and botm in the process. + """Change top from 3d to 2d, removing NaNs in top and botm in the process. This method sets variable `top` to the top of the upper layer (like the definition in MODFLOW). This removes redundant data, as the top of all layers exept the most @@ -1220,8 +1212,7 @@ def remove_layer_dim_from_top( def add_layer_dim_to_top(ds, set_non_existing_layers_to_nan=True, copy=True): - """ - Change top from 2d to 3d, setting top and botm to NaN for non-existent layers. + """Change top from 2d to 3d, setting top and botm to NaN for non-existent layers. Parameters ---------- @@ -1248,28 +1239,23 @@ def add_layer_dim_to_top(ds, set_non_existing_layers_to_nan=True, copy=True): def convert_to_modflow_top_bot(ds, **kwargs): - """ - Removes the layer dimension from top and fills nans in top and botm. + """Removes the layer dimension from top and fills nans in top and botm. Alias to remove_layer_dim_from_top - """ ds = remove_layer_dim_from_top(ds, **kwargs) def convert_to_regis_top_bot(ds, **kwargs): - """ - Adds a layer dimension to top and sets non-existing cells to nan in top and botm. + """Adds a layer dimension to top and sets non-existing cells to nan in top and botm. Alias to add_layer_dim_to_top - """ ds = add_layer_dim_to_top(ds, **kwargs) def remove_inactive_layers(ds): - """ - Remove layers which only contain inactive cells + """Remove layers which only contain inactive cells. Parameters ---------- @@ -1280,7 +1266,6 @@ def remove_inactive_layers(ds): ------- ds : xr.Dataset The model Dataset without inactive layers. - """ idomain = get_idomain(ds) # only keep layers with at least one active cell @@ -1314,8 +1299,9 @@ def get_idomain(ds): idomain.attrs.clear() # set idomain of cells with a positive thickness to 1 thickness = calculate_thickness(ds) - # subtracting floats can result in rounding errors. Mainly anoying for zero thickness layers. - thickness = thickness.where(~np.isclose(thickness, 0.), 0.) + # subtracting floats can result in rounding errors. Mainly anoying for zero + # thickness layers. + thickness = thickness.where(~np.isclose(thickness, 0.0), 0.0) idomain.data[thickness.data > 0.0] = 1 # set idomain above/below the first/last active layer to 0 idomain.data[idomain.where(idomain > 0).ffill(dim="layer").isnull()] = 0 @@ -1347,7 +1333,7 @@ def get_first_active_layer(ds, **kwargs): def get_first_active_layer_from_idomain(idomain, nodata=-999): - """get the first (top) active layer in each cell from the idomain. + """Get the first (top) active layer in each cell from the idomain. Parameters ---------- @@ -1377,7 +1363,7 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999): def get_last_active_layer_from_idomain(idomain, nodata=-999): - """get the last (bottom) active layer in each cell from the idomain. + """Get the last (bottom) active layer in each cell from the idomain. Parameters ---------- @@ -1444,7 +1430,7 @@ def get_layer_of_z(ds, z, above_model=-999, below_model=-999): def update_idomain_from_thickness(idomain, thickness, mask): - """get new idomain from thickness in the cells where mask is 1 (or True). + """Get new idomain from thickness in the cells where mask is 1 (or True). Idomain becomes: - 1: if cell thickness is bigger than 0 @@ -1517,7 +1503,7 @@ def aggregate_by_weighted_mean_to_ds(ds, source_ds, var_name): ValueError if source_ds does not have a layer dimension - See also + See Also -------- nlmod.read.geotop.aggregate_to_ds """ @@ -1559,7 +1545,7 @@ def aggregate_by_weighted_mean_to_ds(ds, source_ds, var_name): def check_elevations_consistency(ds): if "layer" in ds["top"].dims: tops = ds["top"].data - top_ref = np.full(tops.shape[1:], np.NaN) + top_ref = np.full(tops.shape[1:], np.nan) for lay, layer in zip(range(tops.shape[0]), ds.layer.data): top = tops[lay] mask = ~np.isnan(top) @@ -1572,7 +1558,7 @@ def check_elevations_consistency(ds): top_ref[mask] = top[mask] bots = ds["botm"].data - bot_ref = np.full(bots.shape[1:], np.NaN) + bot_ref = np.full(bots.shape[1:], np.nan) for lay, layer in zip(range(bots.shape[0]), ds.layer.data): bot = bots[lay] mask = ~np.isnan(bot) @@ -1591,8 +1577,7 @@ def check_elevations_consistency(ds): def insert_layer(ds, name, top, bot, kh=None, kv=None, copy=True): - """ - Inserts a layer in a model Dataset, burning it in an existing layer model. + """Inserts a layer in a model Dataset, burning it in an existing layer model. This method loops over the existing layers, and checks if (part of) the new layer needs to be inserted above the existing layer, and if the top or bottom of the @@ -1652,7 +1637,6 @@ def insert_layer(ds, name, top, bot, kh=None, kv=None, copy=True): ------- ds : xarray.Dataset xarray Dataset containing the new layer(s) - """ shape = ds["botm"].shape[1:] assert top.shape == shape diff --git a/nlmod/dims/rdp.py b/nlmod/dims/rdp.py index 3b24b819..df996789 100644 --- a/nlmod/dims/rdp.py +++ b/nlmod/dims/rdp.py @@ -1,10 +1,11 @@ """ -rdp -~~~ +rdp +~~~ Python implementation of the Ramer-Douglas-Peucker algorithm. -:copyright: 2014-2016 Fabian Hirschmann +:copyright: 2014-2016 Fabian Hirschmann :license: MIT, see LICENSE.txt for more details. """ + import sys from functools import partial @@ -15,8 +16,8 @@ def pldist(point, start, end): - """Calculates the distance from ``point`` to the line given by the points - ``start`` and ``end``. + """Calculates the distance from ``point`` to the line given by the points ``start`` + and ``end``. :param point: a point :type point: numpy array @@ -113,9 +114,8 @@ def rdp_iter(M, epsilon, dist=pldist, return_mask=False): def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False): - """ - Simplifies a given array of points using the Ramer-Douglas-Peucker - algorithm. + """Simplifies a given array of points using the Ramer-Douglas-Peucker algorithm. + Example: >>> from rdp import rdp >>> rdp([[1, 1], [2, 2], [3, 3], [4, 4]]) @@ -153,7 +153,6 @@ def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False): :param return_mask: return mask instead of simplified array :type return_mask: bool """ - if algo == "iter": algo = partial(rdp_iter, return_mask=return_mask) elif algo == "rec": diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py index 01725ca7..f6fbf843 100644 --- a/nlmod/dims/resample.py +++ b/nlmod/dims/resample.py @@ -16,8 +16,7 @@ def get_xy_mid_structured(extent, delr, delc, descending_y=True): - """Calculates the x and y coordinates of the cell centers of a structured - grid. + """Calculates the x and y coordinates of the cell centers of a structured grid. Parameters ---------- @@ -154,8 +153,8 @@ def ds_to_structured_grid( angrot=0.0, method="nearest", ): - """Resample a dataset (xarray) from a structured grid to a new dataset from - a different structured grid. + """Resample a dataset (xarray) from a structured grid to a new dataset from a + different structured grid. Parameters ---------- @@ -193,7 +192,6 @@ def ds_to_structured_grid( dataset with dimensions (layer, y, x). y and x are from the new grid. """ - assert isinstance(ds_in, xr.core.dataset.Dataset) if hasattr(ds_in, "gridtype"): assert ds_in.attrs["gridtype"] == "structured" @@ -226,8 +224,7 @@ def ds_to_structured_grid( def _set_angrot_attributes(extent, xorigin, yorigin, angrot, attrs): - """Internal method to set the properties of the grid in an attribute - dictionary. + """Internal method to set the properties of the grid in an attribute dictionary. Parameters ---------- @@ -285,7 +282,7 @@ def _set_angrot_attributes(extent, xorigin, yorigin, angrot, attrs): def fillnan_da_structured_grid(xar_in, method="nearest"): - """fill not-a-number values in a structured grid, DataArray. + """Fill not-a-number values in a structured grid, DataArray. The fill values are determined using the 'nearest' method of the scipy.interpolate.griddata function @@ -345,7 +342,7 @@ def fillnan_da_structured_grid(xar_in, method="nearest"): def fillnan_da_vertex_grid(xar_in, ds=None, x=None, y=None, method="nearest"): - """fill not-a-number values in a vertex grid, DataArray. + """Fill not-a-number values in a vertex grid, DataArray. The fill values are determined using the 'nearest' method of the scipy.interpolate.griddata function @@ -375,7 +372,6 @@ def fillnan_da_vertex_grid(xar_in, ds=None, x=None, y=None, method="nearest"): ----- can be slow if the xar_in is a large raster """ - if xar_in.dims != ("icell2d",): raise ValueError( f"expected dataarray with dimensions ('icell2d'), got dimensions -> {xar_in.dims}" @@ -407,7 +403,7 @@ def fillnan_da_vertex_grid(xar_in, ds=None, x=None, y=None, method="nearest"): def fillnan_da(da, ds=None, method="nearest"): - """fill not-a-number values in a DataArray. + """Fill not-a-number values in a DataArray. The fill values are determined using the 'nearest' method of the scipy.interpolate.griddata function @@ -457,7 +453,6 @@ def vertex_da_to_ds(da, ds, method="nearest"): xarray.DataArray A DataArray, with the same gridtype as ds. """ - if "icell2d" not in da.dims: return structured_da_to_ds(da, ds, method=method) points = np.array((da.x.data, da.y.data)).T @@ -515,7 +510,7 @@ def dim_to_regular_dim(da, dims, z): return xr.DataArray(z, dims=dims, coords=coords) -def structured_da_to_ds(da, ds, method="average", nodata=np.NaN): +def structured_da_to_ds(da, ds, method="average", nodata=np.nan): """Resample a DataArray to the coordinates of a model dataset. Parameters diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py index f1b21c17..0af4d927 100644 --- a/nlmod/dims/time.py +++ b/nlmod/dims/time.py @@ -75,7 +75,6 @@ def set_ds_time_deprecated( ds : xarray.Dataset dataset with time variant model data """ - warnings.warn( "this function is deprecated and will eventually be removed, " "please use nlmod.time.set_ds_time() in the future.", @@ -200,7 +199,6 @@ def set_ds_time( ------- ds : xarray.Dataset model dataset with added time coordinate - """ logger.info( "Function set_ds_time() has changed since nlmod version 0.7." @@ -280,7 +278,6 @@ def set_ds_time( def ds_time_idx_from_tdis_settings(start, perlen, nstp=1, tsmult=1.0, time_units="D"): """Get time index from TDIS perioddata: perlen, nstp, tsmult. - Parameters ---------- start : str, pd.Timestamp @@ -356,7 +353,6 @@ def estimate_nstp( if `return_dt_arr` is `True` returns the durations of the timesteps corresponding with the returned nstp. """ - nt = len(forcing) # Scaled linear between min and max. array nstp will be modified along the way @@ -409,8 +405,7 @@ def estimate_nstp( def get_time_step_length(perlen, nstp, tsmult): - """ - Get the length of the timesteps within a singe stress-period. + """Get the length of the timesteps within a singe stress-period. Parameters ---------- @@ -426,7 +421,6 @@ def get_time_step_length(perlen, nstp, tsmult): t : np.ndarray An array with the length of each of the timesteps within the stress period, in the same unit as perlen. - """ t = np.array([tsmult**x for x in range(nstp)]) t = t * perlen / t.sum() @@ -455,7 +449,6 @@ def ds_time_idx_from_model(gwf): IndexVariable time coordinate for xarray data-array or dataset """ - return ds_time_idx_from_modeltime(gwf.modeltime) @@ -481,7 +474,6 @@ def ds_time_idx_from_modeltime(modeltime): IndexVariable time coordinate for xarray data-array or dataset """ - return ds_time_idx( np.cumsum(modeltime.perlen), start_datetime=modeltime.start_datetime, diff --git a/nlmod/epsg28992.py b/nlmod/epsg28992.py index 6429fd9b..665cb1a7 100644 --- a/nlmod/epsg28992.py +++ b/nlmod/epsg28992.py @@ -1,6 +1,6 @@ """ NOTE: this is the correct epsg:28992 definition for plotting backgroundmaps in RD -More info (in Dutch) here: +More info (in Dutch) here: https://qgis.nl/2011/12/05/epsg28992-of-rijksdriehoekstelsel-verschuiving/ This was still a problem in October 2023 """ diff --git a/nlmod/gis.py b/nlmod/gis.py index 17d4709c..7c835e55 100644 --- a/nlmod/gis.py +++ b/nlmod/gis.py @@ -5,8 +5,8 @@ import numpy as np from .dims.grid import polygons_from_model_ds -from .dims.resample import get_affine_mod_to_world from .dims.layers import calculate_thickness +from .dims.resample import get_affine_mod_to_world logger = logging.getLogger(__name__) @@ -220,7 +220,6 @@ def ds_to_vector_file( fnames : str or list of str filename(s) of exported geopackage or shapefiles. """ - # get default combination dictionary if combine_dic is None: combine_dic = { @@ -450,7 +449,6 @@ def _break_down_dimension( Copied and altered from imod-python. """ - keep_vars = [] for var in variables: if dim in ds[var].dims: diff --git a/nlmod/gwf/__init__.py b/nlmod/gwf/__init__.py index 5fa031b0..abdc270b 100644 --- a/nlmod/gwf/__init__.py +++ b/nlmod/gwf/__init__.py @@ -1,3 +1,4 @@ +# ruff: noqa: F401 F403 from . import output, surface_water, wells from .gwf import * from .horizontal_flow_barrier import * diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py index ce45e29c..e234b5a0 100644 --- a/nlmod/gwf/gwf.py +++ b/nlmod/gwf/gwf.py @@ -1,8 +1,3 @@ -# -*- coding: utf-8 -*- -"""Created on Thu Jan 7 17:20:34 2021. - -@author: oebbe -""" import logging import numbers import warnings @@ -11,8 +6,8 @@ import xarray as xr from ..dims import grid -from ..dims.resample import get_delr, get_delc from ..dims.layers import get_idomain +from ..dims.resample import get_delc, get_delr from ..sim import ims, sim, tdis from ..util import _get_value_from_ds_attr, _get_value_from_ds_datavar from . import recharge @@ -21,7 +16,7 @@ def gwf(ds, sim, under_relaxation=False, **kwargs): - """create groundwater flow model from the model dataset. + """Create groundwater flow model from the model dataset. Parameters ---------- @@ -37,7 +32,6 @@ def gwf(ds, sim, under_relaxation=False, **kwargs): gwf : flopy ModflowGwf groundwaterflow object. """ - # start creating model logger.info("creating mf6 GWF") @@ -63,7 +57,7 @@ def gwf(ds, sim, under_relaxation=False, **kwargs): def dis(ds, gwf, length_units="METERS", pname="dis", **kwargs): - """create discretisation package from the model dataset. + """Create discretisation package from the model dataset. Parameters ---------- @@ -85,7 +79,7 @@ def dis(ds, gwf, length_units="METERS", pname="dis", **kwargs): def _dis(ds, model, length_units="METERS", pname="dis", **kwargs): - """create discretisation package from the model dataset. + """Create discretisation package from the model dataset. Parameters ---------- @@ -164,7 +158,7 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs): def disv(ds, gwf, length_units="METERS", pname="disv", **kwargs): - """create discretisation vertices package from the model dataset. + """Create discretisation vertices package from the model dataset. Parameters ---------- @@ -186,7 +180,7 @@ def disv(ds, gwf, length_units="METERS", pname="disv", **kwargs): def _disv(ds, model, length_units="METERS", pname="disv", **kwargs): - """create discretisation vertices package from the model dataset. + """Create discretisation vertices package from the model dataset. Parameters ---------- @@ -267,7 +261,7 @@ def _disv(ds, model, length_units="METERS", pname="disv", **kwargs): def npf( ds, gwf, k="kh", k33="kv", icelltype=0, save_flows=False, pname="npf", **kwargs ): - """create node property flow package from model dataset. + """Create node property flow package from model dataset. Parameters ---------- @@ -331,7 +325,7 @@ def ghb( layer=None, **kwargs, ): - """create general head boundary from model dataset. + """Create general head boundary from model dataset. Parameters ---------- @@ -414,7 +408,7 @@ def drn( layer=None, **kwargs, ): - """create drain from model dataset. + """Create drain from model dataset. Parameters ---------- @@ -487,7 +481,7 @@ def riv( layer=None, **kwargs, ): - """create river package from model dataset. + """Create river package from model dataset. Parameters ---------- @@ -570,7 +564,7 @@ def chd( layer=0, **kwargs, ): - """create constant head package from model dataset. + """Create constant head package from model dataset. Parameters ---------- @@ -645,7 +639,7 @@ def chd( def ic(ds, gwf, starting_head="starting_head", pname="ic", **kwargs): - """create initial condictions package from model dataset. + """Create initial condictions package from model dataset. Parameters ---------- @@ -689,7 +683,7 @@ def sto( pname="sto", **kwargs, ): - """create storage package from model dataset. + """Create storage package from model dataset. Parameters ---------- @@ -741,8 +735,7 @@ def sto( def surface_drain_from_ds(ds, gwf, resistance, elev="ahn", pname="drn", **kwargs): - """create surface level drain (maaivelddrainage in Dutch) from the model - dataset. + """Create surface level drain (maaivelddrainage in Dutch) from the model dataset. Parameters ---------- @@ -765,7 +758,6 @@ def surface_drain_from_ds(ds, gwf, resistance, elev="ahn", pname="drn", **kwargs drn : flopy ModflowGwfdrn drn package """ - ds.attrs["surface_drn_resistance"] = resistance maskarr = _get_value_from_ds_datavar(ds, "elev", elev, return_da=True) @@ -793,7 +785,7 @@ def surface_drain_from_ds(ds, gwf, resistance, elev="ahn", pname="drn", **kwargs def rch(ds, gwf, pname="rch", **kwargs): - """create recharge package from model dataset. + """Create recharge package from model dataset. Parameters ---------- @@ -817,7 +809,7 @@ def rch(ds, gwf, pname="rch", **kwargs): def evt(ds, gwf, pname="evt", **kwargs): - """create evapotranspiration package from model dataset. + """Create evapotranspiration package from model dataset. Parameters ---------- @@ -842,7 +834,7 @@ def evt(ds, gwf, pname="evt", **kwargs): def uzf(ds, gwf, pname="uzf", **kwargs): - """create unsaturated zone flow package from model dataset. + """Create unsaturated zone flow package from model dataset. Parameters ---------- @@ -886,7 +878,8 @@ def _set_record(out, budget, output="head"): def buy(ds, gwf, pname="buy", **kwargs): - """create buoyancy package from model dataset. + """Create buoyancy package from model dataset. + Parameters ---------- ds : xarray.Dataset @@ -946,7 +939,7 @@ def oc( pname="oc", **kwargs, ): - """create output control package from model dataset. + """Create output control package from model dataset. Parameters ---------- @@ -1026,7 +1019,6 @@ def ds_to_gwf(ds, complexity="SIMPLE", icelltype=0, under_relaxation=False): flopy.mf6.ModflowGwf MODFLOW6 GroundwaterFlow model object. """ - # create simulation mf_sim = sim(ds) diff --git a/nlmod/gwf/horizontal_flow_barrier.py b/nlmod/gwf/horizontal_flow_barrier.py index fe97f823..ffd0f053 100644 --- a/nlmod/gwf/horizontal_flow_barrier.py +++ b/nlmod/gwf/horizontal_flow_barrier.py @@ -9,10 +9,10 @@ def get_hfb_spd(gwf, linestrings, hydchr=1 / 100, depth=None, elevation=None): - """Generate a stress period data for horizontal flow barrier between two - cell nodes, with several limitations. The stress period data can be used - directly in the HFB package of flopy. The hfb is placed at the cell - interface; it follows the sides of the cells. + """Generate a stress period data for horizontal flow barrier between two cell nodes, + with several limitations. The stress period data can be used directly in the HFB + package of flopy. The hfb is placed at the cell interface; it follows the sides of + the cells. The estimation of the cross-sectional area at the interface is pretty crude, as the thickness at the cell interface is just the average of the thicknesses of the two @@ -101,8 +101,8 @@ def get_hfb_spd(gwf, linestrings, hydchr=1 / 100, depth=None, elevation=None): def line2hfb(gdf, gwf, prevent_rings=True, plot=False): - """Obtain the cells with a horizontal flow barrier between them from a - geodataframe with line elements. + """Obtain the cells with a horizontal flow barrier between them from a geodataframe + with line elements. Parameters ---------- @@ -278,7 +278,7 @@ def polygon_to_hfb( def plot_hfb(cellids, gwf, ax=None, color="red", **kwargs): - """plots a horizontal flow barrier. + """Plots a horizontal flow barrier. Parameters ---------- diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py index f0d0f356..44523534 100644 --- a/nlmod/gwf/output.py +++ b/nlmod/gwf/output.py @@ -11,10 +11,10 @@ from ..dims.resample import get_affine_world_to_mod from ..mfoutput.mfoutput import ( _get_budget_da, - _get_heads_da, - _get_time_index, _get_flopy_data_object, _get_grb_file, + _get_heads_da, + _get_time_index, ) logger = logging.getLogger(__name__) @@ -56,7 +56,6 @@ def get_heads_da( ): """Read binary heads file. - Parameters ---------- ds : xarray.Dataset @@ -202,9 +201,8 @@ def get_budget_da( def get_gwl_from_wet_cells(head, layer="layer", botm=None): - """Get the groundwater level from a multi-dimensional head array where dry - cells are NaN. This methods finds the most upper non-nan-value of each cell - or timestep. + """Get the groundwater level from a multi-dimensional head array where dry cells are + NaN. This methods finds the most upper non-nan-value of each cell or timestep. Parameters ---------- @@ -249,8 +247,7 @@ def get_gwl_from_wet_cells(head, layer="layer", botm=None): def get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None): - """ - Get the flow residuals of a MODFLOW 6 simulation. + """Get the flow residuals of a MODFLOW 6 simulation. Parameters ---------- @@ -272,7 +269,6 @@ def get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None): ------- da : xr.DataArray The flow residual in each cell, in m3/d. - """ if grb_file is None: grb_file = _get_grb_file(ds) @@ -289,7 +285,7 @@ def get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None): for iflowja in flowja: # residuals.append(flopy.mf6.utils.get_residuals(iflowja, grb_file)) # use our own faster method instead of a for loop: - residual = np.full(grb.shape, np.NaN) + residual = np.full(grb.shape, np.nan) residual.ravel()[mask_active] = iflowja.flatten()[flowja_index] residuals.append(residual) dims = ("time",) + dims @@ -297,7 +293,7 @@ def get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None): else: # residuals = flopy.mf6.utils.get_residuals(flowja[0], grb_file) # use our own faster method instead of a for loop: - residuals = np.full(grb.shape, np.NaN) + residuals = np.full(grb.shape, np.nan) residuals.ravel()[mask_active] = flowja[0].flatten()[flowja_index] da = xr.DataArray(residuals, dims=dims, coords=coords) return da @@ -306,8 +302,7 @@ def get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None): def get_flow_lower_face( ds, gwf=None, fname=None, grb_file=None, kstpkper=None, lays=None ): - """ - Get the flow over the lower face of all model cells + """Get the flow over the lower face of all model cells. The flow Lower Face (flf) used to be written to the budget file in previous versions of MODFLOW. In MODFLOW 6 we determine these flows from the flow-ja-face-records. @@ -335,7 +330,6 @@ def get_flow_lower_face( ------- da : xr.DataArray The flow over the lower face of each cell, in m3/d. - """ if grb_file is None: grb_file = _get_grb_file(ds) @@ -375,7 +369,7 @@ def get_flow_lower_face( flfs = [] for iflowja in flowja: if ds.gridtype == "vertex": - flf = np.full(shape, np.NaN) + flf = np.full(shape, np.nan) mask = flf_index >= 0 flf[mask] = iflowja[0, 0, flf_index[mask]] else: @@ -385,7 +379,7 @@ def get_flow_lower_face( coords = dict(coords) | {"time": _get_time_index(cbf, ds)} else: if ds.gridtype == "vertex": - flfs = np.full(shape, np.NaN) + flfs = np.full(shape, np.nan) mask = flf_index >= 0 flfs[mask] = flowja[0][0, 0, flf_index[mask]] else: diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py index fe06889c..7f54708b 100644 --- a/nlmod/gwf/recharge.py +++ b/nlmod/gwf/recharge.py @@ -20,8 +20,7 @@ def ds_to_rch( gwf, ds, mask=None, pname="rch", recharge="recharge", auxiliary=None, **kwargs ): - """Convert the recharge data in the model dataset to a rch package with - time series. + """Convert the recharge data in the model dataset to a rch package with time series. Parameters ---------- @@ -110,8 +109,8 @@ def ds_to_evt( auxiliary=None, **kwargs, ): - """Convert the evaporation data in the model dataset to a evt package with - time series. + """Convert the evaporation data in the model dataset to a evt package with time + series. Parameters ---------- @@ -141,7 +140,6 @@ def ds_to_evt( Raises ------ - DESCRIPTION. ValueError DESCRIPTION. @@ -542,8 +540,7 @@ def ds_to_uzf( def _get_unique_series(ds, var, pname): - """Get the location and values of unique time series from a variable var in - ds. + """Get the location and values of unique time series from a variable var in ds. Parameters ---------- diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py index c61f646e..7c14d8d5 100644 --- a/nlmod/gwf/surface_water.py +++ b/nlmod/gwf/surface_water.py @@ -10,11 +10,11 @@ from shapely.strtree import STRtree from tqdm import tqdm +from ..cache import cache_pickle from ..dims.grid import gdf_to_grid from ..dims.layers import get_idomain -from ..dims.resample import get_extent_polygon, extent_to_polygon, get_delr, get_delc +from ..dims.resample import extent_to_polygon, get_delc, get_delr, get_extent_polygon from ..read import bgt, waterboard -from ..cache import cache_pickle logger = logging.getLogger(__name__) @@ -42,7 +42,6 @@ def aggregate(gdf, method, ds=None): celldata : pd.DataFrame DataFrame with aggregated surface water parameters per grid cell """ - required_cols = {"stage", "c0", "botm"} missing_cols = required_cols.difference(gdf.columns) if len(missing_cols) > 0: @@ -300,9 +299,10 @@ def estimate_polygon_length(gdf): def distribute_cond_over_lays( cond, cellid, rivbot, laytop, laybot, idomain=None, kh=None, stage=None ): - """Distribute the conductance in a cell over the layers in that cell, based - on the the river-bottom and the layer bottoms, and optionally based on the - stage and the hydraulic conductivity.""" + """Distribute the conductance in a cell over the layers in that cell, based on the + the river-bottom and the layer bottoms, and optionally based on the stage and the + hydraulic conductivity. + """ if isinstance(rivbot, (np.ndarray, xr.DataArray)): rivbot = float(rivbot[cellid]) if len(laybot.shape) == 3: @@ -393,7 +393,6 @@ def build_spd( - DRN: [(cellid), elev, cond] - GHB: [(cellid), elev, cond] """ - spd = [] top = ds.top.data @@ -517,7 +516,9 @@ def add_info_to_gdf( geom_type="Polygon", add_index_from_column=None, ): - """Add information from 'gdf_from' to 'gdf_to', based on the spatial intersection.""" + """Add information from 'gdf_from' to 'gdf_to', based on the spatial + intersection. + """ gdf_to = gdf_to.copy() if columns is None: columns = gdf_from.columns[~gdf_from.columns.isin(gdf_to.columns)] @@ -702,7 +703,7 @@ def download_watercourses( def _get_waterboard_selection(gdf=None, extent=None, config=None): - """Internal method to select waterboards to get data from""" + """Internal method to select waterboards to get data from.""" if config is None: config = waterboard.get_configuration() if gdf is None and extent is None: @@ -759,7 +760,7 @@ def add_stages_from_waterboards( la = download_level_areas(gdf, extent=extent, config=config) if columns is None: columns = ["summer_stage", "winter_stage"] - gdf[columns] = np.NaN + gdf[columns] = np.nan for wb in la.keys(): if len(la[wb]) == 0: continue @@ -813,7 +814,7 @@ def add_bottom_height_from_waterboards( wc = download_watercourses(gdf, extent=extent, config=config) if columns is None: columns = ["bottom_height"] - gdf[columns] = np.NaN + gdf[columns] = np.nan for wb in wc.keys(): if len(wc[wb]) == 0: continue @@ -883,8 +884,8 @@ def get_gdf(ds=None, extent=None, fname_ahn=None, ahn=None, buffer=0.0): def add_min_ahn_to_gdf(gdf, ahn, buffer=0.0, column="ahn_min"): - """Add a column names with the minimum surface level height near surface - water features. + """Add a column names with the minimum surface level height near surface water + features. Parameters ---------- @@ -905,7 +906,6 @@ def add_min_ahn_to_gdf(gdf, ahn, buffer=0.0, column="ahn_min"): A GeoDataFrame with surface water features, with an added column containing the minimum surface level height near the features. """ - from geocube.api.core import make_geocube from geocube.rasterize import rasterize_image @@ -995,7 +995,7 @@ def gdf_to_seasonal_pkg( # make sure we have a bottom height if "rbot" not in gdf: - gdf["rbot"] = np.NaN + gdf["rbot"] = np.nan mask = gdf["rbot"].isna() if mask.any(): logger.info( diff --git a/nlmod/gwf/wells.py b/nlmod/gwf/wells.py index 3c8f0ced..0e22fe10 100644 --- a/nlmod/gwf/wells.py +++ b/nlmod/gwf/wells.py @@ -24,8 +24,7 @@ def wel_from_df( auxmultname="multiplier", **kwargs, ): - """ - Add a Well (WEL) package based on input from a (Geo)DataFrame. + """Add a Well (WEL) package based on input from a (Geo)DataFrame. Parameters ---------- @@ -70,7 +69,6 @@ def wel_from_df( ------- wel : flopy.mf6.ModflowGwfwel wel package. - """ if aux is None: aux = [] @@ -140,8 +138,7 @@ def maw_from_df( ds=None, **kwargs, ): - """ - Add a Multi-AquiferWell (MAW) package based on input from a (Geo)DataFrame. + """Add a Multi-AquiferWell (MAW) package based on input from a (Geo)DataFrame. Parameters ---------- @@ -188,7 +185,6 @@ def maw_from_df( ------- wel : flopy.mf6.ModflowGwfmaw maw package. - """ if aux is None: aux = [] @@ -267,8 +263,7 @@ def maw_from_df( def _add_cellid(df, ds=None, gwf=None, x="x", y="y"): - """ - Intersect a DataFrame of point Data with the model grid, and add cellid-column. + """Intersect a DataFrame of point Data with the model grid, and add cellid-column. Parameters ---------- @@ -290,7 +285,6 @@ def _add_cellid(df, ds=None, gwf=None, x="x", y="y"): df : gpd.GeoDataFrame A GeoDataFrame with a column named cellid that contains the icell2d-number (vertex-grid) or (row, column) (structured grid). - """ if not isinstance(df, gpd.GeoDataFrame): df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df[x], df[y])) @@ -300,8 +294,7 @@ def _add_cellid(df, ds=None, gwf=None, x="x", y="y"): def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): - """ - Get factors (pandas.DataFrame) for each layer that well screens intersects with. + """Get factors (pandas.DataFrame) for each layer that well screens intersects with. Parameters ---------- @@ -322,7 +315,6 @@ def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): multipliers : pd.DataFrame A DataFrame containg the multiplication factors, with the layers as the index and the name of the well screens (the index of df) as columns. - """ # get required data either from gwf or ds if ds is not None: @@ -351,8 +343,7 @@ def _get_layer_multiplier_for_wells(df, top, botm, ds=None, gwf=None): def _get_layer_multiplier_for_well(cid, well_top, well_bot, ml_top, ml_bot, ml_kh): - """ - Get a factor (numpy array) for each layer that a well screen intersects with. + """Get a factor (numpy array) for each layer that a well screen intersects with. Parameters ---------- @@ -373,7 +364,6 @@ def _get_layer_multiplier_for_well(cid, well_top, well_bot, ml_top, ml_bot, ml_k ------- multiplier : numpy array An array with a factor (between 0 and 1) for each of the model layers. - """ # keep the tops and botms of the cell where the well is in ml_top_cid = ml_top[cid].copy() diff --git a/nlmod/gwt/__init__.py b/nlmod/gwt/__init__.py index 78096246..869b583f 100644 --- a/nlmod/gwt/__init__.py +++ b/nlmod/gwt/__init__.py @@ -1,3 +1,4 @@ +# ruff: noqa: F401 F403 from . import output, prepare from .gwt import * from .output import * diff --git a/nlmod/gwt/gwt.py b/nlmod/gwt/gwt.py index 7311715a..5ae02377 100644 --- a/nlmod/gwt/gwt.py +++ b/nlmod/gwt/gwt.py @@ -11,7 +11,7 @@ def gwt(ds, sim, modelname=None, **kwargs): - """create groundwater transport model from the model dataset. + """Create groundwater transport model from the model dataset. Parameters ---------- @@ -29,7 +29,6 @@ def gwt(ds, sim, modelname=None, **kwargs): gwt : flopy ModflowGwt groundwater transport object. """ - # start creating model logger.info("creating mf6 GWT") @@ -46,7 +45,7 @@ def gwt(ds, sim, modelname=None, **kwargs): def dis(ds, gwt, length_units="METERS", pname="dis", **kwargs): - """create discretisation package from the model dataset. + """Create discretisation package from the model dataset. Parameters ---------- @@ -68,7 +67,7 @@ def dis(ds, gwt, length_units="METERS", pname="dis", **kwargs): def disv(ds, gwt, length_units="METERS", pname="disv", **kwargs): - """create discretisation vertices package from the model dataset. + """Create discretisation vertices package from the model dataset. Parameters ---------- @@ -90,7 +89,7 @@ def disv(ds, gwt, length_units="METERS", pname="disv", **kwargs): def adv(ds, gwt, scheme=None, **kwargs): - """create advection package for groundwater transport model. + """Create advection package for groundwater transport model. Parameters ---------- @@ -114,7 +113,7 @@ def adv(ds, gwt, scheme=None, **kwargs): def dsp(ds, gwt, **kwargs): - """create dispersion package for groundwater transport model. + """Create dispersion package for groundwater transport model. Parameters ---------- @@ -139,7 +138,7 @@ def dsp(ds, gwt, **kwargs): def ssm(ds, gwt, sources=None, **kwargs): - """create source-sink mixing package for groundwater transport model. + """Create source-sink mixing package for groundwater transport model. Parameters ---------- @@ -177,7 +176,7 @@ def ssm(ds, gwt, sources=None, **kwargs): def mst(ds, gwt, porosity=None, **kwargs): - """create mass storage transfer package for groundwater transport model. + """Create mass storage transfer package for groundwater transport model. Parameters ---------- @@ -213,7 +212,7 @@ def mst(ds, gwt, porosity=None, **kwargs): def cnc(ds, gwt, da_mask, da_conc, pname="cnc", **kwargs): - """create constant concentration package for groundwater transport model. + """Create constant concentration package for groundwater transport model. Parameters ---------- @@ -251,7 +250,7 @@ def oc( pname="oc", **kwargs, ): - """create output control package for groundwater transport model. + """Create output control package for groundwater transport model. Parameters ---------- @@ -290,7 +289,7 @@ def oc( def ic(ds, gwt, strt, pname="ic", **kwargs): - """create initial condictions package for groundwater transport model. + """Create initial condictions package for groundwater transport model. Parameters ---------- @@ -319,7 +318,7 @@ def ic(ds, gwt, strt, pname="ic", **kwargs): def gwfgwt(ds, sim, exgtype="GWF6-GWT6", **kwargs): - """create GWF-GWT exchange package for mf6 simulation. + """Create GWF-GWT exchange package for mf6 simulation. Parameters ---------- diff --git a/nlmod/gwt/output.py b/nlmod/gwt/output.py index 955c2061..5dfd76b5 100644 --- a/nlmod/gwt/output.py +++ b/nlmod/gwt/output.py @@ -4,7 +4,7 @@ import xarray as xr from ..dims.layers import calculate_thickness -from ..mfoutput.mfoutput import _get_heads_da, _get_time_index, _get_flopy_data_object +from ..mfoutput.mfoutput import _get_flopy_data_object, _get_heads_da, _get_time_index logger = logging.getLogger(__name__) @@ -91,9 +91,9 @@ def get_concentration_da( def get_concentration_at_gw_surface(conc, layer="layer"): - """Get the concentration level from a multi-dimensional concentration array - where dry or inactive cells are NaN. This methods finds the most upper non- - nan-value of each cell or timestep. + """Get the concentration level from a multi-dimensional concentration array where + dry or inactive cells are NaN. This methods finds the most upper non- nan-value of + each cell or timestep. Parameters ---------- @@ -137,8 +137,8 @@ def get_concentration_at_gw_surface(conc, layer="layer"): def freshwater_head(ds, hp, conc, denseref=None, drhodc=None): - """Calculate equivalent freshwater head from point water heads. - Heads file produced by mf6 contains point water heads. + """Calculate equivalent freshwater head from point water heads. Heads file produced + by mf6 contains point water heads. Parameters ---------- @@ -180,8 +180,8 @@ def freshwater_head(ds, hp, conc, denseref=None, drhodc=None): def pointwater_head(ds, hf, conc, denseref=None, drhodc=None): - """Calculate point water head from freshwater heads. - Heads file produced by mf6 contains point water heads. + """Calculate point water head from freshwater heads. Heads file produced by mf6 + contains point water heads. Parameters ---------- diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py index f7b21ed5..a2771f86 100644 --- a/nlmod/mfoutput/mfoutput.py +++ b/nlmod/mfoutput/mfoutput.py @@ -1,11 +1,10 @@ -import os import logging +import os import warnings import dask -import xarray as xr - import flopy +import xarray as xr from ..dims.grid import get_dims_coords_from_modelgrid, modelgrid_from_ds from ..dims.resample import get_affine_mod_to_world @@ -243,7 +242,7 @@ def _get_flopy_data_object( var, ds=None, gwml=None, fname=None, grb_file=None, **kwargs ): """Get modflow HeadFile or CellBudgetFile object, containg heads, budgets or - concentrations + concentrations. Provide one of ds, gwf or fname. diff --git a/nlmod/modpath/__init__.py b/nlmod/modpath/__init__.py index 396e6362..0154abca 100644 --- a/nlmod/modpath/__init__.py +++ b/nlmod/modpath/__init__.py @@ -1 +1,2 @@ +# ruff: noqa: F403 from .modpath import * diff --git a/nlmod/modpath/modpath.py b/nlmod/modpath/modpath.py index 336db050..37c74409 100644 --- a/nlmod/modpath/modpath.py +++ b/nlmod/modpath/modpath.py @@ -16,8 +16,8 @@ def write_and_run(mpf, remove_prev_output=True, script_path=None, silent=False): - """write modpath files and run the model. Extra options include removing - previous output and copying the modelscript to the model workspace. + """Write modpath files and run the model. Extra options include removing previous + output and copying the modelscript to the model workspace. Parameters ---------- @@ -55,9 +55,8 @@ def write_and_run(mpf, remove_prev_output=True, script_path=None, silent=False): def xy_to_nodes(xy_list, mpf, ds, layer=0): - """convert a list of points, defined by x and y coordinates, to a list of - nodes. A node is a unique cell in a model. The icell2d is a unique cell in - a layer. + """Convert a list of points, defined by x and y coordinates, to a list of nodes. A + node is a unique cell in a model. The icell2d is a unique cell in a layer. Parameters ---------- @@ -145,7 +144,7 @@ def package_to_nodes(gwf, package_name, mpf): def layer_to_nodes(mpf, modellayer): - """get the nodes of all cells in one ore more model layer(s). + """Get the nodes of all cells in one ore more model layer(s). Parameters ---------- @@ -264,16 +263,14 @@ def bas(mpf, porosity=0.3, **kwargs): mpfbas : flopy.modpath.mp7bas.Modpath7Bas modpath bas package. """ - mpfbas = flopy.modpath.Modpath7Bas(mpf, porosity=porosity, **kwargs) return mpfbas def remove_output(mpf): - """Remove the output of a previous modpath run. Commonly used before - starting a new modpath run to avoid loading the wrong data when a modpath - run has failed. + """Remove the output of a previous modpath run. Commonly used before starting a new + modpath run to avoid loading the wrong data when a modpath run has failed. Parameters ---------- @@ -456,6 +453,7 @@ def sim( simulationtype="combined", weaksinkoption="pass_through", weaksourceoption="pass_through", + **kwargs, ): """Create a modpath backward simulation from a particle group. @@ -505,6 +503,7 @@ def sim( stoptimeoption=stoptimeoption, stoptime=stoptime, particlegroups=particlegroups, + **kwargs, ) return mpsim diff --git a/nlmod/plot/__init__.py b/nlmod/plot/__init__.py index ffeb6bd0..e77e6c9b 100644 --- a/nlmod/plot/__init__.py +++ b/nlmod/plot/__init__.py @@ -1,3 +1,4 @@ +# ruff: noqa: F401 from . import flopy from .dcs import DatasetCrossSection from .plot import ( @@ -7,6 +8,7 @@ geotop_lithok_in_cross_section, geotop_lithok_on_map, map_array, + # modelextent, modelgrid, surface_water, ) diff --git a/nlmod/plot/dcs.py b/nlmod/plot/dcs.py index f27162b0..cc44ed24 100644 --- a/nlmod/plot/dcs.py +++ b/nlmod/plot/dcs.py @@ -1,14 +1,13 @@ -import flopy import logging -import matplotlib +from functools import partial +import flopy import geopandas as gpd +import matplotlib import matplotlib.pyplot as plt import numpy as np import pandas as pd import xarray as xr - -from functools import partial from matplotlib.animation import FFMpegWriter, FuncAnimation from matplotlib.collections import LineCollection, PatchCollection from matplotlib.patches import Rectangle @@ -19,7 +18,6 @@ from ..dims.resample import get_affine_world_to_mod from .plotutil import get_map - logger = logging.getLogger(__name__) @@ -342,7 +340,7 @@ def plot_map_cs( label="cross section", **kwargs, ): - """Creates a different figure with the map of the cross section + """Creates a different figure with the map of the cross section. Parameters ---------- @@ -365,7 +363,6 @@ def plot_map_cs( matplotlib Axes axes """ - if ax is None: _, ax = get_map( self.ds.extent, background=background, figsize=figsize, **kwargs @@ -376,7 +373,7 @@ def plot_map_cs( return ax def get_patches_array(self, z): - """similar to plot_array function, only computes the array to update an existing plot_array. + """Similar to plot_array function, only computes the array to update an existing plot_array. Parameters ---------- @@ -525,7 +522,7 @@ def animate( cbar_label=None, fname=None, ): - """animate a cross section + """Animate a cross section. Parameters ---------- @@ -553,7 +550,6 @@ def animate( matplotlib.animation.FuncAnimation animation object """ - f = self.ax.get_figure() # plot first timeframe diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py index 7ff8d0fe..66fb315e 100644 --- a/nlmod/plot/plot.py +++ b/nlmod/plot/plot.py @@ -2,10 +2,10 @@ from functools import partial import flopy as fp -import xarray as xr import matplotlib.pyplot as plt import numpy as np import pandas as pd +import xarray as xr from matplotlib.animation import FFMpegWriter, FuncAnimation from matplotlib.collections import PatchCollection from matplotlib.colors import ListedColormap, Normalize @@ -59,7 +59,7 @@ def facet_plot( xlim=None, ylim=None, ): - """make a 2d plot of every modellayer, store them in a grid. + """Make a 2d plot of every modellayer, store them in a grid. Parameters ---------- @@ -92,7 +92,6 @@ def facet_plot( axes : TYPE DESCRIPTION. """ - warnings.warn( "this function is out of date and will probably be removed in a future version", DeprecationWarning, @@ -216,7 +215,14 @@ def data_array(da, ds=None, ax=None, rotated=False, edgecolor=None, **kwargs): def geotop_lithok_in_cross_section( - line, gt=None, ax=None, legend=True, legend_loc=None, lithok_props=None, alpha=None, **kwargs + line, + gt=None, + ax=None, + legend=True, + legend_loc=None, + lithok_props=None, + alpha=None, + **kwargs, ): """PLot the lithoclass-data of GeoTOP in a cross-section. @@ -269,7 +275,7 @@ def geotop_lithok_in_cross_section( cs = DatasetCrossSection(gt, line, layer="z", ax=ax, **kwargs) array, cmap, norm = _get_geotop_cmap_and_norm(gt["lithok"], lithok_props) cs.plot_array(array, norm=norm, cmap=cmap, alpha=alpha) - + if legend: # make a legend with dummy handles _add_geotop_lithok_legend(lithok_props, ax, lithok=gt["lithok"], loc=legend_loc) @@ -306,7 +312,6 @@ def geotop_lithok_on_map( Returns ------- qm : matplotlib.collections.QuadMesh - """ if ax is None: ax = plt.gca() @@ -330,7 +335,7 @@ def geotop_lithok_on_map( def _add_geotop_lithok_legend(lithok_props, ax, lithok=None, **kwargs): - """Add a legend with lithok-data""" + """Add a legend with lithok-data.""" handles = [] if lithok is None: lithoks = lithok_props.index @@ -345,10 +350,10 @@ def _add_geotop_lithok_legend(lithok_props, ax, lithok=None, **kwargs): def _get_geotop_cmap_and_norm(lithok, lithok_props): - """Get an array of lithok-values, with a corresponding colormap and norm""" + """Get an array of lithok-values, with a corresponding colormap and norm.""" lithok_un = np.unique(lithok) lithok_un = lithok_un[~np.isnan(lithok_un)] - array = np.full(lithok.shape, np.NaN) + array = np.full(lithok.shape, np.nan) colors = [] for i, ilithok in enumerate(lithok_un): ilithok = int(ilithok) diff --git a/nlmod/plot/plotutil.py b/nlmod/plot/plotutil.py index 00c134e3..67c6c09d 100644 --- a/nlmod/plot/plotutil.py +++ b/nlmod/plot/plotutil.py @@ -186,8 +186,7 @@ def rotate_yticklabels(ax): def rd_ticks(ax, base=1000.0, fmt_base=1000.0, fmt="{:.0f}"): - """Add ticks every 1000 (base) m, and divide ticklabels by 1000 - (fmt_base)""" + """Add ticks every 1000 (base) m, and divide ticklabels by 1000 (fmt_base).""" def fmt_rd_ticks(x, _): return fmt.format(x / fmt_base) diff --git a/nlmod/read/__init__.py b/nlmod/read/__init__.py index f6c7f2ce..5e605993 100644 --- a/nlmod/read/__init__.py +++ b/nlmod/read/__init__.py @@ -1,3 +1,4 @@ +# ruff: noqa: F401 from . import ( administrative, ahn, @@ -9,11 +10,11 @@ knmi, knmi_data_platform, meteobase, + nhi, regis, rws, waterboard, webservices, - nhi, ) from .geotop import get_geotop from .regis import get_regis diff --git a/nlmod/read/ahn.py b/nlmod/read/ahn.py index 1bd2fda8..859f2ceb 100644 --- a/nlmod/read/ahn.py +++ b/nlmod/read/ahn.py @@ -1,10 +1,10 @@ import datetime as dt import logging +import geopandas as gpd import matplotlib.pyplot as plt import numpy as np import pandas as pd -import geopandas as gpd import rasterio import rioxarray import xarray as xr @@ -90,8 +90,7 @@ def get_ahn_at_point( res=0.5, **kwargs, ): - """ - Get the height of the surface level at a certain point, defined by x and y. + """Get the height of the surface level at a certain point, defined by x and y. Parameters ---------- @@ -117,7 +116,6 @@ def get_ahn_at_point( ------- float The surface level value at the requested point. - """ extent = [x - buffer, x + buffer, y - buffer, y + buffer] ahn = get_latest_ahn_from_wcs(extent, identifier=identifier, res=res, **kwargs) @@ -133,8 +131,7 @@ def get_ahn_at_point( def get_ahn_along_line(line, ahn=None, dx=None, num=None, method="linear", plot=False): - """ - Get the height of the surface level along a line. + """Get the height of the surface level along a line. Parameters ---------- @@ -161,7 +158,6 @@ def get_ahn_along_line(line, ahn=None, dx=None, num=None, method="linear", plot= ------- z : xr.DataArray A DataArray with dimension s, containing surface level values along the line. - """ if ahn is None: bbox = line.bounds @@ -236,7 +232,6 @@ def get_latest_ahn_from_wcs( xr.DataArray or MemoryFile DataArray (if as_data_array is True) or Rasterio MemoryFile of the AHN """ - url = "https://service.pdok.nl/rws/ahn/wcs/v1_0?SERVICE=WCS&request=GetCapabilities" if isinstance(extent, xr.DataArray): @@ -275,10 +270,9 @@ def get_latest_ahn_from_wcs( def get_ahn2_tiles(extent=None): """Get the tiles (kaartbladen) of AHN3 as a GeoDataFrame. - The links in the tiles are cuurently incorrect. Thereore - get_ahn3_tiles is used in get_ahn2 and get_ahn1, as the tiles from - get_ahn3_tiles also contain information about the tiles of ahn1 and - ahn2 + The links in the tiles are cuurently incorrect. Thereore get_ahn3_tiles is used in + get_ahn2 and get_ahn1, as the tiles from get_ahn3_tiles also contain information + about the tiles of ahn1 and ahn2 """ url = "https://services.arcgis.com/nSZVuSZjHpEZZbRo/arcgis/rest/services/Kaartbladen_AHN2/FeatureServer" layer = 0 @@ -299,8 +293,7 @@ def get_ahn3_tiles(extent=None): def get_ahn4_tiles(extent=None): - """Get the tiles (kaartbladen) of AHN4 as a GeoDataFrame with download - links.""" + """Get the tiles (kaartbladen) of AHN4 as a GeoDataFrame with download links.""" url = "https://services.arcgis.com/nSZVuSZjHpEZZbRo/arcgis/rest/services/Kaartbladen_AHN4/FeatureServer" layer = 0 gdf = arcrest(url, layer, extent) diff --git a/nlmod/read/bro.py b/nlmod/read/bro.py index fe19dbbf..352c67aa 100644 --- a/nlmod/read/bro.py +++ b/nlmod/read/bro.py @@ -11,7 +11,7 @@ def add_modelled_head(oc, ml=None, ds=None, method="linear"): - """add modelled heads as seperate observations to the ObsCollection. + """Add modelled heads as seperate observations to the ObsCollection. Parameters ---------- @@ -30,7 +30,6 @@ def add_modelled_head(oc, ml=None, ds=None, method="linear"): ObsCollection combination of observed and modelled groundwater heads. """ - oc["modellayer"] = oc.gwobs.get_modellayers(gwf=ml) if ds is not None and "heads" in ds: heads = ds["heads"] @@ -100,7 +99,7 @@ def get_bro( max_screen_top=None, min_screen_bot=None, ): - """get bro groundwater measurements within an extent. + """Get bro groundwater measurements within an extent. Parameters ---------- @@ -168,10 +167,10 @@ def get_bro( @cache.cache_pickle def get_bro_metadata(extent, max_dx=20000, max_dy=20000): - """wrapper around hpd.read_bro that deals with large extents and only - returns metadata (location, tube top/bot, ground level, ..) of the wells - and no actual measurements. This is useful when the extent is too big - to obtain all measurements at once. + """Wrapper around hpd.read_bro that deals with large extents and only returns + metadata (location, tube top/bot, ground level, ..) of the wells and no actual + measurements. This is useful when the extent is too big to obtain all measurements + at once. Parameters ---------- @@ -188,7 +187,6 @@ def get_bro_metadata(extent, max_dx=20000, max_dy=20000): ------- ObsCollection """ - # check if extent is within limits dx = extent[1] - extent[0] dy = extent[3] - extent[2] diff --git a/nlmod/read/geotop.py b/nlmod/read/geotop.py index 536172c1..9457e53a 100644 --- a/nlmod/read/geotop.py +++ b/nlmod/read/geotop.py @@ -103,7 +103,6 @@ def to_model_layers( ds: xr.DataSet dataset with top and botm (and optionally kh and kv) per geotop layer """ - if strat_props is None: strat_props = get_strat_props() @@ -135,8 +134,8 @@ def to_model_layers( geulen = [] for layer, unit in enumerate(units): mask = strat == unit - top[layer] = np.nanmax(np.where(mask, z, np.NaN), 0) + 0.25 - bot[layer] = np.nanmin(np.where(mask, z, np.NaN), 0) - 0.25 + top[layer] = np.nanmax(np.where(mask, z, np.nan), 0) + 0.25 + bot[layer] = np.nanmin(np.where(mask, z, np.nan), 0) - 0.25 if int(unit) in strat_props.index: layers.append(strat_props.at[unit, "code"]) else: @@ -236,8 +235,8 @@ def to_model_layers( @cache.cache_netcdf() def get_geotop(extent, url=GEOTOP_URL, probabilities=False): """Get a slice of the geotop netcdf url within the extent, set the x and y - coordinates to match the cell centers and keep only the strat and lithok - data variables. + coordinates to match the cell centers and keep only the strat and lithok data + variables. Parameters ---------- @@ -319,12 +318,12 @@ def add_top_and_botm(ds): """ bottom = np.expand_dims(ds.z.values - 0.25, axis=(1, 2)) bottom = np.repeat(np.repeat(bottom, len(ds.y), 1), len(ds.x), 2) - bottom[np.isnan(ds.strat.values)] = np.NaN + bottom[np.isnan(ds.strat.values)] = np.nan ds["botm"] = ("z", "y", "x"), bottom top = np.expand_dims(ds.z.values + 0.25, axis=(1, 2)) top = np.repeat(np.repeat(top, len(ds.y), 1), len(ds.x), 2) - top[np.isnan(ds.strat.values)] = np.NaN + top[np.isnan(ds.strat.values)] = np.nan ds["top"] = ("z", "y", "x"), top return ds @@ -380,7 +379,6 @@ def add_kh_and_kv( Raises ------ - DESCRIPTION. Returns @@ -411,8 +409,8 @@ def add_kh_and_kv( if stochastic is None: # calculate kh and kv from most likely lithoclass lithok = gt["lithok"].values - kh_ar = np.full(lithok.shape, np.NaN) - kv_ar = np.full(lithok.shape, np.NaN) + kh_ar = np.full(lithok.shape, np.nan) + kv_ar = np.full(lithok.shape, np.nan) if "strat" in df: combs = np.column_stack((strat.ravel(), lithok.ravel())) # drop nans @@ -445,7 +443,7 @@ def add_kh_and_kv( probality = gt[f"kans_{ilithok}"].values if "strat" in df: khi, kvi = _handle_nans_in_stochastic_approach( - np.NaN, np.NaN, kh_method, kv_method + np.nan, np.nan, kh_method, kv_method ) khi = np.full(strat.shape, khi) kvi = np.full(strat.shape, kvi) @@ -508,7 +506,7 @@ def _get_kh_kv_from_df(df, ilithok, istrat=None, anisotropy=1.0, mask=None): else: msg = f"{msg}. Setting values of {mask.sum()} voxels to NaN." logger.warning(msg) - return np.NaN, np.NaN + return np.nan, np.nan kh = df.loc[mask_df, "kh"].mean() if "kv" in df: @@ -540,8 +538,8 @@ def _handle_nans_in_stochastic_approach(kh, kv, kh_method, kv_method): def aggregate_to_ds( gt, ds, kh="kh", kv="kv", kd="kD", c="c", kh_gt="kh", kv_gt="kv", add_kd_and_c=False ): - """Aggregate voxels from GeoTOP to layers in a model DataSet with top and - botm, to calculate kh and kv. + """Aggregate voxels from GeoTOP to layers in a model DataSet with top and botm, to + calculate kh and kv. Parameters ---------- diff --git a/nlmod/read/jarkus.py b/nlmod/read/jarkus.py index 5d962973..4ffe921c 100644 --- a/nlmod/read/jarkus.py +++ b/nlmod/read/jarkus.py @@ -1,12 +1,13 @@ -"""module with functions to deal with the northsea by: +"""Module with functions to deal with the northsea. - - identifying model cells with the north sea + - identify model cells with the north sea - add bathymetry of the northsea to the layer model - - extrpolate the layer model below the northsea bed. + - extrapolate the layer model below the northsea bed. Note: if you like jazz please check this out: https://www.northseajazz.com """ + import datetime as dt import logging import os @@ -25,7 +26,7 @@ @cache.cache_netcdf() def get_bathymetry(ds, northsea, kind="jarkus", method="average"): - """get bathymetry of the Northsea from the jarkus dataset. + """Get bathymetry of the Northsea from the jarkus dataset. Parameters ---------- @@ -126,7 +127,6 @@ def get_dataset_jarkus(extent, kind="jarkus", return_tiles=False, time=-1): dataset containing bathymetry data """ - extent = [int(x) for x in extent] netcdf_tile_names = get_jarkus_tilenames(extent, kind) diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py index b142848c..466641bd 100644 --- a/nlmod/read/knmi.py +++ b/nlmod/read/knmi.py @@ -15,7 +15,7 @@ @cache.cache_netcdf(coords_3d=True, coords_time=True) def get_recharge(ds, method="linear", most_common_station=False): - """add multiple recharge packages to the groundwater flow model with knmi + """Add multiple recharge packages to the groundwater flow model with knmi data by following these steps: 1. check for each cell (structured or vertex) which knmi measurement stations (prec and evap) are the closest. @@ -50,7 +50,6 @@ def get_recharge(ds, method="linear", most_common_station=False): ds : xr.DataSet dataset with spatial model data including the rch raster """ - if "time" not in ds: raise ( AttributeError( @@ -159,7 +158,7 @@ def _add_ts_to_ds(timeseries, loc_sel, variable, ds): def get_locations_vertex(ds): - """get dataframe with the locations of the grid cells of a vertex grid. + """Get dataframe with the locations of the grid cells of a vertex grid. Parameters ---------- @@ -193,7 +192,7 @@ def get_locations_vertex(ds): def get_locations_structured(ds): - """get dataframe with the locations of the grid cells of a structured grid. + """Get dataframe with the locations of the grid cells of a structured grid. Parameters ---------- @@ -206,7 +205,6 @@ def get_locations_structured(ds): DataFrame with the locations of all active grid cells. includes the columns: x, y, row, col and layer """ - # store x and y mids in locations of active cells fal = get_first_active_layer(ds) rows, columns = np.where(fal != fal.attrs["nodata"]) @@ -228,7 +226,7 @@ def get_locations_structured(ds): def get_knmi_at_locations(ds, start="2010", end=None, most_common_station=False): - """get knmi data at the locations of the active grid cells in ds. + """Get knmi data at the locations of the active grid cells in ds. Parameters ---------- diff --git a/nlmod/read/knmi_data_platform.py b/nlmod/read/knmi_data_platform.py index 44c1b90c..afda1362 100644 --- a/nlmod/read/knmi_data_platform.py +++ b/nlmod/read/knmi_data_platform.py @@ -58,7 +58,7 @@ def get_list_of_files( start_after_filename: Optional[str] = None, timeout: int = 120, ) -> List[str]: - """Download list of files from KNMI data platform""" + """Download list of files from KNMI data platform.""" if api_key is None: api_key = get_anonymous_api_key() files = [] @@ -89,7 +89,7 @@ def download_file( api_key: Optional[str] = None, timeout: int = 120, ) -> None: - """Download file from KNMI data platform""" + """Download file from KNMI data platform.""" if api_key is None: api_key = get_anonymous_api_key() url = ( @@ -119,7 +119,7 @@ def download_files( api_key: Optional[str] = None, timeout: int = 120, ) -> None: - """Download multiple files from KNMI data platform""" + """Download multiple files from KNMI data platform.""" for fname in tqdm(fnames): download_file( dataset_name=dataset_name, @@ -132,7 +132,7 @@ def download_files( def read_nc(fo: Union[str, FileIO], **kwargs: dict) -> xr.Dataset: - """Read netcdf (.nc) file to xarray Dataset""" + """Read netcdf (.nc) file to xarray Dataset.""" # could help to provide argument: engine="h5netcdf" return xr.open_dataset(fo, **kwargs) @@ -161,7 +161,7 @@ def get_timestamp_from_fname(fname: str) -> Union[Timestamp, None]: def add_h5_meta(meta: Dict[str, Any], h5obj: Any, orig_ky: str = "") -> Dict[str, Any]: - """Read metadata from hdf5 (.h5) file and add to existing metadata dictionary""" + """Read metadata from hdf5 (.h5) file and add to existing metadata dictionary.""" def cleanup(val: Any) -> Any: if isinstance(val, (ndarray, list)): @@ -174,7 +174,7 @@ def cleanup(val: Any) -> Any: return val if hasattr(h5obj, "attrs"): - attrs = getattr(h5obj, "attrs") + attrs = h5obj.attrs submeta = {f"{orig_ky}/{ky}": cleanup(val) for ky, val in attrs.items()} meta.update(submeta) @@ -186,7 +186,7 @@ class MultipleDatasetsFound(Exception): def read_h5_contents(h5fo: FileIO) -> Tuple[ndarray, Dict[str, Any]]: - """Read contents from a hdf5 (.h5) file""" + """Read contents from a hdf5 (.h5) file.""" from h5py import Dataset as h5Dataset data = None @@ -206,7 +206,7 @@ def read_h5_contents(h5fo: FileIO) -> Tuple[ndarray, Dict[str, Any]]: def read_h5(fo: Union[str, FileIO]) -> xr.Dataset: - """Read hdf5 (.h5) file to xarray Dataset""" + """Read hdf5 (.h5) file to xarray Dataset.""" from h5py import File as h5File with h5File(fo) as h5fo: @@ -231,7 +231,7 @@ def read_h5(fo: Union[str, FileIO]) -> xr.Dataset: def read_grib( fo: Union[str, FileIO], filter_by_keys=None, **kwargs: dict ) -> xr.Dataset: - """Read GRIB file to xarray Dataset""" + """Read GRIB file to xarray Dataset.""" if kwargs is None: kwargs = {} @@ -248,7 +248,7 @@ def read_grib( def read_dataset_from_zip( fname: str, hour: Optional[int] = None, **kwargs: dict ) -> xr.Dataset: - """Read KNMI data platfrom .zip file to xarray Dataset""" + """Read KNMI data platfrom .zip file to xarray Dataset.""" if fname.endswith(".zip"): with ZipFile(fname) as zipfo: fnames = sorted([x for x in zipfo.namelist() if not x.endswith("/")]) @@ -276,7 +276,7 @@ def read_dataset( hour: Optional[int] = None, **kwargs: dict, ) -> xr.Dataset: - """Read xarray dataset from different file types; .nc, .h5 or grib file""" + """Read xarray dataset from different file types; .nc, .h5 or grib file.""" if hour is not None: if hour == 24: hour = 0 diff --git a/nlmod/read/meteobase.py b/nlmod/read/meteobase.py index a5df07db..5b3a9e11 100644 --- a/nlmod/read/meteobase.py +++ b/nlmod/read/meteobase.py @@ -11,8 +11,7 @@ class MeteobaseType(Enum): - """Enum class to couple folder names to observation type (from in - LEESMIJ.txt)""" + """Enum class to couple folder names to observation type (from in LEESMIJ.txt)""" NEERSLAG = "Neerslagradargegevens in Arc/Info-formaat." MAKKINK = "Verdampingsgegevens volgens Makkink." @@ -56,8 +55,7 @@ def read_leesmij(fo: FileIO) -> Dict[str, Dict[str, str]]: def get_timestamp_from_fname(fname: str) -> Timestamp: - """Get the Timestamp from a filename (with some assumptions about the - formatting)""" + """Get the Timestamp from a filename (with some assumptions about the formatting)""" datestr = re.search("([0-9]{8})", fname) # assumes YYYYMMDD if datestr is not None: match = datestr.group(0) @@ -130,7 +128,7 @@ def read_ascii(fo: FileIO) -> Union[np.ndarray, dict]: def get_xy_from_ascii_meta( - meta: Dict[str, Union[int, float]] + meta: Dict[str, Union[int, float]], ) -> Tuple[np.ndarray, np.ndarray]: """Get the xy coordinates Esri ASCII raster format header. @@ -268,7 +266,6 @@ def read_meteobase( ------- List[DataArray] """ - with ZipFile(Path(path)) as zfile: with zfile.open("LEESMIJ.TXT") as fo: meta = read_leesmij(fo) diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py index 0ddaf34b..4e721d60 100644 --- a/nlmod/read/nhi.py +++ b/nlmod/read/nhi.py @@ -1,12 +1,11 @@ +import io import logging import os -import io -import requests +import geopandas as gpd import numpy as np import pandas as pd -import geopandas as gpd - +import requests import rioxarray from ..dims.resample import structured_da_to_ds @@ -15,8 +14,7 @@ def download_file(url, pathname, filename=None, overwrite=False, timeout=120.0): - """ - Download a file from the NHI website. + """Download a file from the NHI website. Parameters ---------- @@ -37,7 +35,6 @@ def download_file(url, pathname, filename=None, overwrite=False, timeout=120.0): ------- fname : str The full path of the downloaded file. - """ if filename is None: filename = url.split("/")[-1] @@ -51,8 +48,7 @@ def download_file(url, pathname, filename=None, overwrite=False, timeout=120.0): def download_buisdrainage(pathname, overwrite=False): - """ - Download resistance and depth of buisdrainage from the NHI website + """Download resistance and depth of buisdrainage from the NHI website. Parameters ---------- @@ -67,7 +63,6 @@ def download_buisdrainage(pathname, overwrite=False): The full path of the downloaded file containing the resistance of buisdrainage. fname_d : str The full path of the downloaded file containing the depth of buisdrainage. - """ url_bas = "https://thredds.data.nhi.nu/thredds/fileServer/opendap/models/nhi3_2/25m" @@ -90,8 +85,7 @@ def add_buisdrainage( cond_method="average", depth_method="mode", ): - """ - Add data about the buisdrainage to the model Dataset. + """Add data about the buisdrainage to the model Dataset. This data consists of the conductance of buisdrainage (m2/d) and the depth of buisdrainage (m to surface level). With the default settings for `cond_method` and @@ -129,7 +123,6 @@ def add_buisdrainage( ds : xr.Dataset The model dataset with added variables with the names `cond_var` and `depth_var`. - """ if pathname is None: pathname = ds.cachedir @@ -190,8 +183,7 @@ def get_gwo_wells( timeout=120, **kwargs, ): - """ - Get metadata of extraction wells from the NHI GWO database + """Get metadata of extraction wells from the NHI GWO database. Parameters ---------- @@ -228,7 +220,6 @@ def get_gwo_wells( ------- gdf : geopandas.GeoDataFrame A GeoDataFrame containing the properties of the wells and their filters. - """ # zie https://gwo.nhi.nu/api/v1/download/ url = "https://gwo.nhi.nu/api/v1/well_filters/" @@ -284,8 +275,7 @@ def get_gwo_measurements( timeout=120, **kwargs, ): - """ - Get extraction rates and metadata of wells from the NHI GWO database + """Get extraction rates and metadata of wells from the NHI GWO database. Parameters ---------- @@ -321,7 +311,6 @@ def get_gwo_measurements( A DataFrame containing the extraction rates of the wells in the database. gdf : geopandas.GeoDataFrame A GeoDataFrame containing the properties of the wells and their filters. - """ url = "http://gwo.nhi.nu/api/v1/measurements/" properties = [] diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py index 1f27c6c7..f1ba8683 100644 --- a/nlmod/read/regis.py +++ b/nlmod/read/regis.py @@ -25,7 +25,7 @@ def get_combined_layer_models( geotop_layers="HLc", geotop_k=None, ): - """combine layer models into a single layer model. + """Combine layer models into a single layer model. Possibilities so far include: - use_regis -> full model based on regis @@ -64,7 +64,6 @@ def get_combined_layer_models( ValueError if an invalid combination of layers is used. """ - if use_regis: regis_ds = get_regis( extent, regis_botm_layer, remove_nan_layers=remove_nan_layers @@ -101,7 +100,7 @@ def get_regis( drop_layer_dim_from_top=True, probabilities=False, ): - """get a regis dataset projected on the modelgrid. + """Get a regis dataset projected on the modelgrid. Parameters ---------- @@ -132,7 +131,6 @@ def get_regis( regis_ds : xarray dataset dataset with regis data projected on the modelgrid. """ - ds = xr.open_dataset(REGIS_URL, decode_times=False) # set x and y dimensions to cell center @@ -203,8 +201,8 @@ def add_geotop_to_regis_layers( anisotropy=1.0, gt_layered=None, ): - """Combine geotop and regis in such a way that the one or more layers in - Regis are replaced by the geo_eenheden of geotop. + """Combine geotop and regis in such a way that the one or more layers in Regis are + replaced by the geo_eenheden of geotop. Parameters ---------- @@ -303,14 +301,13 @@ def add_geotop_to_regis_layers( def get_layer_names(): - """get all the available regis layer names. + """Get all the available regis layer names. Returns ------- layer_names : np.array array with names of all the regis layers. """ - layer_names = xr.open_dataset(REGIS_URL).layer.astype(str).values return layer_names diff --git a/nlmod/read/rws.py b/nlmod/read/rws.py index bc2f36e3..28b572b7 100644 --- a/nlmod/read/rws.py +++ b/nlmod/read/rws.py @@ -15,8 +15,8 @@ def get_gdf_surface_water(ds): - """read a shapefile with surface water as a geodataframe, cut by the extent - of the model. + """Read a shapefile with surface water as a geodataframe, cut by the extent of the + model. Parameters ---------- @@ -39,7 +39,7 @@ def get_gdf_surface_water(ds): @cache.cache_netcdf(coords_3d=True) def get_surface_water(ds, da_basename): - """create 3 data-arrays from the shapefile with surface water: + """Create 3 data-arrays from the shapefile with surface water: - area: area of the shape in the cell - cond: conductance based on the area and "bweerstand" column in shapefile @@ -58,7 +58,6 @@ def get_surface_water(ds, da_basename): ds : xarray.Dataset dataset with modelgrid data. """ - modelgrid = dims.modelgrid_from_ds(ds) gdf = get_gdf_surface_water(ds) @@ -93,8 +92,8 @@ def get_surface_water(ds, da_basename): @cache.cache_netcdf(coords_2d=True) def get_northsea(ds, da_name="northsea"): - """Get Dataset which is 1 at the northsea and 0 everywhere else. Sea is - defined by rws surface water shapefile. + """Get Dataset which is 1 at the northsea and 0 everywhere else. Sea is defined by + rws surface water shapefile. Parameters ---------- @@ -109,7 +108,6 @@ def get_northsea(ds, da_name="northsea"): Dataset with a single DataArray, this DataArray is 1 at sea and 0 everywhere else. Grid dimensions according to ds. """ - gdf_surf_water = get_gdf_surface_water(ds) # find grid cells with sea @@ -140,7 +138,6 @@ def add_northsea(ds, cachedir=None): b) fill top, bot, kh and kv add northsea cell by extrapolation c) get bathymetry (northsea depth) from jarkus. """ - logger.info( "Filling NaN values in top/botm and kh/kv in " "North Sea using bathymetry data from jarkus" @@ -181,8 +178,7 @@ def calculate_sea_coverage( nodata=-1, return_filled_dtm=False, ): - """ - Determine where the sea is by interpreting the digital terrain model. + """Determine where the sea is by interpreting the digital terrain model. This method assumes the pixel defined in xy_sea (by default top-left) of the DTM-DataArray is sea. It then determines the height of the sea that is required for @@ -190,7 +186,6 @@ def calculate_sea_coverage( Parameters ---------- - dtm : xr.DataArray The digital terrain data, which can be of higher resolution than ds, Nans are filled by the minial value of dtm. @@ -223,7 +218,6 @@ def calculate_sea_coverage( sea : xr.DataArray A DataArray with value of 1 where the sea is and 0 where it is not. """ - from skimage.morphology import reconstruction if not (dtm < zmax).any(): diff --git a/nlmod/read/waterboard.py b/nlmod/read/waterboard.py index 64de1fda..64eaa5a8 100644 --- a/nlmod/read/waterboard.py +++ b/nlmod/read/waterboard.py @@ -511,7 +511,6 @@ def get_data(wb, data_kind, extent=None, max_record_count=None, config=None, **k Raises ------ - DESCRIPTION. Returns @@ -594,11 +593,12 @@ def get_data(wb, data_kind, extent=None, max_record_count=None, config=None, **k def _set_column_from_columns(gdf, set_column, from_columns, nan_values=None): - """Retrieve values from one or more Geo)DataFrame-columns and set these - values as another column.""" + """Retrieve values from one or more Geo)DataFrame-columns and set these values as + another column. + """ if set_column in gdf.columns: raise (Exception(f"Column {set_column} allready exists")) - gdf[set_column] = np.NaN + gdf[set_column] = np.nan if from_columns is None: return gdf if isinstance(from_columns, str): @@ -634,5 +634,5 @@ def _set_column_from_columns(gdf, set_column, from_columns, nan_values=None): if nan_values is not None: if isinstance(nan_values, (float, int)): nan_values = [nan_values] - gdf.loc[gdf[set_column].isin(nan_values), set_column] = np.NaN + gdf.loc[gdf[set_column].isin(nan_values), set_column] = np.nan return gdf diff --git a/nlmod/read/webservices.py b/nlmod/read/webservices.py index 2e018824..97cc90bf 100644 --- a/nlmod/read/webservices.py +++ b/nlmod/read/webservices.py @@ -183,7 +183,7 @@ def arcrest( def _get_data(url, params, timeout=120, **kwargs): - """get data using a request + """Get data using a request. Parameters ---------- @@ -197,7 +197,6 @@ def _get_data(url, params, timeout=120, **kwargs): Returns ------- data - """ r = requests.get(url, params=params, timeout=timeout, **kwargs) if not r.ok: @@ -432,9 +431,8 @@ def _split_wcs_extent( fmt, crs, ): - """There is a max height and width limit for the wcs server. This function - splits your extent in chunks smaller than the limit. It returns a list of - Memory files. + """There is a max height and width limit for the wcs server. This function splits + your extent in chunks smaller than the limit. It returns a list of Memory files. Parameters ---------- @@ -463,12 +461,12 @@ def _split_wcs_extent( ------- MemoryFile Rasterio MemoryFile of the merged data + Notes ----- 1. The resolution is used to obtain the data from the wcs server. Not sure what kind of interpolation is used to resample the original grid. """ - # write tiles datasets = [] start_x = extent[0] diff --git a/nlmod/sim/__init__.py b/nlmod/sim/__init__.py index 1ca20b53..57b71d59 100644 --- a/nlmod/sim/__init__.py +++ b/nlmod/sim/__init__.py @@ -1 +1,2 @@ +# ruff: noqa: F403 from .sim import * diff --git a/nlmod/sim/sim.py b/nlmod/sim/sim.py index 052d2d9c..c33a3c64 100644 --- a/nlmod/sim/sim.py +++ b/nlmod/sim/sim.py @@ -13,9 +13,9 @@ def write_and_run(sim, ds, write_ds=True, script_path=None, silent=False): - """write modflow files and run the model. Extra options include writing the - model dataset to a netcdf file in the model workspace and copying the - modelscript to the model workspace. + """Write modflow files and run the model. Extra options include writing the model + dataset to a netcdf file in the model workspace and copying the modelscript to the + model workspace. Parameters ---------- @@ -108,7 +108,7 @@ def get_tdis_perioddata(ds, nstp="nstp", tsmult="tsmult"): def sim(ds, exe_name=None, version_tag=None): - """create sim from the model dataset. + """Create sim from the model dataset. Parameters ---------- @@ -131,7 +131,6 @@ def sim(ds, exe_name=None, version_tag=None): sim : flopy MFSimulation simulation object. """ - # start creating model logger.info("creating mf6 SIM") @@ -161,7 +160,7 @@ def sim(ds, exe_name=None, version_tag=None): def tdis(ds, sim, pname="tdis", nstp="nstp", tsmult="tsmult", **kwargs): - """create tdis package from the model dataset. + """Create tdis package from the model dataset. Parameters ---------- @@ -180,7 +179,6 @@ def tdis(ds, sim, pname="tdis", nstp="nstp", tsmult="tsmult", **kwargs): dis : flopy TDis tdis object. """ - # start creating model logger.info("creating mf6 TDIS") @@ -201,7 +199,7 @@ def tdis(ds, sim, pname="tdis", nstp="nstp", tsmult="tsmult", **kwargs): def ims(sim, complexity="MODERATE", pname="ims", **kwargs): - """create IMS package. + """Create IMS package. Parameters ---------- @@ -217,7 +215,6 @@ def ims(sim, complexity="MODERATE", pname="ims", **kwargs): ims : flopy ModflowIms ims object. """ - logger.info("creating mf6 IMS") print_option = kwargs.pop("print_option", "summary") diff --git a/nlmod/util.py b/nlmod/util.py index 0ea38fa8..e2f653cc 100644 --- a/nlmod/util.py +++ b/nlmod/util.py @@ -3,16 +3,16 @@ import os import re import sys -from pathlib import Path import warnings +from pathlib import Path from typing import Dict, Optional -from flopy.utils import get_modflow -from flopy.utils.get_modflow import flopy_appdata_path, get_release import geopandas as gpd import requests import xarray as xr from colorama import Back, Fore, Style +from flopy.utils import get_modflow +from flopy.utils.get_modflow import flopy_appdata_path, get_release from shapely.geometry import box logger = logging.getLogger(__name__) @@ -125,16 +125,16 @@ def get_exe_path( download_if_not_found : bool, optional Download the executables if they are not found, by default True. repo : str, default "executables" - Name of GitHub repository. Choose one of "executables" (default), - "modflow6", or "modflow6-nightly-build". If repo and version_tag are - provided the most recent installation location of MODFLOW is found in flopy metadata - that respects `version_tag` and `repo`. If not found, the executables are downloaded + Name of GitHub repository. Choose one of "executables" (default), "modflow6", + or "modflow6-nightly-build". If repo and version_tag are provided the most + recent installation location of MODFLOW is found in flopy metadata that + respects `version_tag` and `repo`. If not found, the executables are downloaded using repo and version_tag. version_tag : str, default None GitHub release ID: for example "18.0" or "latest". If repo and version_tag are - provided the most recent installation location of MODFLOW is found in flopy metadata - that respects `version_tag` and `repo`. If not found, the executables are downloaded - using repo and version_tag. + provided the most recent installation location of MODFLOW is found in flopy + metadata that respects `version_tag` and `repo`. If not found, the executables + are downloaded using repo and version_tag. Returns ------- @@ -181,8 +181,7 @@ def get_bin_directory( version_tag=None, repo="executables", ) -> Path: - """ - Get the directory where the executables are stored. + """Get the directory where the executables are stored. Searching for the executables is done in the following order: 0. If exe_name is a full path, return the full path of the executable. @@ -195,8 +194,8 @@ def get_bin_directory( Else: 4. Download the executables using `version_tag` and `repo`. - The returned directory is checked to contain exe_name if exe_name is provided. If exe_name - is set to None only the existence of the directory is checked. + The returned directory is checked to contain exe_name if exe_name is provided. If + exe_name is set to None only the existence of the directory is checked. Parameters ---------- @@ -207,17 +206,16 @@ def get_bin_directory( download_if_not_found : bool, optional Download the executables if they are not found, by default True. repo : str, default "executables" - Name of GitHub repository. Choose one of "executables" (default), - "modflow6", or "modflow6-nightly-build". If repo and version_tag are - provided the most recent installation location of MODFLOW is found in flopy metadata - that respects `version_tag` and `repo`. If not found, the executables are downloaded - using repo and version_tag. repo cannot be None. + Name of GitHub repository. Choose one of "executables" (default), "modflow6", + or "modflow6-nightly-build". If repo and version_tag are provided the most + recent installation location of MODFLOW is found in flopy metadata that + respects `version_tag` and `repo`. If not found, the executables are downloaded + using repo and version_tag. version_tag : str, default None GitHub release ID: for example "18.0" or "latest". If repo and version_tag are - provided the most recent installation location of MODFLOW is found in flopy metadata - that respects `version_tag` and `repo`. If not found, the executables are downloaded - using repo and version_tag. If version_tag is None, no version check is performed - on present executables and if no exe is found, the latest version is downloaded. + provided the most recent installation location of MODFLOW is found in flopy + metadata that respects `version_tag` and `repo`. If not found, the executables + are downloaded using repo and version_tag. Returns ------- @@ -248,7 +246,10 @@ def get_bin_directory( # If bindir is provided if bindir is not None and enable_version_check: - msg = "Incompatible arguments. If bindir is provided, unable to check the version." + msg = ( + "Incompatible arguments. If bindir is provided, " + "unable to check the version." + ) raise ValueError(msg) use_bindir = ( @@ -287,7 +288,7 @@ def get_bin_directory( download_mfbinaries( bindir=bindir, version_tag=version_tag if version_tag is not None else "latest", - repo=repo + repo=repo, ) # Rerun this function @@ -300,7 +301,10 @@ def get_bin_directory( ) else: - msg = f"Could not find {exe_name} in {bindir}, {nlmod_bindir} and {flopy_bindirs}." + msg = ( + f"Could not find {exe_name} in {bindir}, " + f"{nlmod_bindir} and {flopy_bindirs}." + ) raise FileNotFoundError(msg) @@ -315,15 +319,15 @@ def get_flopy_bin_directories(version_tag=None, repo="executables"): ---------- repo : str, default "executables" Name of GitHub repository. Choose one of "executables" (default), - "modflow6", or "modflow6-nightly-build". If repo and version_tag are - provided the most recent installation location of MODFLOW is found in flopy metadata - that respects `version_tag` and `repo`. If not found, the executables are downloaded - using repo and version_tag. + "modflow6", or "modflow6-nightly-build". If repo and version_tag are provided + the most recent installation location of MODFLOW is found in flopy metadata + that respects `version_tag` and `repo`. If not found, the executables are + downloaded using repo and version_tag. version_tag : str, default None GitHub release ID: for example "18.0" or "latest". If repo and version_tag are - provided the most recent installation location of MODFLOW is found in flopy metadata - that respects `version_tag` and `repo`. If not found, the executables are downloaded - using repo and version_tag. + provided the most recent installation location of MODFLOW is found in flopy + metadata that respects `version_tag` and `repo`. If not found, the executables + are downloaded using repo and version_tag. Returns ------- @@ -398,7 +402,6 @@ def download_mfbinaries(bindir=None, version_tag="latest", repo="executables"): "modflow6", or "modflow6-nightly-build". version_tag : str, default "latest" GitHub release ID. - """ if bindir is None: # Path objects are immutable so a copy is implied @@ -497,8 +500,10 @@ def get_da_from_da_ds(da_ds, dims=("y", "x"), data=None): def find_most_recent_file(folder, name, extension=".pklz"): - """Find the most recent file in a folder. File must startwith name and end width - extension. If you want to look for the most recent folder use extension = ''. + """Find the most recent file in a folder. + + File must startwith name and end width extension. If you want to look for the most + recent folder use extension = ''. Parameters ---------- @@ -514,7 +519,6 @@ def find_most_recent_file(folder, name, extension=".pklz"): newest_file : str name of the most recent file """ - i = 0 for file in os.listdir(folder): if file.startswith(name) and file.endswith(extension): @@ -551,7 +555,6 @@ def compare_model_extents(extent1, extent2): 1: extent1 is completely within extent2 2: extent2 is completely within extent1 """ - # option1 extent1 is completely within extent2 check_xmin = extent1[0] >= extent2[0] check_xmax = extent1[1] <= extent2[1] @@ -602,7 +605,6 @@ def polygon_from_extent(extent): polygon_ext : shapely.geometry.polygon.Polygon polygon of the extent. """ - bbox = (extent[0], extent[2], extent[1], extent[3]) polygon_ext = box(*tuple(bbox)) @@ -625,7 +627,6 @@ def gdf_from_extent(extent, crs="EPSG:28992"): gdf_extent : GeoDataFrame geodataframe with extent. """ - geom_extent = polygon_from_extent(extent) gdf_extent = gpd.GeoDataFrame(geometry=[geom_extent], crs=crs) @@ -633,8 +634,9 @@ def gdf_from_extent(extent, crs="EPSG:28992"): def gdf_within_extent(gdf, extent): - """Select only parts of the geodataframe within the extent. Only accepts Polygon and - Linestring geometry types. + """Select only parts of the geodataframe within the extent. + + Only accepts Polygon and Linestring geometry types. Parameters ---------- @@ -688,6 +690,7 @@ def get_google_drive_filename(fid, timeout=120): warnings.warn( "this function is no longer supported use the gdown package instead", DeprecationWarning, + stacklevel=1, ) if isinstance(id, requests.Response): @@ -714,6 +717,7 @@ def download_file_from_google_drive(fid, destination=None): warnings.warn( "this function is no longer supported use the gdown package instead", DeprecationWarning, + stacklevel=1, ) def get_confirm_token(response): @@ -805,14 +809,12 @@ def __init__( self, *args, colors: Optional[Dict[str, str]] = None, **kwargs ) -> None: """Initialize the formatter with specified format strings.""" - super().__init__(*args, **kwargs) self.colors = colors if colors else {} def format(self, record) -> str: """Format the specified record as text.""" - record.color = self.colors.get(record.levelname, "") record.reset = Style.RESET_ALL @@ -820,6 +822,18 @@ def format(self, record) -> str: def get_color_logger(level="INFO"): + """Get a logger with colored output. + + Parameters + ---------- + level : str, optional + The logging level to set for the logger. Default is "INFO". + + Returns + ------- + logger : logging.Logger + The configured logger object. + """ if level == "DEBUG": FORMAT = "{color}{levelname}:{name}.{funcName}:{lineno}:{message}{reset}" else: diff --git a/nlmod/version.py b/nlmod/version.py index 5e1c0340..e281423c 100644 --- a/nlmod/version.py +++ b/nlmod/version.py @@ -6,15 +6,12 @@ def show_versions() -> None: """Method to print the version of dependencies.""" - msg = ( - f"Python version: {python_version()}\n" - f"NumPy version: {metadata.version('numpy')}\n" - f"Xarray version: {metadata.version('xarray')}\n" - f"Matplotlib version: {metadata.version('matplotlib')}\n" - f"Flopy version: {metadata.version('flopy')}\n" + f"Python version : {python_version()}\n" + f"NumPy version : {metadata.version('numpy')}\n" + f"Xarray version : {metadata.version('xarray')}\n" + f"Matplotlib version : {metadata.version('matplotlib')}\n" + f"Flopy version : {metadata.version('flopy')}\n\n" + f"nlmod version : {__version__}" ) - - msg += f"\nnlmod version: {__version__}" - - return print(msg) + print(msg) diff --git a/pyproject.toml b/pyproject.toml index e01548c2..4867067e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,13 +57,7 @@ repository = "https://github.com/gwmod/nlmod" documentation = "https://nlmod.readthedocs.io/en/latest/" [project.optional-dependencies] -full = [ - "nlmod[knmi]", - "gdown", - "geocube", - "contextily", - "scikit-image", -] +full = ["nlmod[knmi]", "gdown", "geocube", "contextily", "scikit-image"] knmi = ["h5netcdf", "nlmod[grib]"] grib = ["cfgrib", "ecmwflibs"] test = ["pytest>=7", "pytest-cov", "pytest-dependency"] @@ -103,9 +97,32 @@ line-length = 88 profile = "black" [tool.ruff] -line-length = 88 +line-length = 88 extend-include = ["*.ipynb"] +[tool.ruff.lint] +# See: https://docs.astral.sh/ruff/rules/ +select = [ + "C4", # flake8-comprehensions + "E", # pycodestyle + "F", # pyflakes + "I", # isort + "PT", # pytest-style + "D", # pydocstyle + "B", # flake8-bugbear + "NPY", # numpy +] +ignore = [ + "D401", # Imperative mood for docstring. Be glad we have docstrings at all :P! + "D100", # Missing docstring in module. + "D104", # Missing docstring in public package. +] + +[tool.ruff.format] + +[tool.ruff.lint.pydocstyle] +convention = "numpy" + [tool.pytest.ini_options] addopts = "--strict-markers --durations=0 --cov-report xml:coverage.xml --cov nlmod -v" markers = ["notebooks: run notebooks", "slow: slow tests", "skip: skip tests"] diff --git a/tests/test_001_model.py b/tests/test_001_model.py index 13e2a66e..3f0e4c1d 100644 --- a/tests/test_001_model.py +++ b/tests/test_001_model.py @@ -75,7 +75,7 @@ def test_get_ds_variable_delrc(): ) -@pytest.mark.slow +@pytest.mark.slow() def test_create_small_model_grid_only(tmpdir, model_name="test"): extent = [98700.0, 99000.0, 489500.0, 489700.0] # extent, nrow, ncol = nlmod.read.regis.fit_extent_to_regis(extent, 100, 100) @@ -117,7 +117,7 @@ def test_create_small_model_grid_only(tmpdir, model_name="test"): ds.to_netcdf(os.path.join(tst_model_dir, "small_model.nc")) -@pytest.mark.slow +@pytest.mark.slow() def test_create_sea_model_grid_only(tmpdir, model_name="test"): extent = [95000.0, 105000.0, 494000.0, 500000.0] # extent, nrow, ncol = nlmod.read.regis.fit_extent_to_regis(extent, 100, 100) @@ -143,7 +143,7 @@ def test_create_sea_model_grid_only(tmpdir, model_name="test"): ds.to_netcdf(os.path.join(tst_model_dir, "basic_sea_model.nc")) -@pytest.mark.slow +@pytest.mark.slow() def test_create_sea_model_grid_only_delr_delc_50(tmpdir, model_name="test"): ds = get_ds_time_transient(tmpdir) extent = [95000.0, 105000.0, 494000.0, 500000.0] @@ -160,7 +160,7 @@ def test_create_sea_model_grid_only_delr_delc_50(tmpdir, model_name="test"): ds.to_netcdf(os.path.join(tst_model_dir, "sea_model_grid_50.nc")) -@pytest.mark.slow +@pytest.mark.slow() def test_create_sea_model(tmpdir): ds = xr.open_dataset( os.path.join(tst_model_dir, "basic_sea_model.nc"), mask_and_scale=False @@ -210,7 +210,7 @@ def test_create_sea_model(tmpdir): _ = nlmod.sim.write_and_run(sim, ds) -@pytest.mark.slow +@pytest.mark.slow() def test_create_sea_model_perlen_list(tmpdir): ds = xr.open_dataset(os.path.join(tst_model_dir, "basic_sea_model.nc")) @@ -280,7 +280,7 @@ def test_create_sea_model_perlen_list(tmpdir): nlmod.sim.write_and_run(sim, ds) -@pytest.mark.slow +@pytest.mark.slow() def test_create_sea_model_perlen_14(tmpdir): ds = xr.open_dataset(os.path.join(tst_model_dir, "basic_sea_model.nc")) diff --git a/tests/test_002_regis_geotop.py b/tests/test_002_regis_geotop.py index 8b46145c..dced58ba 100644 --- a/tests/test_002_regis_geotop.py +++ b/tests/test_002_regis_geotop.py @@ -1,4 +1,5 @@ import matplotlib.pyplot as plt + import nlmod diff --git a/tests/test_003_mfpackages.py b/tests/test_003_mfpackages.py index 86b08b1b..6857d9d2 100644 --- a/tests/test_003_mfpackages.py +++ b/tests/test_003_mfpackages.py @@ -94,7 +94,7 @@ def get_value_from_ds_datavar(): }, ) shape = list(ds.sizes.values()) - ds["test_var"] = ("layer", "y", "x"), np.arange(np.product(shape)).reshape(shape) + ds["test_var"] = ("layer", "y", "x"), np.arange(np.prod(shape)).reshape(shape) # get value from ds v0 = nlmod.util._get_value_from_ds_datavar( diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py index 709bdaf0..1d40d0d1 100644 --- a/tests/test_005_external_data.py +++ b/tests/test_005_external_data.py @@ -1,8 +1,8 @@ import pandas as pd -import xarray as xr -from shapely.geometry import LineString import pytest import test_001_model +import xarray as xr +from shapely.geometry import LineString import nlmod diff --git a/tests/test_007_run_notebooks.py b/tests/test_007_run_notebooks.py index 1ed78c4d..8b27e8b8 100644 --- a/tests/test_007_run_notebooks.py +++ b/tests/test_007_run_notebooks.py @@ -1,4 +1,5 @@ """run notebooks in the examples directory.""" +# ruff: noqa: D103 import os import nbformat @@ -19,91 +20,91 @@ def _run_notebook(nbdir, fname): return out -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_00_model_from_scratch(): _run_notebook(nbdir, "00_model_from_scratch.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_01_basic_model(): _run_notebook(nbdir, "01_basic_model.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_02_surface_water(): _run_notebook(nbdir, "02_surface_water.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_03_local_grid_refinement(): _run_notebook(nbdir, "03_local_grid_refinement.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_04_modifying_layermodels(): _run_notebook(nbdir, "04_modifying_layermodels.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_05_caching(): _run_notebook(nbdir, "05_caching.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_06_gridding_vector_data(): _run_notebook(nbdir, "06_gridding_vector_data.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_07_resampling(): _run_notebook(nbdir, "07_resampling.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_08_gis(): _run_notebook(nbdir, "08_gis.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_09_schoonhoven(): _run_notebook(nbdir, "09_schoonhoven.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_10_modpath(): _run_notebook(nbdir, "10_modpath.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_11_grid_rotation(): _run_notebook(nbdir, "11_grid_rotation.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_12_layer_generation(): _run_notebook(nbdir, "12_layer_generation.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_13_plot_methods(): _run_notebook(nbdir, "13_plot_methods.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_14_stromingen_example(): _run_notebook(nbdir, "14_stromingen_example.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_15_geotop(): _run_notebook(nbdir, "15_geotop.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_16_groundwater_transport(): _run_notebook(nbdir, "16_groundwater_transport.ipynb") -@pytest.mark.notebooks +@pytest.mark.notebooks() def test_run_notebook_17_unsaturated_zone_flow(): _run_notebook(nbdir, "17_unsaturated_zone_flow.ipynb") diff --git a/tests/test_008_waterschappen.py b/tests/test_008_waterschappen.py index 6b81bdc2..e0e952f3 100644 --- a/tests/test_008_waterschappen.py +++ b/tests/test_008_waterschappen.py @@ -3,7 +3,6 @@ import nlmod - # def test_download_polygons(): # is tested in test_024_administrative.test_get_waterboards # nlmod.read.waterboard.get_polygons() diff --git a/tests/test_009_layers.py b/tests/test_009_layers.py index ce19af49..3b304872 100644 --- a/tests/test_009_layers.py +++ b/tests/test_009_layers.py @@ -1,7 +1,7 @@ import os -import numpy as np import matplotlib.pyplot as plt +import numpy as np from shapely.geometry import LineString import nlmod diff --git a/tests/test_018_knmi_data_platform.py b/tests/test_018_knmi_data_platform.py index 9c53bbb0..db64b8c0 100644 --- a/tests/test_018_knmi_data_platform.py +++ b/tests/test_018_knmi_data_platform.py @@ -1,11 +1,14 @@ +# ruff: noqa: D103 import os from pathlib import Path +import pytest + from nlmod.read import knmi_data_platform data_path = Path(__file__).parent / "data" - +@pytest.mark.skip(reason="FileNotFoundError: download/INTER_OPER_R___EV24____L3__20240626T000000_20240627T000000_0003.nc not found") def test_download_multiple_nc_files() -> None: dataset_name = "EV24" dataset_version = "2" @@ -30,7 +33,7 @@ def test_download_multiple_nc_files() -> None: # plot the mean evaporation ds["prediction"].mean("time").plot() - +@pytest.mark.skip(reason="KeyError: 'files'") def test_download_read_zip_file() -> None: dataset_name = "rad_nl25_rac_mfbs_24h_netcdf4" dataset_version = "2.0" diff --git a/tests/test_019_attributes_encodings.py b/tests/test_019_attributes_encodings.py index 8e18ce38..03333f9d 100644 --- a/tests/test_019_attributes_encodings.py +++ b/tests/test_019_attributes_encodings.py @@ -1,5 +1,4 @@ import os -import time from tempfile import TemporaryDirectory import numpy as np diff --git a/tests/test_021_nhi.py b/tests/test_021_nhi.py index 6f5e06f8..8be1c5be 100644 --- a/tests/test_021_nhi.py +++ b/tests/test_021_nhi.py @@ -1,16 +1,19 @@ +# ruff: noqa: D103 import os -import numpy as np -import geopandas as gpd import tempfile -import nlmod -import pytest + +import geopandas as gpd import matplotlib.pyplot as plt +import numpy as np +import pytest + +import nlmod tmpdir = tempfile.gettempdir() -@pytest.mark.slow -def test_buidrainage(): +@pytest.mark.slow() +def test_buisdrainage(): model_ws = os.path.join(tmpdir, "buidrain") ds = nlmod.get_ds([110_000, 130_000, 435_000, 445_000], model_ws=model_ws) ds = nlmod.read.nhi.add_buisdrainage(ds) diff --git a/tests/test_022_gwt.py b/tests/test_022_gwt.py index 3f864c4f..5d9fce1e 100644 --- a/tests/test_022_gwt.py +++ b/tests/test_022_gwt.py @@ -1,7 +1,9 @@ -import tempfile import os +import tempfile + import pandas as pd import xarray as xr + import nlmod diff --git a/tests/test_023_hfb.py b/tests/test_023_hfb.py index f7e2a73f..1902f23d 100644 --- a/tests/test_023_hfb.py +++ b/tests/test_023_hfb.py @@ -1,8 +1,10 @@ -from shapely.geometry import LineString, Polygon -import geopandas as gpd +# ruff: noqa: D103 import flopy -import nlmod +import geopandas as gpd import util +from shapely.geometry import LineString, Polygon + +import nlmod def test_get_hfb_spd(): diff --git a/tests/test_025_modpath.py b/tests/test_025_modpath.py index 6e22c1a6..e18fb5f7 100644 --- a/tests/test_025_modpath.py +++ b/tests/test_025_modpath.py @@ -1,6 +1,8 @@ import os -import xarray as xr + import flopy +import xarray as xr + import nlmod diff --git a/tests/test_026_grid.py b/tests/test_026_grid.py index b646e5a0..80b001e7 100644 --- a/tests/test_026_grid.py +++ b/tests/test_026_grid.py @@ -1,9 +1,11 @@ -import tempfile import os -import numpy as np -import xarray as xr +import tempfile + import geopandas as gpd import matplotlib.pyplot as plt +import numpy as np +import xarray as xr + import nlmod model_ws = os.path.join(tempfile.gettempdir(), "test_grid") @@ -119,13 +121,13 @@ def test_vertex_da_to_ds(): def test_fillnan_da(): # for a structured grid ds = get_structured_model_ds() - ds["top"][5, 5] = np.NaN + ds["top"][5, 5] = np.nan top = nlmod.resample.fillnan_da(ds["top"], ds=ds) assert not np.isnan(top[5, 5]) # also for a vertex grid ds = get_vertex_model_ds() - ds["top"][100] = np.NaN + ds["top"][100] = np.nan mask = ds["top"].isnull() assert mask.any() top = nlmod.resample.fillnan_da(ds["top"], ds=ds)