From 2fd9326534d2782428111ae47dfac1af0c21e04d Mon Sep 17 00:00:00 2001
From: mkennard-aquaveo <42077662+mkennard-aquaveo@users.noreply.github.com>
Date: Fri, 7 Aug 2020 06:44:13 -0600
Subject: [PATCH] fix(#958): fix GridIntersect._vtx_grid_to_shape_generator
 rtree for cells with more than 3 vertices (#959)

Also added a test.
---
 autotest/t065_test_gridintersect.py | 44 +++++++++++++++++++++++++++++
 flopy/utils/gridintersect.py        |  2 +-
 2 files changed, 45 insertions(+), 1 deletion(-)

diff --git a/autotest/t065_test_gridintersect.py b/autotest/t065_test_gridintersect.py
index 3f8c922550..efa4518f24 100644
--- a/autotest/t065_test_gridintersect.py
+++ b/autotest/t065_test_gridintersect.py
@@ -59,6 +59,27 @@ def get_rect_grid(angrot=0., xyoffset=0., top=None, botm=None):
     return sgr
 
 
+def get_rect_vertex_grid(angrot=0., xyoffset=0.):
+    cell2d = [[0, 5.0, 5.0, 4, 0, 1, 4, 3],
+              [1, 15.0, 5.0, 4, 1, 2, 5, 4],
+              [2, 5.0, 15.0, 4, 3, 4, 7, 6],
+              [3, 15.0, 15.0, 4, 4, 5, 8, 7]]
+    vertices = [[0, 0.0, 0.0],
+                [1, 10.0, 0.0],
+                [2, 20.0, 0.0],
+                [3, 0.0, 10.0],
+                [4, 10.0, 10.0],
+                [5, 20.0, 10.0],
+                [6, 0.0, 20.0],
+                [7, 10.0, 20.0],
+                [8, 20.0, 20.0]]
+    tgr = fgrid.VertexGrid(vertices, cell2d,
+                           botm=np.atleast_2d(np.zeros(len(cell2d))),
+                           top=np.ones(len(cell2d)), xoff=xyoffset,
+                           yoff=xyoffset, angrot=angrot)
+    return tgr
+
+
 def plot_structured_grid(sgr):
     _, ax = plt.subplots(1, 1, figsize=(8, 8))
     sgr.plot(ax=ax)
@@ -254,6 +275,29 @@ def test_rect_grid_point_on_inner_boundary_shapely(rtree=True):
     return result
 
 
+def test_rect_vertex_grid_point_in_one_cell_shapely(rtree=True):
+    # avoid test fail when shapely not available
+    try:
+        import shapely
+    except:
+        return
+    gr = get_rect_vertex_grid()
+    ix = GridIntersect(gr, method='vertex', rtree=rtree)
+    result = ix.intersect_point(Point(4., 4.))
+    assert len(result) == 1
+    assert result.cellids[0] == 0
+    result = ix.intersect_point(Point(4., 6.))
+    assert len(result) == 1
+    assert result.cellids[0] == 0
+    result = ix.intersect_point(Point(6., 6.))
+    assert len(result) == 1
+    assert result.cellids[0] == 0
+    result = ix.intersect_point(Point(6., 4.))
+    assert len(result) == 1
+    assert result.cellids[0] == 0
+    return result
+
+
 def test_rect_grid_multipoint_in_one_cell_shapely(rtree=True):
     # avoid test fail when shapely not available
     try:
diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index e2b31ad316..8dbab97891 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -227,7 +227,7 @@ def _vtx_grid_to_shape_generator(self):
         elif isinstance(self.mfgrid._cell2d, list):
             for icell in range(len(self.mfgrid._cell2d)):
                 points = []
-                for iv in self.mfgrid._cell2d[icell][-3:]:
+                for iv in self.mfgrid._cell2d[icell][4:]:
                     points.append(
                         (
                             self.mfgrid._vertices[iv][1],