diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py index ff61b229d7..983e3771c4 100644 --- a/autotest/test_gridintersect.py +++ b/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 diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 6d845e8cd6..b02e3a9271 100644 --- a/flopy/discretization/grid.py +++ b/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 -------- diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 465ca23229..8b3c5cda2f 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -142,7 +142,7 @@ class GridIntersect: structured routines, especially for larger grids. """ - def __init__(self, mfgrid, method=None, rtree=True): + def __init__(self, mfgrid, method=None, rtree=True, local=False): """Intersect shapes (Point, Linestring, Polygon) with a modflow grid. Parameters @@ -150,16 +150,21 @@ def __init__(self, mfgrid, method=None, rtree=True): mfgrid : flopy modflowgrid MODFLOW grid as implemented in flopy method : str, optional - default is None, which determines intersection method based on - the grid type. Options are either 'vertex' which uses shapely - interesection operations or 'structured' which uses optimized - methods that only work for structured grids + Options are either 'vertex' which uses shapely interesection operations + or 'structured' which uses optimized methods that only work for structured + grids. The default is None, which determines intersection method based on + the grid type. rtree : bool, optional whether to build an STR-Tree, default is True. If False no STR-tree is built, but intersects will loop through all model gridcells (which is generally slower). Only read when `method='vertex'`. + local : bool, optional + use local model coordinates from model grid to build grid geometries, + default is False and uses real-world coordinates (with offset and rotation), + if specified. """ self.mfgrid = mfgrid + self.local = local if method is None: # determine method from grid_type self.method = self.mfgrid.grid_type @@ -356,8 +361,11 @@ def _rect_grid_to_geoms_cellids(self): ncol = self.mfgrid.ncol ncells = nrow * ncol cellids = np.arange(ncells) - xvertices = self.mfgrid.xvertices - yvertices = self.mfgrid.yvertices + if self.local: + xvertices, yvertices = np.meshgrid(*self.mfgrid.xyedges) + else: + xvertices = self.mfgrid.xvertices + yvertices = self.mfgrid.yvertices # arrays of coordinates for rectangle cells I, J = np.ogrid[0:nrow, 0:ncol] @@ -434,12 +442,17 @@ def _vtx_grid_to_geoms_cellids(self): for i in range(self.mfgrid._cell2d["ncvert"][icell]) ] for iv in self.mfgrid._cell2d[icverts][icell]: - points.append( - ( + if self.local: + xy = ( self.mfgrid._vertices.xv[iv], self.mfgrid._vertices.yv[iv], ) - ) + else: + xy = ( + self.mfgrid.verts[iv, 0], + self.mfgrid.verts[iv, 1], + ) + points.append(xy) # close the polygon, if necessary if points[0] != points[-1]: points.append(points[0]) @@ -450,12 +463,17 @@ def _vtx_grid_to_geoms_cellids(self): for icell in range(len(self.mfgrid._cell2d)): points = [] for iv in self.mfgrid._cell2d[icell][4:]: - points.append( - ( + if self.local: + xy = ( self.mfgrid._vertices[iv][1], self.mfgrid._vertices[iv][2], ) - ) + else: + xy = ( + self.mfgrid.verts[iv, 0], + self.mfgrid.verts[iv, 1], + ) + points.append(xy) # close the polygon, if necessary if points[0] != points[-1]: points.append(points[0]) @@ -1207,7 +1225,7 @@ def _intersect_point_structured(self, shp, return_all_intersections=False): self.mfgrid.angrot != 0.0 or self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 - ): + ) and not self.local: rx, ry = transform( p.x, p.y, @@ -1347,11 +1365,13 @@ def _intersect_linestring_structured( pl = shapely_geo.box(xmin, ymin, xmax, ymax) # rotate and translate linestring to local coords - if self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0: + if ( + self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 + ) and not self.local: shp = affinity_loc.translate( shp, xoff=-self.mfgrid.xoffset, yoff=-self.mfgrid.yoffset ) - if self.mfgrid.angrot != 0.0: + if self.mfgrid.angrot != 0.0 and not self.local: shp = affinity_loc.rotate( shp, -self.mfgrid.angrot, origin=(0.0, 0.0) ) @@ -1380,7 +1400,7 @@ def _intersect_linestring_structured( self.mfgrid.angrot != 0.0 or self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 - ): + ) and not self.local: v_realworld = [] for pt in v: pt = np.array(pt) @@ -1424,7 +1444,7 @@ def _intersect_linestring_structured( self.mfgrid.angrot != 0.0 or self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 - ): + ) and not self.local: v_realworld = [] for pt in vertices: pt = np.array(pt) @@ -1546,7 +1566,7 @@ def _get_nodes_intersecting_linestring( self.mfgrid.angrot != 0.0 or self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 - ): + ) and not self.local: x0, y0 = transform( [x[0]], [y[0]], @@ -1913,11 +1933,13 @@ def _intersect_polygon_structured( ixshapes = [] # transform polygon to local grid coordinates - if self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0: + if ( + self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 + ) and not self.local: shp = affinity_loc.translate( shp, xoff=-self.mfgrid.xoffset, yoff=-self.mfgrid.yoffset ) - if self.mfgrid.angrot != 0.0: + if self.mfgrid.angrot != 0.0 and not self.local: shp = affinity_loc.rotate( shp, -self.mfgrid.angrot, origin=(0.0, 0.0) ) @@ -1982,7 +2004,7 @@ def _intersect_polygon_structured( self.mfgrid.angrot != 0.0 or self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0 - ): + ) and not self.local: v_realworld = [] if intersect.geom_type.startswith("Multi"): for ipoly in intersect.geoms: