update(Grid.intersect): added optional z-coordinate to intersect() method

autotest/t062_test_intersect.py

@@ -1,3 +1,4 @@
+import os
 import flopy
 import numpy as np
 import matplotlib.pyplot as plt
@@ -94,6 +95,27 @@ def get_vlist(i, j, nrow, ncol):
     return gwf
 
 
+# simple functions to load vertices and index lists for unstructured grid
+def load_verts(fname):
+    verts = np.genfromtxt(fname, dtype=[int, float, float],
+                          names=['iv', 'x', 'y'])
+    verts['iv'] -= 1  # zero based
+    return verts
+
+
+def load_iverts(fname):
+    f = open(fname, 'r')
+    iverts = []
+    xc = []
+    yc = []
+    for line in f:
+        ll = line.strip().split()
+        iverts.append([int(i) - 1 for i in ll[4:]])
+        xc.append(float(ll[1]))
+        yc.append(float(ll[2]))
+    return iverts, np.array(xc), np.array(yc)
+
+
 def test_intersection():
     ml_dis = dis_model()
     ml_disv = disv_model()
@@ -140,9 +162,11 @@ def test_intersection():
         else:
             print("In real_world coordinates:")
         try:
-            row, col = ml_dis.modelgrid.intersect(x, y, local, forgive=forgive)
+            row, col = ml_dis.modelgrid.intersect(
+                x, y, local=local, forgive=forgive
+            )
             cell2d_disv = ml_disv.modelgrid.intersect(
-                x, y, local, forgive=forgive
+                x, y, local=local, forgive=forgive
             )
         except Exception as e:
             if not forgive and any(
@@ -162,5 +186,98 @@ def test_intersection():
             assert all(np.isnan([row, col, cell2d_disv]))
 
 
+def test_structured_xyz_intersect():
+    model_ws = os.path.join(
+        "..", "examples", "data", "freyberg_multilayer_transient"
+    )
+    fname = "freyberg.nam"
+
+    ml = flopy.modflow.Modflow.load(fname, model_ws=model_ws)
+    mg = ml.modelgrid
+    top_botm = ml.modelgrid.top_botm
+    xc, yc, zc = mg.xyzcellcenters
+
+    for _ in range(10):
+        k = np.random.randint(0, mg.nlay, 1)[0]
+        i = np.random.randint(0, mg.nrow, 1)[0]
+        j = np.random.randint(0, mg.ncol, 1)[0]
+        x = xc[i, j]
+        y = yc[i, j]
+        z = zc[k, i, j]
+        k2, i2, j2 = ml.modelgrid.intersect(x, y, z)
+        if (k, i, j) != (k2, i2, j2):
+            raise AssertionError("Structured grid intersection failed")
+
+
+def test_vertex_xyz_intersect():
+    sim_ws = os.path.join("..", "examples", "data", "mf6", "test003_gwfs_disv")
+
+    sim = flopy.mf6.MFSimulation.load(sim_ws=sim_ws)
+    ml = sim.get_model(list(sim.model_names)[0])
+    mg = ml.modelgrid
+
+    xc, yc, zc = mg.xyzcellcenters
+    for _ in range(10):
+        icell = np.random.randint(0, mg.ncpl, 1)[0]
+        lay = np.random.randint(0, mg.nlay, 1)[0]
+        x = xc[icell]
+        y = yc[icell]
+        z = zc[lay, icell]
+        lay1, icell1 = mg.intersect(x, y, z)
+
+        if (lay, icell) != (lay1, icell1):
+            raise AssertionError("Vertex grid intersection failed")
+
+
+def test_unstructured_xyz_intersect():
+    ws = os.path.join("..", "examples", "data", "unstructured")
+    # usg example
+    name = os.path.join(ws, "ugrid_verts.dat")
+    verts = load_verts(name)
+
+    name = os.path.join(ws, "ugrid_iverts.dat")
+    iverts, xc, yc = load_iverts(name)
+
+    # create a 3 layer model grid
+    ncpl = np.array(3 * [len(iverts)])
+    nnodes = np.sum(ncpl)
+
+    top = np.ones((nnodes), )
+    botm = np.ones((nnodes), )
+
+    # set top and botm elevations
+    i0 = 0
+    i1 = ncpl[0]
+    elevs = [100, 0, -100, -200]
+    for ix, cpl in enumerate(ncpl):
+        top[i0:i1] *= elevs[ix]
+        botm[i0:i1] *= elevs[ix + 1]
+        i0 += cpl
+        i1 += cpl
+
+    # create the modelgrid
+    mg = flopy.discretization.UnstructuredGrid(vertices=verts,
+                                               iverts=iverts,
+                                               xcenters=xc,
+                                               ycenters=yc, top=top,
+                                               botm=botm, ncpl=ncpl)
+
+    xc, yc, zc = mg.xyzcellcenters
+    zc = zc[0].reshape(mg.nlay, mg.ncpl[0])
+    for _ in range(10):
+        icell = np.random.randint(0, mg.ncpl[0], 1)[0]
+        lay = np.random.randint(0, mg.nlay, 1)[0]
+        x = xc[icell]
+        y = yc[icell]
+        z = zc[lay, icell]
+        icell1 = mg.intersect(x, y, z)
+        icell = icell + (mg.ncpl[0] * lay)
+        if icell != icell1:
+            raise AssertionError("Unstructured grid intersection failed")
+
+
 if __name__ == "__main__":
     test_intersection()
+    test_structured_xyz_intersect()
+    test_vertex_xyz_intersect()
+    test_unstructured_xyz_intersect()

flopy/discretization/grid.py

@@ -671,6 +671,32 @@ def intersect(self, x, y, local=False, forgive=False):
         else:
             return x, y
 
+    def _warn_intersect(self, module, lineno):
+        """
+        Warning for modelgrid intersect() interface change.
+
+        Should be be removed after a couple of releases. Added in 3.3.5
+
+        Parameters
+        ----------
+        module : str
+            module name path
+        lineno : int
+            line number where warning is called from
+
+        Returns
+        -------
+            None
+        """
+        module = os.path.split(module)[-1]
+        warning = (
+            "The interface 'intersect(self, x, y, local=False, "
+            "forgive=False)' has been deprecated. Use the "
+            "intersect(self, x, y, z=None, local=False, "
+            "forgive=False) interface instead."
+        )
+        warnings.warn_explicit(warning, UserWarning, module, lineno)
+
     def set_coord_info(
         self,
         xoff=None,

flopy/discretization/structuredgrid.py

@@ -1,5 +1,6 @@
 import copy
 import os.path
+import inspect
 
 import numpy as np
 from .grid import Grid, CachedData
@@ -741,7 +742,7 @@ def map_polygons(self):
     ###############
     ### Methods ###
     ###############
-    def intersect(self, x, y, local=False, forgive=False):
+    def intersect(self, x, y, z=None, local=False, forgive=False):
         """
         Get the row and column of a point with coordinates x and y
 
@@ -754,6 +755,9 @@ def intersect(self, x, y, local=False, forgive=False):
             The x-coordinate of the requested point
         y : float
             The y-coordinate of the requested point
+        z : float
+            Optional z-coordinate of the requested point (will return layer,
+            row, column) if supplied
         local: bool (optional)
             If True, x and y are in local coordinates (defaults to False)
         forgive: bool (optional)
@@ -768,6 +772,10 @@ def intersect(self, x, y, local=False, forgive=False):
             The column number
 
         """
+        # trigger interface change warning
+        frame_info = inspect.getframeinfo(inspect.currentframe())
+        self._warn_intersect(frame_info.filename, frame_info.lineno)
+
         # transform x and y to local coordinates
         x, y = super().intersect(x, y, local, forgive)
 
@@ -797,7 +805,28 @@ def intersect(self, x, y, local=False, forgive=False):
             row = np.where(ycomp)[0][-1]
         if np.any(np.isnan([row, col])):
             row = col = np.nan
-        return row, col
+            if z is not None:
+                return None, row, col
+
+        if z is None:
+            return row, col
+
+        lay = np.nan
+        for layer in range(self.__nlay):
+            if (
+                self.top_botm[layer, row, col]
+                >= z
+                >= self.top_botm[layer + 1, row, col]
+            ):
+                lay = layer
+                break
+
+        if np.any(np.isnan([lay, row, col])):
+            lay = row = col = np.nan
+            if not forgive:
+                raise Exception("point given is outside the model area")
+
+        return lay, row, col
 
     def _cell_vert_list(self, i, j):
         """Get vertices for a single cell or sequence of i, j locations."""

flopy/discretization/unstructuredgrid.py

@@ -1,7 +1,12 @@
 import os
 import copy
 import numpy as np
+import inspect
+
+from matplotlib.path import Path
+
 from .grid import Grid, CachedData
+from ..utils.geometry import is_clockwise
 
 
 class UnstructuredGrid(Grid):
@@ -506,9 +511,77 @@ def map_polygons(self):
 
         return copy.copy(self._polygons)
 
-    def intersect(self, x, y, local=False, forgive=False):
-        x, y = super().intersect(x, y, local, forgive)
-        raise Exception("Not implemented yet")
+    def intersect(self, x, y, z=None, local=False, forgive=False):
+        """
+        Get the CELL2D number of a point with coordinates x and y
+
+        When the point is on the edge of two cells, the cell with the lowest
+        CELL2D number is returned.
+
+        Parameters
+        ----------
+        x : float
+            The x-coordinate of the requested point
+        y : float
+            The y-coordinate of the requested point
+        z : float, None
+            optional, z-coordiante of the requested point
+        local: bool (optional)
+            If True, x and y are in local coordinates (defaults to False)
+        forgive: bool (optional)
+            Forgive x,y arguments that fall outside the model grid and
+            return NaNs instead (defaults to False - will throw exception)
+
+        Returns
+        -------
+        icell2d : int
+            The CELL2D number
+
+        """
+        frame_info = inspect.getframeinfo(inspect.currentframe())
+        self._warn_intersect(frame_info.filename, frame_info.lineno)
+
+        if local:
+            # transform x and y to real-world coordinates
+            x, y = super().get_coords(x, y)
+        xv, yv, zv = self.xyzvertices
+
+        if self.grid_varies_by_layer:
+            ncpl = self.nnodes
+        else:
+            ncpl = self.ncpl[0]
+
+        for icell2d in range(ncpl):
+            xa = np.array(xv[icell2d])
+            ya = np.array(yv[icell2d])
+            # x and y at least have to be within the bounding box of the cell
+            if (
+                np.any(x <= xa)
+                and np.any(x >= xa)
+                and np.any(y <= ya)
+                and np.any(y >= ya)
+            ):
+                if is_clockwise(xa, ya):
+                    radius = -1e-9
+                else:
+                    radius = 1e-9
+                path = Path(np.stack((xa, ya)).transpose())
+                # use a small radius, so that the edge of the cell is included
+                if path.contains_point((x, y), radius=radius):
+                    if z is None:
+                        return icell2d
+
+                    for lay in range(self.nlay):
+                        if lay != 0 and not self.grid_varies_by_layer:
+                            icell2d += self.ncpl[lay - 1]
+                        if zv[0, icell2d] >= z >= zv[1, icell2d]:
+                            return icell2d
+
+        if forgive:
+            icell2d = np.nan
+            return icell2d
+
+        raise Exception("point given is outside of the model area")
 
     @property
     def top_botm(self):