fix(gridintersect): gridintersect does not work for rotated vertex grids

autotest/test_gridintersect.py

@@ -135,9 +135,7 @@ def test_rect_grid_3d_point_outside():
 
 @requires_pkg("shapely")
 def test_rect_grid_3d_point_inside():
-    botm = np.concatenate([np.ones(4), 0.5 * np.ones(4), np.zeros(4)]).reshape(
-        3, 2, 2
-    )
+    botm = np.concatenate([np.ones(4), 0.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(2.0, 2.0, 0.2))
@@ -456,9 +454,7 @@ def test_rect_grid_linestring_in_and_out_of_cell():
 def test_rect_grid_linestring_in_and_out_of_cell2():
     gr = get_rect_grid()
     ix = GridIntersect(gr, method="structured")
-    result = ix.intersect(
-        LineString([(5, 15), (5.0, 9), (15.0, 5.0), (5.0, 1.0)])
-    )
+    result = ix.intersect(LineString([(5, 15), (5.0, 9), (15.0, 5.0), (5.0, 1.0)]))
     assert len(result) == 3
     # assert result.cellids[0] == (1, 0)
     # assert result.cellids[1] == (1, 1)
@@ -755,9 +751,7 @@ def test_rect_grid_polygon_outside():
 def test_rect_grid_polygon_in_2cells():
     gr = get_rect_grid()
     ix = GridIntersect(gr, method="structured")
-    result = ix.intersect(
-        Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.0), (2.5, 15.0)])
-    )
+    result = ix.intersect(Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.0), (2.5, 15.0)]))
     assert len(result) == 2
     assert result.areas.sum() == 50.0
 
@@ -794,9 +788,7 @@ def test_rect_grid_polygon_running_along_boundary():
 def test_rect_grid_polygon_on_inner_boundary():
     gr = get_rect_grid()
     ix = GridIntersect(gr, method="structured")
-    result = ix.intersect(
-        Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)])
-    )
+    result = ix.intersect(Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]))
     assert len(result) == 2
     assert result.areas.sum() == 50.0
 
@@ -977,9 +969,7 @@ def test_rect_grid_polygon_outside_shapely(rtree):
 def test_rect_grid_polygon_in_2cells_shapely(rtree):
     gr = get_rect_grid()
     ix = GridIntersect(gr, method="vertex", rtree=rtree)
-    result = ix.intersect(
-        Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.0), (2.5, 15.0)])
-    )
+    result = ix.intersect(Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.0), (2.5, 15.0)]))
     assert len(result) == 2
     assert result.areas.sum() == 50.0
 
@@ -1000,9 +990,7 @@ def test_rect_grid_polygon_on_outer_boundary_shapely(rtree):
 def test_rect_grid_polygon_on_inner_boundary_shapely(rtree):
     gr = get_rect_grid()
     ix = GridIntersect(gr, method="vertex", rtree=rtree)
-    result = ix.intersect(
-        Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)])
-    )
+    result = ix.intersect(Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]))
     assert len(result) == 2
     assert result.areas.sum() == 50.0
 
@@ -1085,9 +1073,7 @@ def test_tri_grid_polygon_in_2cells(rtree):
     if gr == -1:
         return
     ix = GridIntersect(gr, rtree=rtree)
-    result = ix.intersect(
-        Polygon([(2.5, 5.0), (5.0, 5.0), (5.0, 15.0), (2.5, 15.0)])
-    )
+    result = ix.intersect(Polygon([(2.5, 5.0), (5.0, 5.0), (5.0, 15.0), (2.5, 15.0)]))
     assert len(result) == 2
     assert result.areas.sum() == 25.0
 
@@ -1112,9 +1098,7 @@ def test_tri_grid_polygon_on_inner_boundary(rtree):
     if gr == -1:
         return
     ix = GridIntersect(gr, rtree=rtree)
-    result = ix.intersect(
-        Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)])
-    )
+    result = ix.intersect(Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]))
     assert len(result) == 4
     assert result.areas.sum() == 50.0
 
@@ -1205,6 +1189,10 @@ def test_point_offset_rot_structured_grid():
     ix = GridIntersect(sgr, method="structured")
     result = ix.intersect(p)
     assert len(result) == 1
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="structured", local=True)
+    result = ix.intersect(p)
+    assert len(result) == 0
 
 
 @requires_pkg("shapely")
@@ -1213,8 +1201,11 @@ def test_linestring_offset_rot_structured_grid():
     ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))])
     ix = GridIntersect(sgr, method="structured")
     result = ix.intersect(ls)
-    # NOTE: in shapely 2.0, this returns a Linestring with length 10^-15 in cell (0, 1)
-    assert len(result) == 2 or len(result) == 3
+    assert len(result) == 2
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="structured", local=True)
+    result = ix.intersect(ls)
+    assert len(result) == 0
 
 
 @requires_pkg("shapely")
@@ -1231,6 +1222,10 @@ def test_polygon_offset_rot_structured_grid():
     ix = GridIntersect(sgr, method="structured")
     result = ix.intersect(p)
     assert len(result) == 3
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="structured", local=True)
+    result = ix.intersect(p)
+    assert len(result) == 0
 
 
 @requires_pkg("shapely")
@@ -1241,6 +1236,10 @@ def test_point_offset_rot_structured_grid_shapely(rtree):
     ix = GridIntersect(sgr, method="vertex", rtree=rtree)
     result = ix.intersect(p)
     assert len(result) == 1
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
+    result = ix.intersect(p)
+    assert len(result) == 0
 
 
 @requires_pkg("shapely")
@@ -1251,6 +1250,10 @@ def test_linestring_offset_rot_structured_grid_shapely(rtree):
     ix = GridIntersect(sgr, method="vertex", rtree=rtree)
     result = ix.intersect(ls)
     assert len(result) == 2
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
+    result = ix.intersect(ls)
+    assert len(result) == 0
 
 
 @requires_pkg("shapely")
@@ -1268,6 +1271,59 @@ def test_polygon_offset_rot_structured_grid_shapely(rtree):
     ix = GridIntersect(sgr, method="vertex", rtree=rtree)
     result = ix.intersect(p)
     assert len(result) == 3
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
+    result = ix.intersect(p)
+    assert len(result) == 0
+
+
+@requires_pkg("shapely")
+@rtree_toggle
+def test_point_offset_rot_vertex_grid_shapely(rtree):
+    sgr = get_rect_vertex_grid(angrot=45.0, xyoffset=10.0)
+    p = Point(10.0, 10 + np.sqrt(200.0))
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree)
+    result = ix.intersect(p)
+    assert len(result) == 1
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
+    result = ix.intersect(p)
+    assert len(result) == 0
+
+
+@requires_pkg("shapely")
+@rtree_toggle
+def test_linestring_offset_rot_vertex_grid_shapely(rtree):
+    sgr = get_rect_vertex_grid(angrot=45.0, xyoffset=10.0)
+    ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))])
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree)
+    result = ix.intersect(ls)
+    assert len(result) == 2
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
+    result = ix.intersect(ls)
+    assert len(result) == 0
+
+
+@requires_pkg("shapely")
+@rtree_toggle
+def test_polygon_offset_rot_vertex_grid_shapely(rtree):
+    sgr = get_rect_vertex_grid(angrot=45.0, xyoffset=10.0)
+    p = Polygon(
+        [
+            (5, 10.0 + np.sqrt(200.0)),
+            (15, 10.0 + np.sqrt(200.0)),
+            (15, 10.0 + 1.5 * np.sqrt(200.0)),
+            (5, 10.0 + 1.5 * np.sqrt(200.0)),
+        ]
+    )
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree)
+    result = ix.intersect(p)
+    assert len(result) == 3
+    # check empty result when using local model coords
+    ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
+    result = ix.intersect(p)
+    assert len(result) == 0
 
 
 # %% test rasters
@@ -1318,9 +1374,7 @@ def test_rasters(example_data_path):
     if (np.max(data) - 2608.557) > 1e-4:
         raise AssertionError
 
-    data = rio.resample_to_grid(
-        ml.modelgrid, band=rio.bands[0], method="nearest"
-    )
+    data = rio.resample_to_grid(ml.modelgrid, band=rio.bands[0], method="nearest")
     if data.size != 5913:
         raise AssertionError
     if abs(np.min(data) - 1942.1735) > 1e-4:
@@ -1369,12 +1423,19 @@ def test_raster_sampling_methods(example_data_path):
     }
 
     for method, value in methods.items():
-        data = rio.resample_to_grid(
-            ml.modelgrid, band=rio.bands[0], method=method
-        )
+        data = rio.resample_to_grid(ml.modelgrid, band=rio.bands[0], method=method)
 
         print(data[30, 37])
         if np.abs(data[30, 37] - value) > 1e-05:
-            raise AssertionError(
-                f"{method} resampling returning incorrect values"
-            )
+            raise AssertionError(f"{method} resampling returning incorrect values")
+
+
+if __name__ == "__main__":
+    sgr = get_rect_grid(angrot=45.0, xyoffset=10.0)
+    ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))])
+    ix = GridIntersect(sgr, method="structured")
+    result = ix.intersect(ls)
+    assert len(result) == 2
+    ix = GridIntersect(sgr, method="structured", local=True)
+    result = ix.intersect(ls)
+    assert len(result) == 0

flopy/discretization/grid.py

@@ -113,26 +113,12 @@ class Grid:
         rotation angle of model grid, as it is rotated around the origin point
     angrot_radians : float
         rotation angle of model grid, in radians
-    xgrid : ndarray
-        returns numpy meshgrid of x edges in reference frame defined by
-        point_type
-    ygrid : ndarray
-        returns numpy meshgrid of y edges in reference frame defined by
-        point_type
-    zgrid : ndarray
-        returns numpy meshgrid of z edges in reference frame defined by
-        point_type
     xcenters : ndarray
         returns x coordinate of cell centers
     ycenters : ndarray
         returns y coordinate of cell centers
     ycenters : ndarray
         returns z coordinate of cell centers
-    xyzgrid : [ndarray, ndarray, ndarray]
-        returns the location of grid edges of all model cells. if the model
-        grid contains spatial reference information, the grid edges are in the
-        coordinate system provided by the spatial reference information.
-        returns a list of three ndarrays for the x, y, and z coordinates
     xyzcellcenters : [ndarray, ndarray, ndarray]
         returns the cell centers of all model cells in the model grid.  if
         the model grid contains spatial reference information, the cell centers
@@ -150,10 +136,6 @@ class Grid:
         returns the model grid lines in a list.  each line is returned as a
         list containing two tuples in the format [(x1,y1), (x2,y2)] where
         x1,y1 and x2,y2 are the endpoints of the line.
-    xyvertices : (point_type) : ndarray
-        1D array of x and y coordinates of cell vertices for whole grid
-        (single layer) in C-style (row-major) order
-        (same as np.ravel())
 
     See Also
     --------