Solve issue #442

nlmod/dims/grid.py

@@ -36,65 +36,19 @@
     remove_inactive_layers,
 )
 from .rdp import rdp
-from .shared import GridTypeDims, get_area, get_delc, get_delr  # noqa: F401
+from .shared import (
+    GridTypeDims,
+    get_area,
+    get_delc,
+    get_delr,
+    is_structured,
+    is_vertex,
+    is_layered,
+)  # noqa: F401
 
 logger = logging.getLogger(__name__)
 
 
-def is_structured(ds):
-    """Check if a dataset is structured.
-
-    Parameters
-    ----------
-    ds : xr.Dataset or xr.Dataarray
-        dataset or dataarray
-
-    Returns
-    -------
-    bool
-        True if the dataset is structured.
-    """
-    return GridTypeDims.parse_dims(ds) in (
-        GridTypeDims.STRUCTURED,
-        GridTypeDims.STRUCTURED_LAYERED,
-    )
-
-
-def is_vertex(ds):
-    """Check if a dataset is vertex.
-
-    Parameters
-    ----------
-    ds : xr.Dataset or xr.Dataarray
-        dataset or dataarray
-
-    Returns
-    -------
-    bool
-        True if the dataset is structured.
-    """
-    return GridTypeDims.parse_dims(ds) in (
-        GridTypeDims.VERTEX,
-        GridTypeDims.VERTEX_LAYERED,
-    )
-
-
-def is_layered(ds):
-    """Check if a dataset is layered.
-
-    Parameters
-    ----------
-    ds : xr.Dataset or xr.Dataarray
-        dataset or dataarray
-
-    Returns
-    -------
-    bool
-        True if the dataset is layered.
-    """
-    return "layer" in ds.dims
-
-
 def snap_extent(extent, delr, delc):
     """Snap the extent in such a way that an integer number of columns and rows fit in
     the extent. The new extent is always equal to, or bigger than the original extent.

nlmod/dims/shared.py

@@ -38,6 +38,60 @@ def parse_dims(cls, ds):
         return cls(ds.dims)
 
 
+def is_structured(ds):
+    """Check if a dataset is structured.
+
+    Parameters
+    ----------
+    ds : xr.Dataset or xr.Dataarray
+        dataset or dataarray
+
+    Returns
+    -------
+    bool
+        True if the dataset is structured.
+    """
+    return GridTypeDims.parse_dims(ds) in (
+        GridTypeDims.STRUCTURED,
+        GridTypeDims.STRUCTURED_LAYERED,
+    )
+
+
+def is_vertex(ds):
+    """Check if a dataset is vertex.
+
+    Parameters
+    ----------
+    ds : xr.Dataset or xr.Dataarray
+        dataset or dataarray
+
+    Returns
+    -------
+    bool
+        True if the dataset is structured.
+    """
+    return GridTypeDims.parse_dims(ds) in (
+        GridTypeDims.VERTEX,
+        GridTypeDims.VERTEX_LAYERED,
+    )
+
+
+def is_layered(ds):
+    """Check if a dataset is layered.
+
+    Parameters
+    ----------
+    ds : xr.Dataset or xr.Dataarray
+        dataset or dataarray
+
+    Returns
+    -------
+    bool
+        True if the dataset is layered.
+    """
+    return "layer" in ds.dims
+
+
 def get_delr(ds):
     """
     Get the distance along rows (delr) from the x-coordinate of a structured model ds.
@@ -53,8 +107,12 @@ def get_delr(ds):
         The cell-size along rows (of length ncol).
 
     """
-    assert ds.gridtype == "structured"
-    x = (ds.x - ds.extent[0]).values
+    assert is_structured(ds)
+    if "extent" in ds.attrs:
+        west_model = ds.extent[0]
+    else:
+        west_model = float(ds.x[0] - (ds.x[1] - ds.x[0]) / 2)
+    x = (ds.x - west_model).values
     delr = _get_delta_along_axis(x)
     return delr
 
@@ -75,8 +133,12 @@ def get_delc(ds):
         The cell-size along columns (of length nrow).
 
     """
-    assert ds.gridtype == "structured"
-    y = (ds.extent[3] - ds.y).values
+    assert is_structured(ds)
+    if "extent" in ds.attrs:
+        north_model = ds.extent[3]
+    else:
+        north_model = float(ds.y[0] + (ds.y[0] - ds.y[1]) / 2)
+    y = (north_model - ds.y).values
     delc = _get_delta_along_axis(y)
     return delc
 

nlmod/plot/dcs.py

@@ -1,6 +1,4 @@
 import logging
-from functools import partial
-
 import flopy
 import geopandas as gpd
 import matplotlib
@@ -14,7 +12,7 @@
 from shapely.affinity import affine_transform
 from shapely.geometry import LineString, MultiLineString, Point, Polygon
 
-from ..dims.grid import get_affine_world_to_mod, modelgrid_from_ds
+from ..dims.grid import get_affine_world_to_mod, modelgrid_from_ds, get_delr, get_delc
 from .plotutil import get_map
 
 logger = logging.getLogger(__name__)
@@ -125,9 +123,12 @@ def get_grid_edges(self):
             x and y edges of the dataset.
         """
         x = self.ds[self.x].values
-        x = np.hstack((x[:-1] - np.diff(x) / 2, x[-2:] + np.diff(x[-3:]) / 2))
+        dx = get_delr(self.ds)
+        x = np.hstack((x - dx / 2, x[-1] + dx[-1] / 2))
+
         y = self.ds[self.y].values
-        y = np.hstack((y[:-1] - np.diff(y) / 2, y[-2:] + np.diff(y[-3:]) / 2))
+        dy = get_delc(self.ds)
+        y = np.hstack((y + dy / 2, y[-1] - dy[-1] / 2))
         return x, y
 
     def coordinates_in_dataset(self, xy):
@@ -381,7 +382,7 @@ def plot_map_cs(
 
     def iterate_active_cells(self, zcs):
         """Iterate over the cell indices of the cells in an array that are visible in the cross section and active in the model.
-        
+
         Parameters
         ----------
         zcs : np.ndarray
@@ -390,8 +391,8 @@ def iterate_active_cells(self, zcs):
         Yields
         ------
         tuple
-            i, j indices of the cells that are active and visible in cross section. 
-        
+            i, j indices of the cells that are active and visible in cross section.
+
         """
         for i in range(zcs.shape[0]):
             for j in range(zcs.shape[1]):
@@ -431,7 +432,6 @@ def _get_rect(self, i, j, hcs=None):
         rect = Rectangle(xy, width, height)
         return rect
 
-
     def array_on_cs(self, z):
         """Select cells in an array that are in the cross section.
 
@@ -473,15 +473,16 @@ def plot_array(self, z, head=None, **kwargs):
         else:
             hcs = None
 
-        patches = [self._get_rect(i,j, hcs=hcs) for i, j in self.iterate_active_cells(zcs)]
+        patches = [
+            self._get_rect(i, j, hcs=hcs) for i, j in self.iterate_active_cells(zcs)
+        ]
         array = [zcs[i, j] for i, j in self.iterate_active_cells(zcs)]
 
         patch_collection = PatchCollection(patches, **kwargs)
         patch_collection.set_array(np.array(array))
         self.ax.add_collection(patch_collection)
         return patch_collection
 
-
     def plot_surface(self, z, **kwargs):
         if isinstance(z, xr.DataArray):
             z = z.data
@@ -577,7 +578,12 @@ def animate(
         iper = 0
         if head is not None:
             plot_head = head
-            self.pc = self.plot_array(da[iper].squeeze(), cmap=cmap, norm=norm, head=head.values[iper].squeeze())
+            self.pc = self.plot_array(
+                da[iper].squeeze(),
+                cmap=cmap,
+                norm=norm,
+                head=head.values[iper].squeeze(),
+            )
         else:
             self.pc = self.plot_array(da[iper].squeeze(), cmap=cmap, norm=norm)
             plot_head = None
@@ -607,9 +613,12 @@ def update(iper):
             if plot_head is not None:
                 # create new patches
                 hcs = self.array_on_cs(plot_head[iper].squeeze())
-                patches = [self._get_rect(i,j, hcs=hcs) for i, j in self.iterate_active_cells(zcs)]
+                patches = [
+                    self._get_rect(i, j, hcs=hcs)
+                    for i, j in self.iterate_active_cells(zcs)
+                ]
                 array = [zcs[i, j] for i, j in self.iterate_active_cells(zcs)]
-                self.pc.remove() # remove previous patches
+                self.pc.remove()  # remove previous patches
                 self.pc = PatchCollection(patches, cmap=cmap, norm=norm)
                 self.pc.set_array(np.array(array))
                 self.ax.add_collection(self.pc)
@@ -618,7 +627,6 @@ def update(iper):
                 array = [zcs[i, j] for i, j in self.iterate_active_cells(zcs)]
                 self.pc.set_array(np.array(array))
 
-
             # update title
             if da.time.dtype.kind == "M":
                 t = pd.Timestamp(da.time.values[iper]).strftime(date_fmt)

tests/test_011_dcs.py

@@ -7,6 +7,18 @@
 def test_dcs_structured():
     ds = util.get_ds_structured()
     line = [(0, 0), (1000, 1000)]
+    plt.subplots()
+    dcs = nlmod.plot.DatasetCrossSection(ds, line)
+    dcs.plot_layers()
+    dcs.label_layers()
+    dcs.plot_array(ds["kh"], alpha=0.5)
+    dcs.plot_grid()
+
+
+def test_dcs_cs_model():
+    ds = util.get_ds_structured([0, 1000, 0, 100], delr=100)
+    line = [(0, 0), (1000, 100)]
+    plt.subplots()
     dcs = nlmod.plot.DatasetCrossSection(ds, line)
     dcs.plot_layers()
     dcs.label_layers()
@@ -17,6 +29,7 @@ def test_dcs_structured():
 def test_dcs_vertex():
     ds = util.get_ds_vertex()
     line = [(0, 0), (1000, 1000)]
+    plt.subplots()
     dcs = nlmod.plot.DatasetCrossSection(ds, line)
     dcs.plot_layers()
     dcs.label_layers()