fix(GridIntersect): fix crash intersecting 3D points that are outside… (

#957)

Fixes #956

autotest/t065_test_gridintersect.py

@@ -50,11 +50,11 @@ def get_tri_grid(angrot=0., xyoffset=0., triangle_exe=None):
     return tgr
 
 
-def get_rect_grid(angrot=0., xyoffset=0.):
+def get_rect_grid(angrot=0., xyoffset=0., top=None, botm=None):
     delc = 10 * np.ones(2, dtype=np.float)
     delr = 10 * np.ones(2, dtype=np.float)
     sgr = fgrid.StructuredGrid(
-        delc, delr, top=None, botm=None, xoff=xyoffset, yoff=xyoffset,
+        delc, delr, top=top, botm=botm, xoff=xyoffset, yoff=xyoffset,
         angrot=angrot)
     return sgr
 
@@ -98,6 +98,48 @@ def plot_ix_point_result(rec, ax):
 # %% test point structured
 
 
+def test_rect_grid_3d_point_outside():
+    # avoid test fail when shapely not available
+    try:
+        import shapely
+    except:
+        return
+    botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape(2, 2, 2)
+    gr = get_rect_grid(top=np.ones(4), botm=botm)
+    ix = GridIntersect(gr, method="structured")
+    result = ix.intersect_point(Point(25., 25., .5))
+    assert len(result) == 0
+    return result
+
+
+def test_rect_grid_3d_point_inside():
+    # avoid test fail when shapely not available
+    try:
+        import shapely
+    except:
+        return
+    botm = np.concatenate([np.ones(4), .5 * np.ones(4), np.zeros(4)]).reshape(3, 2, 2)
+    gr = get_rect_grid(top=np.ones(4), botm=botm)
+    ix = GridIntersect(gr, method="structured")
+    result = ix.intersect_point(Point(2., 2., .2))
+    assert result.cellids[0] == (1, 1, 0)
+    return result
+
+
+def test_rect_grid_3d_point_above():
+    # avoid test fail when shapely not available
+    try:
+        import shapely
+    except:
+        return
+    botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape(2, 2, 2)
+    gr = get_rect_grid(top=np.ones(4), botm=botm)
+    ix = GridIntersect(gr, method="structured")
+    result = ix.intersect_point(Point(2., 2., 2))
+    assert len(result) == 0
+    return result
+
+
 def test_rect_grid_point_outside():
     # avoid test fail when shapely not available
     try:

flopy/utils/gridintersect.py

@@ -624,17 +624,18 @@ def _intersect_point_structured(self, shp):
             ipos = ModflowGridIndices.find_position_in_array(Ye, ry)
 
             if jpos is not None and ipos is not None:
-                nodelist.append((ipos, jpos))
-                ixshapes.append(p)
-
-            # three dimensional point
-            if p._ndim == 3:
-                # find k
-                kpos = ModflowGridIndices.find_position_in_array(
-                    self.mfgrid.botm[:, ipos, jpos], p.z
-                )
-                if kpos is not None:
-                    nodelist.append((kpos, ipos, jpos))
+                # three dimensional point
+                if p._ndim == 3:
+                    # find k
+                    kpos = ModflowGridIndices.find_position_in_array(
+                        self.mfgrid.botm[:, ipos, jpos], p.z
+                    )
+                    if kpos is not None:
+                        nodelist.append((kpos, ipos, jpos))
+                        ixshapes.append(p)
+                else:
+                    nodelist.append((ipos, jpos))
+                    ixshapes.append(p)
 
         # remove duplicates
         tempnodes = []