update_data(contour_array): added routine to mask errant triangles (#…

…1562)

* update_data(contour_array): added routine to mask errant triangles

* added `tri_mask` **kwarg to contour_array to check triangle elements against grid cell nearest neighbors and mask triangles that do not conform to nearest neighbor requirements. `tri_mask` is only needed when holes are present or a modelgrid has concave boundaries. 
* added neighbors() method to Grid, StructuredGrid, and UnstructuredGrid that returns a cells nearest neighbors 
* updated UnstructuredGrid to optionally include iac, ja arrays
* updated MFModel and Modflow to provide iac and ja arrays during UnstructuredGrid construction
* added tests for neighbors routine

* Fill out docstrings

autotest/test_grid.py

@@ -339,6 +339,53 @@ def test_unstructured_xyz_intersect(example_data_path):
             raise AssertionError("Unstructured grid intersection failed")
 
 
+def test_structured_neighbors(example_data_path):
+    ws = str(example_data_path / "freyberg")
+    ml = Modflow.load("freyberg.nam", model_ws=ws)
+    modelgrid = ml.modelgrid
+    k, i, j = 0, 5, 5
+    neighbors = modelgrid.neighbors(k, i, j)
+    for neighbor in neighbors:
+        if (
+                neighbor != (k, i + 1, j)
+                and neighbor != (k, i - 1, j)
+                and neighbor != (k, i, j + 1)
+                and neighbor != (k, i, j - 1)
+        ):
+            raise AssertionError(
+                "modelgid.neighbors not returning proper values"
+            )
+
+def test_vertex_neighbors(example_data_path):
+    ws = str(example_data_path / "mf6" / "test003_gwfs_disv")
+    sim = MFSimulation.load(sim_ws=ws)
+    gwf = sim.get_model("gwf_1")
+    modelgrid = gwf.modelgrid
+    node = 63
+    neighbors = modelgrid.neighbors(node)
+    for neighbor in neighbors:
+        if (
+                neighbor != node + 1
+                and neighbor != node - 1
+                and neighbor != node + 10
+                and neighbor != node - 10
+        ):
+            raise AssertionError(
+                "modelgid.neighbors not returning proper values"
+            )
+
+
+def test_unstructured_neighbors(example_data_path):
+    ws = str(example_data_path / "mf6" / "test006_gwf3")
+    sim = MFSimulation.load(sim_ws=ws)
+    gwf = sim.get_model("gwf_1")
+    modelgrid = gwf.modelgrid
+    truth = [3, 5, 11]
+    neighbors = modelgrid.neighbors(4)
+    if not truth == neighbors:
+        raise AssertionError("modelgid.neighbors not returning proper values")
+
+
 @pytest.mark.parametrize("spc_file", ["grd.spc", "grdrot.spc"])
 def test_structured_from_gridspec(example_data_path, spc_file):
     fn = str(example_data_path / "specfile" / spc_file)
@@ -959,7 +1006,7 @@ def test_get_lni_unstructured(grid):
                     list(grid.ncpl)
                     if not isinstance(grid.ncpl, int)
                     else [grid.ncpl for _ in range(grid.nlay)]
-                )
+        )
             )
         )
         assert csum[layer] + i == nn

flopy/discretization/grid.py

@@ -3,6 +3,7 @@
 
 import numpy as np
 
+from ..plot.plotutil import UnstructuredPlotUtilities
 from ..utils import geometry
 from ..utils.gridutil import get_lni
 
@@ -184,6 +185,7 @@ def __init__(
         self._iverts = None
         self._verts = None
         self._laycbd = None
+        self._neighbors = None
 
     ###################################
     # access to basic grid properties
@@ -467,6 +469,60 @@ def xyzvertices(self):
     def cross_section_vertices(self):
         return self.xyzvertices[0], self.xyzvertices[1]
 
+    def neighbors(self, node, **kwargs):
+        """
+        Method to get nearest neighbors of a cell
+
+        Parameters
+        ----------
+        node : int
+            model grid node number
+
+        Returns
+        -------
+            list : list of cell node numbers
+        """
+        ncpl = self.ncpl
+        if isinstance(ncpl, list):
+            ncpl = self.nnodes
+
+        lay = 0
+        while node >= ncpl:
+            node -= ncpl
+            lay += 1
+
+        if self._neighbors is None:
+            neigh = {}
+            xverts, yverts = self.cross_section_vertices
+            xverts, yverts = UnstructuredPlotUtilities.irregular_shape_patch(
+                xverts, yverts
+            )
+            for nn in range(ncpl):
+                conn = []
+                verts = self.get_cell_vertices(nn)
+                for ix in range(0, len(verts)):
+                    xv0, yv0 = verts[ix - 1]
+                    xv1, yv1 = verts[ix]
+                    idx0 = np.unique(
+                        np.where((xverts == xv0) & (yverts == yv0))[0]
+                    )
+                    idx1 = np.unique(
+                        np.where((xverts == xv1) & (yverts == yv1))[0]
+                    )
+                    mask = np.isin(idx0, idx1)
+                    conn += [i for i in idx0[mask]]
+
+                conn = list(set(conn))
+                conn.pop(conn.index(nn))
+                neigh[nn] = conn
+            self._neighbors = neigh
+
+        neighbors = self._neighbors[node]
+        if lay > 0:
+            neighbors = [i + (ncpl * lay) for i in neighbors]
+
+        return neighbors
+
     def remove_confining_beds(self, array):
         """
         Method to remove confining bed layers from an array

flopy/discretization/structuredgrid.py

@@ -743,6 +743,59 @@ def map_polygons(self):
     ###############
     ### Methods ###
     ###############
+    def neighbors(self, *args, **kwargs):
+        """
+        Method to get nearest neighbors for a cell
+
+        Parameters
+        ----------
+        *args
+            lay (int), row (int), column (int)
+            or
+            node (int)
+
+        **kwargs
+            k : int
+                layer number
+            i : int
+                row number
+            j : int
+                column number
+            as_node : bool
+                flag to return neighbors as node numbers
+
+        Returns
+        -------
+            list of neighboring cells
+        """
+        nn = None
+        if kwargs:
+            if "node" in kwargs:
+                nn = kwargs.pop("node")
+            else:
+                k = kwargs.pop("k", 0)
+                i = kwargs.pop("i")
+                j = kwargs.pop("j")
+
+        if len(args) > 0:
+            if len(args) == 1:
+                nn = args[0]
+            elif len(args) == 2:
+                k = 0
+                i, j = args[0:2]
+            else:
+                k, i, j = args[0:3]
+
+        if nn is None:
+            nn = self.get_node([(k, i, j)])[0]
+
+        as_nodes = kwargs.pop("as_nodes", False)
+
+        neighbors = super().neighbors(nn)
+        if not as_nodes:
+            neighbors = self.get_lrc(neighbors)
+        return neighbors
+
     def intersect(self, x, y, z=None, local=False, forgive=False):
         """
         Get the row and column of a point with coordinates x and y

flopy/discretization/unstructuredgrid.py

@@ -48,6 +48,10 @@ class UnstructuredGrid(Grid):
         top elevations for all cells in the grid.
     botm : list or ndarray
         bottom elevations for all cells in the grid.
+    iac : list or ndarray
+        optional number of connections per node array
+    ja : list or ndarray
+        optional jagged connection array
 
     Properties
     ----------
@@ -97,6 +101,8 @@ def __init__(
         xoff=0.0,
         yoff=0.0,
         angrot=0.0,
+        iac=None,
+        ja=None,
     ):
         super().__init__(
             "unstructured",
@@ -143,6 +149,9 @@ def __init__(
                 msg = f"Length of iverts must equal ncpl ({len(iverts)} {self.ncpl})"
                 assert np.all([cpl == len(iverts) for cpl in self.ncpl]), msg
 
+        self._iac = iac
+        self._ja = ja
+
         return
 
     def set_ncpl(self, ncpl):
@@ -219,15 +228,11 @@ def verts(self):
             return np.array([list(t)[1:] for t in self._vertices], dtype=float)
 
     @property
-    def ia(self):
-        if self._ia is None:
-            self._set_unstructured_iaja()
-        return self._ia
+    def iac(self):
+        return self._iac
 
     @property
     def ja(self):
-        if self._ja is None:
-            self._set_unstructured_iaja()
         return self._ja
 
     @property
@@ -350,6 +355,25 @@ def cross_section_vertices(self):
             yv *= self.nlay
         return xv, yv
 
+    def neighbors(self, node, **kwargs):
+        """
+        Method to get a list of nearest neighbors
+
+        Parameters
+        ----------
+        node : int
+            node number
+
+        Returns
+        -------
+            list of nearest neighbors
+        """
+        nrec = self.iac[node]
+        idx0 = np.sum(self.iac[:node])
+        idx1 = idx0 + nrec
+        neighbors = self.ja[idx0:idx1][1:]
+        return list(neighbors)
+
     def cross_section_lay_ncpl_ncb(self, ncb):
         """
         Get PlotCrossSection compatible layers, ncpl, and ncb

flopy/mf6/mfmodel.py

@@ -459,6 +459,8 @@ def modelgrid(self):
                 xoff=self._modelgrid.xoffset,
                 yoff=self._modelgrid.yoffset,
                 angrot=self._modelgrid.angrot,
+                iac=dis.iac,
+                ja=dis.ja,
             )
         elif self.get_grid_type() == DiscretizationType.DISL:
             dis = self.get_package("disl")

flopy/modflow/mf.py

@@ -294,6 +294,8 @@ def modelgrid(self):
                 xoff=self._modelgrid.xoffset,
                 yoff=self._modelgrid.yoffset,
                 angrot=self._modelgrid.angrot,
+                iac=self.disu.iac,
+                ja=self.disu.ja,
             )
             print(
                 "WARNING: Model grid functionality limited for unstructured "

flopy/plot/map.py

@@ -207,6 +207,7 @@ def contour_array(self, a, masked_values=None, **kwargs):
 
         filled = kwargs.pop("filled", False)
         plot_triplot = kwargs.pop("plot_triplot", False)
+        tri_mask = kwargs.pop("tri_mask", False)
 
         # Get vertices for the selected layer
         xcentergrid = self.mg.get_xcellcenters_for_layer(self.layer)
@@ -229,13 +230,27 @@ def contour_array(self, a, masked_values=None, **kwargs):
         xcentergrid = xcentergrid.flatten()
         ycentergrid = ycentergrid.flatten()
         triang = tri.Triangulation(xcentergrid, ycentergrid)
+        analyze = tri.TriAnalyzer(triang)
+        mask = analyze.get_flat_tri_mask(rescale=False)
+
+        # mask out holes, optional???
+        if tri_mask:
+            triangles = triang.triangles
+            for i in range(2):
+                for ix, nodes in enumerate(triangles):
+                    neighbors = self.mg.neighbors(nodes[i], as_nodes=True)
+                    isin = np.isin(nodes[i + 1 :], neighbors)
+                    if not np.alltrue(isin):
+                        mask[ix] = True
 
         if ismasked is not None:
             ismasked = ismasked.flatten()
-            mask = np.any(
+            mask2 = np.any(
                 np.where(ismasked[triang.triangles], True, False), axis=1
             )
-            triang.set_mask(mask)
+            mask[mask2] = True
+
+        triang.set_mask(mask)
 
         if filled:
             contour_set = ax.tricontourf(triang, plotarray, **kwargs)