From 60ead3d4fcedc2efa19268638ca8b803173971c8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= <d.brakenhoff@artesia-water.nl>
Date: Thu, 27 Jun 2024 13:39:39 +0200
Subject: [PATCH 01/10] add modelextent plot

---
 nlmod/plot/plot.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py
index 7ff8d0fe..2134d556 100644
--- a/nlmod/plot/plot.py
+++ b/nlmod/plot/plot.py
@@ -11,6 +11,8 @@
 from matplotlib.colors import ListedColormap, Normalize
 from matplotlib.patches import Patch
 from mpl_toolkits.axes_grid1 import make_axes_locatable
+from geopandas import GeoDataFrame
+from shapely.geometry import Polygon
 
 from ..dims.grid import modelgrid_from_ds
 from ..dims.resample import get_affine_mod_to_world, get_extent
@@ -48,6 +50,47 @@ def modelgrid(ds, ax=None, **kwargs):
     return ax
 
 
+def modelextent(ds, ax=None, **kwargs):
+    """Plot model extent.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        The dataset containing the data.
+    ax : matplotlib.axes.Axes, optional
+        The axes object to plot on. If not provided, a new figure and axes will be
+        created.
+    **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=True)
+    dx = 0.05 * (xmax - xmin)
+    dy = 0.05 * (ymax - ymin)
+    if ax is None:
+        _, ax = plt.subplots(figsize=(10, 10))
+        ax.axis("scaled")
+
+        ax.axis([xmin - dx, xmax + dx, ymin - dy, ymax + dy])
+    xy = [
+        (xmin, ymin),
+        (xmax, ymin),
+        (xmax, ymax),
+        (xmin, ymax),
+        (xmin, ymin),
+    ]
+    gdf = GeoDataFrame(geometry=[Polygon(xy)])
+    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
+
+
 def facet_plot(
     gwf,
     ds,

From 0c295be3cd50cdd6e02b694d0277fdfe2ffe7635 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 2 Jul 2024 13:47:48 +0200
Subject: [PATCH 02/10] Move methods from dims.reample to dims.grid and util

---
 nlmod/dims/base.py         |   6 +-
 nlmod/dims/grid.py         | 159 +++++++++++++++++++++++++++++++--
 nlmod/dims/resample.py     | 177 +++++++++++--------------------------
 nlmod/gis.py               |   3 +-
 nlmod/gwf/gwf.py           |   2 +-
 nlmod/gwf/output.py        |   3 +-
 nlmod/gwf/surface_water.py |   9 +-
 nlmod/mfoutput/mfoutput.py |   7 +-
 nlmod/plot/dcs.py          |   3 +-
 nlmod/plot/plot.py         |  40 +++++----
 nlmod/plot/plotutil.py     |   2 +-
 nlmod/read/ahn.py          |   3 +-
 nlmod/read/bgt.py          |   2 +-
 nlmod/read/jarkus.py       |   4 +-
 nlmod/read/knmi.py         |   2 +-
 nlmod/util.py              |  59 ++++++++++++-
 16 files changed, 306 insertions(+), 175 deletions(-)

diff --git a/nlmod/dims/base.py b/nlmod/dims/base.py
index c2b729b6..43c17363 100644
--- a/nlmod/dims/base.py
+++ b/nlmod/dims/base.py
@@ -8,7 +8,7 @@
 
 from .. import util
 from ..epsg28992 import EPSG_28992
-from . import resample
+from . import resample, grid
 from .layers import fill_nan_top_botm_kh_kv
 
 logger = logging.getLogger(__name__)
@@ -364,7 +364,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)
@@ -645,7 +645,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)
diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 2aadff70..50245f61 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -22,6 +22,7 @@
 from flopy.utils.gridintersect import GridIntersect
 from packaging import version
 from scipy.interpolate import griddata
+from affine import Affine
 from shapely.affinity import affine_transform
 from shapely.geometry import Point, Polygon
 from tqdm import tqdm
@@ -35,14 +36,6 @@
     remove_inactive_layers,
 )
 from .rdp import rdp
-from .resample import (
-    affine_transform_gdf,
-    get_affine_mod_to_world,
-    get_affine_world_to_mod,
-    structured_da_to_ds,
-    get_delr,
-    get_delc,
-)
 
 logger = logging.getLogger(__name__)
 
@@ -229,6 +222,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."""
 
@@ -575,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
@@ -683,6 +730,8 @@ 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)
     ds = extrapolate_ds(ds)
@@ -1978,3 +2027,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)
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index f3ff7fe1..ff88491a 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -4,14 +4,13 @@
 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
 
+
 logger = logging.getLogger(__name__)
 
 
@@ -92,58 +91,6 @@ def get_xy_mid_structured(extent, delr, delc, descending_y=True):
         raise TypeError("unexpected type for delr and/or delc")
 
 
-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 ds_to_structured_grid(
     ds_in,
     extent,
@@ -576,6 +523,8 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
             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)),
@@ -595,6 +544,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)
@@ -635,98 +586,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 = 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)
+    logger.warning(
+        "nlmod.resample.get_affine is deprecated. " "Use nlmod.grid.get_affine instead"
+    )
+    from .grid import get_affine
+
+    return get_affine(ds)
diff --git a/nlmod/gis.py b/nlmod/gis.py
index 17d4709c..21941092 100644
--- a/nlmod/gis.py
+++ b/nlmod/gis.py
@@ -4,8 +4,7 @@
 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 .dims.grid import polygons_from_model_ds, get_affine_mod_to_world
 from .dims.layers import calculate_thickness
 
 logger = logging.getLogger(__name__)
diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py
index ce45e29c..65e17cdc 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -11,7 +11,7 @@
 import xarray as xr
 
 from ..dims import grid
-from ..dims.resample import get_delr, get_delc
+from ..dims.grid import get_delr, get_delc
 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
diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py
index f0d0f356..a296cd0f 100644
--- a/nlmod/gwf/output.py
+++ b/nlmod/gwf/output.py
@@ -7,8 +7,7 @@
 import xarray as xr
 from shapely.geometry import Point
 
-from ..dims.grid import modelgrid_from_ds
-from ..dims.resample import get_affine_world_to_mod
+from ..dims.grid import modelgrid_from_ds, get_affine_world_to_mod
 from ..mfoutput.mfoutput import (
     _get_budget_da,
     _get_heads_da,
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index c61f646e..9a3851c6 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -10,9 +10,14 @@
 from shapely.strtree import STRtree
 from tqdm import tqdm
 
-from ..dims.grid import gdf_to_grid
+from ..util import extent_to_polygon
+from ..dims.grid import (
+    gdf_to_grid,
+    get_extent_polygon,
+    get_delr,
+    get_delc,
+)
 from ..dims.layers import get_idomain
-from ..dims.resample import get_extent_polygon, extent_to_polygon, get_delr, get_delc
 from ..read import bgt, waterboard
 from ..cache import cache_pickle
 
diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py
index f7b21ed5..e764b718 100644
--- a/nlmod/mfoutput/mfoutput.py
+++ b/nlmod/mfoutput/mfoutput.py
@@ -7,8 +7,11 @@
 
 import flopy
 
-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_dims_coords_from_modelgrid,
+    modelgrid_from_ds,
+    get_affine_mod_to_world,
+)
 from ..dims.time import ds_time_idx
 from .binaryfile import _get_binary_budget_data, _get_binary_head_data
 
diff --git a/nlmod/plot/dcs.py b/nlmod/plot/dcs.py
index f27162b0..b89d08f9 100644
--- a/nlmod/plot/dcs.py
+++ b/nlmod/plot/dcs.py
@@ -15,8 +15,7 @@
 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 modelgrid_from_ds, get_affine_world_to_mod
 from .plotutil import get_map
 
 
diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py
index 2134d556..b1c4f0b2 100644
--- a/nlmod/plot/plot.py
+++ b/nlmod/plot/plot.py
@@ -11,11 +11,13 @@
 from matplotlib.colors import ListedColormap, Normalize
 from matplotlib.patches import Patch
 from mpl_toolkits.axes_grid1 import make_axes_locatable
-from geopandas import GeoDataFrame
-from shapely.geometry import Polygon
 
-from ..dims.grid import modelgrid_from_ds
-from ..dims.resample import get_affine_mod_to_world, get_extent
+from ..dims.grid import (
+    modelgrid_from_ds,
+    get_affine_mod_to_world,
+    get_extent,
+    get_extent_gdf,
+)
 from ..read import geotop, rws
 from .dcs import DatasetCrossSection
 from .plotutil import (
@@ -50,7 +52,7 @@ def modelgrid(ds, ax=None, **kwargs):
     return ax
 
 
-def modelextent(ds, ax=None, **kwargs):
+def modelextent(ds, dx=None, ax=None, rotated=True, **kwargs):
     """Plot model extent.
 
     Parameters
@@ -68,22 +70,15 @@ def modelextent(ds, ax=None, **kwargs):
     ax : matplotlib.axes.Axes
         axes object
     """
-    extent = xmin, xmax, ymin, ymax = get_extent(ds, rotated=True)
-    dx = 0.05 * (xmax - xmin)
-    dy = 0.05 * (ymax - ymin)
+    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 - dy, ymax + dy])
-    xy = [
-        (xmin, ymin),
-        (xmax, ymin),
-        (xmax, ymax),
-        (xmin, ymax),
-        (xmin, ymin),
-    ]
-    gdf = GeoDataFrame(geometry=[Polygon(xy)])
+        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:
@@ -259,7 +254,14 @@ def data_array(da, ds=None, ax=None, rotated=False, edgecolor=None, **kwargs):
 
 
 def geotop_lithok_in_cross_section(
-    line, gt=None, ax=None, legend=True, legend_loc=None, lithok_props=None, alpha=None, **kwargs
+    line,
+    gt=None,
+    ax=None,
+    legend=True,
+    legend_loc=None,
+    lithok_props=None,
+    alpha=None,
+    **kwargs,
 ):
     """PLot the lithoclass-data of GeoTOP in a cross-section.
 
@@ -312,7 +314,7 @@ def geotop_lithok_in_cross_section(
     cs = DatasetCrossSection(gt, line, layer="z", ax=ax, **kwargs)
     array, cmap, norm = _get_geotop_cmap_and_norm(gt["lithok"], lithok_props)
     cs.plot_array(array, norm=norm, cmap=cmap, alpha=alpha)
-    
+
     if legend:
         # make a legend with dummy handles
         _add_geotop_lithok_legend(lithok_props, ax, lithok=gt["lithok"], loc=legend_loc)
diff --git a/nlmod/plot/plotutil.py b/nlmod/plot/plotutil.py
index 00c134e3..07908b31 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
 
 
diff --git a/nlmod/read/ahn.py b/nlmod/read/ahn.py
index 1bd2fda8..0c74eba1 100644
--- a/nlmod/read/ahn.py
+++ b/nlmod/read/ahn.py
@@ -13,7 +13,8 @@
 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
 
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/jarkus.py b/nlmod/read/jarkus.py
index 5d962973..1f342882 100644
--- a/nlmod/read/jarkus.py
+++ b/nlmod/read/jarkus.py
@@ -7,6 +7,7 @@
 
 Note: if you like jazz please check this out: https://www.northseajazz.com
 """
+
 import datetime as dt
 import logging
 import os
@@ -17,7 +18,8 @@
 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__)
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index b142848c..ddfc317d 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -8,7 +8,7 @@
 
 from .. import cache, util
 from ..dims.layers import get_first_active_layer
-from ..dims.resample import get_affine_mod_to_world
+from ..dims.grid import get_affine_mod_to_world
 
 logger = logging.getLogger(__name__)
 
diff --git a/nlmod/util.py b/nlmod/util.py
index 16c65a6d..e3b996ff 100644
--- a/nlmod/util.py
+++ b/nlmod/util.py
@@ -13,7 +13,7 @@
 import requests
 import xarray as xr
 from colorama import Back, Fore, Style
-from shapely.geometry import box
+from shapely.geometry import box, Polygon
 
 logger = logging.getLogger(__name__)
 
@@ -561,6 +561,51 @@ 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.
 
@@ -574,7 +619,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))
 
@@ -597,7 +645,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)
 
@@ -621,7 +672,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()

From 6e12f995e2ffa28e499c307660d99172a273a321 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 2 Jul 2024 14:17:50 +0200
Subject: [PATCH 03/10] Remove circular import

---
 nlmod/dims/grid.py     | 6 +++++-
 nlmod/dims/resample.py | 3 ---
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 50245f61..f312e48e 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -28,7 +28,6 @@
 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,
@@ -312,6 +311,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,
@@ -326,6 +326,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,
@@ -734,6 +736,8 @@ def update_ds_from_layer_ds(ds, layer_ds, method="nearest", **kwargs):
 
         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
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index ff88491a..7311a6a0 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -6,11 +6,8 @@
 import xarray as xr
 from scipy.interpolate import griddata
 from scipy.spatial import cKDTree
-
-
 from ..util import get_da_from_da_ds
 
-
 logger = logging.getLogger(__name__)
 
 

From a04b0c65a51672416439de6cbe8f7e5f0e4f648c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= <d.brakenhoff@artesia-water.nl>
Date: Tue, 2 Jul 2024 15:35:26 +0200
Subject: [PATCH 04/10] sort imports

---
 nlmod/dims/base.py         | 2 +-
 nlmod/dims/grid.py         | 2 +-
 nlmod/dims/rdp.py          | 6 +++---
 nlmod/dims/resample.py     | 1 +
 nlmod/gis.py               | 2 +-
 nlmod/gwf/gwf.py           | 2 +-
 nlmod/gwf/output.py        | 2 +-
 nlmod/gwf/surface_water.py | 6 +++---
 nlmod/mfoutput/mfoutput.py | 2 +-
 nlmod/plot/dcs.py          | 2 +-
 nlmod/plot/plot.py         | 2 +-
 nlmod/read/knmi.py         | 2 +-
 nlmod/util.py              | 4 +---
 13 files changed, 17 insertions(+), 18 deletions(-)

diff --git a/nlmod/dims/base.py b/nlmod/dims/base.py
index c1d3279d..5c32aacd 100644
--- a/nlmod/dims/base.py
+++ b/nlmod/dims/base.py
@@ -8,7 +8,7 @@
 
 from .. import util
 from ..epsg28992 import EPSG_28992
-from . import resample, grid
+from . import grid, resample
 from .layers import fill_nan_top_botm_kh_kv
 
 logger = logging.getLogger(__name__)
diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 74da2dfa..4e762276 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -16,13 +16,13 @@
 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
 from flopy.utils.gridintersect import GridIntersect
 from packaging import version
 from scipy.interpolate import griddata
-from affine import Affine
 from shapely.affinity import affine_transform
 from shapely.geometry import Point, Polygon
 from tqdm import tqdm
diff --git a/nlmod/dims/rdp.py b/nlmod/dims/rdp.py
index df996789..bd688d38 100644
--- a/nlmod/dims/rdp.py
+++ b/nlmod/dims/rdp.py
@@ -1,8 +1,8 @@
 """
-rdp 
-~~~ 
+rdp
+~~~
 Python implementation of the Ramer-Douglas-Peucker algorithm.
-:copyright: 2014-2016 Fabian Hirschmann <fabian@hirschmann.email> 
+:copyright: 2014-2016 Fabian Hirschmann <fabian@hirschmann.email>
 :license: MIT, see LICENSE.txt for more details.
 """
 
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index ac8acaee..dcd5ed01 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -6,6 +6,7 @@
 import xarray as xr
 from scipy.interpolate import griddata
 from scipy.spatial import cKDTree
+
 from ..util import get_da_from_da_ds
 
 logger = logging.getLogger(__name__)
diff --git a/nlmod/gis.py b/nlmod/gis.py
index d70d5e17..3ef06361 100644
--- a/nlmod/gis.py
+++ b/nlmod/gis.py
@@ -4,7 +4,7 @@
 import geopandas as gpd
 import numpy as np
 
-from .dims.grid import polygons_from_model_ds, get_affine_mod_to_world
+from .dims.grid import get_affine_mod_to_world, polygons_from_model_ds
 from .dims.layers import calculate_thickness
 
 logger = logging.getLogger(__name__)
diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py
index 1471a768..2afc8586 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -6,7 +6,7 @@
 import xarray as xr
 
 from ..dims import grid
-from ..dims.grid import get_delr, get_delc
+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
diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py
index e4fd590c..ac32f29e 100644
--- a/nlmod/gwf/output.py
+++ b/nlmod/gwf/output.py
@@ -7,7 +7,7 @@
 import xarray as xr
 from shapely.geometry import Point
 
-from ..dims.grid import modelgrid_from_ds, get_affine_world_to_mod
+from ..dims.grid import get_affine_world_to_mod, modelgrid_from_ds
 from ..mfoutput.mfoutput import (
     _get_budget_da,
     _get_flopy_data_object,
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index c3e6e7dc..b912490e 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -11,15 +11,15 @@
 from tqdm import tqdm
 
 from ..cache import cache_pickle
-from ..util import extent_to_polygon
 from ..dims.grid import (
     gdf_to_grid,
-    get_extent_polygon,
-    get_delr,
     get_delc,
+    get_delr,
+    get_extent_polygon,
 )
 from ..dims.layers import get_idomain
 from ..read import bgt, waterboard
+from ..util import extent_to_polygon
 
 logger = logging.getLogger(__name__)
 
diff --git a/nlmod/mfoutput/mfoutput.py b/nlmod/mfoutput/mfoutput.py
index d6600607..f1d06eb9 100644
--- a/nlmod/mfoutput/mfoutput.py
+++ b/nlmod/mfoutput/mfoutput.py
@@ -7,9 +7,9 @@
 import xarray as xr
 
 from ..dims.grid import (
+    get_affine_mod_to_world,
     get_dims_coords_from_modelgrid,
     modelgrid_from_ds,
-    get_affine_mod_to_world,
 )
 from ..dims.time import ds_time_idx
 from .binaryfile import _get_binary_budget_data, _get_binary_head_data
diff --git a/nlmod/plot/dcs.py b/nlmod/plot/dcs.py
index 9430e90e..0a879445 100644
--- a/nlmod/plot/dcs.py
+++ b/nlmod/plot/dcs.py
@@ -14,7 +14,7 @@
 from shapely.affinity import affine_transform
 from shapely.geometry import LineString, MultiLineString, Point, Polygon
 
-from ..dims.grid import modelgrid_from_ds, 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__)
diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py
index dc78bcc0..fd3c39f9 100644
--- a/nlmod/plot/plot.py
+++ b/nlmod/plot/plot.py
@@ -13,10 +13,10 @@
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 
 from ..dims.grid import (
-    modelgrid_from_ds,
     get_affine_mod_to_world,
     get_extent,
     get_extent_gdf,
+    modelgrid_from_ds,
 )
 from ..read import geotop, rws
 from .dcs import DatasetCrossSection
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index 86455489..9344324b 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -7,8 +7,8 @@
 from hydropandas.io import knmi as hpd_knmi
 
 from .. import cache, util
-from ..dims.layers import get_first_active_layer
 from ..dims.grid import get_affine_mod_to_world
+from ..dims.layers import get_first_active_layer
 
 logger = logging.getLogger(__name__)
 
diff --git a/nlmod/util.py b/nlmod/util.py
index ac6ba69b..ed052967 100644
--- a/nlmod/util.py
+++ b/nlmod/util.py
@@ -13,7 +13,7 @@
 from colorama import Back, Fore, Style
 from flopy.utils import get_modflow
 from flopy.utils.get_modflow import flopy_appdata_path, get_release
-from shapely.geometry import box, Polygon
+from shapely.geometry import Polygon, box
 
 logger = logging.getLogger(__name__)
 
@@ -595,7 +595,6 @@ def compare_model_extents(extent1, extent2):
 def extent_to_polygon(extent):
     """Generate a shapely Polygon from an extent ([xmin, xmax, ymin, ymax])
 
-
     Parameters
     ----------
     extent : tuple, list or array
@@ -630,7 +629,6 @@ def extent_to_gdf(extent, crs="EPSG:28992"):
     gdf_extent : geopandas.GeoDataFrame
         geodataframe with extent.
     """
-
     geom_extent = extent_to_polygon(extent)
     gdf_extent = gpd.GeoDataFrame(geometry=[geom_extent], crs=crs)
 

From d2910f39de03411686913c60bccc134a736d4beb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 2 Jul 2024 16:24:58 +0200
Subject: [PATCH 05/10] codacy stuff and import order of dims

---
 nlmod/dims/__init__.py | 4 ++--
 nlmod/dims/grid.py     | 1 +
 nlmod/dims/resample.py | 1 +
 3 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/nlmod/dims/__init__.py b/nlmod/dims/__init__.py
index 18c147c6..14c2d185 100644
--- a/nlmod/dims/__init__.py
+++ b/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 *
diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 4e762276..440df8d9 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -6,6 +6,7 @@
 -   fill, interpolate and resample grid data
 """
 
+# ruff: noqa: E402
 import logging
 import os
 import warnings
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index dcd5ed01..8952f496 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -1,3 +1,4 @@
+# ruff: noqa: E402
 import logging
 import numbers
 

From f46101c2f9c4a31fb458e553b8e550a70c08c442 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 2 Jul 2024 16:37:47 +0200
Subject: [PATCH 06/10] Small improvements to notebooks

---
 docs/examples/01_basic_model.ipynb           | 1 -
 docs/examples/04_modifying_layermodels.ipynb | 4 ++--
 docs/examples/11_grid_rotation.ipynb         | 6 +++---
 nlmod/dims/grid.py                           | 1 -
 nlmod/dims/resample.py                       | 1 -
 5 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/docs/examples/01_basic_model.ipynb b/docs/examples/01_basic_model.ipynb
index ed358fb3..a8bd308e 100644
--- a/docs/examples/01_basic_model.ipynb
+++ b/docs/examples/01_basic_model.ipynb
@@ -18,7 +18,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
     "import nlmod"
    ]
   },
diff --git a/docs/examples/04_modifying_layermodels.ipynb b/docs/examples/04_modifying_layermodels.ipynb
index f7fe8687..f8ef9abd 100644
--- a/docs/examples/04_modifying_layermodels.ipynb
+++ b/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\")"
    ]
diff --git a/docs/examples/11_grid_rotation.ipynb b/docs/examples/11_grid_rotation.ipynb
index 61f3d8d0..8b2dc97f 100644
--- a/docs/examples/11_grid_rotation.ipynb
+++ b/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\")"
    ]
   },
diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 440df8d9..4e762276 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -6,7 +6,6 @@
 -   fill, interpolate and resample grid data
 """
 
-# ruff: noqa: E402
 import logging
 import os
 import warnings
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index 8952f496..dcd5ed01 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -1,4 +1,3 @@
-# ruff: noqa: E402
 import logging
 import numbers
 

From 2c228b689f9da6b56ea644e4778c9bb7bdae1429 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 2 Jul 2024 17:06:31 +0200
Subject: [PATCH 07/10] Add codacy rules and replace one use of
 polygon_from_extent

---
 .prospector.yaml   | 2 ++
 nlmod/dims/grid.py | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/.prospector.yaml b/.prospector.yaml
index 57c8ec7c..4656c947 100644
--- a/.prospector.yaml
+++ b/.prospector.yaml
@@ -29,6 +29,8 @@ pylint:
     - too-many-branches
     - too-many-statements
     - logging-fstring-interpolation
+    - import-outside-toplevel
+    - implicit-str-concat
 
 mccabe:
   disable:
diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 4e762276..2b8a3625 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -1956,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

From 0224b305187f317d108cdad449c44944d2fb061b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 2 Jul 2024 22:26:37 +0200
Subject: [PATCH 08/10] Update resample.py

---
 nlmod/dims/resample.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index dcd5ed01..dc5594a5 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -615,7 +615,7 @@ def affine_transform_gdf(gdf, affine):
 def get_extent(ds, rotated=True):
     """Get the model extent, corrected for angrot if necessary."""
     logger.warning(
-        "nlmod.resample.get_extent is deprecated. " "Use nlmod.grid.get_extent instead"
+        "nlmod.resample.get_extent is deprecated. Use nlmod.grid.get_extent instead"
     )
     from .grid import get_extent
 
@@ -647,8 +647,8 @@ def get_affine_world_to_mod(ds):
 def get_affine(ds, sx=None, sy=None):
     """Get the affine-transformation, from pixel to real-world coordinates."""
     logger.warning(
-        "nlmod.resample.get_affine is deprecated. " "Use nlmod.grid.get_affine instead"
+        "nlmod.resample.get_affine is deprecated. Use nlmod.grid.get_affine instead"
     )
     from .grid import get_affine
 
-    return get_affine(ds)
+    return get_affine(ds, sx=sx, sy=sy)

From fe8d2cae116337f7655880a9763e5ee2e91fac31 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= <d.brakenhoff@artesia-water.nl>
Date: Wed, 3 Jul 2024 14:09:39 +0200
Subject: [PATCH 09/10] docstring comments @OnnoEbbens

---
 nlmod/plot/plot.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py
index fd3c39f9..c0ab4e7d 100644
--- a/nlmod/plot/plot.py
+++ b/nlmod/plot/plot.py
@@ -59,9 +59,14 @@ def modelextent(ds, dx=None, ax=None, rotated=True, **kwargs):
     ----------
     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
+        Plot the model extent in real-world coordinates for rotated grids. Default is
+        True. Set to False to plot model extent in local coordinates.
     **kwargs
         Additional keyword arguments to pass to the boundary plot.
 

From 0366f8d6a5c13c05177f3dcf6bb10e71b280d627 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Wed, 3 Jul 2024 15:37:00 +0200
Subject: [PATCH 10/10] Set rotated to False in plot-methods

---
 nlmod/plot/plot.py | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/nlmod/plot/plot.py b/nlmod/plot/plot.py
index c0ab4e7d..30407da2 100644
--- a/nlmod/plot/plot.py
+++ b/nlmod/plot/plot.py
@@ -1,3 +1,4 @@
+import logging
 import warnings
 from functools import partial
 
@@ -28,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)
@@ -52,7 +55,7 @@ def modelgrid(ds, ax=None, **kwargs):
     return ax
 
 
-def modelextent(ds, dx=None, ax=None, rotated=True, **kwargs):
+def modelextent(ds, dx=None, ax=None, rotated=False, **kwargs):
     """Plot model extent.
 
     Parameters
@@ -65,8 +68,8 @@ def modelextent(ds, dx=None, ax=None, rotated=True, **kwargs):
         The axes object to plot on. If not provided, a new figure and axes will be
         created.
     rotated : bool, optional
-        Plot the model extent in real-world coordinates for rotated grids. Default is
-        True. Set to False to plot model extent in local coordinates.
+        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.
 
@@ -212,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).
 
@@ -407,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
@@ -458,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,
@@ -519,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:
@@ -578,7 +585,7 @@ def animate_map(
     colorbar=True,
     colorbar_label="",
     plot_grid=True,
-    rotated=True,
+    rotated=False,
     background=False,
     figsize=None,
     ax=None,
@@ -625,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