: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 df309016..dc5594a5 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -4,11 +4,8 @@
import numpy as np
import rasterio
import xarray as xr
-from affine import Affine
from scipy.interpolate import griddata
from scipy.spatial import cKDTree
-from shapely.affinity import affine_transform
-from shapely.geometry import Polygon
from ..util import get_da_from_da_ds
@@ -16,8 +13,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
----------
@@ -102,8 +98,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
----------
@@ -117,13 +113,20 @@ def ds_to_structured_grid(
delc : int or float
cell size along columns of the desired grid (dy).
xorigin : int or float, optional
- lower left x coordinate of the model grid only used if angrot != 0.
- Default is 0.0.
+ lower left x coordinate of the model grid. When angrot == 0, xorigin is added to
+ the first two values of extent. Otherwise it is the x-coordinate of the point
+ the grid is rotated around, and xorigin is added to the Dataset-attributes.
+ The default is 0.0.
yorigin : int or float, optional
- lower left y coordinate of the model grid only used if angrot != 0.
- Default is 0.0.
+ lower left y coordinate of the model grid. When angrot == 0, yorigin is added to
+ the last two values of extent. Otherwise it is the y-coordinate of the point
+ the grid is rotated around, and yorigin is added to the Dataset-attributes.
+ The default is 0.0.
angrot : int or float, optinal
- the rotation of the grid in counter clockwise degrees, default is 0.0
+ the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid
+ is rotated, and all coordinates of the model are in model coordinates. See
+ https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more
+ infomation. The default is 0.0.
method : str, optional
type of interpolation used to resample. Sea structured_da_to_ds for
possible values of method. The default is 'nearest'.
@@ -134,21 +137,23 @@ 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"
if delc is None:
delc = delr
- x, y = get_xy_mid_structured(extent, delr, delc)
-
attrs = ds_in.attrs.copy()
_set_angrot_attributes(extent, xorigin, yorigin, angrot, attrs)
+ x, y = get_xy_mid_structured(attrs["extent"], delr, delc)
+
# add new attributes
attrs["gridtype"] = "structured"
+
if isinstance(delr, numbers.Number) and isinstance(delc, numbers.Number):
- attrs["delr"] = delr
- attrs["delc"] = delc
+ delr = np.full_like(x, delr)
+ delc = np.full_like(y, delc)
if method in ["nearest", "linear"] and angrot == 0.0:
ds_out = ds_in.interp(
@@ -164,22 +169,27 @@ 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
----------
extent : list, tuple or np.array of length 4
extent (xmin, xmax, ymin, ymax) of the desired grid.
- xorigin : float
- x-position of the lower-left corner of the model grid. Only used when angrot is
- not 0.
- yorigin : float
- y-position of the lower-left corner of the model grid. Only used when angrot is
- not 0.
- angrot : float
- counter-clockwise rotation angle (in degrees) of the lower-left corner of the
- model grid.
+ xorigin : int or float, optional
+ lower left x coordinate of the model grid. When angrot == 0, xorigin is added to
+ the first two values of extent. Otherwise it is the x-coordinate of the point
+ the grid is rotated around, and xorigin is added to the Dataset-attributes.
+ The default is 0.0.
+ yorigin : int or float, optional
+ lower left y coordinate of the model grid. When angrot == 0, yorigin is added to
+ the last two values of extent. Otherwise it is the y-coordinate of the point
+ the grid is rotated around, and yorigin is added to the Dataset-attributes.
+ The default is 0.0.
+ angrot : int or float, optinal
+ the rotation of the grid in counter clockwise degrees. When angrot != 0 the grid
+ is rotated, and all coordinates of the model are in model coordinates. See
+ https://nlmod.readthedocs.io/en/stable/examples/11_grid_rotation.html for more
+ infomation. The default is 0.0.
attrs : dict
Attributes of a model dataset.
@@ -217,7 +227,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
@@ -277,7 +287,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
@@ -307,6 +317,10 @@ 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}"
+ )
# get list of coordinates from all points in raster
if x is None:
@@ -334,7 +348,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
@@ -384,7 +398,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
@@ -431,7 +444,9 @@ def dim_to_regular_dim(da, dims, z):
coords = dict(da.coords)
coords["x"] = ds.x
coords["y"] = ds.y
- coords.pop("icell2d")
+ for key in list(coords):
+ if "icell2d" in coords[key].dims:
+ coords.pop(key)
else:
# just use griddata
z = griddata(points, da.data, xi, method=method)
@@ -440,7 +455,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
@@ -499,13 +514,12 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
# xmin, xmax, ymin, ymax
dx = (ds.attrs["extent"][1] - ds.attrs["extent"][0]) / len(ds.x)
dy = (ds.attrs["extent"][3] - ds.attrs["extent"][2]) / len(ds.y)
- elif "delr" in ds.attrs and "delc" in ds.attrs:
- dx = ds.attrs["delr"]
- dy = ds.attrs["delc"]
else:
raise ValueError(
"No extent or delr and delc in ds. Cannot determine affine."
)
+ from .grid import get_affine
+
da_out = da.rio.reproject(
dst_crs=ds.rio.crs,
shape=(len(ds.y), len(ds.x)),
@@ -525,6 +539,8 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
dims.remove("x")
dims.append("icell2d")
da_out = get_da_from_da_ds(ds, dims=tuple(dims), data=nodata)
+ from .grid import get_affine
+
for area in np.unique(ds["area"]):
dx = dy = np.sqrt(area)
x, y = get_xy_mid_structured(ds.extent, dx, dy)
@@ -565,94 +581,74 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
def extent_to_polygon(extent):
- """Generate a shapely Polygon from an extent ([xmin, xmax, ymin, ymax])"""
- nw = (extent[0], extent[2])
- no = (extent[1], extent[2])
- zo = (extent[1], extent[3])
- zw = (extent[0], extent[3])
- return Polygon([nw, no, zo, zw])
-
+ logger.warning(
+ "nlmod.resample.extent_to_polygon is deprecated. "
+ "Use nlmod.util.extent_to_polygon instead"
+ )
+ from ..util import extent_to_polygon
-def _get_attrs(ds):
- if isinstance(ds, dict):
- return ds
- else:
- return ds.attrs
+ return extent_to_polygon(extent)
def get_extent_polygon(ds, rotated=True):
"""Get the model extent, as a shapely Polygon."""
- attrs = _get_attrs(ds)
- polygon = extent_to_polygon(attrs["extent"])
- if rotated and "angrot" in ds.attrs and attrs["angrot"] != 0.0:
- affine = get_affine_mod_to_world(ds)
- polygon = affine_transform(polygon, affine.to_shapely())
- return polygon
+ logger.warning(
+ "nlmod.resample.get_extent_polygon is deprecated. "
+ "Use nlmod.grid.get_extent_polygon instead"
+ )
+ from .grid import get_extent_polygon
+
+ return get_extent_polygon(ds, rotated=rotated)
def affine_transform_gdf(gdf, affine):
"""Apply an affine transformation to a geopandas GeoDataFrame."""
- if isinstance(affine, Affine):
- affine = affine.to_shapely()
- gdfm = gdf.copy()
- gdfm.geometry = gdf.affine_transform(affine)
- return gdfm
+ logger.warning(
+ "nlmod.resample.affine_transform_gdf is deprecated. "
+ "Use nlmod.grid.affine_transform_gdf instead"
+ )
+ from .grid import affine_transform_gdf
+
+ return affine_transform_gdf(gdf, affine)
def get_extent(ds, rotated=True):
"""Get the model extent, corrected for angrot if necessary."""
- attrs = _get_attrs(ds)
- extent = attrs["extent"]
- if rotated and "angrot" in attrs and attrs["angrot"] != 0.0:
- affine = get_affine_mod_to_world(ds)
- xc = np.array([extent[0], extent[1], extent[1], extent[0]])
- yc = np.array([extent[2], extent[2], extent[3], extent[3]])
- xc, yc = affine * (xc, yc)
- extent = [xc.min(), xc.max(), yc.min(), yc.max()]
- return extent
+ logger.warning(
+ "nlmod.resample.get_extent is deprecated. Use nlmod.grid.get_extent instead"
+ )
+ from .grid import get_extent
+
+ return get_extent(ds, rotated=rotated)
def get_affine_mod_to_world(ds):
"""Get the affine-transformation from model to real-world coordinates."""
- attrs = _get_attrs(ds)
- xorigin = attrs["xorigin"]
- yorigin = attrs["yorigin"]
- angrot = attrs["angrot"]
- return Affine.translation(xorigin, yorigin) * Affine.rotation(angrot)
+ logger.warning(
+ "nlmod.resample.get_affine_mod_to_world is deprecated. "
+ "Use nlmod.grid.get_affine_mod_to_world instead"
+ )
+ from .grid import get_affine_mod_to_world
+
+ return get_affine_mod_to_world(ds)
def get_affine_world_to_mod(ds):
"""Get the affine-transformation from real-world to model coordinates."""
- attrs = _get_attrs(ds)
- xorigin = attrs["xorigin"]
- yorigin = attrs["yorigin"]
- angrot = attrs["angrot"]
- return Affine.rotation(-angrot) * Affine.translation(-xorigin, -yorigin)
+ logger.warning(
+ "nlmod.resample.get_affine_world_to_mod is deprecated. "
+ "Use nlmod.grid.get_affine_world_to_mod instead"
+ )
+ from .grid import get_affine_world_to_mod
+
+ return get_affine_world_to_mod(ds)
def get_affine(ds, sx=None, sy=None):
"""Get the affine-transformation, from pixel to real-world coordinates."""
- attrs = _get_attrs(ds)
- if sx is None:
- sx = attrs["delr"]
- if sy is None:
- sy = -attrs["delc"]
-
- if "angrot" in attrs:
- xorigin = attrs["xorigin"]
- yorigin = attrs["yorigin"]
- angrot = -attrs["angrot"]
- # xorigin and yorigin represent the lower left corner, while for the transform we
- # need the upper left
- dy = attrs["extent"][3] - attrs["extent"][2]
- xoff = xorigin + dy * np.sin(angrot * np.pi / 180)
- yoff = yorigin + dy * np.cos(angrot * np.pi / 180)
- return (
- Affine.translation(xoff, yoff)
- * Affine.scale(sx, sy)
- * Affine.rotation(angrot)
- )
- else:
- xoff = attrs["extent"][0]
- yoff = attrs["extent"][3]
- return Affine.translation(xoff, yoff) * Affine.scale(sx, sy)
+ logger.warning(
+ "nlmod.resample.get_affine is deprecated. Use nlmod.grid.get_affine instead"
+ )
+ from .grid import get_affine
+
+ return get_affine(ds, sx=sx, sy=sy)
diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py
index f766354a..9f481d63 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,13 +199,7 @@ 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."
- " For the old behavior, use `nlmod.time.set_ds_time_deprecated()`."
- )
-
if time is None and perlen is None:
raise (ValueError("Please specify either time or perlen in set_ds_time"))
elif perlen is not None:
@@ -280,7 +273,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 +348,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
@@ -408,6 +399,29 @@ def estimate_nstp(
return nstp_ceiled
+def get_time_step_length(perlen, nstp, tsmult):
+ """Get the length of the timesteps within a singe stress-period.
+
+ Parameters
+ ----------
+ perlen : float
+ The length of the stress period, in the time unit of the model (generally days).
+ nstp : int
+ The numer of timesteps within the stress period.
+ tsmult : float
+ THe time step multiplier, generally equal or lager than 1.
+
+ Returns
+ -------
+ 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()
+ return t
+
+
def ds_time_from_model(gwf):
warnings.warn(
"this function was renamed to `ds_time_idx_from_model`. "
@@ -430,7 +444,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)
@@ -456,7 +469,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,
@@ -530,3 +542,24 @@ def dataframe_to_flopy_timeseries(
time_series_namerecord=time_series_namerecord,
interpolation_methodrecord=interpolation_methodrecord,
)
+
+
+def ds_time_to_pandas_index(ds, include_start=True):
+ """Convert xarray time index to pandas datetime index.
+
+ Parameters
+ ----------
+ ds : xarray.Dataset
+ dataset with time index
+ include_start : bool, optional
+ include the start time in the index, by default True
+
+ Returns
+ -------
+ pd.DatetimeIndex
+ pandas datetime index
+ """
+ if include_start:
+ return ds.time.to_index().insert(0, pd.Timestamp(ds.time.start))
+ else:
+ return ds.time.to_index()
diff --git a/nlmod/epsg28992.py b/nlmod/epsg28992.py
index 6429fd9b..2c65a088 100644
--- a/nlmod/epsg28992.py
+++ b/nlmod/epsg28992.py
@@ -1,8 +1,10 @@
"""
-NOTE: this is the correct epsg:28992 definition for plotting backgroundmaps in RD
-More info (in Dutch) here:
-https://qgis.nl/2011/12/05/epsg28992-of-rijksdriehoekstelsel-verschuiving/
-This was still a problem in October 2023
+NOTE: this is the correct epsg:28992 definition for plotting backgroundmaps in RD.
+
+Related information (in Dutch):
+https://geoforum.nl/t/betrouwbare-bron-voor-proj4-definitie-van-rd-new-epsg-28992/5144/15
+This was still a problem in July 2024.
+
"""
EPSG_28992 = (
diff --git a/nlmod/gis.py b/nlmod/gis.py
index d9ee9ce3..c6739749 100644
--- a/nlmod/gis.py
+++ b/nlmod/gis.py
@@ -4,13 +4,16 @@
import geopandas as gpd
import numpy as np
-from .dims.grid import polygons_from_model_ds
-from .dims.resample import get_affine_mod_to_world
+from nlmod.dims.grid import get_affine_mod_to_world, polygons_from_model_ds
+from nlmod.dims.layers import calculate_thickness
+from nlmod.epsg28992 import EPSG_28992
logger = logging.getLogger(__name__)
-def vertex_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time="mean"):
+def vertex_da_to_gdf(
+ model_ds, data_variables, polygons=None, dealing_with_time="mean", crs=EPSG_28992
+):
"""Convert one or more DataArrays from a vertex model dataset to a Geodataframe.
Parameters
@@ -28,6 +31,9 @@ def vertex_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time=
becomes very slow. For now only the time averaged data will be
saved in the geodataframe. Later this can be extended with multiple
possibilities. The default is 'mean'.
+ crs : str, optional
+ coordinate reference system for the geodataframe. The default
+ is EPSG:28992 (RD).
Raises
------
@@ -65,7 +71,8 @@ def vertex_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time=
dv_dic[f"{da_name}_mean"] = da_mean.values
else:
raise NotImplementedError(
- "Can only use the mean of a DataArray with dimension time, use dealing_with_time='mean'"
+ "Can only use the mean of a DataArray with dimension time, "
+ "use dealing_with_time='mean'"
)
else:
raise ValueError(
@@ -73,7 +80,8 @@ def vertex_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time=
)
else:
raise NotImplementedError(
- f"expected one or two dimensions got {no_dims} for data variable {da_name}"
+ f"expected one or two dimensions got {no_dims} for "
+ f"data variable {da_name}"
)
# create geometries
@@ -81,12 +89,14 @@ def vertex_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time=
polygons = polygons_from_model_ds(model_ds)
# construct geodataframe
- gdf = gpd.GeoDataFrame(dv_dic, geometry=polygons)
+ gdf = gpd.GeoDataFrame(dv_dic, geometry=polygons, crs=crs)
return gdf
-def struc_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time="mean"):
+def struc_da_to_gdf(
+ model_ds, data_variables, polygons=None, dealing_with_time="mean", crs=EPSG_28992
+):
"""Convert one or more DataArrays from a structured model dataset to a Geodataframe.
Parameters
@@ -99,6 +109,9 @@ def struc_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time="
polygons : list of shapely Polygons, optional
geometries used for the GeoDataframe, if None the polygons are created
from the data variable 'vertices' in model_ds. The default is None.
+ crs : str, optional
+ coordinate reference system for the geodataframe. The default
+ is EPSG:28992 (RD).
Raises
------
@@ -152,7 +165,7 @@ def struc_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time="
polygons = polygons_from_model_ds(model_ds)
# construct geodataframe
- gdf = gpd.GeoDataFrame(dv_dic, geometry=polygons)
+ gdf = gpd.GeoDataFrame(dv_dic, geometry=polygons, crs=crs)
return gdf
@@ -173,7 +186,6 @@ def dataarray_to_shapefile(model_ds, data_variables, fname, polygons=None):
geometries used for the GeoDataframe, if None the polygons are created
from the data variable 'vertices' in model_ds. The default is None.
-
Returns
-------
None.
@@ -191,6 +203,7 @@ def ds_to_vector_file(
driver="GPKG",
combine_dic=None,
exclude=("x", "y", "time_steps", "area", "vertices", "rch_name", "icvert"),
+ crs=EPSG_28992,
):
"""Save all data variables in a model dataset to multiple shapefiles.
@@ -213,13 +226,15 @@ def ds_to_vector_file(
exclude : tuple of str, optional
data variables that are not exported to shapefiles. The default is
('x', 'y', 'time_steps', 'area', 'vertices').
+ crs : str, optional
+ coordinate reference system for the vector file. The default
+ is EPSG:28992 (RD).
Returns
-------
fnames : str or list of str
filename(s) of exported geopackage or shapefiles.
"""
-
# get default combination dictionary
if combine_dic is None:
combine_dic = {
@@ -265,9 +280,9 @@ def ds_to_vector_file(
for key, item in combine_dic.items():
if set(item).issubset(da_names):
if model_ds.gridtype == "structured":
- gdf = struc_da_to_gdf(model_ds, item, polygons=polygons)
+ gdf = struc_da_to_gdf(model_ds, item, polygons=polygons, crs=crs)
elif model_ds.gridtype == "vertex":
- gdf = vertex_da_to_gdf(model_ds, item, polygons=polygons)
+ gdf = vertex_da_to_gdf(model_ds, item, polygons=polygons, crs=crs)
if driver == "GPKG":
gdf.to_file(fname_gpkg, layer=key, driver=driver)
else:
@@ -283,9 +298,9 @@ def ds_to_vector_file(
# create unique shapefiles for the other data variables
for da_name in da_names:
if model_ds.gridtype == "structured":
- gdf = struc_da_to_gdf(model_ds, (da_name,), polygons=polygons)
+ gdf = struc_da_to_gdf(model_ds, (da_name,), polygons=polygons, crs=crs)
elif model_ds.gridtype == "vertex":
- gdf = vertex_da_to_gdf(model_ds, (da_name,), polygons=polygons)
+ gdf = vertex_da_to_gdf(model_ds, (da_name,), polygons=polygons, crs=crs)
if driver == "GPKG":
gdf.to_file(fname_gpkg, layer=da_name, driver=driver)
else:
@@ -308,6 +323,9 @@ def ds_to_ugrid_nc_file(
xv="xv",
yv="yv",
face_node_connectivity="icvert",
+ split_layer_dimension=True,
+ split_time_dimension=False,
+ for_imod_qgis_plugin=False,
):
"""Save a model dataset to a UGRID NetCDF file, so it can be opened as a Mesh Layer
in qgis.
@@ -335,13 +353,24 @@ def ds_to_ugrid_nc_file(
face_node_connectivity : str, optional
The name of the variable that contains the indexes of the vertices for
each face. The default is 'icvert'.
+ split_layer_dimension : bool, optional
+ Splits the layer dimension into seperate variables when True. The defaults is
+ True.
+ split_time_dimension : bool, optional
+ Splits the time dimension into seperate variables when True. The defaults is
+ False.
+ for_imod_qgis_plugin : bool, optional
+ When True, set some properties of the netcdf file to improve compatibility with
+ the iMOD-QGIS plugin. Layers are renamed to 'layer_i' until 'layer_n', a
+ variable 'top' is added for each layer, and the variable 'botm' is renamed to
+ 'bottom'. The default is False.
Returns
-------
ds : xr.DataSet
The dataset that was saved to a NetCDF-file. Can be used for debugging.
"""
- assert model_ds.gridtype == "vertex", "Only vertex grids are supported"
+ assert model_ds.gridtype == "vertex", "Only vertex grids are supported for now"
# copy the dataset, so we do not alter the original one
ds = model_ds.copy()
@@ -377,6 +406,10 @@ def ds_to_ugrid_nc_file(
ds[face_node_connectivity].attrs["cf_role"] = "face_node_connectivity"
ds[face_node_connectivity].attrs["start_index"] = 0
+ if for_imod_qgis_plugin and "botm" in ds:
+ ds["top"] = ds["botm"] + calculate_thickness(ds)
+ ds = ds.rename({"botm": "bottom"})
+
# set for each of the variables that they describe the faces
if variables is None:
variables = list(ds.keys())
@@ -405,9 +438,16 @@ def ds_to_ugrid_nc_file(
ds[var].encoding["dtype"] = np.int32
# Breaks down variables with a layer dimension into separate variables.
- ds, variables = _break_down_dimension(ds, variables, "layer")
- # Breaks down variables with a time dimension into separate variables.
- ds, variables = _break_down_dimension(ds, variables, "time")
+ if split_layer_dimension:
+ if for_imod_qgis_plugin:
+ ds, variables = _break_down_dimension(
+ ds, variables, "layer", add_dim_name=True, add_one_based_index=True
+ )
+ else:
+ ds, variables = _break_down_dimension(ds, variables, "layer")
+ if split_time_dimension:
+ # Breaks down variables with a time dimension into separate variables.
+ ds, variables = _break_down_dimension(ds, variables, "time")
# only keep the selected variables
ds = ds[variables + [dummy_var, xv, yv, face_node_connectivity]]
@@ -417,14 +457,26 @@ def ds_to_ugrid_nc_file(
return ds
-def _break_down_dimension(ds, variables, dim):
- # Copied and altered from imod-python.
+def _break_down_dimension(
+ ds, variables, dim, add_dim_name=False, add_one_based_index=False
+):
+ """Internal method to split a dimension of a variable into multiple variables.
+
+ Copied and altered from imod-python.
+ """
keep_vars = []
for var in variables:
if dim in ds[var].dims:
stacked = ds[var]
- for value in stacked[dim].values:
- name = f"{var}_{value}"
+ for i, value in enumerate(stacked[dim].values):
+ name = var
+ if add_dim_name:
+ name = f"{name}_{dim}"
+ if add_one_based_index:
+ name = f"{name}_{i+1}"
+ else:
+ name = f"{name}_{value}"
+
ds[name] = stacked.sel({dim: value}, drop=True)
if "long_name" in ds[name].attrs:
long_name = ds[name].attrs["long_name"]
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 393bb4a7..2afc8586 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -1,17 +1,12 @@
-# -*- coding: utf-8 -*-
-"""Created on Thu Jan 7 17:20:34 2021.
-
-@author: oebbe
-"""
import logging
import numbers
import warnings
import flopy
-import numpy as np
import xarray as xr
from ..dims import grid
+from ..dims.grid import get_delc, get_delr
from ..dims.layers import get_idomain
from ..sim import ims, sim, tdis
from ..util import _get_value_from_ds_attr, _get_value_from_ds_datavar
@@ -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
----------
@@ -109,11 +103,6 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs):
return disv(ds, model, length_units=length_units)
# check attributes
- for att in ["delr", "delc"]:
- if att in ds.attrs:
- if isinstance(ds.attrs[att], np.float32):
- ds.attrs[att] = float(ds.attrs[att])
-
if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:
xorigin = ds.attrs["xorigin"]
yorigin = ds.attrs["yorigin"]
@@ -135,8 +124,8 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs):
nlay=ds.sizes["layer"],
nrow=ds.sizes["y"],
ncol=ds.sizes["x"],
- delr=ds["delr"].values if "delr" in ds else ds.delr,
- delc=ds["delc"].values if "delc" in ds else ds.delc,
+ delr=get_delr(ds),
+ delc=get_delc(ds),
top=ds["top"].data,
botm=ds["botm"].data,
idomain=idomain,
@@ -154,8 +143,8 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs):
nlay=ds.sizes["layer"],
nrow=ds.sizes["y"],
ncol=ds.sizes["x"],
- delr=ds["delr"].values if "delr" in ds else ds.delr,
- delc=ds["delc"].values if "delc" in ds else ds.delc,
+ delr=get_delr(ds),
+ delc=get_delc(ds),
top=ds["top"].data,
botm=ds["botm"].data,
idomain=idomain,
@@ -169,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
----------
@@ -191,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
----------
@@ -272,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
----------
@@ -336,7 +325,7 @@ def ghb(
layer=None,
**kwargs,
):
- """create general head boundary from model dataset.
+ """Create general head boundary from model dataset.
Parameters
----------
@@ -419,7 +408,7 @@ def drn(
layer=None,
**kwargs,
):
- """create drain from model dataset.
+ """Create drain from model dataset.
Parameters
----------
@@ -492,7 +481,7 @@ def riv(
layer=None,
**kwargs,
):
- """create river package from model dataset.
+ """Create river package from model dataset.
Parameters
----------
@@ -575,7 +564,7 @@ def chd(
layer=0,
**kwargs,
):
- """create constant head package from model dataset.
+ """Create constant head package from model dataset.
Parameters
----------
@@ -650,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
----------
@@ -694,7 +683,7 @@ def sto(
pname="sto",
**kwargs,
):
- """create storage package from model dataset.
+ """Create storage package from model dataset.
Parameters
----------
@@ -746,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
----------
@@ -770,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)
@@ -798,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
----------
@@ -822,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
----------
@@ -847,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
----------
@@ -891,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
@@ -951,7 +939,7 @@ def oc(
pname="oc",
**kwargs,
):
- """create output control package from model dataset.
+ """Create output control package from model dataset.
Parameters
----------
@@ -1031,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 0c6cc07a..ffd0f053 100644
--- a/nlmod/gwf/horizontal_flow_barrier.py
+++ b/nlmod/gwf/horizontal_flow_barrier.py
@@ -2,16 +2,17 @@
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
+import xarray as xr
from shapely.geometry import Point, Polygon
from ..dims.grid import gdf_to_da, gdf_to_grid
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
@@ -100,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
----------
@@ -217,7 +218,9 @@ def line2hfb(gdf, gwf, prevent_rings=True, plot=False):
def polygon_to_hfb(
gdf, ds, column=None, gwf=None, lay=0, hydchr=1 / 100, add_data=False
):
- if isinstance(gdf, str):
+ if isinstance(gdf, xr.DataArray):
+ da = gdf
+ elif isinstance(gdf, str):
da = ds[gdf]
else:
if column is None:
@@ -275,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 ed4ed748..ac32f29e 100644
--- a/nlmod/gwf/output.py
+++ b/nlmod/gwf/output.py
@@ -7,18 +7,19 @@
import xarray as xr
from shapely.geometry import Point
-from ..dims.grid import modelgrid_from_ds
+from ..dims.grid import get_affine_world_to_mod, modelgrid_from_ds
from ..mfoutput.mfoutput import (
_get_budget_da,
+ _get_flopy_data_object,
+ _get_grb_file,
_get_heads_da,
_get_time_index,
- _get_flopy_data_object,
)
logger = logging.getLogger(__name__)
-def get_headfile(ds=None, gwf=None, fname=None, grbfile=None):
+def get_headfile(ds=None, gwf=None, fname=None, grb_file=None, **kwargs):
"""Get flopy HeadFile object.
Provide one of ds, gwf or fname.
@@ -31,7 +32,7 @@ def get_headfile(ds=None, gwf=None, fname=None, grbfile=None):
groundwater flow model, by default None
fname : str, optional
path to heads file, by default None
- grbfile : str
+ grb_file : str
path to file containing binary grid information
Returns
@@ -39,21 +40,21 @@ def get_headfile(ds=None, gwf=None, fname=None, grbfile=None):
flopy.utils.HeadFile
HeadFile object handle
"""
- return _get_flopy_data_object("head", ds, gwf, fname, grbfile)
+ return _get_flopy_data_object("head", ds, gwf, fname, grb_file, **kwargs)
def get_heads_da(
ds=None,
gwf=None,
fname=None,
- grbfile=None,
+ grb_file=None,
delayed=False,
chunked=False,
+ precision="auto",
**kwargs,
):
"""Read binary heads file.
-
Parameters
----------
ds : xarray.Dataset
@@ -62,20 +63,26 @@ def get_heads_da(
Flopy groundwaterflow object.
fname : path, optional
path to a binary heads file
- grbfile : str, optional
+ grb_file : str, optional
path to file containing binary grid information, only needed if reading
output from file using fname
delayed : bool, optional
if delayed is True, do not load output data into memory, default is False.
chunked : bool, optional
chunk data array containing output, default is False.
+ precision : str, optional
+ precision of floating point data in the head-file. Accepted values are 'auto',
+ 'single' or 'double'. When precision is 'auto', it is determined from the
+ head-file. Default is 'auto'.
Returns
-------
da : xarray.DataArray
heads data array.
"""
- hobj = get_headfile(ds=ds, gwf=gwf, fname=fname, grbfile=grbfile)
+ hobj = get_headfile(
+ ds=ds, gwf=gwf, fname=fname, grb_file=grb_file, precision=precision
+ )
# gwf.output.head() defaults to a structured grid
if gwf is not None and ds is None and fname is None:
kwargs["modelgrid"] = gwf.modelgrid
@@ -99,7 +106,7 @@ def get_heads_da(
return da
-def get_cellbudgetfile(ds=None, gwf=None, fname=None, grbfile=None):
+def get_cellbudgetfile(ds=None, gwf=None, fname=None, grb_file=None, **kwargs):
"""Get flopy CellBudgetFile object.
Provide one of ds, gwf or fname.
@@ -111,9 +118,9 @@ def get_cellbudgetfile(ds=None, gwf=None, fname=None, grbfile=None):
gwf : flopy.mf6.ModflowGwf, optional
groundwater flow model, by default None
fname_cbc : str, optional
- path to cell budget file, by default None\
- grbfile : str, optional
- path to file containing binary grid information, only needed if
+ path to cell budget file, by default None
+ grb_file : str, optional
+ path to file containing binary grid information, only needed if
fname_cbc is passed as only argument.
Returns
@@ -121,7 +128,7 @@ def get_cellbudgetfile(ds=None, gwf=None, fname=None, grbfile=None):
flopy.utils.CellBudgetFile
CellBudgetFile object handle
"""
- return _get_flopy_data_object("budget", ds, gwf, fname, grbfile)
+ return _get_flopy_data_object("budget", ds, gwf, fname, grb_file, **kwargs)
def get_budget_da(
@@ -129,10 +136,11 @@ def get_budget_da(
ds=None,
gwf=None,
fname=None,
- grbfile=None,
+ grb_file=None,
column="q",
delayed=False,
chunked=False,
+ precision="auto",
**kwargs,
):
"""Read binary budget file.
@@ -148,7 +156,7 @@ def get_budget_da(
fname : path, optional
specify the budget file to load, if not provided budget file will
be obtained from ds or gwf.
- grbfile : str
+ grb_file : str
path to file containing binary grid information, only needed if reading
output from file using fname
column : str
@@ -158,13 +166,19 @@ def get_budget_da(
if delayed is True, do not load output data into memory, default is False.
chunked : bool, optional
chunk data array containing output, default is False.
+ precision : str, optional
+ precision of floating point data in the budget-file. Accepted values are 'auto',
+ 'single' or 'double'. When precision is 'auto', it is determined from the
+ budget-file. Default is 'auto'.
Returns
-------
da : xarray.DataArray
budget data array.
"""
- cbcobj = get_cellbudgetfile(ds=ds, gwf=gwf, fname=fname, grbfile=grbfile)
+ cbcobj = get_cellbudgetfile(
+ ds=ds, gwf=gwf, fname=fname, grb_file=grb_file, precision=precision
+ )
da = _get_budget_da(cbcobj, text, column=column, **kwargs)
da.attrs["units"] = "m3/d"
@@ -186,9 +200,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
----------
@@ -232,7 +245,153 @@ def get_gwl_from_wet_cells(head, layer="layer", botm=None):
return gwl
-def get_head_at_point(head, x, y, ds=None, gi=None, drop_nan_layers=True):
+def get_flow_residuals(ds, gwf=None, fname=None, grb_file=None, kstpkper=None):
+ """Get the flow residuals of a MODFLOW 6 simulation.
+
+ Parameters
+ ----------
+ ds : xarray.Dataset
+ Xarray dataset with model data.
+ gwf : flopy ModflowGwf, optional
+ Flopy groundwaterflow object. One of ds or gwf must be provided.
+ fname : path, optional
+ specify the budget file to load, if not provided budget file will
+ be obtained from ds or gwf.
+ grb_file : str
+ The location of the grb-file. grb_file is determied from ds when None. The
+ default is None.
+ kstpkper : tuple of 2 ints, optional
+ The index of the timestep and the stress period to include in the result. Include
+ all data in the budget-file when None. The default is None.
+
+ Returns
+ -------
+ da : xr.DataArray
+ The flow residual in each cell, in m3/d.
+ """
+ if grb_file is None:
+ grb_file = _get_grb_file(ds)
+ grb = flopy.mf6.utils.MfGrdFile(grb_file)
+ cbf = get_cellbudgetfile(ds=ds, gwf=gwf, fname=fname, grb_file=grb_file)
+ dims = ds["botm"].dims
+ coords = ds["botm"].coords
+ flowja = cbf.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)
+ mask_active = np.diff(grb.ia) > 0
+ flowja_index = grb.ia[:-1][mask_active]
+ if kstpkper is None:
+ # loop over all timesteps/stress-periods
+ residuals = []
+ 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.ravel()[mask_active] = iflowja.flatten()[flowja_index]
+ residuals.append(residual)
+ dims = ("time",) + dims
+ coords = dict(coords) | {"time": _get_time_index(cbf, ds)}
+ 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.ravel()[mask_active] = flowja[0].flatten()[flowja_index]
+ da = xr.DataArray(residuals, dims=dims, coords=coords)
+ return da
+
+
+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.
+
+ 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.
+
+ Parameters
+ ----------
+ ds : xarray.Dataset
+ Xarray dataset with model data.
+ gwf : flopy ModflowGwf, optional
+ Flopy groundwaterflow object. One of ds or gwf must be provided.
+ fname : path, optional
+ specify the budget file to load, if not provided budget file will
+ be obtained from ds or gwf.
+ grb_file : str, optional
+ The location of the grb-file. grb_file is determied from ds when None. The
+ default is None.
+ kstpkper : tuple of 2 ints, optional
+ The index of the timestep and the stress period to include in the result. Include
+ all data in the budget-file when None. The default is None.
+ lays : int or list of ints, optional
+ The layers to include in the result. When lays is None, all layers are included.
+ The default is None.
+
+ Returns
+ -------
+ 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)
+ cbf = get_cellbudgetfile(ds=ds, gwf=gwf, fname=fname, grb_file=grb_file)
+ flowja = cbf.get_data(text="FLOW-JA-FACE", kstpkper=kstpkper)
+
+ if ds.gridtype == "vertex":
+ # determine flf_index first
+ grb = flopy.mf6.utils.MfGrdFile(grb_file)
+
+ if lays is None:
+ lays = range(grb.nlay)
+ if isinstance(lays, int):
+ lays = [lays]
+ shape = (len(lays), len(ds.icell2d))
+
+ flf_index = np.full(shape, -1)
+ # get these properties outside of the for loop to increase speed
+ grb_ia = grb.ia
+ grb_ja = grb.ja
+ for ilay, lay in enumerate(lays):
+ ja_start_next_layer = (lay + 1) * grb.ncpl
+ for icell2d in range(grb.ncpl):
+ node = lay * grb.ncpl + icell2d
+ ia = np.arange(grb_ia[node], grb_ia[node + 1])
+ mask = grb_ja[ia] >= ja_start_next_layer
+ if mask.any():
+ # assert mask.sum() == 1
+ flf_index[ilay, icell2d] = int(ia[mask][0])
+ coords = ds["botm"][lays].coords
+ else:
+ coords = ds["botm"].coords
+ dims = ds["botm"].dims
+
+ if kstpkper is None:
+ # loop over all tiesteps/stress-periods
+ flfs = []
+ for iflowja in flowja:
+ if ds.gridtype == "vertex":
+ flf = np.full(shape, np.nan)
+ mask = flf_index >= 0
+ flf[mask] = iflowja[0, 0, flf_index[mask]]
+ else:
+ _, _, flf = flopy.mf6.utils.get_structured_faceflows(iflowja, grb_file)
+ flfs.append(flf)
+ dims = ("time",) + dims
+ coords = dict(coords) | {"time": _get_time_index(cbf, ds)}
+ else:
+ if ds.gridtype == "vertex":
+ flfs = np.full(shape, np.nan)
+ mask = flf_index >= 0
+ flfs[mask] = flowja[0][0, 0, flf_index[mask]]
+ else:
+ _, _, flfs = flopy.mf6.utils.get_structured_faceflows(flowja[0], grb_file)
+ da = xr.DataArray(flfs, dims=dims, coords=coords)
+ if ds.gridtype != "vertex" and lays is not None:
+ da = da.isel(layer=lays)
+ return da
+
+
+def get_head_at_point(
+ head, x, y, ds=None, gi=None, drop_nan_layers=True, rotated=False
+):
"""Get the head at a certain point from a head DataArray for all cells.
Parameters
@@ -253,12 +412,20 @@ def get_head_at_point(head, x, y, ds=None, gi=None, drop_nan_layers=True):
None.
drop_nan_layers : bool, optional
Drop layers that are NaN at all timesteps. The default is True.
+ rotated : bool, optional
+ If the model grid has a rotation, and rotated is False, x and y are in model
+ coordinates. Otherwise x and y are in real world coordinates. The defaults is
+ False.
Returns
-------
head_point : xarray.DataArray
A DataArray with dimensions (time, layer).
"""
+ if rotated and "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:
+ # calculate model coordinates from the specified real-world coordinates
+ x, y = get_affine_world_to_mod(ds) * (x, y)
+
if "icell2d" in head.dims:
if gi is None:
if ds is None:
diff --git a/nlmod/gwf/recharge.py b/nlmod/gwf/recharge.py
index f8f29e55..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.
@@ -387,10 +385,11 @@ def ds_to_uzf(
if landflag is None:
landflag = xr.full_like(ds["botm"], 0, dtype=int)
# set the landflag in the top layer to 1
- fal = get_first_active_layer_from_idomain(idomain, nodata=0)
+ fal = get_first_active_layer_from_idomain(idomain)
+ # for the inactive domain set fal to 0 (setting nodata to 0 gives problems)
+ fal.data[fal == fal.nodata] = 0
landflag[fal] = 1
-
- # set landflag to 0 in inactivate domain
+ # set landflag to 0 in inactivate domain (where we set fal to 0 before)
landflag = xr.where(idomain > 0, landflag, 0)
# determine ivertcon, by setting its value to iuzno of the layer below
@@ -541,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 e5f6eca1..4d7acd4e 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -10,11 +10,16 @@
from shapely.strtree import STRtree
from tqdm import tqdm
-from ..dims.grid import gdf_to_grid
+from ..cache import cache_pickle
+from ..dims.grid import (
+ gdf_to_grid,
+ get_delc,
+ get_delr,
+ get_extent_polygon,
+)
from ..dims.layers import get_idomain
-from ..dims.resample import get_extent_polygon
from ..read import bgt, waterboard
-from ..cache import cache_pickle
+from ..util import extent_to_polygon
logger = logging.getLogger(__name__)
@@ -42,7 +47,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:
@@ -147,7 +151,11 @@ def agg_de_lange(group, cid, ds, c1=0.0, c0=1.0, N=1e-3, crad_positive=True):
# correction if group contains multiple shapes
# but covers whole cell
if group.area.sum() == A:
- li = A / np.max([ds.delr, ds.delc])
+ delr = get_delr(ds)
+ assert len(np.unique(delr)) == 1, "Variable grid size is not yet supported"
+ delc = get_delc(ds)
+ assert len(np.unique(delc)) == 1, "Variable grid size is not yet supported"
+ li = A / np.max([delr[0], delc[0]])
# width
B = group.area.sum(skipna=True) / li
@@ -296,9 +304,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:
@@ -366,6 +375,7 @@ def build_spd(
celldata : geopandas.GeoDataFrame
GeoDataFrame containing data. Cellid must be the index,
and must have columns "rbot", "stage" and "cond".
+ Optional columns are 'boundname' and 'aux'. These columns should have type str.
pkg : str
Modflow package: RIV, DRN or GHB
ds : xarray.DataSet
@@ -388,7 +398,6 @@ def build_spd(
- DRN: [(cellid), elev, cond]
- GHB: [(cellid), elev, cond]
"""
-
spd = []
top = ds.top.data
@@ -512,7 +521,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)]
@@ -571,15 +582,20 @@ def get_gdf_stage(gdf, season="winter"):
def download_level_areas(
- gdf, extent=None, config=None, raise_exceptions=True, **kwargs
+ gdf=None,
+ extent=None,
+ config=None,
+ raise_exceptions=True,
+ drop_duplicates=True,
+ **kwargs,
):
"""Download level areas (peilgebieden) of bronhouders.
Parameters
----------
- gdf : geopandas.GeoDataFrame
+ gdf : geopandas.GeoDataFrame, optional
A GeoDataFrame with surface water features, containing the column "bronhouder".
- extent : list, tuple or np.array
+ extent : list, tuple or np.array, optional
Model extent (xmin, xmax, ymin, ymax). When extent is None, all data of the
water boards in gdf are downloaded downloaded.
config : dict, optional
@@ -590,6 +606,9 @@ def download_level_areas(
Raises exceptions, mostly caused by a webservice that is offline. When
raise_exceptions is False, the error is raised as a warning. The default is
True.
+ drop_duplicates : bool, optional
+ Drop features with a duplicate index, keeping the first occurence. The default
+ is True.
Returns
-------
@@ -599,48 +618,57 @@ def download_level_areas(
"""
if config is None:
config = waterboard.get_configuration()
- bronhouders = gdf["bronhouder"].unique()
+ wbs = _get_waterboard_selection(gdf=gdf, extent=extent, config=config)
+
la = {}
data_kind = "level_areas"
- for wb in config.keys():
- if config[wb]["bgt_code"] in bronhouders:
- logger.info(f"Downloading {data_kind} for {wb}")
- try:
- lawb = waterboard.get_data(wb, data_kind, extent, **kwargs)
- if len(lawb) == 0:
- logger.info(f"No {data_kind} for {wb} found within model area")
- continue
- la[wb] = lawb
- mask = ~la[wb].is_valid
+ for wb in wbs:
+ logger.info(f"Downloading {data_kind} for {wb}")
+ try:
+ lawb = waterboard.get_data(wb, data_kind, extent, **kwargs)
+ if len(lawb) == 0:
+ logger.info(f"No {data_kind} for {wb} found within model area")
+ continue
+ if drop_duplicates:
+ mask = lawb.index.duplicated()
if mask.any():
- logger.warning(
- f"{mask.sum()} geometries of level areas of {wb} are invalid. Thet are made valid by adding a buffer of 0.0."
- )
- # first copy to prevent ValueError: assignment destination is read-only
- la[wb] = la[wb].copy()
- la[wb].loc[mask, "geometry"] = la[wb][mask].buffer(0.0)
- except Exception as e:
- if str(e) == f"{data_kind} not available for {wb}":
- logger.warning(e)
- elif raise_exceptions:
- raise
- else:
- logger.warning(e)
+ msg = "Dropping {} level area(s) of {} with duplicate indexes"
+ logger.warning(msg.format(mask.sum(), wb))
+ lawb = lawb.loc[~mask]
+
+ la[wb] = lawb
+ mask = ~la[wb].is_valid
+ if mask.any():
+ logger.warning(
+ f"{mask.sum()} geometries of level areas of {wb} are invalid. Thet are made valid by adding a buffer of 0.0."
+ )
+ # first copy to prevent ValueError: assignment destination is read-only
+ la[wb] = la[wb].copy()
+ la[wb].loc[mask, "geometry"] = la[wb][mask].buffer(0.0)
+ except Exception as e:
+ if str(e) == f"{data_kind} not available for {wb}":
+ logger.warning(e)
+ elif raise_exceptions:
+ raise
+ else:
+ logger.warning(e)
return la
def download_watercourses(
- gdf, extent=None, config=None, raise_exceptions=True, **kwargs
+ gdf=None, extent=None, config=None, raise_exceptions=True, **kwargs
):
"""Download watercourses of bronhouders.
Parameters
----------
- gdf : geopandas.GeoDataFrame
+ gdf : geopandas.GeoDataFrame, optional
A GeoDataFrame with surface water features, containing the column "bronhouder".
- extent : list, tuple or np.array
+ Determine the required waterboards for this gdf, when not None. The default is
+ None.
+ extent : list, tuple or np.array, optional
Model extent (xmin, xmax, ymin, ymax). When extent is None, all data of the
- water boards in gdf are downloaded downloaded.
+ water boards in gdf are downloaded downloaded. The default is None.
config : dict, optional
A dictionary with information about the webservices of the water boards. When
config is None, it is created with nlmod.read.waterboard.get_configuration().
@@ -658,28 +686,46 @@ def download_watercourses(
"""
if config is None:
config = waterboard.get_configuration()
- bronhouders = gdf["bronhouder"].unique()
+ wbs = _get_waterboard_selection(gdf=gdf, extent=extent, config=config)
wc = {}
data_kind = "watercourses"
- for wb in config.keys():
- if config[wb]["bgt_code"] in bronhouders:
- logger.info(f"Downloading {data_kind} for {wb}")
- try:
- wcwb = waterboard.get_data(wb, data_kind, extent, **kwargs)
- if len(wcwb) == 0:
- logger.info(f"No {data_kind} for {wb} found within model area")
- continue
- wc[wb] = wcwb
- except Exception as e:
- if str(e) == f"{data_kind} not available for {wb}":
- logger.warning(e)
- elif raise_exceptions:
- raise
- else:
- logger.warning(e)
+ for wb in wbs:
+ logger.info(f"Downloading {data_kind} for {wb}")
+ try:
+ wcwb = waterboard.get_data(wb, data_kind, extent, **kwargs)
+ if len(wcwb) == 0:
+ logger.info(f"No {data_kind} for {wb} found within model area")
+ continue
+ wc[wb] = wcwb
+ except Exception as e:
+ if str(e) == f"{data_kind} not available for {wb}":
+ logger.warning(e)
+ elif raise_exceptions:
+ raise
+ else:
+ logger.warning(e)
return wc
+def _get_waterboard_selection(gdf=None, extent=None, config=None):
+ """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:
+ raise (ValueError("Please specify either gdf or extent"))
+
+ if gdf is not None:
+ bronhouders = gdf["bronhouder"].unique()
+ wbs = []
+ for wb in config.keys():
+ if config[wb]["bgt_code"] in bronhouders:
+ wbs.append(wb)
+ elif extent is not None:
+ wb_gdf = waterboard.get_polygons()
+ wbs = wb_gdf.index[wb_gdf.intersects(extent_to_polygon(extent))]
+ return wbs
+
+
def add_stages_from_waterboards(
gdf, la=None, extent=None, columns=None, config=None, min_total_overlap=0.0
):
@@ -719,7 +765,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
@@ -773,7 +819,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
@@ -843,8 +889,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
----------
@@ -865,7 +911,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
@@ -877,7 +922,7 @@ def add_min_ahn_to_gdf(gdf, ahn, buffer=0.0, column="ahn_min"):
rasterize_function=partial(rasterize_image, all_touched=True),
)
gc["ahn"] = ahn
-
+ gc = gc.set_coords("index")
ahn_min = gc.groupby("index").min()["ahn"].to_pandas()
ahn_min.index = ahn_min.index.astype(int)
gdf[column] = ahn_min
@@ -955,7 +1000,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(
@@ -999,10 +1044,11 @@ def gdf_to_seasonal_pkg(
# the value to scale is also represented with a time series
# So we switch the conductance (column 2) and the multiplier (column 3/4)
spd = np.array(spd, dtype=object)
- if pkg == "RIV":
- spd[:, [2, 4]] = spd[:, [4, 2]]
- else:
- spd[:, [2, 3]] = spd[:, [3, 2]]
+ if len(spd) > 0:
+ if pkg == "RIV":
+ spd[:, [2, 4]] = spd[:, [4, 2]]
+ else:
+ spd[:, [2, 3]] = spd[:, [3, 2]]
spd = spd.tolist()
if boundname_column is None:
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/binaryfile.py b/nlmod/mfoutput/binaryfile.py
index a3197f0c..9d228a1b 100644
--- a/nlmod/mfoutput/binaryfile.py
+++ b/nlmod/mfoutput/binaryfile.py
@@ -138,13 +138,24 @@ def _get_binary_budget_data(kstpkper, fobj, text, column="q"):
idx = np.array([idx])
header = fobj.recordarray[idx]
- ipos = fobj.iposarray[idx].item()
- imeth = header["imeth"][0]
t = header["text"][0]
if isinstance(t, bytes):
t = t.decode("utf-8")
+ data = []
+ for ipos in fobj.iposarray[idx]:
+ data.append(_get_binary_budget_record(fobj, ipos, header, column))
+
+ if len(data) == 1:
+ return data[0]
+ else:
+ return np.ma.sum(data, axis=0)
+
+
+def _get_binary_budget_record(fobj, ipos, header, column):
+ """Get a single data record from the budget file."""
+ imeth = header["imeth"][0]
nlay = abs(header["nlay"][0])
nrow = header["nrow"][0]
ncol = header["ncol"][0]
diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py
index 64158e3d..f1d06eb9 100644
--- a/nlmod/mfoutput/mfoutput.py
+++ b/nlmod/mfoutput/mfoutput.py
@@ -1,14 +1,16 @@
-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
+from ..dims.grid import (
+ get_affine_mod_to_world,
+ get_dims_coords_from_modelgrid,
+ modelgrid_from_ds,
+)
from ..dims.time import ds_time_idx
from .binaryfile import _get_binary_budget_data, _get_binary_head_data
@@ -71,8 +73,8 @@ def _get_time_index(fobj, ds=None, gwf_or_gwt=None):
elif ds is not None:
tindex = ds_time_idx(
fobj.get_times(),
- start_datetime=ds.time.attrs["start"],
- time_units=ds.time.attrs["time_units"],
+ start_datetime=(ds.time.attrs["start"] if "time" in ds else None),
+ time_units=(ds.time.attrs["time_units"] if "time" in ds else None),
)
return tindex
@@ -239,9 +241,11 @@ def _get_budget_da(
return da
-def _get_flopy_data_object(var, ds=None, gwml=None, fname=None, grbfile=None):
+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.
@@ -255,7 +259,7 @@ def _get_flopy_data_object(var, ds=None, gwml=None, fname=None, grbfile=None):
groundwater flow or transport model, by default None
fname : str, optional
path to Head- or CellBudgetFile, by default None
- grbfile : str, optional
+ grb_file : str, optional
path to file containing binary grid information, if None modelgrid
information is obtained from ds. By default None
@@ -283,14 +287,11 @@ def _get_flopy_data_object(var, ds=None, gwml=None, fname=None, grbfile=None):
# return gwf.output.head(), gwf.output.budget() or gwt.output.concentration()
return getattr(gwml.output, var)()
fname = os.path.join(ds.model_ws, ds.model_name + extension)
- if grbfile is None and ds is not None:
+ if grb_file is None and ds is not None:
# get grb file
- if ds.gridtype == "vertex":
- grbfile = os.path.join(ds.model_ws, ds.model_name + ".disv.grb")
- elif ds.gridtype == "structured":
- grbfile = os.path.join(ds.model_ws, ds.model_name + ".dis.grb")
- if grbfile is not None and os.path.exists(grbfile):
- modelgrid = flopy.mf6.utils.MfGrdFile(grbfile).modelgrid
+ grb_file = _get_grb_file(ds)
+ if grb_file is not None and os.path.exists(grb_file):
+ modelgrid = flopy.mf6.utils.MfGrdFile(grb_file).modelgrid
elif ds is not None:
modelgrid = modelgrid_from_ds(ds)
else:
@@ -301,9 +302,17 @@ def _get_flopy_data_object(var, ds=None, gwml=None, fname=None, grbfile=None):
if modelgrid is None:
logger.error(msg)
raise ValueError(msg)
- return flopy.utils.CellBudgetFile(fname, modelgrid=modelgrid)
+ return flopy.utils.CellBudgetFile(fname, modelgrid=modelgrid, **kwargs)
else:
if modelgrid is None:
logger.warning(msg)
warnings.warn(msg)
- return flopy.utils.HeadFile(fname, text=var, modelgrid=modelgrid)
+ return flopy.utils.HeadFile(fname, text=var, modelgrid=modelgrid, **kwargs)
+
+
+def _get_grb_file(ds):
+ if ds.gridtype == "vertex":
+ grb_file = os.path.join(ds.model_ws, ds.model_name + ".disv.grb")
+ elif ds.gridtype == "structured":
+ grb_file = os.path.join(ds.model_ws, ds.model_name + ".dis.grb")
+ return grb_file
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 125c3aa3..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
----------
@@ -215,9 +214,11 @@ def mpf(gwf, exe_name=None, modelname=None, model_ws=None):
"the save_flows option of the npf package should be True not None"
)
- # get executable
+ # get executable. version_tag not supported yet
if exe_name is None:
- exe_name = util.get_exe_path("mp7_2_002_provisional")
+ exe_name = util.get_exe_path(exe_name="mp7_2_002_provisional")
+ else:
+ exe_name = util.get_exe_path(exe_name=exe_name)
# create mpf model
mpf = flopy.modpath.Modpath7(
@@ -262,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
----------
@@ -454,6 +453,7 @@ def sim(
simulationtype="combined",
weaksinkoption="pass_through",
weaksourceoption="pass_through",
+ **kwargs,
):
"""Create a modpath backward simulation from a particle group.
@@ -503,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..1563924f 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 52fb7fde..0a879445 100644
--- a/nlmod/plot/dcs.py
+++ b/nlmod/plot/dcs.py
@@ -1,16 +1,23 @@
+import logging
+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 matplotlib.animation import FFMpegWriter, FuncAnimation
from matplotlib.collections import LineCollection, PatchCollection
from matplotlib.patches import Rectangle
from shapely.affinity import affine_transform
from shapely.geometry import LineString, MultiLineString, Point, Polygon
-from ..dims.grid import modelgrid_from_ds
-from ..dims.resample import get_affine_world_to_mod
+from ..dims.grid import get_affine_world_to_mod, modelgrid_from_ds
+from .plotutil import get_map
+
+logger = logging.getLogger(__name__)
class DatasetCrossSection:
@@ -322,6 +329,93 @@ def plot_grid(
self.ax.add_collection(patch_collection)
return patch_collection
+ def plot_map_cs(
+ self,
+ ax=None,
+ figsize=5,
+ background=True,
+ lw=5,
+ ls="--",
+ label="cross section",
+ **kwargs,
+ ):
+ """Creates a different figure with the map of the cross section.
+
+ Parameters
+ ----------
+ ax : None or matplotlib.Axes, optional
+ if None a new axis object is created using nlmod.plot.get_map()
+ figsize : int, optional
+ size of the figure, only used if ax is None, by default 5
+ background : bool, optional
+ add a backgroun map, only used if ax is None, by default True
+ lw : int, optional
+ linewidth of the cross section, by default 10
+ ls : str, optional
+ linestyle of the cross section, by default "--"
+ label : str, optional
+ label of the cross section, by default "cross section"
+ **kwargs are passed to the nlmod.plot.get_map() function. Only if ax is None
+
+ Returns
+ -------
+ matplotlib Axes
+ axes
+ """
+ if ax is None:
+ _, ax = get_map(
+ self.ds.extent, background=background, figsize=figsize, **kwargs
+ )
+ gpd.GeoDataFrame(geometry=[self.line]).plot(ax=ax, ls=ls, lw=lw, label=label)
+ ax.legend()
+
+ return ax
+
+ def get_patches_array(self, z):
+ """Similar to plot_array function, only computes the array to update an existing plot_array.
+
+ Parameters
+ ----------
+ z : DataArray
+ data to plot on the patches.
+
+ Returns
+ -------
+ list
+ plot data.
+ """
+ if isinstance(z, xr.DataArray):
+ z = z.data
+
+ if self.icell2d in self.ds.dims:
+ assert len(z.shape) == 2
+ assert z.shape[0] == len(self.layer)
+ assert z.shape[1] == len(self.ds[self.icell2d])
+
+ zcs = z[:, self.icell2ds]
+ else:
+ assert len(z.shape) == 3
+ assert z.shape[0] == len(self.layer)
+ assert z.shape[1] == len(self.ds[self.y])
+ assert z.shape[2] == len(self.ds[self.x])
+
+ zcs = z[:, self.rows, self.cols]
+
+ array = []
+ for i in range(zcs.shape[0]):
+ for j in range(zcs.shape[1]):
+ if not (
+ np.isnan(self.top[i, j])
+ or np.isnan(self.bot[i, j])
+ or np.isnan(zcs[i, j])
+ ):
+ if self.bot[i, j] == self.zmax or self.top[i, j] == self.zmin:
+ continue
+
+ array.append(zcs[i, j])
+
+ return array
+
def plot_array(self, z, head=None, **kwargs):
if isinstance(z, xr.DataArray):
z = z.data
@@ -415,3 +509,99 @@ def get_top_and_bot(self, top, bot):
top[top > self.zmax] = self.zmax
bot[bot > self.zmax] = self.zmax
return top, bot
+
+ def animate(
+ self,
+ da,
+ cmap="Spectral_r",
+ norm=None,
+ head=None,
+ plot_title="",
+ date_fmt="%Y-%m-%d",
+ cbar_label=None,
+ fname=None,
+ ):
+ """Animate a cross section.
+
+ Parameters
+ ----------
+ da : DataArray
+ should have dimensions structured: time, y, x or vertex: time, icell2d
+ cmap : str, optional
+ passed to plot_array function, by default "Spectral_r"
+ norm : , optional
+ norm for the colorbar of the datarray, by default None
+ head : DataArray, optional
+ If not given the top cell is completely filled, by default None
+ plot_title : str or None, optional
+ if not None a title is added which is updated with every timestep (using
+ date_fmt for the date format), by default ""
+ date_fmt : str, optional
+ date format for plot title, by default "%Y-%m-%d"
+ cbar_label : str, optional
+ label for the colorbar, by default None
+ fname : str or Path, optional
+ filename if not None this is where the aniation is saved as mp4, by
+ default None
+
+ Returns
+ -------
+ matplotlib.animation.FuncAnimation
+ animation object
+ """
+ f = self.ax.get_figure()
+
+ # plot first timeframe
+ iper = 0
+ if head is not None:
+ plot_head = head[iper].values
+ logger.info("varying head not supported for animation yet")
+
+ pc = self.plot_array(da[iper].squeeze(), cmap=cmap, norm=norm, head=plot_head)
+ cbar = f.colorbar(pc, ax=self.ax, shrink=1.0)
+
+ if cbar_label is not None:
+ cbar.set_label(cbar_label)
+ elif "units" in da.attrs:
+ cbar.set_label(da.units)
+
+ t = pd.Timestamp(da.time.values[iper])
+ if plot_title is None:
+ title = None
+ else:
+ title = self.ax.set_title(f"{plot_title}, t = {t.strftime(date_fmt)}")
+
+ # update func
+ def update(iper, pc, title):
+ array = self.get_patches_array(da[iper].squeeze())
+ pc.set_array(array)
+
+ # update title
+ t = pd.Timestamp(da.time.values[iper])
+ if title is not None:
+ title.set_text(f"{plot_title}, t = {t.strftime(date_fmt)}")
+
+ return pc, title
+
+ # create animation
+ anim = FuncAnimation(
+ f,
+ partial(update, pc=pc, title=title),
+ frames=da["time"].shape[0],
+ blit=False,
+ interval=100,
+ )
+
+ # save animation
+ if fname is None:
+ return anim
+ else:
+ # save animation as mp4
+ writer = FFMpegWriter(
+ fps=10,
+ bitrate=-1,
+ extra_args=["-pix_fmt", "yuv420p"],
+ codec="libx264",
+ )
+ anim.save(fname, writer=writer)
+ return anim
diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py
index 176dc81b..30407da2 100644
--- a/nlmod/plot/plot.py
+++ b/nlmod/plot/plot.py
@@ -1,19 +1,24 @@
+import logging
import warnings
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
from matplotlib.patches import Patch
from mpl_toolkits.axes_grid1 import make_axes_locatable
-from ..dims.grid import modelgrid_from_ds
-from ..dims.resample import get_affine_mod_to_world, get_extent
+from ..dims.grid import (
+ get_affine_mod_to_world,
+ get_extent,
+ get_extent_gdf,
+ modelgrid_from_ds,
+)
from ..read import geotop, rws
from .dcs import DatasetCrossSection
from .plotutil import (
@@ -24,6 +29,8 @@
title_inside,
)
+logger = logging.getLogger(__name__)
+
def surface_water(model_ds, ax=None, **kwargs):
surf_water = rws.get_gdf_surface_water(model_ds)
@@ -38,9 +45,52 @@ def surface_water(model_ds, ax=None, **kwargs):
def modelgrid(ds, ax=None, **kwargs):
if ax is None:
_, ax = plt.subplots(figsize=(10, 10))
- ax.axis("scaled")
+ ax.set_aspect("auto")
modelgrid = modelgrid_from_ds(ds)
+ extent = None if ax.get_autoscale_on() else ax.axis()
modelgrid.plot(ax=ax, **kwargs)
+ if extent is not None:
+ ax.axis(extent)
+
+ return ax
+
+
+def modelextent(ds, dx=None, ax=None, rotated=False, **kwargs):
+ """Plot model extent.
+
+ Parameters
+ ----------
+ ds : xarray.Dataset
+ The dataset containing the data.
+ dx : float, optional
+ The buffer around the model extent. Default is 5% of the longest model edge.
+ ax : matplotlib.axes.Axes, optional
+ The axes object to plot on. If not provided, a new figure and axes will be
+ created.
+ rotated : bool, optional
+ When True, plot the model extent in real-world coordinates for rotated grids.
+ The default is False, which plots the model extent in local coordinates.
+ **kwargs
+ Additional keyword arguments to pass to the boundary plot.
+
+ Returns
+ -------
+ ax : matplotlib.axes.Axes
+ axes object
+ """
+ extent = xmin, xmax, ymin, ymax = get_extent(ds, rotated=rotated)
+ if dx is None:
+ dx = max(0.05 * (xmax - xmin), 0.05 * (ymax - ymin))
+ if ax is None:
+ _, ax = plt.subplots(figsize=(10, 10))
+ ax.axis("scaled")
+
+ ax.axis([xmin - dx, xmax + dx, ymin - dx, ymax + dx])
+ gdf = get_extent_gdf(ds, rotated=rotated)
+ extent = None if ax.get_autoscale_on() else ax.axis()
+ gdf.boundary.plot(ax=ax, **kwargs)
+ if extent is not None:
+ ax.axis(extent)
return ax
@@ -55,7 +105,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
----------
@@ -88,7 +138,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,
@@ -166,7 +215,8 @@ def data_array(da, ds=None, ax=None, rotated=False, edgecolor=None, **kwargs):
ax : matplotlib.Axes, optional
The axes used for plotting. Set to current axes when None. The default is None.
rotated : bool, optional
- Plot the data-array in rotated coordinates
+ When True, plot the data-array in real-world coordinates for rotated grids.
+ The default is False, which plots the data-array in local coordinates.
**kwargs : cit
Kwargs are passed to PatchCollection (vertex) or pcolormesh (structured).
@@ -212,7 +262,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, **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.
@@ -234,6 +291,9 @@ def geotop_lithok_in_cross_section(
lithok_props : pd.DataFrame, optional
A DataFrame containing the properties of the lithoclasses.
Will call nlmod.read.geotop.get_lithok_props() when None. The default is None.
+ alpha : float, optional
+ Opacity for plot_array function, The default is None.
+
**kwargs : dict
kwargs are passed onto DatasetCrossSection.
@@ -261,7 +321,8 @@ 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)
+ 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)
@@ -298,7 +359,6 @@ def geotop_lithok_on_map(
Returns
-------
qm : matplotlib.collections.QuadMesh
-
"""
if ax is None:
ax = plt.gca()
@@ -322,7 +382,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
@@ -337,10 +397,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)
@@ -351,7 +411,7 @@ def _get_geotop_cmap_and_norm(lithok, lithok_props):
return array, cmap, norm
-def _get_figure(ax=None, da=None, ds=None, figsize=None, rotated=True, extent=None):
+def _get_figure(ax=None, da=None, ds=None, figsize=None, rotated=False, extent=None):
# figure
if ax is not None:
f = ax.figure
@@ -402,7 +462,7 @@ def map_array(
colorbar=True,
colorbar_label="",
plot_grid=True,
- rotated=True,
+ rotated=False,
add_to_plot=None,
background=False,
figsize=None,
@@ -463,7 +523,10 @@ def map_array(
# bgmap
if background:
- add_background_map(ax, map_provider="nlmaps.water", alpha=0.5)
+ if not rotated and "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:
+ logger.warning("Background map not supported in in model coordinates")
+ else:
+ add_background_map(ax, map_provider="nlmaps.water", alpha=0.5)
# add other info to plot
if add_to_plot is not None:
@@ -522,7 +585,7 @@ def animate_map(
colorbar=True,
colorbar_label="",
plot_grid=True,
- rotated=True,
+ rotated=False,
background=False,
figsize=None,
ax=None,
@@ -569,7 +632,7 @@ def animate_map(
plot_grid : bool, optional
Whether to plot the model grid. Default is True.
rotated : bool, optional
- Whether to plot rotated model, if applicable. Default is True.
+ Whether to plot rotated model, if applicable. Default is False.
background : bool, optional
Whether to add a background map. Default is False.
figsize : tuple, optional
diff --git a/nlmod/plot/plotutil.py b/nlmod/plot/plotutil.py
index 00c134e3..3ebe5756 100644
--- a/nlmod/plot/plotutil.py
+++ b/nlmod/plot/plotutil.py
@@ -3,7 +3,7 @@
from matplotlib.patches import Polygon
from matplotlib.ticker import FuncFormatter, MultipleLocator
-from ..dims.resample import get_affine_mod_to_world
+from ..dims.grid import get_affine_mod_to_world
from ..epsg28992 import EPSG_28992
@@ -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 c2444e79..25f0370f 100644
--- a/nlmod/read/ahn.py
+++ b/nlmod/read/ahn.py
@@ -1,6 +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 rasterio
import rioxarray
import xarray as xr
@@ -9,14 +13,15 @@
from tqdm import tqdm
from .. import cache
-from ..dims.resample import get_extent, structured_da_to_ds
+from ..dims.grid import get_extent
+from ..dims.resample import structured_da_to_ds
from ..util import get_ds_empty
from .webservices import arcrest, wcs
logger = logging.getLogger(__name__)
-@cache.cache_netcdf
+@cache.cache_netcdf(coords_2d=True)
def get_ahn(ds=None, identifier="AHN4_DTM_5m", method="average", extent=None):
"""Get a model dataset with ahn variable.
@@ -86,6 +91,33 @@ 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.
+
+ Parameters
+ ----------
+ x : float
+ The x-coordinate fo the point.
+ y : float
+ The y-coordinate fo the point..
+ buffer : float, optional
+ The buffer around x and y that is downloaded. The default is 0.75.
+ return_da : bool, optional
+ Return the downloaded DataArray when True. The default is False.
+ return_mean : bool, optional
+ Resturn the mean of all non-nan pixels within buffer. Return the center pixel
+ when False. The default is False.
+ identifier : str, optional
+ The identifier passed onto get_latest_ahn_from_wcs. The default is "dsm_05m".
+ res : float, optional
+ The resolution that is passed onto get_latest_ahn_from_wcs. The default is 0.5.
+ **kwargs : dict
+ kwargs are passed onto the method get_latest_ahn_from_wcs.
+
+ Returns
+ -------
+ 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)
if return_da:
@@ -99,7 +131,66 @@ def get_ahn_at_point(
return ahn.data[int((ahn.shape[0] - 1) / 2), int((ahn.shape[1] - 1) / 2)]
-@cache.cache_netcdf
+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.
+
+ Parameters
+ ----------
+ line : shapely.LineString
+ The line along which the surface level is calculated.
+ ahn : xr.DataArray, optional
+ The 2d DataArray containing surface level values. If None, ahn4-values are
+ downloaded from the web. The default is None.
+ dx : float, optional
+ The distance between the points along the line at which the surface level is
+ calculated. Only used when num is None. When dx is None, it is set to the
+ resolution of ahn. The default is None.
+ num : int, optional
+ If not None, the surface level is calculated at num equally spaced points along
+ the line. The default is None.
+ method : string, optional
+ The method to interpolate the 2d surface level values to the points along the
+ line. The default is "linear".
+ plot : bool, optional
+ if True, plot the 2d surface level, the line and the calculated heights. The
+ default is False.
+
+ Returns
+ -------
+ z : xr.DataArray
+ A DataArray with dimension s, containing surface level values along the line.
+ """
+ if ahn is None:
+ bbox = line.bounds
+ extent = [bbox[0], bbox[2], bbox[1], bbox[3]]
+ ahn = get_ahn4(extent)
+ if num is not None:
+ s = np.linspace(0.0, line.length, num)
+ else:
+ if dx is None:
+ dx = float(ahn.x[1] - ahn.x[0])
+ s = np.arange(0.0, line.length, dx)
+
+ x, y = zip(*[p.xy for p in line.interpolate(s)])
+
+ x = np.array(x)[:, 0]
+ y = np.array(y)[:, 0]
+
+ x = xr.DataArray(x, dims="s", coords={"s": s})
+ y = xr.DataArray(y, dims="s", coords={"s": s})
+ z = ahn.interp(x=x, y=y, method=method)
+
+ if plot:
+ _, ax = plt.subplots(figsize=(10, 10))
+ ahn.plot(ax=ax)
+ gpd.GeoDataFrame(geometry=[line]).plot(ax=ax)
+
+ _, ax = plt.subplots(figsize=(10, 10))
+ z.plot(ax=ax)
+ return z
+
+
+@cache.cache_netcdf()
def get_latest_ahn_from_wcs(
extent=None,
identifier="dsm_05m",
@@ -142,7 +233,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):
@@ -181,10 +271,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
@@ -205,8 +294,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)
@@ -215,7 +303,7 @@ def get_ahn4_tiles(extent=None):
return gdf
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_ahn1(extent, identifier="ahn1_5m", as_data_array=True):
"""Download AHN1.
@@ -242,7 +330,7 @@ def get_ahn1(extent, identifier="ahn1_5m", as_data_array=True):
return da
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_ahn2(extent, identifier="ahn2_5m", as_data_array=True):
"""Download AHN2.
@@ -266,7 +354,7 @@ def get_ahn2(extent, identifier="ahn2_5m", as_data_array=True):
return _download_and_combine_tiles(tiles, identifier, extent, as_data_array)
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_ahn3(extent, identifier="AHN3_5m_DTM", as_data_array=True):
"""Download AHN3.
@@ -289,7 +377,7 @@ def get_ahn3(extent, identifier="AHN3_5m_DTM", as_data_array=True):
return _download_and_combine_tiles(tiles, identifier, extent, as_data_array)
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_ahn4(extent, identifier="AHN4_DTM_5m", as_data_array=True):
"""Download AHN4.
@@ -319,6 +407,11 @@ def _download_and_combine_tiles(tiles, identifier, extent, as_data_array):
datasets = []
for name in tqdm(tiles.index, desc=f"Downloading tiles of {identifier}"):
url = tiles.at[name, identifier]
+ if isinstance(url, pd.Series):
+ logger.warning(
+ f"Multiple tiles with the same name: {name}. Choosing the first one."
+ )
+ url = url.iloc[0]
path = url.split("/")[-1].replace(".zip", ".TIF")
if path.lower().endswith(".tif.tif"):
path = path[:-4]
diff --git a/nlmod/read/bgt.py b/nlmod/read/bgt.py
index db7b2ef9..6ac8d2cd 100644
--- a/nlmod/read/bgt.py
+++ b/nlmod/read/bgt.py
@@ -10,7 +10,7 @@
import requests
from shapely.geometry import LineString, MultiPolygon, Point, Polygon
-from ..dims.resample import extent_to_polygon
+from ..util import extent_to_polygon
def get_bgt(
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 5ded82f0..9457e53a 100644
--- a/nlmod/read/geotop.py
+++ b/nlmod/read/geotop.py
@@ -13,7 +13,7 @@
logger = logging.getLogger(__name__)
-GEOTOP_URL = r"http://www.dinodata.nl/opendap/GeoTOP/geotop.nc"
+GEOTOP_URL = "https://dinodata.nl/opendap/GeoTOP/geotop.nc"
def get_lithok_props(rgb_colors=True):
@@ -59,7 +59,7 @@ def get_kh_kv_table(kind="Brabant"):
return df
-@cache.cache_netcdf
+@cache.cache_netcdf()
def to_model_layers(
geotop_ds,
strat_props=None,
@@ -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:
@@ -233,11 +232,11 @@ def to_model_layers(
return ds
-@cache.cache_netcdf
+@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
----------
@@ -591,7 +589,7 @@ def aggregate_to_ds(
if "layer" in top.dims:
top = top[0].drop_vars("layer")
else:
- if "layer" in ds["top"].dims:
+ if "layer" in ds["top"].dims:
top = ds["top"][ilay].drop_vars("layer")
else:
top = ds["botm"][ilay - 1].drop_vars("layer")
diff --git a/nlmod/read/jarkus.py b/nlmod/read/jarkus.py
index e0672a5f..257fe7ca 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
@@ -17,15 +18,16 @@
import xarray as xr
from .. import cache
-from ..dims.resample import fillnan_da, get_extent, structured_da_to_ds
+from ..dims.grid import get_extent
+from ..dims.resample import fillnan_da, structured_da_to_ds
from ..util import get_da_from_da_ds, get_ds_empty
logger = logging.getLogger(__name__)
-@cache.cache_netcdf
+@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
----------
@@ -75,7 +77,7 @@ def get_bathymetry(ds, northsea, kind="jarkus", method="average"):
# fill nan values in bathymetry
da_bathymetry_filled = fillnan_da(da_bathymetry_raw)
- # bathymetrie mag nooit groter zijn dan NAP 0.0
+ # bathymetry can never be larger than NAP 0.0
da_bathymetry_filled = xr.where(da_bathymetry_filled > 0, 0, da_bathymetry_filled)
# bathymetry projected on model grid
@@ -92,7 +94,7 @@ def get_bathymetry(ds, northsea, kind="jarkus", method="average"):
return ds_out
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_dataset_jarkus(extent, kind="jarkus", return_tiles=False, time=-1):
"""Get bathymetry from Jarkus within a certain extent. If return_tiles is False, the
following actions are performed:
@@ -126,7 +128,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)
@@ -272,7 +273,7 @@ def add_bathymetry_to_top_bot_kh_kv(ds, bathymetry, fill_mask, kh_sea=10, kv_sea
ds["kv"][lay] = xr.where(fill_mask, kv_sea, ds["kv"][lay])
- # reset bot for all layers based on bathymetrie
+ # reset bot for all layers based on bathymetry
for lay in range(1, ds.sizes["layer"]):
ds["botm"][lay] = np.where(
ds["botm"][lay] > ds["botm"][lay - 1],
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index 512f34a1..9344324b 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -7,15 +7,15 @@
from hydropandas.io import knmi as hpd_knmi
from .. import cache, util
+from ..dims.grid import get_affine_mod_to_world
from ..dims.layers import get_first_active_layer
-from ..dims.resample import get_affine_mod_to_world
logger = logging.getLogger(__name__)
-@cache.cache_netcdf
+@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
----------
@@ -277,20 +275,38 @@ def get_knmi_at_locations(ds, start="2010", end=None, most_common_station=False)
stns_ev24 = locations["stn_ev24"].unique()
# get knmi data stations closest to any grid cell
- oc_knmi_prec = hpd.ObsCollection.from_knmi(
- stns=stns_rd,
- starts=[start],
- ends=[end],
- meteo_vars=["RD"],
- fill_missing_obs=True,
- )
+ olist_rd, new_stns_rd = [], []
+ for stnrd in stns_rd:
+ o = hpd.PrecipitationObs.from_knmi(
+ meteo_var="RD", stn=stnrd, start=start, end=end, fill_missing_obs=True
+ )
- oc_knmi_evap = hpd.ObsCollection.from_knmi(
- stns=stns_ev24,
- starts=[start],
- ends=[end],
- meteo_vars=["EV24"],
- fill_missing_obs=True,
- )
+ # if a station has no data in the given period another station is selected
+ if o.station != stnrd:
+ locations["stn_rd"] = locations["stn_rd"].replace(stnrd, o.station)
+
+ # only add the station if it does not exist yet
+ if o.station not in new_stns_rd:
+ olist_rd.append(o)
+ new_stns_rd.append(o.station)
+
+ oc_knmi_prec = hpd.ObsCollection(olist_rd)
+
+ olist_ev24, new_stns_ev24 = [], []
+ for stnev24 in stns_ev24:
+ o = hpd.EvaporationObs.from_knmi(
+ meteo_var="EV24", stn=stnev24, start=start, end=end, fill_missing_obs=True
+ )
+
+ # if a station has no data in the given period another station is selected
+ if o.station != stnev24:
+ locations["stn_ev24"] = locations["stn_rd"].replace(stnev24, o.station)
+
+ # only add the station if it does not exist yet
+ if o.station not in new_stns_ev24:
+ olist_ev24.append(o)
+ new_stns_ev24.append(o.station)
+
+ oc_knmi_evap = hpd.ObsCollection(olist_ev24)
return locations, oc_knmi_prec, oc_knmi_evap
diff --git a/nlmod/read/knmi_data_platform.py b/nlmod/read/knmi_data_platform.py
index f6efadf9..a2b57413 100644
--- a/nlmod/read/knmi_data_platform.py
+++ b/nlmod/read/knmi_data_platform.py
@@ -16,30 +16,31 @@
logger = logging.getLogger(__name__)
# base_url = "https://api.dataplatform.knmi.nl/dataset-content/v1/datasets"
-base_url = "https://api.dataplatform.knmi.nl/open-data"
+base_url = "https://api.dataplatform.knmi.nl/open-data/v1"
def get_anonymous_api_key() -> Union[str, None]:
try:
- url = "https://developer.dataplatform.knmi.nl/get-started"
- tables = read_html(url) # get all tables from url
- for table in tables:
- for coln in table.columns:
- if "KEY" in coln.upper(): # look for columns with key
- api_key_str = table.iloc[0].loc[
- coln
- ] # get entry with key (first row)
- api_key = max(
- api_key_str.split(), key=len
- ) # get key base on str length
- logger.info(f"Retrieved anonymous API Key from {url}")
- return api_key
+ url = "https://developer.dataplatform.knmi.nl/open-data-api#token"
+ webpage = requests.get(url, timeout=120) # get webpage
+ api_key = (
+ webpage.text.split("")[0].split("")[-1].strip()
+ ) # obtain apikey from codeblock on webpage
+ if len(api_key) != 120:
+ msg = f"Could not obtain API Key from {url}, trying API Key from memory. Found API Key = {api_key}"
+ logger.error(msg)
+ raise ValueError(msg)
+ logger.info(f"Retrieved anonymous API Key from {url}")
+ return api_key
except Exception as exc:
- if Timestamp.today() < Timestamp("2024-07-01"):
- logger.info("Retrieved anonymous API Key from memory")
+ api_key_memory_date = "2025-07-01"
+ if Timestamp.today() < Timestamp(api_key_memory_date):
+ logger.info(
+ f"Retrieved anonymous API Key (available till {api_key_memory_date}) from memory"
+ )
api_key = (
- "eyJvcmciOiI1ZTU1NGUxOTI3NGE5NjAwMDEyYTNlYjEiLCJpZCI6ImE1OGI5"
- "NGZmMDY5NDRhZDNhZjFkMDBmNDBmNTQyNjBkIiwiaCI6Im11cm11cjEyOCJ9"
+ "eyJvcmciOiI1ZTU1NGUxOTI3NGE5NjAwMDEyYTNlYjEiLCJpZCI6ImE1OGI5N"
+ "GZmMDY5NDRhZDNhZjFkMDBmNDBmNTQyNjBkIiwiaCI6Im11cm11cjEyOCJ9"
)
return api_key
else:
@@ -58,7 +59,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 = []
@@ -69,6 +70,7 @@ def get_list_of_files(
params = {"maxKeys": f"{max_keys}"}
if start_after_filename is not None:
params["startAfterFilename"] = start_after_filename
+ logger.debug(f"Request to {url=} with {params=}")
r = requests.get(
url, params=params, headers={"Authorization": api_key}, timeout=timeout
)
@@ -88,7 +90,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 = (
@@ -118,7 +120,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,
@@ -131,7 +133,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)
@@ -160,7 +162,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)):
@@ -173,7 +175,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)
@@ -185,7 +187,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
@@ -205,7 +207,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:
@@ -230,7 +232,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 = {}
@@ -247,7 +249,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("/")])
@@ -275,7 +277,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 858e3a16..4e721d60 100644
--- a/nlmod/read/nhi.py
+++ b/nlmod/read/nhi.py
@@ -1,7 +1,10 @@
+import io
import logging
import os
+import geopandas as gpd
import numpy as np
+import pandas as pd
import requests
import rioxarray
@@ -11,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
----------
@@ -33,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]
@@ -47,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
----------
@@ -63,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"
@@ -86,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
@@ -125,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
@@ -173,3 +170,208 @@ def add_buisdrainage(
ds[depth_var] = ds[depth_var] / 100.0
return ds
+
+
+def get_gwo_wells(
+ username,
+ password,
+ n_well_filters=1_000,
+ well_site=None,
+ organisation=None,
+ status=None,
+ well_index="Name",
+ timeout=120,
+ **kwargs,
+):
+ """Get metadata of extraction wells from the NHI GWO database.
+
+ Parameters
+ ----------
+ username : str
+ The username of the NHI GWO database. To retrieve a username and password visit
+ https://gwo.nhi.nu/register/.
+ password : str
+ The password of the NHI GWO database. To retrieve a username and password visit
+ https://gwo.nhi.nu/register/.
+ n_well_filters : int, optional
+ The number of wells that are requested per page. This number determines in how
+ many pieces the request is split. The default is 1000.
+ organisation : str, optional
+ The organisation that manages the wells. If not None, the organisation will be
+ used to filter the wells. The default is None.
+ well_site : str, optional
+ The name of well site the wells belong to. If not None, the well site will be
+ used to filter the wells. The default is None.
+ status : str, optional
+ The status of the wells. If not None, the status will be used to filter the
+ wells. Possible values are "Active", "Inactive" or "Abandoned". The default is
+ None.
+ well_index : str, tuple or list, optional
+ The column(s) in the resulting GeoDataFrame that is/are used as the index of
+ this GeoDataFrame. The default is "Name".
+ timeout : int, optional
+ The timeout time (in seconds) for requests to the database. The default is
+ 120 seconds.
+ **kwargs : dict
+ Kwargs are passed as additional parameters in the request to the database. For
+ available parameters see https://gwo.nhi.nu/api/v1/download/.
+
+ Returns
+ -------
+ 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/"
+
+ page = 1
+ properties = []
+ while page is not None:
+ params = {"format": "csv", "n_well_filters": n_well_filters, "page": page}
+ if status is not None:
+ params["well__status"] = status
+ if organisation is not None:
+ params["well__organization"] = organisation
+ if well_site is not None:
+ params["well__site"] = well_site
+ params.update(kwargs)
+
+ r = requests.get(url, auth=(username, password), params=params, timeout=timeout)
+ content = r.content.decode("utf-8")
+ if len(content) == 0:
+ if page == 1:
+ msg = "No extraction wells found for the requested parameters"
+ raise ValueError(msg)
+ else:
+ # the number of wells is exactly a multiple of n_well_filters
+ page = None
+ continue
+ lines = content.split("\n")
+ empty_lines = np.where([set(line) == set(";") for line in lines])[0]
+ assert len(empty_lines) == 1, "Returned extraction wells cannot be interpreted"
+ skiprows = list(range(empty_lines[0] + 1)) + [empty_lines[0] + 2]
+ df = pd.read_csv(io.StringIO(content), skiprows=skiprows, sep=";")
+ properties.append(df)
+
+ if len(df) == n_well_filters:
+ page += 1
+ else:
+ page = None
+ df = pd.concat(properties)
+ geometry = gpd.points_from_xy(df.XCoordinate, df.YCoordinate)
+ gdf = gpd.GeoDataFrame(df, geometry=geometry, crs=28992)
+ if well_index is not None:
+ gdf = gdf.set_index(well_index)
+ return gdf
+
+
+def get_gwo_measurements(
+ username,
+ password,
+ n_measurements=10_000,
+ well_site=None,
+ well_index="Name",
+ measurement_index=("Name", "DateTime"),
+ timeout=120,
+ **kwargs,
+):
+ """Get extraction rates and metadata of wells from the NHI GWO database.
+
+ Parameters
+ ----------
+ username : str
+ The username of the NHI GWO database. To retrieve a username and password visit
+ https://gwo.nhi.nu/register/.
+ password : str
+ The password of the NHI GWO database. To retrieve a username and password visit
+ https://gwo.nhi.nu/register/.
+ n_measurements : int, optional
+ The number of measurements that are requested per page, with a maximum of
+ 200,000. This number determines in how many pieces the request is split. The
+ default is 10,000.
+ well_site : str, optional
+ The name of well site the wells belong to. If not None, the well site will be
+ used to filter the wells. The default is None.
+ well_index : str, tuple or list, optional
+ The column(s) in the resulting GeoDataFrame that is/are used as the index of
+ this GeoDataFrame. The default is "Name".
+ measurement_index : str, tuple or list, optional, optional
+ The column(s) in the resulting measurement-DataFrame that is/are used as the
+ index of this DataFrame. The default is ("Name", "DateTime").
+ timeout : int, optional
+ The timeout time (in seconds) of requests to the database. The default is
+ 120 seconds.
+ **kwargs : dict
+ Kwargs are passed as additional parameters in the request to the database. For
+ available parameters see https://gwo.nhi.nu/api/v1/download/.
+
+ Returns
+ -------
+ measurements : pandas.DataFrame
+ 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 = []
+ measurements = []
+ page = 1
+ while page is not None:
+ params = {
+ "format": "csv",
+ "n_measurements": n_measurements,
+ "page": page,
+ }
+ if well_site is not None:
+ params["filter__well__site"] = well_site
+ params.update(kwargs)
+ r = requests.get(url, auth=(username, password), params=params, timeout=timeout)
+
+ content = r.content.decode("utf-8")
+ if len(content) == 0:
+ if page == 1:
+ msg = "No extraction rates found for the requested parameters"
+ raise (ValueError(msg))
+ else:
+ # the number of measurements is exactly a multiple of n_measurements
+ page = None
+ continue
+ lines = content.split("\n")
+ empty_lines = np.where([set(line) == set(";") for line in lines])[0]
+ assert len(empty_lines) == 2, "Returned extraction rates cannot be interpreted"
+
+ # read properties
+ skiprows = list(range(empty_lines[0] + 1)) + [empty_lines[0] + 2]
+ nrows = empty_lines[1] - empty_lines[0] - 3
+ df = pd.read_csv(io.StringIO(content), sep=";", skiprows=skiprows, nrows=nrows)
+ properties.append(df)
+
+ # read measurements
+ skiprows = list(range(empty_lines[1] + 1)) + [empty_lines[1] + 2]
+ df = pd.read_csv(
+ io.StringIO(content),
+ skiprows=skiprows,
+ sep=";",
+ parse_dates=["DateTime"],
+ dayfirst=True,
+ )
+ measurements.append(df)
+ if len(df) == n_measurements:
+ page += 1
+ else:
+ page = None
+ measurements = pd.concat(measurements)
+ # drop columns without measurements
+ measurements = measurements.loc[:, ~measurements.isna().all()]
+ if measurement_index is not None:
+ if isinstance(measurement_index, tuple):
+ measurement_index = list(measurement_index)
+ measurements = measurements.set_index(["Name", "DateTime"])
+ df = pd.concat(properties)
+ geometry = gpd.points_from_xy(df.XCoordinate, df.YCoordinate)
+ gdf = gpd.GeoDataFrame(df, geometry=geometry)
+ if well_index is not None:
+ gdf = gdf.set_index(well_index)
+ # drop duplicate properties from multiple pages
+ gdf = gdf[~gdf.index.duplicated()]
+ return measurements, gdf
diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py
index 729d7b44..f1ba8683 100644
--- a/nlmod/read/regis.py
+++ b/nlmod/read/regis.py
@@ -12,11 +12,10 @@
logger = logging.getLogger(__name__)
-REGIS_URL = "http://www.dinodata.nl:80/opendap/REGIS/REGIS.nc"
-# REGIS_URL = 'https://www.dinodata.nl/opendap/hyrax/REGIS/REGIS.nc'
+REGIS_URL = "https://dinodata.nl/opendap/REGIS/REGIS.nc"
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_combined_layer_models(
extent,
regis_botm_layer="AKc",
@@ -26,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
@@ -65,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
@@ -93,7 +91,7 @@ def get_combined_layer_models(
return combined_ds
-@cache.cache_netcdf
+@cache.cache_netcdf()
def get_regis(
extent,
botm_layer="AKc",
@@ -102,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
----------
@@ -133,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
@@ -196,10 +193,16 @@ def get_regis(
def add_geotop_to_regis_layers(
- rg, gt, layers="HLc", geotop_k=None, remove_nan_layers=True, anisotropy=1.0
+ rg,
+ gt,
+ layers="HLc",
+ geotop_k=None,
+ remove_nan_layers=True,
+ 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
----------
@@ -217,6 +220,10 @@ def add_geotop_to_regis_layers(
anisotropy : float, optional
The anisotropy value (kh/kv) used when there are no kv values in df. The
default is 1.0.
+ gt_layered : xarray.Dataset
+ A layered representation of the geotop-dataset. By supplying this parameter, the
+ user can change the GeoTOP-layering, which is usueally defined by
+ nlmod.read.geotop.to_model_layers(gt).
Returns
-------
@@ -254,8 +261,11 @@ def add_geotop_to_regis_layers(
rg["top"] = rg["botm"] + calculate_thickness(rg)
for layer in layers:
- # transform geotop data into layers
- gtl = geotop.to_model_layers(gt)
+ if gt_layered is not None:
+ gtl = gt_layered.copy(deep=True)
+ else:
+ # transform geotop data into layers
+ gtl = geotop.to_model_layers(gt)
# temporarily add layer dimension to top in gtl
gtl["top"] = gtl["botm"] + calculate_thickness(gtl)
@@ -291,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 7af2a991..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
----------
@@ -37,9 +37,9 @@ def get_gdf_surface_water(ds):
return gdf_swater
-@cache.cache_netcdf
+@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)
@@ -91,10 +90,10 @@ def get_surface_water(ds, da_basename):
return ds_out
-@cache.cache_netcdf
+@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 b44d44c2..64eaa5a8 100644
--- a/nlmod/read/waterboard.py
+++ b/nlmod/read/waterboard.py
@@ -49,15 +49,16 @@ def get_configuration():
config["Amstel, Gooi en Vecht"] = {
"bgt_code": "W0155",
"watercourses": {
- "url": "https://maps.waternet.nl/arcgis/rest/services/AGV_Legger/AGV_Onderh_Secundaire_Watergangen/MapServer",
- "layer": 40,
- "bottom_width": "BODEMBREEDTE",
- "bottom_height": "BODEMHOOGTE",
- "water_depth": "WATERDIEPTE",
+ "url": "https://maps.waternet.nl/arcgis/rest/services/Publiek/WNET_GEO_LEGGER_WL_2021/MapServer",
+ "layer": 0, # Primaire Waterloop Legger
+ "bottom_width": "AVVBODDR",
+ "bottom_height": "AVVBODH",
+ "water_depth": "AVVDIEPT",
+ "index": "OVKIDENT",
},
"level_areas": {
- "url": "https://maps.waternet.nl/arcgis/rest/services/AGV_Legger/Vastgestelde_Waterpeilen/MapServer",
- "layer": 0,
+ "url": "https://maps.waternet.nl/arcgis/rest/services/Publiek/GW_GPG/MapServer",
+ "layer": 5, # Vigerende peilgebieden
"index": "GPGIDENT",
"summer_stage": [
"GPGZMRPL",
@@ -201,18 +202,6 @@ def get_configuration():
},
"level_areas": {
"url": "https://kaarten.hhnk.nl/arcgis/rest/services/ws/ws_peilgebieden_vigerend/MapServer",
- "layer": 4,
- "table": {
- "id": 6,
- "SOORTSTREEFPEIL": {
- 901: "STREEFPEIL_JAARROND", # vast peilbeheer
- 902: "STREEFPEIL_WINTER",
- 903: "STREEFPEIL_ZOMER",
- 904: "STREEFPEIL_JAARROND", # dynamisch peilbeheer
- 905: "ONDERGRENS_JAARROND",
- 906: "BOVENGRENS_JAARROND",
- },
- },
"summer_stage": [
"ZOMER",
"STREEFPEIL_ZOMER",
@@ -522,7 +511,6 @@ def get_data(wb, data_kind, extent=None, max_record_count=None, config=None, **k
Raises
------
-
DESCRIPTION.
Returns
@@ -605,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):
@@ -645,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 64c72742..97cc90bf 100644
--- a/nlmod/read/webservices.py
+++ b/nlmod/read/webservices.py
@@ -150,31 +150,40 @@ def arcrest(
else:
gdf = gpd.GeoDataFrame.from_features(features, crs=sr)
if table is not None:
- url_query = f"{url}/{table.pop('id')}/query"
- pgbids = ",".join([str(v) for v in gdf["OBJECTID"].values])
- params["where"] = f"PEILGEBIEDVIGERENDID IN ({pgbids})"
params["f"] = "json"
- data = _get_data(url_query, params, timeout=timeout)
+ url_query = f"{url}/{table.pop('id')}/query"
+
+ # loop over chunks of 100 pgbids. Long where clauses can cause
+ # the request to fail. 1300 pgbids fails but 130 works
+ chunk_size = 100
+ ids_chunks = [
+ gdf["OBJECTID"].values[i : i + chunk_size]
+ for i in range(0, len(gdf), chunk_size)
+ ]
+ data = {}
+ features = []
+
+ for ids_chunk in ids_chunks:
+ pgbids = ",".join([str(v) for v in ids_chunk])
+ where = f"PEILGEBIEDVIGERENDID IN ({pgbids})"
+ params["where"] = where
+ _data = _get_data(url_query, params, timeout=timeout, **kwargs)
+
+ data.update(_data)
+ features.extend(_data["features"])
+
+ assert "exceededTransferLimit" not in data, "exceededTransferLimit"
+ data["features"] = features
+
df = pd.DataFrame(
[feature["attributes"] for feature in data["features"]]
)
- # add peilen to gdf
- for col, convert_dic in table.items():
- df[col].replace(convert_dic, inplace=True)
- df.set_index(col, inplace=True)
- for oid in gdf["OBJECTID"]:
- insert_s = df.loc[
- df["PEILGEBIEDVIGERENDID"] == oid, "WATERHOOGTE"
- ]
- gdf.loc[
- gdf["OBJECTID"] == oid, insert_s.index
- ] = insert_s.values
return gdf
def _get_data(url, params, timeout=120, **kwargs):
- """get data using a request
+ """Get data using a request.
Parameters
----------
@@ -188,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:
@@ -423,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
----------
@@ -454,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 119a97bf..4d5e2438 100644
--- a/nlmod/sim/sim.py
+++ b/nlmod/sim/sim.py
@@ -1,6 +1,7 @@
import datetime as dt
import logging
import os
+import pathlib
from shutil import copyfile
import flopy
@@ -13,9 +14,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
----------
@@ -51,7 +52,10 @@ def write_and_run(sim, ds, write_ds=True, script_path=None, silent=False):
ds.attrs["model_dataset_written_to_disk_on"] = dt.datetime.now().strftime(
"%Y%m%d_%H:%M:%S"
)
- ds.to_netcdf(os.path.join(ds.attrs["model_ws"], f"{ds.model_name}.nc"))
+ if isinstance(ds.attrs["model_ws"], pathlib.PurePath):
+ ds.to_netcdf(ds.attrs["model_ws"] / f"{ds.model_name}.nc")
+ else:
+ ds.to_netcdf(os.path.join(ds.attrs["model_ws"], f"{ds.model_name}.nc"))
logger.info("write modflow files to model workspace")
sim.write_simulation(silent=silent)
@@ -107,8 +111,8 @@ def get_tdis_perioddata(ds, nstp="nstp", tsmult="tsmult"):
return tdis_perioddata
-def sim(ds, exe_name=None):
- """create sim from the model dataset.
+def sim(ds, exe_name=None, version_tag=None):
+ """Create sim from the model dataset.
Parameters
----------
@@ -117,21 +121,36 @@ def sim(ds, exe_name=None):
attributes: model_name, mfversion, model_ws, time_units, start,
perlen, nstp, tsmult
exe_name: str, optional
- path to modflow executable, default is None, which assumes binaries
- are available in nlmod/bin directory. Binaries can be downloaded
- using `nlmod.util.download_mfbinaries()`.
+ path to modflow executable, default is None. If None, the path is
+ obtained from the flopy metadata that respects `version_tag`. If not
+ found, the executables are downloaded. Not compatible with version_tag.
+ version_tag : str, default None
+ GitHub release ID: for example "18.0" or "latest". If version_tag is provided,
+ the most recent installation location of MODFLOW is found in flopy metadata
+ that respects `version_tag`. If not found, the executables are downloaded.
+ Not compatible with exe_name.
Returns
-------
sim : flopy MFSimulation
simulation object.
"""
-
# start creating model
logger.info("creating mf6 SIM")
- if exe_name is None:
- exe_name = util.get_exe_path(ds.mfversion)
+ # Most likely exe_name was previously set with to_model_ds()
+ if exe_name is not None:
+ exe_name = util.get_exe_path(exe_name=exe_name, version_tag=version_tag)
+ elif "exe_name" in ds.attrs:
+ exe_name = util.get_exe_path(
+ exe_name=ds.attrs["exe_name"], version_tag=version_tag
+ )
+ elif "mfversion" in ds.attrs:
+ exe_name = util.get_exe_path(
+ exe_name=ds.attrs["mfversion"], version_tag=version_tag
+ )
+ else:
+ raise ValueError("No exe_name provided and no exe_name found in ds.attrs")
# Create the Flopy simulation object
sim = flopy.mf6.MFSimulation(
@@ -145,7 +164,7 @@ def sim(ds, exe_name=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
----------
@@ -164,7 +183,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")
@@ -185,7 +203,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
----------
@@ -201,7 +219,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 a2c1c72b..cfa273e7 100644
--- a/nlmod/util.py
+++ b/nlmod/util.py
@@ -1,19 +1,24 @@
+import json
import logging
import os
import re
import sys
import warnings
+from pathlib import Path
from typing import Dict, Optional
-import flopy
import geopandas as gpd
import requests
import xarray as xr
from colorama import Back, Fore, Style
-from shapely.geometry import box
+from flopy.utils import get_modflow
+from flopy.utils.get_modflow import flopy_appdata_path, get_release
+from shapely.geometry import Polygon, box
logger = logging.getLogger(__name__)
+nlmod_bindir = Path(__file__).parent / "bin"
+
class LayerError(Exception):
"""Generic error when modifying layers."""
@@ -89,29 +94,349 @@ def get_model_dirs(model_ws):
return figdir, cachedir
-def get_exe_path(exe_name="mf6"):
- """Get the full path of the executable. Uses the bin directory in the nlmod package.
+def get_exe_path(
+ exe_name="mf6",
+ bindir=None,
+ download_if_not_found=True,
+ version_tag=None,
+ repo="executables",
+):
+ """Get the full path of the executable.
+
+ 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.
+ 1. The directory specified with `bindir`. Raises error if exe_name is provided
+ and not found.
+ 2. The directory used by nlmod installed in this environment.
+ 3. If the executables were downloaded with flopy/nlmod from an other env,
+ most recent installation location of MODFLOW is found in flopy metadata
+
+ Else:
+ 4. Download the executables using `version_tag` and `repo`.
+
+ The returned directory is checked to contain exe_name if it is provided.
Parameters
----------
exe_name : str, optional
- name of the executable. The default is 'mf6'.
+ The name of the executable, by default "mf6".
+ bindir : Path, optional
+ The directory where the executables are stored, by default None
+ 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.
+ 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.
Returns
-------
- exe_path : str
+ exe_full_path : str
full path of the executable.
"""
- exe_path = os.path.join(os.path.dirname(__file__), "bin", exe_name)
- if sys.platform.startswith("win"):
- exe_path += ".exe"
+ if sys.platform.startswith("win") and not exe_name.endswith(".exe"):
+ exe_name += ".exe"
- if not os.path.exists(exe_path):
- logger.warning(
- f"executable {exe_path} not found, download the binaries using nlmod.util.download_mfbinaries"
+ # If exe_name is a full path
+ if Path(exe_name).exists():
+ enable_version_check = version_tag is not None and repo is not None
+
+ if enable_version_check:
+ msg = (
+ "Incompatible arguments. If exe_name is provided, unable to check "
+ "the version."
+ )
+ raise ValueError(msg)
+ exe_full_path = exe_name
+
+ else:
+ exe_full_path = str(
+ get_bin_directory(
+ exe_name=exe_name,
+ bindir=bindir,
+ download_if_not_found=download_if_not_found,
+ version_tag=version_tag,
+ repo=repo,
+ )
+ / exe_name
)
- return exe_path
+ msg = f"Executable path: {exe_full_path}"
+ logger.debug(msg)
+
+ return exe_full_path
+
+
+def get_bin_directory(
+ exe_name="mf6",
+ bindir=None,
+ download_if_not_found=True,
+ version_tag=None,
+ repo="executables",
+) -> Path:
+ """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.
+ 1. The directory specified with `bindir`. Raises error if exe_name is provided
+ and not found. Requires enable_version_check to be False.
+ 2. The directory used by nlmod installed in this environment.
+ 3. If the executables were downloaded with flopy/nlmod from an other env,
+ most recent installation location of MODFLOW is found in flopy metadata
+
+ 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.
+
+ Parameters
+ ----------
+ exe_name : str, optional
+ The name of the executable, by default mf6.
+ bindir : Path, optional
+ The directory where the executables are stored, by default "mf6".
+ 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.
+ 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.
+
+ Returns
+ -------
+ Path
+ The directory where the executables are stored.
+
+ Raises
+ ------
+ FileNotFoundError
+ If the executables are not found in the specified directories.
+ """
+ bindir = Path(bindir) if bindir is not None else None
+
+ if sys.platform.startswith("win") and not exe_name.endswith(".exe"):
+ exe_name += ".exe"
+
+ enable_version_check = version_tag is not None
+
+ # If exe_name is a full path
+ if Path(exe_name).exists():
+ if enable_version_check:
+ msg = (
+ "Incompatible arguments. If exe_name is provided, unable to check "
+ "the version."
+ )
+ raise ValueError(msg)
+ return Path(exe_name).parent
+
+ # 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."
+ )
+ raise ValueError(msg)
+
+ use_bindir = (
+ bindir is not None and exe_name is not None and (bindir / exe_name).exists()
+ )
+ use_bindir |= bindir is not None and exe_name is None and bindir.exists()
+
+ if use_bindir:
+ return bindir
+
+ # If the executables are in the flopy directory
+ flopy_bindirs = get_flopy_bin_directories(version_tag=version_tag, repo=repo)
+
+ if exe_name is not None:
+ flopy_bindirs = [
+ flopy_bindir
+ for flopy_bindir in flopy_bindirs
+ if Path(flopy_bindir / exe_name).exists()
+ ]
+ else:
+ flopy_bindirs = [
+ flopy_bindir
+ for flopy_bindir in flopy_bindirs
+ if Path(flopy_bindir).exists()
+ ]
+
+ if nlmod_bindir in flopy_bindirs:
+ return nlmod_bindir
+
+ if flopy_bindirs:
+ # Get most recent directory
+ return flopy_bindirs[-1]
+
+ # Else download the executables
+ if download_if_not_found:
+ download_mfbinaries(
+ bindir=bindir,
+ version_tag=version_tag if version_tag is not None else "latest",
+ repo=repo,
+ )
+
+ # Rerun this function
+ return get_bin_directory(
+ exe_name=exe_name,
+ bindir=bindir,
+ download_if_not_found=False,
+ version_tag=version_tag,
+ repo=repo,
+ )
+
+ else:
+ msg = (
+ f"Could not find {exe_name} in {bindir}, "
+ f"{nlmod_bindir} and {flopy_bindirs}."
+ )
+ raise FileNotFoundError(msg)
+
+
+def get_flopy_bin_directories(version_tag=None, repo="executables"):
+ """Get the directories where the executables are stored.
+
+ Obtain the bin directory installed with flopy. If enable_version_check is True,
+ all installation location of MODFLOW are found in flopy metadata that respects
+ `version_tag` and `repo`.
+
+ Parameters
+ ----------
+ 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.
+ 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.
+
+ Returns
+ -------
+ list
+ list of directories where the executables are stored.
+ """
+ flopy_metadata_fp = flopy_appdata_path / "get_modflow.json"
+
+ if not flopy_metadata_fp.exists():
+ return []
+
+ meta_raw = flopy_metadata_fp.read_text()
+
+ # Remove trailing characters that are not part of the JSON.
+ while meta_raw[-3:] != "}\n]":
+ meta_raw = meta_raw[:-1]
+
+ # Get metadata of all flopy installations
+ meta_list = json.loads(meta_raw)
+
+ enable_version_check = version_tag is not None and repo is not None
+
+ if enable_version_check:
+ msg = (
+ "The version of the executables will be checked, because the "
+ f"`version_tag={version_tag}` is passed to `get_flopy_bin_directories()`."
+ )
+
+ # To convert latest into an explicit tag
+ if version_tag == "latest":
+ version_tag_pin = get_release(tag=version_tag, repo=repo, quiet=True)[
+ "tag_name"
+ ]
+ else:
+ version_tag_pin = version_tag
+
+ # get path to the most recent installation. Appended to end of get_modflow.json
+ meta_list_validversion = [
+ meta
+ for meta in meta_list
+ if (meta["release_id"] == version_tag_pin) and (meta["repo"] == repo)
+ ]
+
+ else:
+ msg = (
+ "The version of the executables will not be checked, because the "
+ "`version_tag` is not passed to `get_flopy_bin_directories()`."
+ )
+ meta_list_validversion = meta_list
+ logger.debug(msg)
+
+ path_list = [
+ Path(meta["bindir"])
+ for meta in meta_list_validversion
+ if Path(meta["bindir"]).exists()
+ ]
+ return path_list
+
+
+def download_mfbinaries(bindir=None, version_tag="latest", repo="executables"):
+ """Download and unpack platform-specific modflow binaries.
+
+ Source: USGS
+
+ Parameters
+ ----------
+ binpath : str, optional
+ path to directory to download binaries to, if it doesnt exist it
+ is created. Default is None which sets dir to nlmod/bin.
+ repo : str, default "executables"
+ Name of GitHub repository. Choose one of "executables" (default),
+ "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
+ bindir = nlmod_bindir
+
+ if not os.path.isdir(bindir):
+ os.makedirs(bindir)
+
+ get_modflow(bindir=str(bindir), release_id=version_tag, repo=repo)
+
+ # Ensure metadata is saved.
+ # https://github.com/modflowpy/flopy/blob/
+ # 0748dcb9e4641b5ad9616af115dd3be906f98f50/flopy/utils/get_modflow.py#L623
+ flopy_metadata_fp = flopy_appdata_path / "get_modflow.json"
+
+ if not flopy_metadata_fp.exists():
+ if "pytest" not in str(bindir) and "pytest" not in sys.modules:
+ logger.warning(
+ f"flopy metadata file not found at {flopy_metadata_fp}. "
+ "After downloading and installing the executables. "
+ "Creating a new metadata file."
+ )
+
+ release_metadata = get_release(tag=version_tag, repo=repo, quiet=True)
+ install_metadata = {
+ "release_id": release_metadata["tag_name"],
+ "repo": repo,
+ "bindir": str(bindir),
+ }
+
+ with open(flopy_metadata_fp, "w", encoding="UTF-8") as f:
+ json.dump([install_metadata], f, indent=4)
+
+ # download the provisional version of modpath from Github
+ download_modpath_provisional_exe(bindir=bindir, timeout=120)
def get_ds_empty(ds, keep_coords=None):
@@ -175,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
----------
@@ -192,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):
@@ -229,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]
@@ -267,6 +592,49 @@ def compare_model_extents(extent1, extent2):
raise NotImplementedError("other options are not yet implemented")
+def extent_to_polygon(extent):
+ """Generate a shapely Polygon from an extent ([xmin, xmax, ymin, ymax])
+
+ Parameters
+ ----------
+ extent : tuple, list or array
+ extent (xmin, xmax, ymin, ymax).
+
+ Returns
+ -------
+ shapely.geometry.Polygon
+ polygon of the extent.
+
+ """
+ nw = (extent[0], extent[2])
+ no = (extent[1], extent[2])
+ zo = (extent[1], extent[3])
+ zw = (extent[0], extent[3])
+ return Polygon([nw, no, zo, zw])
+
+
+def extent_to_gdf(extent, crs="EPSG:28992"):
+ """Create a geodataframe with a single polygon with the extent given.
+
+ Parameters
+ ----------
+ extent : tuple, list or array
+ extent.
+ crs : str, optional
+ coördinate reference system of the extent, default is EPSG:28992
+ (RD new)
+
+ Returns
+ -------
+ gdf_extent : geopandas.GeoDataFrame
+ geodataframe with extent.
+ """
+ geom_extent = extent_to_polygon(extent)
+ gdf_extent = gpd.GeoDataFrame(geometry=[geom_extent], crs=crs)
+
+ return gdf_extent
+
+
def polygon_from_extent(extent):
"""Create a shapely polygon from a given extent.
@@ -280,7 +648,10 @@ def polygon_from_extent(extent):
polygon_ext : shapely.geometry.polygon.Polygon
polygon of the extent.
"""
-
+ logger.warning(
+ "nlmod.util.polygon_from_extent is deprecated. "
+ "Use nlmod.util.extent_to_polygon instead"
+ )
bbox = (extent[0], extent[2], extent[1], extent[3])
polygon_ext = box(*tuple(bbox))
@@ -303,7 +674,10 @@ def gdf_from_extent(extent, crs="EPSG:28992"):
gdf_extent : GeoDataFrame
geodataframe with extent.
"""
-
+ logger.warning(
+ "nlmod.util.gdf_from_extent is deprecated. "
+ "Use nlmod.util.extent_to_gdf instead"
+ )
geom_extent = polygon_from_extent(extent)
gdf_extent = gpd.GeoDataFrame(geometry=[geom_extent], crs=crs)
@@ -311,8 +685,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
----------
@@ -327,7 +702,7 @@ def gdf_within_extent(gdf, extent):
dataframe with only polygon features within the extent.
"""
# create geodataframe from the extent
- gdf_extent = gdf_from_extent(extent, crs=gdf.crs)
+ gdf_extent = extent_to_gdf(extent, crs=gdf.crs)
# check type
geom_types = gdf.geom_type.unique()
@@ -366,6 +741,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):
@@ -392,6 +768,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):
@@ -430,30 +807,6 @@ def save_response_content(response, destination):
save_response_content(response, destination)
-def download_mfbinaries(bindir=None):
- """Download and unpack platform-specific modflow binaries.
-
- Source: USGS
-
- Parameters
- ----------
- binpath : str, optional
- path to directory to download binaries to, if it doesnt exist it
- is created. Default is None which sets dir to nlmod/bin.
- version : str, optional
- version string, by default 8.0
- """
-
- if bindir is None:
- bindir = os.path.join(os.path.dirname(__file__), "bin")
- if not os.path.isdir(bindir):
- os.makedirs(bindir)
- flopy.utils.get_modflow(bindir)
- if sys.platform.startswith("win"):
- # download the provisional version of modpath from Github
- download_modpath_provisional_exe(bindir)
-
-
def download_modpath_provisional_exe(bindir=None, timeout=120):
"""Download the provisional version of modpath to the folder with binaries."""
if bindir is None:
@@ -507,14 +860,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
@@ -522,6 +873,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 b6f9c6c1..213fe8f0 100644
--- a/nlmod/version.py
+++ b/nlmod/version.py
@@ -1,20 +1,17 @@
from importlib import metadata
from platform import python_version
-__version__ = "0.7.2"
+__version__ = "0.8.0"
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 6729fa58..4867067e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,7 +35,8 @@ dependencies = [
"matplotlib",
"dask",
"colorama",
- "joblib"
+ "joblib",
+ "bottleneck",
]
keywords = ["hydrology", "groundwater", "modeling", "Modflow 6", "flopy"]
classifiers = [
@@ -56,20 +57,13 @@ repository = "https://github.com/gwmod/nlmod"
documentation = "https://nlmod.readthedocs.io/en/latest/"
[project.optional-dependencies]
-full = [
- "nlmod[knmi]",
- "gdown",
- "geocube",
- "bottleneck",
- "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"]
nbtest = ["nbformat", "nbconvert>6.4.5"]
lint = ["flake8", "isort", "black[jupyter]"]
-ci = ["nlmod[full,lint,test,nbtest]", "netCDF4>=1.6.3", "pandas<2.1.0"]
+ci = ["nlmod[full,lint,test,nbtest]", "netCDF4<1.7.0", "pandas<2.1.0"]
rtd = [
"nlmod[full]",
"ipython",
@@ -78,7 +72,7 @@ rtd = [
"nbsphinx",
"sphinx_rtd_theme==1.0.0",
"nbconvert==7.13.0",
- "netCDF4>=1.6.3",
+ "netCDF4<1.7.0",
]
[tool.setuptools.dynamic]
@@ -92,7 +86,7 @@ include-package-data = true
[tool.setuptools.package-data]
"nlmod.data" = ["*.gleg"]
-"nlmod.data.geotop" = ["*.csv"]
+"nlmod.data.geotop" = ["*"]
"nlmod.data.shapes" = ["*"]
"nlmod.bin" = ["mp7_2_002_provisional"]
@@ -102,6 +96,33 @@ line-length = 88
[tool.isort]
profile = "black"
+[tool.ruff]
+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_004_northsea.py b/tests/test_004_northsea.py
index 2e49171d..d71ccea4 100644
--- a/tests/test_004_northsea.py
+++ b/tests/test_004_northsea.py
@@ -60,7 +60,7 @@ def test_get_bathymetry_seamodel():
assert (~ds_bathymetry.bathymetry.isnull()).sum() > 0
-def test_get_bathymetrie_nosea():
+def test_get_bathymetry_nosea():
# model without sea
ds = test_001_model.get_ds_from_cache("small_model")
ds.update(nlmod.read.rws.get_northsea(ds))
@@ -69,7 +69,7 @@ def test_get_bathymetrie_nosea():
assert (~ds_bathymetry.bathymetry.isnull()).sum() == 0
-def test_add_bathymetrie_to_top_bot_kh_kv_seamodel():
+def test_add_bathymetry_to_top_bot_kh_kv_seamodel():
# model with sea
ds = test_001_model.get_ds_from_cache("basic_sea_model")
ds.update(nlmod.read.rws.get_northsea(ds))
diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py
index f558b796..1d40d0d1 100644
--- a/tests/test_005_external_data.py
+++ b/tests/test_005_external_data.py
@@ -1,6 +1,8 @@
import pandas as pd
import pytest
import test_001_model
+import xarray as xr
+from shapely.geometry import LineString
import nlmod
@@ -65,9 +67,13 @@ def test_get_ahn3():
def test_get_ahn4():
extent = [98000.0, 100000.0, 494000.0, 496000.0]
- da = nlmod.read.ahn.get_ahn4(extent)
+ ahn = nlmod.read.ahn.get_ahn4(extent)
+ assert isinstance(ahn, xr.DataArray)
+ assert not ahn.isnull().all(), "AHN only has nan values"
- assert not da.isnull().all(), "AHN only has nan values"
+ line = LineString([(99000, 495000), (100000, 496000)])
+ ahn_line = nlmod.read.ahn.get_ahn_along_line(line, ahn=ahn)
+ assert isinstance(ahn_line, xr.DataArray)
def test_get_ahn():
@@ -80,6 +86,10 @@ def test_get_ahn():
assert not ahn_ds["ahn"].isnull().all(), "AHN only has nan values"
+def test_get_ahn_at_point():
+ nlmod.read.ahn.get_ahn_at_point(100010, 400010)
+
+
def test_get_surface_water_ghb():
# model with sea
ds = test_001_model.get_ds_from_cache("basic_sea_model")
@@ -88,13 +98,13 @@ def test_get_surface_water_ghb():
sim = nlmod.sim.sim(ds)
# create time discretisation
- tdis = nlmod.sim.tdis(ds, sim)
+ nlmod.sim.tdis(ds, sim)
# create groundwater flow model
gwf = nlmod.gwf.gwf(ds, sim)
# create ims
- ims = nlmod.sim.ims(sim)
+ nlmod.sim.ims(sim)
nlmod.gwf.dis(ds, gwf)
diff --git a/tests/test_006_caching.py b/tests/test_006_caching.py
index 741c1ffd..5bdfb3e0 100644
--- a/tests/test_006_caching.py
+++ b/tests/test_006_caching.py
@@ -1,96 +1,86 @@
+import os
import tempfile
-import pytest
-import test_001_model
-
import nlmod
-tmpdir = tempfile.gettempdir()
-
-
-def test_ds_check_true():
- # two models with the same grid and time dicretisation
- ds = test_001_model.get_ds_from_cache("small_model")
- ds2 = ds.copy()
-
- check = nlmod.cache._check_ds(ds, ds2)
-
- assert check
-
-
-def test_ds_check_time_false():
- # two models with a different time discretisation
- ds = test_001_model.get_ds_from_cache("small_model")
- ds2 = test_001_model.get_ds_time_steady(tmpdir)
-
- check = nlmod.cache._check_ds(ds, ds2)
-
- assert not check
-
-
-def test_ds_check_time_attributes_false():
- # two models with a different time discretisation
- ds = test_001_model.get_ds_from_cache("small_model")
- ds2 = ds.copy()
-
- ds2.time.attrs["time_units"] = "MONTHS"
-
- check = nlmod.cache._check_ds(ds, ds2)
-
- assert not check
-
-def test_cache_data_array():
+def test_cache_ahn_data_array():
+ """Test caching of AHN data array. Does not have dataset as argument."""
extent = [119_900, 120_000, 441_900, 442_000]
- ahn_no_cache = nlmod.read.ahn.get_ahn4(extent)
- ahn_cached = nlmod.read.ahn.get_ahn4(extent, cachedir=tmpdir, cachename="ahn4.nc")
- ahn_cache = nlmod.read.ahn.get_ahn4(extent, cachedir=tmpdir, cachename="ahn4.nc")
- assert ahn_cached.equals(ahn_no_cache)
- assert ahn_cache.equals(ahn_no_cache)
-
-
-@pytest.mark.slow
-def test_ds_check_grid_false(tmpdir):
- # two models with a different grid and same time dicretisation
- ds = test_001_model.get_ds_from_cache("small_model")
- ds2 = test_001_model.get_ds_time_transient(tmpdir)
- extent = [99100.0, 99400.0, 489100.0, 489400.0]
- regis_ds = nlmod.read.regis.get_combined_layer_models(
- extent,
- use_regis=True,
- use_geotop=False,
- cachedir=tmpdir,
- cachename="comb.nc",
+ cache_name = "ahn4.nc"
+
+ with tempfile.TemporaryDirectory() as tmpdir:
+ assert not os.path.exists(os.path.join(tmpdir, cache_name)), "Cache should not exist yet1"
+ ahn_no_cache = nlmod.read.ahn.get_ahn4(extent)
+ assert not os.path.exists(os.path.join(tmpdir, cache_name)), "Cache should not exist yet2"
+
+ ahn_cached = nlmod.read.ahn.get_ahn4(extent, cachedir=tmpdir, cachename=cache_name)
+ assert os.path.exists(os.path.join(tmpdir, cache_name)), "Cache should have existed by now"
+ assert ahn_cached.equals(ahn_no_cache)
+ modification_time1 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+
+ # Check if the cache is used. If not, cache is rewritten and modification time changes
+ ahn_cache = nlmod.read.ahn.get_ahn4(extent, cachedir=tmpdir, cachename=cache_name)
+ assert ahn_cache.equals(ahn_no_cache)
+ modification_time2 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+ assert modification_time1 == modification_time2, "Cache should not be rewritten"
+
+ # Different extent should not lead to using the cache
+ extent = [119_800, 120_000, 441_900, 442_000]
+ ahn_cache = nlmod.read.ahn.get_ahn4(extent, cachedir=tmpdir, cachename=cache_name)
+ modification_time3 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+ assert modification_time1 != modification_time3, "Cache should have been rewritten"
+
+
+def test_cache_northsea_data_array():
+ """Test caching of AHN data array. Does have dataset as argument."""
+ from nlmod.read.rws import get_northsea
+ ds1 = nlmod.get_ds(
+ [119_700, 120_000, 441_900, 442_000],
+ delr=100.,
+ delc=100.,
+ top=0.,
+ botm=[-1., -2.],
+ kh=10.,
+ kv=1.,
)
- ds2 = nlmod.base.to_model_ds(regis_ds, delr=50.0, delc=50.0)
-
- check = nlmod.cache._check_ds(ds, ds2)
-
- assert not check
-
-
-@pytest.mark.skip("too slow")
-def test_use_cached_regis(tmpdir):
- extent = [98700.0, 99000.0, 489500.0, 489700.0]
- regis_ds1 = nlmod.read.regis.get_regis(extent, cachedir=tmpdir, cachename="reg.nc")
-
- regis_ds2 = nlmod.read.regis.get_regis(extent, cachedir=tmpdir, cachename="reg.nc")
-
- assert regis_ds1.equals(regis_ds2)
-
-
-@pytest.mark.skip("too slow")
-def test_do_not_use_cached_regis(tmpdir):
- # cache regis
- extent = [98700.0, 99000.0, 489500.0, 489700.0]
- regis_ds1 = nlmod.read.regis.get_regis(
- extent, cachedir=tmpdir, cachename="regis.nc"
- )
-
- # do not use cache because extent is different
- extent = [99100.0, 99400.0, 489100.0, 489400.0]
- regis_ds2 = nlmod.read.regis.get_regis(
- extent, cachedir=tmpdir, cachename="regis.nc"
+ ds2 = nlmod.get_ds(
+ [119_800, 120_000, 441_900, 444_000],
+ delr=100.,
+ delc=100.,
+ top=0.,
+ botm=[-1., -3.],
+ kh=10.,
+ kv=1.,
)
- assert not regis_ds1.equals(regis_ds2)
+ cache_name = "northsea.nc"
+
+ with tempfile.TemporaryDirectory() as tmpdir:
+ assert not os.path.exists(os.path.join(tmpdir, cache_name)), "Cache should not exist yet1"
+ out1_no_cache = get_northsea(ds1)
+ assert not os.path.exists(os.path.join(tmpdir, cache_name)), "Cache should not exist yet2"
+
+ out1_cached = get_northsea(ds1, cachedir=tmpdir, cachename=cache_name)
+ assert os.path.exists(os.path.join(tmpdir, cache_name)), "Cache should exist by now"
+ assert out1_cached.equals(out1_no_cache)
+ modification_time1 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+
+ # Check if the cache is used. If not, cache is rewritten and modification time changes
+ out1_cache = get_northsea(ds1, cachedir=tmpdir, cachename=cache_name)
+ assert out1_cache.equals(out1_no_cache)
+ modification_time2 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+ assert modification_time1 == modification_time2, "Cache should not be rewritten"
+
+ # Only properties of `coords_2d` determine if the cache is used. Cache should still be used.
+ ds1["toppertje"] = ds1.top + 1
+ out1_cache = get_northsea(ds1, cachedir=tmpdir, cachename=cache_name)
+ assert out1_cache.equals(out1_no_cache)
+ modification_time2 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+ assert modification_time1 == modification_time2, "Cache should not be rewritten"
+
+ # Different extent should not lead to using the cache
+ out2_cache = get_northsea(ds2, cachedir=tmpdir, cachename=cache_name)
+ modification_time3 = os.path.getmtime(os.path.join(tmpdir, cache_name))
+ assert modification_time1 != modification_time3, "Cache should have been rewritten"
+ assert not out2_cache.equals(out1_no_cache)
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_014_gis.py b/tests/test_014_gis.py
index b944a7b5..0ae4231f 100644
--- a/tests/test_014_gis.py
+++ b/tests/test_014_gis.py
@@ -17,4 +17,8 @@ def test_vertex_da_to_gdf():
def test_ds_to_ugrid_nc_file():
ds = util.get_ds_vertex()
- nlmod.gis.ds_to_ugrid_nc_file(ds, os.path.join("data", "ugrid_test.nc"))
+ fname = os.path.join("data", "ugrid_test.nc")
+ nlmod.gis.ds_to_ugrid_nc_file(ds, fname)
+
+ fname = os.path.join("data", "ugrid_test_qgis.nc")
+ nlmod.gis.ds_to_ugrid_nc_file(ds, fname, for_imod_qgis_plugin=True)
diff --git a/tests/test_015_gwf_output.py b/tests/test_015_gwf_output.py
index 074e920b..fc3608ee 100644
--- a/tests/test_015_gwf_output.py
+++ b/tests/test_015_gwf_output.py
@@ -67,17 +67,17 @@ def test_create_small_model_grid_only(tmpdir, model_name="test"):
assert np.array_equal(da.values, heads_correct, equal_nan=True)
fname_hds = os.path.join(ds.model_ws, ds.model_name + ".hds")
- grbfile = os.path.join(ds.model_ws, ds.model_name + ".dis.grb")
- da = get_heads_da(ds=None, gwf=None, fname=fname_hds, grbfile=grbfile) # fname
+ grb_file = os.path.join(ds.model_ws, ds.model_name + ".dis.grb")
+ da = get_heads_da(ds=None, gwf=None, fname=fname_hds, grb_file=grb_file) # fname
assert np.array_equal(da.values, heads_correct, equal_nan=True)
# budget
da = get_budget_da("CHD", ds=ds, gwf=None, fname=None) # ds
da = get_budget_da("CHD", ds=None, gwf=gwf, fname=None) # gwf
fname_cbc = os.path.join(ds.model_ws, ds.model_name + ".cbc")
- get_budget_da("CHD", ds=None, gwf=None, fname=fname_cbc, grbfile=grbfile) # fname
+ get_budget_da("CHD", ds=None, gwf=None, fname=fname_cbc, grb_file=grb_file) # fname
get_budget_da(
- "DATA-SPDIS", column="qz", ds=None, gwf=None, fname=fname_cbc, grbfile=grbfile
+ "DATA-SPDIS", column="qz", ds=None, gwf=None, fname=fname_cbc, grb_file=grb_file
) # fname
# unstructured
@@ -127,18 +127,18 @@ def test_create_small_model_grid_only(tmpdir, model_name="test"):
assert np.array_equal(da.values, heads_correct, equal_nan=True)
fname_hds = os.path.join(ds.model_ws, ds.model_name + ".hds")
- grbfile = os.path.join(ds.model_ws, ds.model_name + ".disv.grb")
- da = get_heads_da(ds=None, gwf=None, fname=fname_hds, grbfile=grbfile) # fname
+ grb_file = os.path.join(ds.model_ws, ds.model_name + ".disv.grb")
+ da = get_heads_da(ds=None, gwf=None, fname=fname_hds, grb_file=grb_file) # fname
assert np.array_equal(da.values, heads_correct, equal_nan=True)
# budget
da = get_budget_da("CHD", ds=ds_unstr, gwf=None, fname=None) # ds
da = get_budget_da("CHD", ds=None, gwf=gwf_unstr, fname=None) # gwf
da = get_budget_da(
- "CHD", ds=None, gwf=None, fname=fname_cbc, grbfile=grbfile
+ "CHD", ds=None, gwf=None, fname=fname_cbc, grb_file=grb_file
) # fname
_ = get_budget_da(
- "DATA-SPDIS", column="qz", ds=None, gwf=None, fname=fname_cbc, grbfile=grbfile
+ "DATA-SPDIS", column="qz", ds=None, gwf=None, fname=fname_cbc, grb_file=grb_file
) # fname
@@ -150,8 +150,8 @@ def test_get_heads_da_from_file_structured_no_grb():
def test_get_heads_da_from_file_structured_with_grb():
fname_hds = "./tests/data/mf6output/structured/test.hds"
- grbfile = "./tests/data/mf6output/structured/test.dis.grb"
- nlmod.gwf.output.get_heads_da(fname=fname_hds, grbfile=grbfile)
+ grb_file = "./tests/data/mf6output/structured/test.dis.grb"
+ nlmod.gwf.output.get_heads_da(fname=fname_hds, grb_file=grb_file)
def test_get_budget_da_from_file_structured_no_grb():
@@ -162,8 +162,8 @@ def test_get_budget_da_from_file_structured_no_grb():
def test_get_budget_da_from_file_structured_with_grb():
fname_cbc = "./tests/data/mf6output/structured/test.cbc"
- grbfile = "./tests/data/mf6output/structured/test.dis.grb"
- nlmod.gwf.output.get_budget_da("CHD", fname=fname_cbc, grbfile=grbfile)
+ grb_file = "./tests/data/mf6output/structured/test.dis.grb"
+ nlmod.gwf.output.get_budget_da("CHD", fname=fname_cbc, grb_file=grb_file)
def test_get_heads_da_from_file_vertex_no_grb():
@@ -174,8 +174,8 @@ def test_get_heads_da_from_file_vertex_no_grb():
def test_get_heads_da_from_file_vertex_with_grb():
fname_hds = "./tests/data/mf6output/vertex/test.hds"
- grbfile = "./tests/data/mf6output/vertex/test.disv.grb"
- nlmod.gwf.output.get_heads_da(fname=fname_hds, grbfile=grbfile)
+ grb_file = "./tests/data/mf6output/vertex/test.disv.grb"
+ nlmod.gwf.output.get_heads_da(fname=fname_hds, grb_file=grb_file)
def test_get_budget_da_from_file_vertex_no_grb():
@@ -186,8 +186,8 @@ def test_get_budget_da_from_file_vertex_no_grb():
def test_get_budget_da_from_file_vertex_with_grb():
fname_cbc = "./tests/data/mf6output/vertex/test.cbc"
- grbfile = "./tests/data/mf6output/vertex/test.disv.grb"
- nlmod.gwf.output.get_budget_da("CHD", fname=fname_cbc, grbfile=grbfile)
+ grb_file = "./tests/data/mf6output/vertex/test.disv.grb"
+ nlmod.gwf.output.get_budget_da("CHD", fname=fname_cbc, grb_file=grb_file)
def test_postprocess_head():
@@ -199,3 +199,19 @@ def test_postprocess_head():
nlmod.gwf.get_gwl_from_wet_cells(head, botm=ds["botm"])
nlmod.gwf.get_head_at_point(head, float(ds.x.mean()), float(ds.y.mean()), ds=ds)
+
+
+def test_get_flow_residuals():
+ ds = test_001_model.get_ds_from_cache("basic_sea_model")
+ da = nlmod.gwf.output.get_flow_residuals(ds)
+ assert "time" in da.dims
+ da = nlmod.gwf.output.get_flow_residuals(ds, kstpkper=(0, 0))
+ assert "time" not in da.dims
+
+
+def test_get_flow_lower_face():
+ ds = test_001_model.get_ds_from_cache("basic_sea_model")
+ da = nlmod.gwf.output.get_flow_lower_face(ds)
+ assert "time" in da.dims
+ da = nlmod.gwf.output.get_flow_lower_face(ds, kstpkper=(0, 0))
+ assert "time" not in da.dims
diff --git a/tests/test_016_time.py b/tests/test_016_time.py
index 1e3339c5..056f7a75 100644
--- a/tests/test_016_time.py
+++ b/tests/test_016_time.py
@@ -28,3 +28,7 @@ def test_ds_time_from_tdis_settings():
elapsed = (tidx.to_numpy() - np.datetime64("2000")) / np.timedelta64(1, "D")
assert np.allclose(elapsed, [100, 150, 200, 233.33333333, 300.0])
+
+
+def test_get_time_step_length():
+ assert (nlmod.time.get_time_step_length(100, 2, 1.5) == np.array([40, 60])).all()
diff --git a/tests/test_018_knmi_data_platform.py b/tests/test_018_knmi_data_platform.py
index 9c53bbb0..05da19c9 100644
--- a/tests/test_018_knmi_data_platform.py
+++ b/tests/test_018_knmi_data_platform.py
@@ -1,3 +1,4 @@
+# ruff: noqa: D103
import os
from pathlib import Path
@@ -19,10 +20,10 @@ def test_download_multiple_nc_files() -> None:
)
# download the last 10 files
- fnames = files[-10:]
+ fnames = files[:2]
dirname = "download"
knmi_data_platform.download_files(
- dataset_name, dataset_version, files[-10:], dirname=dirname
+ dataset_name, dataset_version, fnames, dirname=dirname
)
ds = knmi_data_platform.read_nc(os.path.join(dirname, fnames[0]))
@@ -40,7 +41,7 @@ def test_download_read_zip_file() -> None:
# download the last file
dirname = "download"
- fname = files[-1]
+ fname = files[1]
knmi_data_platform.download_file(
dataset_name, dataset_version, fname=fname, dirname=dirname
)
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 af768339..8be1c5be 100644
--- a/tests/test_021_nhi.py
+++ b/tests/test_021_nhi.py
@@ -1,14 +1,19 @@
+# ruff: noqa: D103
import os
-import numpy as np
import tempfile
-import nlmod
+
+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)
@@ -20,3 +25,44 @@ def test_buidrainage():
# assert that all locations with a positive conductance also have a specified depth
mask = ds["buisdrain_cond"] > 0
assert np.all(~np.isnan(ds["buisdrain_depth"].data[mask]))
+
+
+def test_gwo():
+ username = os.environ["NHI_GWO_USERNAME"]
+ password = os.environ["NHI_GWO_PASSWORD"]
+
+ # download all wells from Brabant Water
+ wells = nlmod.read.nhi.get_gwo_wells(
+ username=username, password=password, organisation="Brabant Water"
+ )
+ assert isinstance(wells, gpd.GeoDataFrame)
+
+ # download extractions from well "13-PP016" of pomping station Veghel
+ measurements, gdf = nlmod.read.nhi.get_gwo_measurements(
+ username, password, well_site="veghel", filter__well__name="13-PP016"
+ )
+ assert measurements.reset_index()["Name"].isin(gdf.index).all()
+
+
+@pytest.mark.skip("too slow")
+def test_gwo_entire_pumping_station():
+ username = os.environ["NHI_GWO_USERNAME"]
+ password = os.environ["NHI_GWO_PASSWORD"]
+ measurements, gdf = nlmod.read.nhi.get_gwo_measurements(
+ username,
+ password,
+ well_site="veghel",
+ )
+ assert measurements.reset_index()["Name"].isin(gdf.index).all()
+
+ ncols = 3
+ nrows = int(np.ceil(len(gdf.index) / ncols))
+ f, axes = plt.subplots(
+ nrows=nrows, ncols=ncols, figsize=(10, 10), sharex=True, sharey=True
+ )
+ axes = axes.ravel()
+ for name, ax in zip(gdf.index, axes):
+ measurements.loc[name, "Volume"].plot(ax=ax)
+ ax.set_xlabel("")
+ ax.set_title(name)
+ f.tight_layout(pad=0.0)
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
new file mode 100644
index 00000000..80b001e7
--- /dev/null
+++ b/tests/test_026_grid.py
@@ -0,0 +1,217 @@
+import os
+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")
+extent = [98000.0, 99000.0, 489000.0, 490000.0]
+
+
+def get_bgt():
+ fname = os.path.join(model_ws, "bgt.gpkg")
+ if not os.path.isfile(fname):
+ if not os.path.isdir(model_ws):
+ os.makedirs(model_ws)
+ bgt = nlmod.read.bgt.get_bgt(extent)
+ bgt.to_file(fname)
+ return gpd.read_file(fname)
+
+
+def get_regis():
+ fname = os.path.join(model_ws, "regis.nc")
+ if not os.path.isfile(fname):
+ if not os.path.isdir(model_ws):
+ os.makedirs(model_ws)
+ regis = nlmod.read.regis.get_regis(extent)
+ regis.to_netcdf(fname)
+ return xr.open_dataset(fname)
+
+
+def get_structured_model_ds():
+ model_ws = os.path.join(tempfile.gettempdir(), "test_grid_structured")
+ fname = os.path.join(model_ws, "ds.nc")
+ if not os.path.isfile(fname):
+ if not os.path.isdir(model_ws):
+ os.makedirs(model_ws)
+ ds = nlmod.get_ds(extent, model_name="test_grid", model_ws=model_ws)
+ ds.to_netcdf(fname)
+ return xr.open_dataset(fname)
+
+
+def get_structured_model_ds_rotated():
+ model_ws = os.path.join(tempfile.gettempdir(), "test_grid_structured_rotated")
+ fname = os.path.join(model_ws, "ds.nc")
+ if not os.path.isfile(fname):
+ if not os.path.isdir(model_ws):
+ os.makedirs(model_ws)
+ ds = nlmod.get_ds(extent, model_name="test_grid", model_ws=model_ws, angrot=15)
+ ds.to_netcdf(fname)
+ return xr.open_dataset(fname)
+
+
+def get_vertex_model_ds(bgt=None):
+ model_ws = os.path.join(tempfile.gettempdir(), "test_grid_vertex")
+ fname = os.path.join(model_ws, "ds.nc")
+ if not os.path.isfile(fname):
+ if not os.path.isdir(model_ws):
+ os.makedirs(model_ws)
+ ds = get_structured_model_ds()
+ if bgt is None:
+ bgt = get_bgt()
+ ds = nlmod.grid.refine(ds, model_ws=model_ws, refinement_features=[(bgt, 1)])
+ ds.to_netcdf(fname)
+ return xr.open_dataset(fname)
+
+
+def get_vertex_model_ds_rotated(bgt=None):
+ model_ws = os.path.join(tempfile.gettempdir(), "test_grid_vertex_rotated")
+ fname = os.path.join(model_ws, "ds.nc")
+ if not os.path.isfile(fname):
+ if not os.path.isdir(model_ws):
+ os.makedirs(model_ws)
+ ds = get_structured_model_ds_rotated()
+ if bgt is None:
+ bgt = get_bgt()
+ ds = nlmod.grid.refine(ds, model_ws=model_ws, refinement_features=[(bgt, 1)])
+ ds.to_netcdf(fname)
+ return xr.open_dataset(fname)
+
+
+def test_get_ds_rotated():
+ ds0 = get_structured_model_ds_rotated()
+ assert ds0.extent[0] == 0 and ds0.extent[2] == 0
+ assert ds0.xorigin == extent[0] and ds0.yorigin == extent[2]
+
+ # test refine method, by refining in all cells that contain surface water polygons
+ ds = get_vertex_model_ds_rotated()
+ assert len(ds.area) > np.prod(ds0.area.shape)
+ assert ds.extent[0] == 0 and ds.extent[2] == 0
+ assert ds.xorigin == extent[0] and ds.yorigin == extent[2]
+
+ f0, ax0 = plt.subplots()
+ nlmod.plot.modelgrid(ds0, ax=ax0)
+ f, ax = plt.subplots()
+ nlmod.plot.modelgrid(ds, ax=ax)
+ assert (np.array(ax.axis()) == np.array(ax0.axis())).all()
+
+
+def test_vertex_da_to_ds():
+ # for a normal grid
+ ds0 = get_structured_model_ds()
+ ds = get_vertex_model_ds()
+ da = nlmod.resample.vertex_da_to_ds(ds["top"], ds0, method="linear")
+ assert not da.isnull().all()
+ da = nlmod.resample.vertex_da_to_ds(ds["botm"], ds0, method="linear")
+ assert not da.isnull().all()
+
+ # for a rotated grid
+ ds0 = get_structured_model_ds_rotated()
+ ds = get_vertex_model_ds_rotated()
+ da = nlmod.resample.vertex_da_to_ds(ds["top"], ds0, method="linear")
+ assert not da.isnull().all()
+ da = nlmod.resample.vertex_da_to_ds(ds["botm"], ds0, method="linear")
+ assert not da.isnull().all()
+
+
+def test_fillnan_da():
+ # for a structured grid
+ ds = get_structured_model_ds()
+ 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
+ mask = ds["top"].isnull()
+ assert mask.any()
+ top = nlmod.resample.fillnan_da(ds["top"], ds=ds)
+ assert not top[mask].isnull().any()
+
+
+def test_gdf_to_bool_da():
+ bgt = get_bgt()
+
+ # test for a structured grid
+ ds = get_structured_model_ds()
+ da = nlmod.grid.gdf_to_bool_da(bgt, ds)
+ assert da.any()
+
+ # test for a vertex grid
+ ds = get_vertex_model_ds()
+ da = nlmod.grid.gdf_to_bool_da(bgt, ds)
+ assert da.any()
+
+ # tets for a slightly rotated structured grid
+ ds = get_structured_model_ds_rotated()
+ da = nlmod.grid.gdf_to_bool_da(bgt, ds)
+ assert da.any()
+
+ # test for a rotated vertex grid
+ ds = get_vertex_model_ds_rotated()
+ da = nlmod.grid.gdf_to_bool_da(bgt, ds)
+ assert da.any()
+
+
+def test_gdf_to_da():
+ bgt = get_bgt()
+
+ # test for a structured grid
+ ds = get_structured_model_ds()
+ da = nlmod.grid.gdf_to_da(bgt, ds, "relatieveHoogteligging", agg_method="max_area")
+ assert not da.isnull().all()
+
+ # test for a vertex grid
+ ds = get_vertex_model_ds()
+ da = nlmod.grid.gdf_to_da(bgt, ds, "relatieveHoogteligging", agg_method="max_area")
+ assert not da.isnull().all()
+
+ # tets for a slightly rotated structured grid
+ ds = get_structured_model_ds_rotated()
+ da = nlmod.grid.gdf_to_da(bgt, ds, "relatieveHoogteligging", agg_method="max_area")
+ assert not da.isnull().all()
+
+ # test for a rotated vertex grid
+ ds = get_vertex_model_ds_rotated()
+ da = nlmod.grid.gdf_to_da(bgt, ds, "relatieveHoogteligging", agg_method="max_area")
+ assert not da.isnull().all()
+
+
+def test_update_ds_from_layer_ds():
+ bgt = get_bgt()
+ regis = get_regis()
+
+ # test for a structured grid
+ ds = nlmod.get_ds(extent, delr=200)
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="nearest")
+ assert len(np.unique(ds["top"])) > 1
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="average")
+ assert len(np.unique(ds["top"])) > 1
+
+ # test for a vertex grid
+ model_ws = os.path.join(tempfile.gettempdir(), "test_grid_vertex_200")
+ ds = nlmod.grid.refine(ds, model_ws=model_ws, refinement_features=[(bgt, 1)])
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="nearest")
+ assert len(np.unique(ds["top"])) > 1
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="average")
+ assert len(np.unique(ds["top"])) > 1
+
+ # tets for a slightly rotated structured grid
+ ds = nlmod.get_ds(extent, delr=200, angrot=15)
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="nearest")
+ assert len(np.unique(ds["top"])) > 1
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="average")
+ assert len(np.unique(ds["top"])) > 1
+
+ # test for a rotated vertex grid
+ model_ws = os.path.join(tempfile.gettempdir(), "test_grid_vertex_200_rotated")
+ ds = nlmod.grid.refine(ds, model_ws=model_ws, refinement_features=[(bgt, 2)])
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="nearest")
+ assert len(np.unique(ds["top"])) > 1
+ ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="average")
+ assert len(np.unique(ds["top"])) > 1