add modelextent plot

.prospector.yaml

@@ -29,6 +29,8 @@ pylint:
     - too-many-branches
     - too-many-statements
     - logging-fstring-interpolation
+    - import-outside-toplevel
+    - implicit-str-concat
 
 mccabe:
   disable:

docs/examples/01_basic_model.ipynb

@@ -18,7 +18,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
     "import nlmod"
    ]
   },

docs/examples/04_modifying_layermodels.ipynb

@@ -455,7 +455,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Set mimimum layer thickness\n",
+    "## Set minimum layer thickness\n",
     "nlmod.layers.set_minimum layer_thickness increases the thickness of a layer if the thickness is less than a specified value.  With a parameter called 'change' you can specify in which direction the layer is changed. The only supported option for now is 'botm', which changes the layer botm. \n",
     "\n",
     "This method only changes the shape of the layer, and does not check if all hydrological properties are defined for cells that had a thickness of 0 before."
@@ -467,7 +467,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# set the mimimum thickness of 'PZWAz2' to 20 m\n",
+    "# set the minimum thickness of 'PZWAz2' to 20 m\n",
     "ds_new = nlmod.layers.set_minimum_layer_thickness(ds.copy(deep=True), \"PZWAz2\", 20.0)\n",
     "compare_layer_models(ds, line, colors, ds2=ds_new, title2=\"Modified\")"
    ]

docs/examples/11_grid_rotation.ipynb

@@ -131,7 +131,7 @@
    "outputs": [],
    "source": [
     "# Download AHN\n",
-    "extent = nlmod.resample.get_extent(ds)\n",
+    "extent = nlmod.grid.get_extent(ds)\n",
     "ahn = nlmod.read.ahn.get_ahn3(extent)\n",
     "\n",
     "# Resample to the grid\n",
@@ -320,8 +320,8 @@
     "cbar = nlmod.plot.colorbar_inside(pc)\n",
     "# as the surface water shapes are in model coordinates, we need to transform them\n",
     "# to real-world coordinates before plotting\n",
-    "affine = nlmod.resample.get_affine_mod_to_world(ds)\n",
-    "bgt_rw = nlmod.resample.affine_transform_gdf(bgt, affine)\n",
+    "affine = nlmod.grid.get_affine_mod_to_world(ds)\n",
+    "bgt_rw = nlmod.grid.affine_transform_gdf(bgt, affine)\n",
     "bgt_rw.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
    ]
   },

nlmod/dims/__init__.py

@@ -2,7 +2,7 @@
 from . import base, grid, layers, resample, time
 from .attributes_encodings import *
 from .base import *
-from .grid import *
-from .layers import *
 from .resample import *
+from .grid import *  # import from grid after resample, to ignore deprecated methods
+from .layers import *
 from .time import *

nlmod/dims/base.py

@@ -8,7 +8,7 @@
 
 from .. import util
 from ..epsg28992 import EPSG_28992
-from . import resample
+from . import grid, resample
 from .layers import fill_nan_top_botm_kh_kv
 
 logger = logging.getLogger(__name__)
@@ -362,7 +362,7 @@ def _get_structured_grid_ds(
     }
 
     if angrot != 0.0:
-        affine = resample.get_affine_mod_to_world(attrs)
+        affine = grid.get_affine_mod_to_world(attrs)
         xc, yc = affine * np.meshgrid(xcenters, ycenters)
         coords["xc"] = (("y", "x"), xc)
         coords["yc"] = (("y", "x"), yc)
@@ -643,7 +643,7 @@ def get_ds(
     x, y = resample.get_xy_mid_structured(attrs["extent"], delr, delc)
     coords = {"x": x, "y": y, "layer": layer}
     if angrot != 0.0:
-        affine = resample.get_affine_mod_to_world(attrs)
+        affine = grid.get_affine_mod_to_world(attrs)
         xc, yc = affine * np.meshgrid(x, y)
         coords["xc"] = (("y", "x"), xc)
         coords["yc"] = (("y", "x"), yc)

nlmod/dims/grid.py

@@ -16,6 +16,7 @@
 import pandas as pd
 import shapely
 import xarray as xr
+from affine import Affine
 from flopy.discretization.structuredgrid import StructuredGrid
 from flopy.discretization.vertexgrid import VertexGrid
 from flopy.utils.gridgen import Gridgen
@@ -27,22 +28,13 @@
 from tqdm import tqdm
 
 from .. import cache, util
-from .base import _get_structured_grid_ds, _get_vertex_grid_ds, extrapolate_ds
 from .layers import (
     fill_nan_top_botm_kh_kv,
     get_first_active_layer,
     get_idomain,
     remove_inactive_layers,
 )
 from .rdp import rdp
-from .resample import (
-    affine_transform_gdf,
-    get_affine_mod_to_world,
-    get_affine_world_to_mod,
-    get_delc,
-    get_delr,
-    structured_da_to_ds,
-)
 
 logger = logging.getLogger(__name__)
 
@@ -227,6 +219,58 @@ def modelgrid_from_ds(ds, rotated=True, nlay=None, top=None, botm=None, **kwargs
     return modelgrid
 
 
+def get_delr(ds):
+    """
+    Get the distance along rows (delr) from the x-coordinate of a structured model
+    dataset.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        A model dataset containing an x-coordinate and an attribute 'extent'.
+
+    Returns
+    -------
+    delr : np.ndarray
+        The cell-size along rows (of length ncol).
+
+    """
+    assert ds.gridtype == "structured"
+    x = (ds.x - ds.extent[0]).values
+    delr = _get_delta_along_axis(x)
+    return delr
+
+
+def get_delc(ds):
+    """
+    Get the distance along columns (delc) from the y-coordinate of a structured model
+    dataset.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        A model dataset containing an y-coordinate and an attribute 'extent'.
+
+    Returns
+    -------
+    delc : np.ndarray
+        The cell-size along columns (of length nrow).
+
+    """
+    assert ds.gridtype == "structured"
+    y = (ds.extent[3] - ds.y).values
+    delc = _get_delta_along_axis(y)
+    return delc
+
+
+def _get_delta_along_axis(x):
+    """Internal method to determine delr or delc from x or y relative to xmin or ymax"""
+    delr = [x[0] * 2]
+    for xi in x[1:]:
+        delr.append((xi - np.sum(delr)) * 2)
+    return np.array(delr)
+
+
 def modelgrid_to_vertex_ds(mg, ds, nodata=-1):
     """Add information about the calculation-grid to a model dataset."""
     warnings.warn(
@@ -264,6 +308,7 @@ def modelgrid_to_ds(mg):
     """
     if mg.grid_type == "structured":
         x, y = mg.xyedges
+        from .base import _get_structured_grid_ds
 
         ds = _get_structured_grid_ds(
             xedges=x,
@@ -278,6 +323,8 @@ def modelgrid_to_ds(mg):
             crs=None,
         )
     elif mg.grid_type == "vertex":
+        from .base import _get_vertex_grid_ds
+
         ds = _get_vertex_grid_ds(
             x=mg.xcellcenters,
             y=mg.ycellcenters,
@@ -573,6 +620,8 @@ def ds_to_gridprops(ds_in, gridprops, method="nearest", icvert_nodata=-1):
         ds_out.coords.update({"layer": ds_in.layer, "x": x, "y": y})
 
         # add other variables
+        from .resample import structured_da_to_ds
+
         for not_interp_var in not_interp_vars:
             ds_out[not_interp_var] = structured_da_to_ds(
                 da=ds_in[not_interp_var], ds=ds_out, method=method, nodata=np.nan
@@ -680,8 +729,12 @@ def update_ds_from_layer_ds(ds, layer_ds, method="nearest", **kwargs):
         for var in layer_ds.data_vars:
             ds[var] = layer_ds[var]
     else:
+        from .resample import structured_da_to_ds
+
         for var in layer_ds.data_vars:
             ds[var] = structured_da_to_ds(layer_ds[var], ds, method=method)
+    from .base import extrapolate_ds
+
     ds = extrapolate_ds(ds)
     ds = fill_nan_top_botm_kh_kv(ds, **kwargs)
     return ds
@@ -1903,7 +1956,7 @@ def mask_model_edge(ds, idomain=None):
             ds["vertices"] = get_vertices(ds)
         polygons_grid = polygons_from_model_ds(ds)
         gdf_grid = gpd.GeoDataFrame(geometry=polygons_grid)
-        extent_edge = util.polygon_from_extent(ds.extent).exterior
+        extent_edge = get_extent_polygon(ds).exterior
         cids_edge = gdf_grid.loc[gdf_grid.touches(extent_edge)].index
         ds_out["edge_mask"] = util.get_da_from_da_ds(
             ds, dims=("layer", "icell2d"), data=0
@@ -1970,3 +2023,97 @@ def polygons_from_model_ds(model_ds):
         polygons = [affine_transform(polygon, affine) for polygon in polygons]
 
     return polygons
+
+
+def _get_attrs(ds):
+    if isinstance(ds, dict):
+        return ds
+    else:
+        return ds.attrs
+
+
+def get_extent_polygon(ds, rotated=True):
+    """Get the model extent, as a shapely Polygon."""
+    attrs = _get_attrs(ds)
+    polygon = util.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
+
+
+def get_extent_gdf(ds, rotated=True, crs="EPSG:28992"):
+    polygon = get_extent_polygon(ds, rotated=rotated)
+    return gpd.GeoDataFrame(geometry=[polygon], crs=crs)
+
+
+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
+
+
+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
+
+
+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)
+
+
+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)
+
+
+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 = get_delr(ds)
+        assert len(np.unique(sx)) == 1, "Affine-transformation needs a constant delr"
+        sx = sx[0]
+    if sy is None:
+        sy = get_delc(ds)
+        assert len(np.unique(sy)) == 1, "Affine-transformation needs a constant delc"
+        sy = sy[0]
+
+    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)

nlmod/dims/rdp.py

@@ -1,8 +1,8 @@
 """
-rdp 
-~~~ 
+rdp
+~~~
 Python implementation of the Ramer-Douglas-Peucker algorithm.
-:copyright: 2014-2016 Fabian Hirschmann <[email protected]> 
+:copyright: 2014-2016 Fabian Hirschmann <[email protected]>
 :license: MIT, see LICENSE.txt for more details.
 """