From af7b92d8cbc589298f2c91a5caf6b5b9d6076edc Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 09:54:35 +0100 Subject: [PATCH 01/16] refactor(GridIntersect): deprecations - deprecate method kwarg - remove keepzerolengths from intersect() - clean up deprecated methods --- flopy/utils/gridintersect.py | 1210 +--------------------------------- 1 file changed, 30 insertions(+), 1180 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 9effe75dc1..577a688db5 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -1,6 +1,8 @@ import warnings import numpy as np +from matplotlib.patches import PathPatch +from matplotlib.path import Path from pandas import DataFrame from .geometry import transform @@ -9,19 +11,6 @@ shapely = import_optional_dependency("shapely", errors="silent") -# TODO: remove the following methods and classes in version 3.10.0 -# - ModflowGridIndices -# - GridIntersect: -# - remove method kwarg from __init__ -# - remove structured methods from intersect() and intersects() -# - _intersect_point_structured() -# - _intersect_linestring_structured() -# - _get_nodes_intersecting_linestring() -# - _check_adjacent_cells_intersecting_line() -# - _intersect_rectangle_structured() -# - _intersect_polygon_structured() -# - _transform_geo_interface_polygon() - def parse_shapely_ix_result(collection, ix_result, shptyps=None): """Recursive function for parsing shapely intersection results. Returns a @@ -78,29 +67,19 @@ class GridIntersect: in slower performance. Therefore, it can sometimes be faster to intersect each individual shape in a collection than it is to intersect with the whole collection at once. - - Building the STR-tree can take a while for large grids. Once built the - intersect routines (for individual shapes) should be pretty fast. """ - def __init__(self, mfgrid, method=None, rtree=True, local=False): + def __init__(self, mfgrid, rtree=True, local=False): """Intersect shapes (Point, Linestring, Polygon) with a modflow grid. Parameters ---------- mfgrid : flopy modflowgrid MODFLOW grid as implemented in flopy - method : str, optional - Method to use for intersection shapes with the grid. Method 'vertex' - will be the only option in the future. Method 'structured' is deprecated. - This keyword argument will be removed in a future release. - - .. deprecated:: 3.9.0 - method="vertex" will be the only option from 3.10.0 - 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). + build an STR-Tree if True (default). If False no STR-tree + is built, but intersects will filter and loop through candidate model + gridcells (which is generally slower and not recommended). 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). @@ -110,67 +89,30 @@ def __init__(self, mfgrid, method=None, rtree=True, local=False): ) self.mfgrid = mfgrid self.local = local - # TODO: remove method kwarg in version v3.10.0 - # keep default behavior for v3.9.0, but warn if method is not vertex - # allow silencing of warning with method="vertex" in v3.9.0 - if method is None: - # determine method from grid_type - self.method = self.mfgrid.grid_type - else: - # set method - self.method = method - if self.method != "vertex": - warnings.warn( - ( - 'Note `method="structured"` is deprecated. ' - 'Pass `method="vertex"` to silence this warning. ' - "This will be the new default in a future release and this " - "keyword argument will be removed." - ), - category=DeprecationWarning, - ) self.rtree = rtree - # really only necessary for method=='vertex' as structured methods - # do not require a full list of shapely geometries, but useful to be - # able to obtain the grid shapes nonetheless self._set_method_get_gridshapes() - if self.method == "vertex": - # build arrays of geoms and cellids - self.geoms, self.cellids = self._get_gridshapes() - - # build STR-tree if specified - if self.rtree: - strtree = import_optional_dependency( - "shapely.strtree", - error_message="STRTree requires shapely", - ) - self.strtree = strtree.STRtree(self.geoms) + # build arrays of geoms and cellids + self.geoms, self.cellids = self._get_gridshapes() - elif self.method == "structured" and mfgrid.grid_type == "structured": - # geoms and cellids do not need to be assigned for structured - # methods - self.geoms = None - self.cellids = None - - else: - raise ValueError( - f"Method '{self.method}' not recognized or not supported " - f"for grid_type '{self.mfgrid.grid_type}'!" + # build STR-tree if specified + if self.rtree: + strtree = import_optional_dependency( + "shapely.strtree", + error_message="STRTree requires shapely", ) + self.strtree = strtree.STRtree(self.geoms) def intersect( self, shp, shapetype=None, sort_by_cellid=True, - keepzerolengths=False, return_all_intersections=False, contains_centroid=False, min_area_fraction=None, geo_dataframe=False, - shapely2=None, ): """Method to intersect a shape with a model grid. @@ -185,9 +127,6 @@ def intersect( sort_by_cellid : bool sort results by cellid, ensures cell with lowest cellid is returned for boundary cases when using vertex methods, default is True - keepzerolengths : bool - boolean method to keep zero length intersections for linestring - intersection, only used if shape type is "linestring" return_all_intersections : bool, optional if True, return multiple intersection results for points or linestrings on grid cell boundaries (e.g. returns 2 intersection @@ -211,53 +150,28 @@ def intersect( a record array containing information about the intersection or a geopandas.GeoDataFrame if geo_dataframe=True """ - if shapely2 is not None: - warnings.warn( - "The shapely2 keyword argument is deprecated. " - "Shapely<2 support was dropped in flopy version 3.9.0." - ) gu = GeoSpatialUtil(shp, shapetype=shapetype) shp = gu.shapely if gu.shapetype in {"Point", "MultiPoint"}: - if self.method == "structured" and self.mfgrid.grid_type == "structured": - rec = self._intersect_point_structured( - shp, return_all_intersections=return_all_intersections - ) - else: - rec = self._intersect_point_shapely( - shp, - sort_by_cellid=sort_by_cellid, - return_all_intersections=return_all_intersections, - ) + rec = self._intersect_point_shapely( + shp, + sort_by_cellid=sort_by_cellid, + return_all_intersections=return_all_intersections, + ) elif gu.shapetype in {"LineString", "MultiLineString"}: - if self.method == "structured" and self.mfgrid.grid_type == "structured": - rec = self._intersect_linestring_structured( - shp, - keepzerolengths, - return_all_intersections=return_all_intersections, - ) - else: - rec = self._intersect_linestring_shapely( - shp, - keepzerolengths, - sort_by_cellid=sort_by_cellid, - return_all_intersections=return_all_intersections, - ) + rec = self._intersect_linestring_shapely( + shp, + sort_by_cellid=sort_by_cellid, + return_all_intersections=return_all_intersections, + ) elif gu.shapetype in {"Polygon", "MultiPolygon"}: - if self.method == "structured" and self.mfgrid.grid_type == "structured": - rec = self._intersect_polygon_structured( - shp, - contains_centroid=contains_centroid, - min_area_fraction=min_area_fraction, - ) - else: - rec = self._intersect_polygon_shapely( - shp, - sort_by_cellid=sort_by_cellid, - contains_centroid=contains_centroid, - min_area_fraction=min_area_fraction, - ) + rec = self._intersect_polygon_shapely( + shp, + sort_by_cellid=sort_by_cellid, + contains_centroid=contains_centroid, + min_area_fraction=min_area_fraction, + ) else: raise TypeError(f"Shapetype {gu.shapetype} is not supported") @@ -492,17 +406,9 @@ def _intersect_point_shapely( def _intersect_linestring_shapely( self, shp, - keepzerolengths=False, sort_by_cellid=True, return_all_intersections=False, ): - if keepzerolengths: - warnings.warn( - "`keepzerolengths` is deprecated. For obtaining all cellids that " - "intersect with a LineString, use `intersects()`.", - DeprecationWarning, - ) - if self.rtree: qcellids = self.strtree.query(shp, predicate="intersects") else: @@ -723,919 +629,6 @@ def intersects(self, shp, shapetype=None, dataframe=False): return DataFrame(rec) return rec - def _intersect_point_structured(self, shp, return_all_intersections=False): - """intersection method for intersecting points with structured grids. - - .. deprecated:: 3.9.0 - use _intersect_point_shapely() or set method="vertex" in GridIntersect. - - Parameters - ---------- - shp : shapely.geometry.Point or MultiPoint - point shape to intersect with grid - return_all_intersections : bool, optional - if True, return multiple intersection results for points on grid - cell boundaries (e.g. returns 2 intersection results if a point - lies on the boundary between two grid cells). The default is False, - which will return a single intersection result for boundary cases. - - Returns - ------- - numpy.recarray - a record array containing information about the intersection - """ - shapely_geo = import_optional_dependency("shapely.geometry") - - nodelist = [] - - Xe, Ye = self.mfgrid.xyedges - - if isinstance(shp, shapely_geo.Point): - shp = [shp] - elif isinstance(shp, shapely_geo.MultiPoint): - shp = list(shp.geoms) - else: - raise ValueError("expected Point or MultiPoint") - - ixshapes = [] - for p in shp: - # if grid is rotated or offset transform point to local coords - if ( - 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, - self.mfgrid.xoffset, - self.mfgrid.yoffset, - self.mfgrid.angrot_radians, - inverse=True, - ) - else: - rx = p.x - ry = p.y - - # two dimensional point - jpos = ModflowGridIndices.find_position_in_array(Xe, rx) - ipos = ModflowGridIndices.find_position_in_array(Ye, ry) - - if jpos is not None and ipos is not None: - # use only first idx if return_all_intersections is False - if not return_all_intersections: - if isinstance(jpos, list): - jpos = jpos[0] - if isinstance(ipos, list): - ipos = ipos[0] - # three dimensional point - if p._ndim == 3: - # find k, if ipos or jpos on boundary, use first entry - if isinstance(jpos, list): - jj = jpos[0] - else: - jj = jpos - if isinstance(ipos, list): - ii = ipos[0] - else: - ii = ipos - kpos = ModflowGridIndices.find_position_in_array( - self.mfgrid.botm[:, ii, jj], p.z - ) - # if z-position on boundary, use first k - if isinstance(kpos, list): - kpos = kpos[0] - if kpos is not None: - # point on boundary, either jpos or ipos has len > 1 - if isinstance(ipos, list) or isinstance(jpos, list): - # convert to list if needed for loop - if not isinstance(ipos, list): - ipos = [ipos] - if not isinstance(jpos, list): - jpos = [jpos] - for ii in ipos: - for jj in jpos: - nodelist.append((kpos, ii, jj)) - ixshapes.append(p) - # point not on boundary - else: - nodelist.append((kpos, ipos, jpos)) - ixshapes.append(p) - else: - # point on boundary, either jpos or ipos has len > 1 - if isinstance(ipos, list) or isinstance(jpos, list): - # convert to list if needed for loop - if not isinstance(ipos, list): - ipos = [ipos] - if not isinstance(jpos, list): - jpos = [jpos] - for ii in ipos: - for jj in jpos: - nodelist.append((ii, jj)) - ixshapes.append(p) - else: - nodelist.append((ipos, jpos)) - ixshapes.append(p) - - # remove duplicates - if not return_all_intersections: - tempnodes = [] - tempshapes = [] - for node, ixs in zip(nodelist, ixshapes): - if node not in tempnodes: - tempnodes.append(node) - tempshapes.append(ixs) - else: - tempshapes[-1] = shapely_geo.MultiPoint([tempshapes[-1], ixs]) - - ixshapes = tempshapes - nodelist = tempnodes - - rec = np.recarray( - len(nodelist), names=["cellids", "ixshapes"], formats=["O", "O"] - ) - rec.cellids = nodelist - rec.ixshapes = ixshapes - return rec - - def _intersect_linestring_structured( - self, shp, keepzerolengths=False, return_all_intersections=False - ): - """method for intersecting linestrings with structured grids. - - .. deprecated:: 3.9.0 - use _intersect_point_shapely() or set method="vertex" in GridIntersect. - - Parameters - ---------- - shp : shapely.geometry.Linestring or MultiLineString - linestring to intersect with grid - keepzerolengths : bool, optional - if True keep intersection results with length=0, in - other words, grid cells the linestring does not cross - but does touch, by default False - return_all_intersections : bool, optional - if True, return multiple intersection results for linestrings on - grid cell boundaries (e.g. returns 2 intersection results if a - linestring lies on the boundary between two grid cells). The - default is False, which will return a single intersection result - for boundary cases. - - Returns - ------- - numpy.recarray - a record array containing information about the intersection - """ - shapely = import_optional_dependency("shapely") - shapely_geo = import_optional_dependency("shapely.geometry") - affinity_loc = import_optional_dependency("shapely.affinity") - - # get local extent of grid - if ( - self.mfgrid.angrot != 0.0 - or self.mfgrid.xoffset != 0.0 - or self.mfgrid.yoffset != 0.0 - ): - xmin = np.min(self.mfgrid.xyedges[0]) - xmax = np.max(self.mfgrid.xyedges[0]) - ymin = np.min(self.mfgrid.xyedges[1]) - ymax = np.max(self.mfgrid.xyedges[1]) - else: - xmin, xmax, ymin, ymax = self.mfgrid.extent - 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 - ) and not self.local: - shp = affinity_loc.translate( - shp, xoff=-self.mfgrid.xoffset, yoff=-self.mfgrid.yoffset - ) - if self.mfgrid.angrot != 0.0 and not self.local: - shp = affinity_loc.rotate(shp, -self.mfgrid.angrot, origin=(0.0, 0.0)) - - # clip line to mfgrid bbox - lineclip = shp.intersection(pl) - - if lineclip.length == 0.0: # linestring does not intersect modelgrid - return np.recarray( - 0, - names=["cellids", "vertices", "lengths", "ixshapes"], - formats=["O", "O", "f8", "O"], - ) - if lineclip.geom_type == "MultiLineString": # there are multiple lines - nodelist, lengths, vertices = [], [], [] - ixshapes = [] - for ls in lineclip.geoms: - n, l, v, ixs = self._get_nodes_intersecting_linestring( - ls, return_all_intersections=return_all_intersections - ) - nodelist += n - lengths += l - # if necessary, transform coordinates back to real - # world coordinates - if ( - 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) - rx, ry = transform( - pt[:, 0], - pt[:, 1], - self.mfgrid.xoffset, - self.mfgrid.yoffset, - self.mfgrid.angrot_radians, - inverse=False, - ) - v_realworld.append(list(zip(rx, ry))) - ixs_realworld = [] - for ix in ixs: - ix_realworld = affinity_loc.rotate( - ix, self.mfgrid.angrot, origin=(0.0, 0.0) - ) - ix_realworld = affinity_loc.translate( - ix_realworld, self.mfgrid.xoffset, self.mfgrid.yoffset - ) - ixs_realworld.append(ix_realworld) - else: - v_realworld = v - ixs_realworld = ixs - vertices += v_realworld - ixshapes += ixs_realworld - else: # linestring is fully within grid - (nodelist, lengths, vertices, ixshapes) = ( - self._get_nodes_intersecting_linestring( - lineclip, return_all_intersections=return_all_intersections - ) - ) - # if necessary, transform coordinates back to real - # world coordinates - if ( - 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) - rx, ry = transform( - pt[:, 0], - pt[:, 1], - self.mfgrid.xoffset, - self.mfgrid.yoffset, - self.mfgrid.angrot_radians, - inverse=False, - ) - v_realworld.append(list(zip(rx, ry))) - vertices = v_realworld - - ix_shapes_realworld = [] - for ixs in ixshapes: - ixs = affinity_loc.rotate( - ixs, self.mfgrid.angrot, origin=(0.0, 0.0) - ) - ixs = affinity_loc.translate( - ixs, self.mfgrid.xoffset, self.mfgrid.yoffset - ) - ix_shapes_realworld.append(ixs) - ixshapes = ix_shapes_realworld - - # bundle linestrings in same cell - tempnodes = [] - templengths = [] - tempverts = [] - tempshapes = [] - unique_nodes = list(set(nodelist)) - parsed_nodes = [] - if len(unique_nodes) < len(nodelist): - for inode in nodelist: - # maintain order of nodes by keeping track of parsed nodes - if inode in parsed_nodes: - continue - templengths.append( - sum(l for l, i in zip(lengths, nodelist) if i == inode) - ) - tempverts.append([v for v, i in zip(vertices, nodelist) if i == inode]) - tempshapes.append( - [ix for ix, i in zip(ixshapes, nodelist) if i == inode] - ) - parsed_nodes.append(inode) - - nodelist = parsed_nodes - lengths = templengths - vertices = tempverts - ixshapes = tempshapes - - # eliminate any nodes that have a zero length - if not keepzerolengths: - tempnodes = [] - templengths = [] - tempverts = [] - tempshapes = [] - for i, _ in enumerate(nodelist): - if lengths[i] > 0: - tempnodes.append(nodelist[i]) - templengths.append(lengths[i]) - tempverts.append(vertices[i]) - ishp = ixshapes[i] - if isinstance(ishp, list): - ishp = shapely.unary_union(ishp) - tempshapes.append(ishp) - nodelist = tempnodes - lengths = templengths - vertices = tempverts - ixshapes = tempshapes - - rec = np.recarray( - len(nodelist), - names=["cellids", "vertices", "lengths", "ixshapes"], - formats=["O", "O", "f8", "O"], - ) - rec.vertices = vertices - rec.lengths = lengths - rec.cellids = nodelist - rec.ixshapes = ixshapes - - return rec - - def _get_nodes_intersecting_linestring( - self, linestring, return_all_intersections=False - ): - """helper function, intersect the linestring with the a structured grid - and return a list of node indices and the length of the line in that - node. - - .. deprecated:: 3.9.0 - method="structured" is deprecated. - - Parameters - ---------- - linestring: shapely.geometry.LineString or MultiLineString - shape to intersect with the grid - - Returns - ------- - nodelist, lengths, vertices: lists - lists containing node ids, lengths of intersects and the - start and end points of the intersects - """ - shapely_geo = import_optional_dependency("shapely.geometry") - - nodelist = [] - lengths = [] - vertices = [] - ixshapes = [] - - # start at the beginning of the line - x, y = linestring.xy - - # linestring already in local coords but - # because intersect_point does transform again - # we transform back to real world here if necessary - if ( - 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]], - self.mfgrid.xoffset, - self.mfgrid.yoffset, - self.mfgrid.angrot_radians, - inverse=False, - ) - else: - x0 = [x[0]] - y0 = [y[0]] - - (i, j) = self.intersect(shapely_geo.Point(x0[0], y0[0])).cellids[0] - Xe, Ye = self.mfgrid.xyedges - xmin = Xe[j] - xmax = Xe[j + 1] - ymax = Ye[i] - ymin = Ye[i + 1] - pl = shapely_geo.box(xmin, ymin, xmax, ymax) - intersect = linestring.intersection(pl) - # if linestring starts in cell, exits, and re-enters - # a MultiLineString is returned. - ixshapes.append(intersect) - length = intersect.length - lengths.append(length) - if hasattr(intersect, "geoms"): - x, y = [], [] - for igeom in intersect.geoms: - x.append(igeom.xy[0]) - y.append(igeom.xy[1]) - x = np.concatenate(x) - y = np.concatenate(y) - else: - x = intersect.xy[0] - y = intersect.xy[1] - verts = [(ixy[0], ixy[1]) for ixy in zip(x, y)] - vertices.append(verts) - nodelist.append((i, j)) - - n = 0 - while True: - (i, j) = nodelist[n] - (node, length, verts, ixshape) = ( - self._check_adjacent_cells_intersecting_line( - linestring, (i, j), nodelist - ) - ) - - for inode, ilength, ivert, ix in zip(node, length, verts, ixshape): - if inode is not None: - if not return_all_intersections: - if ivert not in vertices: - nodelist.append(inode) - lengths.append(ilength) - vertices.append(ivert) - ixshapes.append(ix) - else: - nodelist.append(inode) - lengths.append(ilength) - vertices.append(ivert) - ixshapes.append(ix) - - if n == len(nodelist) - 1: - break - n += 1 - - return nodelist, lengths, vertices, ixshapes - - def _check_adjacent_cells_intersecting_line(self, linestring, i_j, nodelist): - """helper method that follows a line through a structured grid. - - .. deprecated:: 3.9.0 - method="structured" is deprecated. - - Parameters - ---------- - linestring : shapely.geometry.LineString - shape to intersect with the grid - i_j : tuple - tuple containing (nrow, ncol) - nodelist : list of tuples - list of node ids that have already been added - as intersections - - Returns - ------- - node, length, verts: lists - lists containing nodes, lengths and vertices of - intersections with adjacent cells relative to the - current cell (i, j) - """ - shapely_geo = import_optional_dependency("shapely.geometry") - - i, j = i_j - - Xe, Ye = self.mfgrid.xyedges - - node = [] - length = [] - verts = [] - ixshape = [] - - # check to left - if j > 0: - ii = i - jj = j - 1 - if (ii, jj) not in nodelist: - xmin = Xe[jj] - xmax = Xe[jj + 1] - ymax = Ye[ii] - ymin = Ye[ii + 1] - pl = shapely_geo.box(xmin, ymin, xmax, ymax) - if linestring.intersects(pl): - intersect = linestring.intersection(pl) - ixshape.append(intersect) - length.append(intersect.length) - if hasattr(intersect, "geoms"): - x, y = [], [] - for igeom in intersect.geoms: - x.append(igeom.xy[0]) - y.append(igeom.xy[1]) - x = np.concatenate(x) - y = np.concatenate(y) - else: - x = intersect.xy[0] - y = intersect.xy[1] - verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)]) - node.append((ii, jj)) - - # check to right - if j < self.mfgrid.ncol - 1: - ii = i - jj = j + 1 - if (ii, jj) not in nodelist: - xmin = Xe[jj] - xmax = Xe[jj + 1] - ymax = Ye[ii] - ymin = Ye[ii + 1] - pl = shapely_geo.box(xmin, ymin, xmax, ymax) - if linestring.intersects(pl): - intersect = linestring.intersection(pl) - ixshape.append(intersect) - length.append(intersect.length) - if hasattr(intersect, "geoms"): - x, y = [], [] - for igeom in intersect.geoms: - x.append(igeom.xy[0]) - y.append(igeom.xy[1]) - x = np.concatenate(x) - y = np.concatenate(y) - else: - x = intersect.xy[0] - y = intersect.xy[1] - verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)]) - node.append((ii, jj)) - - # check to back - if i > 0: - ii = i - 1 - jj = j - if (ii, jj) not in nodelist: - xmin = Xe[jj] - xmax = Xe[jj + 1] - ymax = Ye[ii] - ymin = Ye[ii + 1] - pl = shapely_geo.box(xmin, ymin, xmax, ymax) - if linestring.intersects(pl): - intersect = linestring.intersection(pl) - ixshape.append(intersect) - length.append(intersect.length) - if hasattr(intersect, "geoms"): - x, y = [], [] - for igeom in intersect.geoms: - x.append(igeom.xy[0]) - y.append(igeom.xy[1]) - x = np.concatenate(x) - y = np.concatenate(y) - else: - x = intersect.xy[0] - y = intersect.xy[1] - verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)]) - node.append((ii, jj)) - - # check to front - if i < self.mfgrid.nrow - 1: - ii = i + 1 - jj = j - if (ii, jj) not in nodelist: - xmin = Xe[jj] - xmax = Xe[jj + 1] - ymax = Ye[ii] - ymin = Ye[ii + 1] - pl = shapely_geo.box(xmin, ymin, xmax, ymax) - if linestring.intersects(pl): - intersect = linestring.intersection(pl) - ixshape.append(intersect) - length.append(intersect.length) - if hasattr(intersect, "geoms"): - x, y = [], [] - for igeom in intersect.geoms: - x.append(igeom.xy[0]) - y.append(igeom.xy[1]) - x = np.concatenate(x) - y = np.concatenate(y) - else: - x = intersect.xy[0] - y = intersect.xy[1] - verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)]) - node.append((ii, jj)) - - # special case for linestrings intersecting in vertex and continuing - # towards bottom right, check diagonally to front-right, if no other - # neighbours found - if np.sum(length) == 0: - if (i < self.mfgrid.nrow - 1) and (j < self.mfgrid.ncol - 1): - ii = i + 1 - jj = j + 1 - if (ii, jj) not in nodelist: - xmin = Xe[jj] - xmax = Xe[jj + 1] - ymax = Ye[ii] - ymin = Ye[ii + 1] - pl = shapely_geo.box(xmin, ymin, xmax, ymax) - if linestring.intersects(pl): - intersect = linestring.intersection(pl) - ixshape.append(intersect) - length.append(intersect.length) - if hasattr(intersect, "geoms"): - x, y = [], [] - for igeom in intersect.geoms: - x.append(igeom.xy[0]) - y.append(igeom.xy[1]) - x = np.concatenate(x) - y = np.concatenate(y) - else: - x = intersect.xy[0] - y = intersect.xy[1] - verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)]) - node.append((ii, jj)) - - return node, length, verts, ixshape - - def _intersect_rectangle_structured(self, rectangle): - """intersect a rectangle with a structured grid to retrieve node ids of - intersecting grid cells. - - .. deprecated:: 3.9.0 - method="structured" is deprecated. - - Note: only works in local coordinates (i.e. non-rotated grid - with origin at (0, 0)) - - Parameters - ---------- - rectangle : list of tuples - list of lower-left coordinate and upper-right - coordinate: [(xmin, ymin), (xmax, ymax)] - - Returns - ------- - nodelist: list of tuples - list of tuples containing node ids with which - the rectangle intersects - """ - - shapely_geo = import_optional_dependency("shapely.geometry") - - nodelist = [] - - # return if rectangle does not contain any cells - if ( - self.mfgrid.angrot != 0.0 - or self.mfgrid.xoffset != 0.0 - or self.mfgrid.yoffset != 0.0 - ): - minx = np.min(self.mfgrid.xyedges[0]) - maxx = np.max(self.mfgrid.xyedges[0]) - miny = np.min(self.mfgrid.xyedges[1]) - maxy = np.max(self.mfgrid.xyedges[1]) - local_extent = [minx, maxx, miny, maxy] - else: - local_extent = self.mfgrid.extent - - xmin, xmax, ymin, ymax = local_extent - bgrid = shapely_geo.box(xmin, ymin, xmax, ymax) - (rxmin, rymin), (rxmax, rymax) = rectangle - b = shapely_geo.box(rxmin, rymin, rxmax, rymax) - - if not b.intersects(bgrid): - # return with nodelist as an empty list - return [] - - Xe, Ye = self.mfgrid.xyedges - - jmin = ModflowGridIndices.find_position_in_array(Xe, xmin) - if jmin is None: - if xmin <= Xe[0]: - jmin = 0 - elif xmin >= Xe[-1]: - jmin = self.mfgrid.ncol - 1 - - jmax = ModflowGridIndices.find_position_in_array(Xe, xmax) - if jmax is None: - if xmax <= Xe[0]: - jmax = 0 - elif xmax >= Xe[-1]: - jmax = self.mfgrid.ncol - 1 - - imin = ModflowGridIndices.find_position_in_array(Ye, ymax) - if imin is None: - if ymax >= Ye[0]: - imin = 0 - elif ymax <= Ye[-1]: - imin = self.mfgrid.nrow - 1 - - imax = ModflowGridIndices.find_position_in_array(Ye, ymin) - if imax is None: - if ymin >= Ye[0]: - imax = 0 - elif ymin <= Ye[-1]: - imax = self.mfgrid.nrow - 1 - - for i in range(imin, imax + 1): - for j in range(jmin, jmax + 1): - nodelist.append((i, j)) - - return nodelist - - def _intersect_polygon_structured( - self, shp, contains_centroid=False, min_area_fraction=None - ): - """intersect polygon with a structured grid. Uses bounding box of the - Polygon to limit search space. - - .. deprecated:: 3.9.0 - method="structured" is deprecated. Use `_intersect_polygon_shapely()`. - - - Notes - ----- - If performance is slow, try setting the method to 'vertex' - in the GridIntersect object. For polygons this is often - faster. - - Parameters - ---------- - shp : shapely.geometry.Polygon - polygon to intersect with the grid - contains_centroid : bool, optional - if True, only store intersection result if cell centroid is - contained within intersection shape - min_area_fraction : float, optional - float defining minimum intersection area threshold, if - intersection area is smaller than min_frac_area * cell_area, do - not store intersection result - - Returns - ------- - numpy.recarray - a record array containing information about the intersection - """ - shapely_geo = import_optional_dependency("shapely.geometry") - affinity_loc = import_optional_dependency("shapely.affinity") - - # initialize the result lists - nodelist = [] - areas = [] - vertices = [] - ixshapes = [] - - # transform polygon to local grid coordinates - 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 and not self.local: - shp = affinity_loc.rotate(shp, -self.mfgrid.angrot, origin=(0.0, 0.0)) - - # use the bounds of the polygon to restrict the cell search - minx, miny, maxx, maxy = shp.bounds - rectangle = ((minx, miny), (maxx, maxy)) - nodes = self._intersect_rectangle_structured(rectangle) - - for i, j in nodes: - if ( - self.mfgrid.angrot != 0.0 - or self.mfgrid.xoffset != 0.0 - or self.mfgrid.yoffset != 0.0 - ): - cell_coords = [ - (self.mfgrid.xyedges[0][j], self.mfgrid.xyedges[1][i]), - (self.mfgrid.xyedges[0][j + 1], self.mfgrid.xyedges[1][i]), - (self.mfgrid.xyedges[0][j + 1], self.mfgrid.xyedges[1][i + 1]), - (self.mfgrid.xyedges[0][j], self.mfgrid.xyedges[1][i + 1]), - ] - else: - cell_coords = self.mfgrid.get_cell_vertices(i, j) - cell_polygon = shapely_geo.Polygon(cell_coords) - if shp.intersects(cell_polygon): - intersect = shp.intersection(cell_polygon) - collection = parse_shapely_ix_result( - [], intersect, shptyps=["Polygon", "MultiPolygon"] - ) - if len(collection) == 0: - continue - if len(collection) > 1: - intersect = shapely_geo.MultiPolygon(collection) - else: - intersect = collection[0] - - # only store results if area > 0.0 - if intersect.area == 0.0: - continue - # option: only store result if cell centroid is contained - # within intersection result - if contains_centroid: - if not intersect.intersects(cell_polygon.centroid): - continue - # option: min_area_fraction, only store if intersected area - # is larger than fraction * cell_area - if min_area_fraction: - if intersect.area < (min_area_fraction * cell_polygon.area): - continue - - nodelist.append((i, j)) - areas.append(intersect.area) - - # if necessary, transform coordinates back to real - # world coordinates - if ( - 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: - v_realworld += self._transform_geo_interface_polygon(ipoly) - else: - v_realworld += self._transform_geo_interface_polygon(intersect) - intersect_realworld = affinity_loc.rotate( - intersect, self.mfgrid.angrot, origin=(0.0, 0.0) - ) - intersect_realworld = affinity_loc.translate( - intersect_realworld, self.mfgrid.xoffset, self.mfgrid.yoffset - ) - else: - v_realworld = intersect.__geo_interface__["coordinates"] - intersect_realworld = intersect - ixshapes.append(intersect_realworld) - vertices.append(v_realworld) - - rec = np.recarray( - len(nodelist), - names=["cellids", "vertices", "areas", "ixshapes"], - formats=["O", "O", "f8", "O"], - ) - rec.vertices = vertices - rec.areas = areas - rec.cellids = nodelist - rec.ixshapes = ixshapes - - return rec - - def _transform_geo_interface_polygon(self, polygon): - """Internal method, helper function to transform geometry - __geo_interface__. - - .. deprecated:: 3.9.0 - method="structured" is deprecated. Only used by - `_intersect_polygon_structured()` - - Used for translating intersection result coordinates back into - real-world coordinates. - - Parameters - ---------- - polygon : shapely.geometry.Polygon - polygon to transform coordinates for - - Returns - ------- - geom_list : list - list containing transformed coordinates in same structure as - the original __geo_interface__. - """ - - if polygon.geom_type.startswith("Multi"): - raise TypeError("Does not support Multi geometries!") - - geom_list = [] - for coords in polygon.__geo_interface__["coordinates"]: - geoms = [] - try: - # test depth of list/tuple - _ = coords[0][0][0] - if len(coords) == 2: - shell, holes = coords - else: - raise ValueError("Cannot parse __geo_interface__") - except TypeError: - shell = coords - holes = None - except Exception as e: - raise e - # transform shell coordinates - shell_pts = [] - for pt in shell: - rx, ry = transform( - [pt[0]], - [pt[1]], - self.mfgrid.xoffset, - self.mfgrid.yoffset, - self.mfgrid.angrot_radians, - inverse=False, - ) - shell_pts.append((rx, ry)) - geoms.append(shell_pts) - # transform holes coordinates if necessary - if holes: - holes_pts = [] - for pt in holes: - rx, ry = transform( - [pt[0]], - [pt[1]], - self.mfgrid.xoffset, - self.mfgrid.yoffset, - self.mfgrid.angrot_radians, - inverse=False, - ) - # append (shells, holes) to transformed coordinates list - geom_list.append(tuple(geoms)) - return geom_list - @staticmethod def plot_polygon(result, ax=None, **kwargs): """method to plot the polygon intersection results from the resulting @@ -1850,150 +843,7 @@ def plot_intersection_result(self, result, plot_grid=True, ax=None, **kwargs): return ax -class ModflowGridIndices: - """Collection of methods that can be used to find cell indices for a - structured, but irregularly spaced MODFLOW grid. - - .. deprecated:: 3.9.0 - This class is deprecated and will be removed in version 3.10.0. - """ - - @staticmethod - def find_position_in_array(arr, x): - """If arr has x positions for the left edge of a cell, then return the - cell index containing x. - - Parameters - ---------- - arr : A one dimensional array (such as Xe) that contains - coordinates for the left cell edge. - - x : float - The x position to find in arr. - """ - jpos = [] - - if np.isclose(x, arr[-1]): - return len(arr) - 2 - - xmin = min(arr[0], arr[-1]) - xmax = max(arr[0], arr[-1]) - - if np.isclose(x, xmin): - x = xmin - if np.isclose(x, xmax): - x = xmax - if not (xmin <= x <= xmax): - return None - - # go through each position - for j in range(len(arr) - 1): - xl = arr[j] - xr = arr[j + 1] - frac = (x - xl) / (xr - xl) - if 0.0 <= frac <= 1.0: - jpos.append(j) - if len(jpos) == 0: - return None - elif len(jpos) == 1: - return jpos[0] - else: - return jpos - - @staticmethod - def kij_from_nodenumber(nodenumber, nlay, nrow, ncol): - """Convert the modflow node number to a zero-based layer, row and - column format. Return (k0, i0, j0). - - Parameters - ---------- - nodenumber: int - The cell nodenumber, ranging from 1 to number of - nodes. - nlay: int - The number of layers. - nrow: int - The number of rows. - ncol: int - The number of columns. - """ - if nodenumber > nlay * nrow * ncol: - raise Exception("Error in function kij_from_nodenumber...") - n = nodenumber - 1 - k = int(n / nrow / ncol) - i = int((n - k * nrow * ncol) / ncol) - j = n - k * nrow * ncol - i * ncol - return (k, i, j) - - @staticmethod - def nodenumber_from_kij(k, i, j, nrow, ncol): - """Calculate the nodenumber using the zero-based layer, row, and column - values. The first node has a value of 1. - - Parameters - ---------- - k : int - The model layer number as a zero-based value. - i : int - The model row number as a zero-based value. - j : int - The model column number as a zero-based value. - nrow : int - The number of model rows. - ncol : int - The number of model columns. - """ - return k * nrow * ncol + i * ncol + j + 1 - - @staticmethod - def nn0_from_kij(k, i, j, nrow, ncol): - """Calculate the zero-based nodenumber using the zero-based layer, row, - and column values. The first node has a value of 0. - - Parameters - ---------- - k : int - The model layer number as a zero-based value. - i : int - The model row number as a zero-based value. - j : int - The model column number as a zero-based value. - nrow : int - The number of model rows. - ncol : int - The number of model columns. - """ - return k * nrow * ncol + i * ncol + j - - @staticmethod - def kij_from_nn0(n, nlay, nrow, ncol): - """Convert the node number to a zero-based layer, row and column - format. Return (k0, i0, j0). - - Parameters - ---------- - nodenumber : int - The cell nodenumber, ranging from 0 to number of - nodes - 1. - nlay : int - The number of layers. - nrow : int - The number of rows. - ncol : int - The number of columns. - """ - if n > nlay * nrow * ncol: - raise Exception("Error in function kij_from_nodenumber...") - k = int(n / nrow / ncol) - i = int((n - k * nrow * ncol) / ncol) - j = n - k * nrow * ncol - i * ncol - return (k, i, j) - - def _polygon_patch(polygon, **kwargs): - from matplotlib.patches import PathPatch - from matplotlib.path import Path - patch = PathPatch( Path.make_compound_path( Path(np.asarray(polygon.exterior.coords)[:, :2]), From 9483c74b92bc34f95e278539b6fde1eebe3cdd22 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 10:25:56 +0100 Subject: [PATCH 02/16] feat(GridIntersect): support numpy arrays with geometries in intersects() - allow arrays of numpy arrays in intersects(), returns all candidate cells for intersection per geometry - add predicate to query_grid() - simplify names ofintersection methods, deprecate old internal methods - move logic to translate node number to rowcol to separate method --- flopy/utils/gridintersect.py | 208 +++++++++++++++++++++++++++++++---- 1 file changed, 185 insertions(+), 23 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 577a688db5..c93e46ce5e 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -1,7 +1,8 @@ import warnings +import matplotlib.pyplot as plt import numpy as np -from matplotlib.patches import PathPatch +from matplotlib.collections import PatchCollection from matplotlib.path import Path from pandas import DataFrame @@ -104,6 +105,40 @@ def __init__(self, mfgrid, rtree=True, local=False): ) self.strtree = strtree.STRtree(self.geoms) + def _parse_input_shape(self, shp, shapetype=None): + """Internal method to parse input shape. + + Allows numpy arrays containing shapely geometries, otherwise delegates to + GeoSpatialUtil. + + Parameters + ---------- + shp : shapely.geometry, geojson object, shapefile.Shape, np.ndarray, + or flopy geometry object + shape to intersect with the grid + shapetype : str, optional + type of shape (i.e. "point", "linestring", "polygon" or their + multi-variants), used by GeoSpatialUtil if shp is passed as a list + of vertices, default is None + + Returns + ------- + shp : shapely.geometry or np.ndarray + shapely geometry or array of shapely geometries + """ + if isinstance(shp, np.ndarray) and isinstance(shp[0], shapely.Geometry): + shapetypes = shapely.get_type_id(shp) + assert len(np.unique(shapetypes)) == 1, ( + "If passing an array of shapely geometries, all geometries must be " + "of the same type." + ) + shapetype = shapely.GeometryType(shapetypes[0]) + else: + gu = GeoSpatialUtil(shp, shapetype=shapetype) + shp = gu.shapely + shapetype = gu.shapetype + return shp, shapetype + def intersect( self, shp, @@ -150,30 +185,52 @@ def intersect( a record array containing information about the intersection or a geopandas.GeoDataFrame if geo_dataframe=True """ - gu = GeoSpatialUtil(shp, shapetype=shapetype) - shp = gu.shapely + shp, shapetype = self._parse_input_shape(shp, shapetype=shapetype) - if gu.shapetype in {"Point", "MultiPoint"}: - rec = self._intersect_point_shapely( + # if array, only accept length 1 + if isinstance(shp, np.ndarray) and len(shp) > 1: + raise ValueError( + "intersect() only accepts arrays containing one " + f"{shapetype.name.lower()} at a time." + ) + + if shapetype in { + "Point", + "MultiPoint", + shapely.GeometryType.POINT, + shapely.GeometryType.MULTIPOINT, + }: + rec = self._intersect_point( shp, sort_by_cellid=sort_by_cellid, return_all_intersections=return_all_intersections, ) - elif gu.shapetype in {"LineString", "MultiLineString"}: - rec = self._intersect_linestring_shapely( + elif shapetype in { + "LineString", + "MultiLineString", + shapely.GeometryType.LINESTRING, + shapely.GeometryType.MULTILINESTRING, + shapely.GeometryType.LINEARRING, + }: + rec = self._intersect_linestring( shp, sort_by_cellid=sort_by_cellid, return_all_intersections=return_all_intersections, ) - elif gu.shapetype in {"Polygon", "MultiPolygon"}: - rec = self._intersect_polygon_shapely( + elif shapetype in { + "Polygon", + "MultiPolygon", + shapely.GeometryType.POLYGON, + shapely.GeometryType.MULTIPOLYGON, + }: + rec = self._intersect_polygon( shp, sort_by_cellid=sort_by_cellid, contains_centroid=contains_centroid, min_area_fraction=min_area_fraction, ) else: - raise TypeError(f"Shapetype {gu.shapetype} is not supported") + raise TypeError(f"Shapetype {shapetype} is not supported") if geo_dataframe: gpd = import_optional_dependency("geopandas") @@ -297,7 +354,7 @@ def _vtx_grid_to_geoms_cellids(self): ] return np.array(geoms), np.arange(self.mfgrid.ncpl) - def query_grid(self, shp): + def query_grid(self, shp, predicate=None): """Perform spatial query on grid with shapely geometry. If no spatial query is possible returns all grid cells. @@ -305,6 +362,9 @@ def query_grid(self, shp): ---------- shp : shapely.geometry shapely geometry + predicate : str, optional + spatial predicate to use for query, default is None. See + documentation of self.strtree.query for options. Returns ------- @@ -312,7 +372,7 @@ def query_grid(self, shp): array containing cellids of grid cells in query result """ if self.rtree: - result = self.strtree.query(shp) + result = self.strtree.query(shp, predicate=predicate) else: # no spatial query result = self.cellids @@ -344,7 +404,18 @@ def filter_query_result(self, cellids, shp): qcellids = cellids[shapely.intersects(self.geoms[cellids], shp)] return qcellids - def _intersect_point_shapely( + def _intersect_point_shapely(self, *args, **kwargs): + """Deprecated method, use _intersect_point instead.""" + import warnings + + warnings.warn( + "_intersect_point_shapely is deprecated, " + "use _intersect_point instead.", + DeprecationWarning, + ) + return self._intersect_point(*args, **kwargs) + + def _intersect_point( self, shp, sort_by_cellid=True, @@ -403,7 +474,18 @@ def _intersect_point_shapely( return rec - def _intersect_linestring_shapely( + def _intersect_linestring_shapely(self, *args, **kwargs): + """Deprecated method, use _intersect_linestring instead.""" + import warnings + + warnings.warn( + "_intersect_linestring_shapely is deprecated, " + "use _intersect_linestring instead.", + DeprecationWarning, + ) + return self._intersect_linestring(*args, **kwargs) + + def _intersect_linestring( self, shp, sort_by_cellid=True, @@ -517,7 +599,17 @@ def parse_linestrings_in_geom_collection(gc): return rec - def _intersect_polygon_shapely( + def _intersect_polygon_shapely(self, *args, **kwargs): + """Deprecated method, use _intersect_polygon instead.""" + import warnings + + warnings.warn( + "_intersect_polygon_shapely is deprecated, use _intersect_polygon instead.", + DeprecationWarning, + ) + return self._intersect_polygon(*args, **kwargs) + + def _intersect_polygon( self, shp, sort_by_cellid=True, @@ -594,7 +686,13 @@ def parse_polygons_in_geom_collection(gc): return rec - def intersects(self, shp, shapetype=None, dataframe=False): + def intersects( + self, + shp, + shapetype=None, + dataframe=False, + return_cellids=True, + ): """Return cellids for grid cells that intersect with shape. Parameters @@ -608,6 +706,14 @@ def intersects(self, shp, shapetype=None, dataframe=False): passed as a list of vertices, default is None dataframe : bool, optional if True, return a pandas.DataFrame, default is False + return_all_intersections : bool, optional + if True (default), return multiple intersection results for points on grid + cell boundaries (e.g. returns 2 intersection results if a point lies on the + boundary between two grid cells). + return_cellids : bool, optional + if True (default), return cellids of intersected grid cells. + If False, only return grid node numbers, i.e. index of entry in + ``GridIntersect.geoms``. Returns ------- @@ -615,18 +721,74 @@ def intersects(self, shp, shapetype=None, dataframe=False): a record array or pandas.DataFrame containing cell IDs of the gridcells the shape intersects with. """ - shp = GeoSpatialUtil(shp, shapetype=shapetype).shapely - qfiltered = self.strtree.query(shp, predicate="intersects") + shp, shapetype = self._parse_input_shape(shp, shapetype=shapetype) + + # query grid or strtree + qcellids = self.query_grid(shp, predicate="intersects") + if not self.rtree: + if isinstance(shp, np.ndarray) and len(shp) > 1: + raise ValueError( + "intersects() only accepts arrays containing one " + f"{shapetype.name.lower()} at a time when rtree=False." + ) + qfiltered = self.filter_query_result(shp, qcellids) + else: + qfiltered = qcellids - # build rec-array - rec = np.recarray(len(qfiltered), names=["cellids"], formats=["O"]) - if self.mfgrid.grid_type == "structured": - rec.cellids = list(zip(*self.mfgrid.get_lrc([qfiltered])[0][1:])) + # sort cellids + if qfiltered.ndim == 1: + qfiltered = np.sort(qfiltered) else: + qfiltered = qfiltered[:, np.lexsort((qfiltered[1], qfiltered[0]))] + + # determine size of output array + nr = len(qfiltered) if qfiltered.ndim == 1 else qfiltered.shape[1] + + # build rec-array + rec = np.recarray( + nr, + names=["shp_ids", "cellids"], + formats=[ + int, + "O" + if (return_cellids and self.mfgrid.grid_type == "structured") + else float, + ], + ) + # shp was passed as single geometry + if qfiltered.ndim == 1: + rec.shp_ids[:] = 0 rec.cellids = qfiltered + # shape passed as array of geometries + else: + rec.shp_ids = qfiltered[0] + rec.cellids = qfiltered[1] + + if self.mfgrid.grid_type == "structured" and return_cellids: + rec.cellids = self._nodenumber_to_rowcol(rec.cellids) if dataframe: - return DataFrame(rec) + return DataFrame(rec).set_index("shp_ids") + return rec + + def _nodenumber_to_rowcol(self, nodes): + """Convert node number to (row, col) tuple. + + Parameters + ---------- + nodes : array_like + array of cellids to convert + + Returns + ------- + array_like + array of (row, col) tuples + """ + # cast to float and allow nans + idx = np.nonzero(~np.isnan(nodes.astype(float))) + rc = np.full_like(nodes, np.nan, dtype=object) + rc[idx] = list(zip(*self.mfgrid.get_lrc([nodes[idx].astype(int)])[0][1:])) + return rc return rec @staticmethod From 3dc27c475449eac62d7132f6ccd9d0a09817c917 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 10:31:38 +0100 Subject: [PATCH 03/16] refactor(GridIntersect): consistency and cleanup - flip arguments in filter_query_result to be consistent with other class methods - remove get_gridshapes method() - move imports to top of file --- flopy/utils/gridintersect.py | 63 ++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 28 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index c93e46ce5e..7c2b5fbcd7 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -3,6 +3,7 @@ import matplotlib.pyplot as plt import numpy as np from matplotlib.collections import PatchCollection +from matplotlib.patches import PathPatch from matplotlib.path import Path from pandas import DataFrame @@ -92,10 +93,18 @@ def __init__(self, mfgrid, rtree=True, local=False): self.local = local self.rtree = rtree - self._set_method_get_gridshapes() - # build arrays of geoms and cellids - self.geoms, self.cellids = self._get_gridshapes() + if self.mfgrid.grid_type == "structured": + self.geoms, self.cellids = self._rect_grid_to_geoms_cellids() + elif self.mfgrid.grid_type == "vertex": + self.geoms, self.cellids = self._vtx_grid_to_geoms_cellids() + elif self.mfgrid.grid_type == "unstructured": + raise NotImplementedError() + self.geoms, self.cellids = self._usg_grid_to_geoms_cellids() + else: + raise NotImplementedError( + f"Grid type {self.mfgrid.grid_type} not supported" + ) # build STR-tree if specified if self.rtree: @@ -245,17 +254,6 @@ def intersect( return rec - def _set_method_get_gridshapes(self): - """internal method, set self._get_gridshapes to the appropriate method - for obtaining grid cell geometries.""" - # Set method for obtaining grid shapes - if self.mfgrid.grid_type == "structured": - self._get_gridshapes = self._rect_grid_to_geoms_cellids - elif self.mfgrid.grid_type == "vertex": - self._get_gridshapes = self._vtx_grid_to_geoms_cellids - elif self.mfgrid.grid_type == "unstructured": - raise NotImplementedError() - def _rect_grid_to_geoms_cellids(self): """internal method, return shapely polygons and cellids for structured grid cells. @@ -306,7 +304,6 @@ def _rect_grid_to_geoms_cellids(self): indices=np.repeat(cellids, 4), ) ) - return geoms, cellids def _usg_grid_to_geoms_cellids(self): @@ -378,7 +375,7 @@ def query_grid(self, shp, predicate=None): result = self.cellids return result - def filter_query_result(self, cellids, shp): + def filter_query_result(self, shp, cellids): """Filter array of geometries to obtain grid cells that intersect with shape. @@ -387,19 +384,32 @@ def filter_query_result(self, cellids, shp): Parameters ---------- - cellids : iterable - iterable of cellids, query result shp : shapely.geometry shapely geometry that is prepared and used to filter query result + cellids : iterable + iterable of cellids, query result Returns ------- array_like filter or generator containing polygons that intersect with shape """ + # flipped arguments to be consistent with all other methods in class + msg = ( + "the cellids and shp arguments were flipped, please" + " pass them as filter_query_result(shp, cellids)" + ) + if isinstance(cellids, np.ndarray): + if isinstance(cellids[0], shapely.Geometry): + warnings.warn(msg) + cellids, shp = shp, cellids + elif isinstance(cellids, shapely.Geometry): + warnings.warn(msg) + cellids, shp = shp, cellids + # get only gridcells that intersect - if not shapely.is_prepared(shp): + if not shapely.is_prepared(shp).all(): shapely.prepare(shp) qcellids = cellids[shapely.intersects(self.geoms[cellids], shp)] return qcellids @@ -409,8 +419,7 @@ def _intersect_point_shapely(self, *args, **kwargs): import warnings warnings.warn( - "_intersect_point_shapely is deprecated, " - "use _intersect_point instead.", + "_intersect_point_shapely is deprecated, use _intersect_point instead.", DeprecationWarning, ) return self._intersect_point(*args, **kwargs) @@ -476,7 +485,6 @@ def _intersect_point( def _intersect_linestring_shapely(self, *args, **kwargs): """Deprecated method, use _intersect_linestring instead.""" - import warnings warnings.warn( "_intersect_linestring_shapely is deprecated, " @@ -494,7 +502,7 @@ def _intersect_linestring( if self.rtree: qcellids = self.strtree.query(shp, predicate="intersects") else: - qcellids = self.filter_query_result(self.cellids, shp) + qcellids = self.filter_query_result(shp, self.cellids) if sort_by_cellid: qcellids = np.sort(qcellids) @@ -619,7 +627,7 @@ def _intersect_polygon( if self.rtree: qcellids = self.strtree.query(shp, predicate="intersects") else: - qcellids = self.filter_query_result(self.cellids, shp) + qcellids = self.filter_query_result(shp, self.cellids) if sort_by_cellid: qcellids = np.sort(qcellids) @@ -812,10 +820,6 @@ def plot_polygon(result, ax=None, **kwargs): matplotlib.pyplot.axes returns the axes handle """ - - import matplotlib.pyplot as plt - from matplotlib.collections import PatchCollection - if ax is None: _, ax = plt.subplots() ax.set_aspect("equal", adjustable="box") @@ -1006,6 +1010,9 @@ def plot_intersection_result(self, result, plot_grid=True, ax=None, **kwargs): def _polygon_patch(polygon, **kwargs): + from matplotlib.patches import PathPatch + from matplotlib.path import Path + patch = PathPatch( Path.make_compound_path( Path(np.asarray(polygon.exterior.coords)[:, :2]), From 79faf9f4b0802dfbc9a8c30d8b2a46bd536728e0 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 10:35:03 +0100 Subject: [PATCH 04/16] feat(GridIntersect): fast point location method - add points_to_cellids() for fast 1:1 mapping between points and grid cells, return nan when outside of grid --- flopy/utils/gridintersect.py | 80 +++++++++++++++++++++++++++++++++++- 1 file changed, 79 insertions(+), 1 deletion(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 7c2b5fbcd7..1748b1a99d 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -779,7 +779,7 @@ def intersects( return DataFrame(rec).set_index("shp_ids") return rec - def _nodenumber_to_rowcol(self, nodes): + def _nodenumbers_to_rowcol(self, nodes): """Convert node number to (row, col) tuple. Parameters @@ -797,6 +797,84 @@ def _nodenumber_to_rowcol(self, nodes): rc = np.full_like(nodes, np.nan, dtype=object) rc[idx] = list(zip(*self.mfgrid.get_lrc([nodes[idx].astype(int)])[0][1:])) return rc + + def points_to_cellids( + self, + pts, + dataframe=False, + return_nodenumbers=False, + ): + """Return cellids of grid cells that intersect with shape. + + Parameters + ---------- + pts : shapely.geometry, geojson geometry, shapefile.shape, + or flopy geometry object + points shape to intersect with the grid + dataframe : bool, optional + if True, return a pandas.DataFrame, default is False + return_nodenumbers : bool, optional + if False (default), return cellids of intersected grid cells. + If True, return grid node numbers, i.e. index of entry in + ``GridIntersect.geoms``. + + Returns + ------- + numpy.recarray or pandas.DataFrame + a record array or pandas.DataFrame containing cell IDs of the gridcells + the shape intersects with. + """ + if not self.rtree: + raise ValueError( + "points_to_cellids() requires rtree=True when" + " initializing GridIntersect" + ) + if not isinstance(pts, np.ndarray): + if shapely.get_type_id(pts) == shapely.GeometryType.MULTIPOINT: + pts = np.array(shapely.get_parts(pts)) + else: + pts = np.array([pts]) + + # query grid or strtree + qfiltered = self.query_grid(pts, predicate="intersects") + + # sort cellids + if qfiltered.ndim == 1: + qfiltered = np.sort(qfiltered) + else: + qfiltered = qfiltered[:, np.lexsort((qfiltered[1], qfiltered[0]))] + + # determine size of output array + if isinstance(pts, np.ndarray): + # one result per point + nr = len(pts) + else: + # single point + nr = 1 if len(qfiltered) > 0 else 0 # 1 if intersects, else 0 + + # build rec-array + rec = np.recarray( + nr, + names=["shp_ids", "cellids"], + formats=[ + int, + float # to allow nans + if (return_nodenumbers or self.mfgrid.grid_type != "structured") + else "O", + ], + ) + + rec.cellids = np.nan + # return only 1 gr id cell intersection result for each point + uniques, idx = np.unique(qfiltered[0], return_index=True) + rec.shp_ids = np.arange(nr) + rec.cellids[uniques] = qfiltered[1, idx] + + if self.mfgrid.grid_type == "structured" and not return_nodenumbers: + rec.cellids = self._nodenumbers_to_rowcol(rec.cellids) + + if dataframe: + return DataFrame(rec).set_index("shp_ids") return rec @staticmethod From 4107bffa5152c23ab8ef24a2bdf32f805fb7afc6 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 11:31:32 +0100 Subject: [PATCH 05/16] refactor(GridIntersect): improve naming and plotting - rename return_cellids to return_nodenumbers for clarity (node numbers always one number, cellid is somewhat ambiguous) - use shapely plotting for linestrings and points --- flopy/utils/gridintersect.py | 46 ++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 26 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 1748b1a99d..d9acaf2b00 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -699,7 +699,7 @@ def intersects( shp, shapetype=None, dataframe=False, - return_cellids=True, + return_nodenumbers=False, ): """Return cellids for grid cells that intersect with shape. @@ -718,9 +718,9 @@ def intersects( if True (default), return multiple intersection results for points on grid cell boundaries (e.g. returns 2 intersection results if a point lies on the boundary between two grid cells). - return_cellids : bool, optional - if True (default), return cellids of intersected grid cells. - If False, only return grid node numbers, i.e. index of entry in + return_nodenumbers : bool, optional + if False (default), return cellids of intersected grid cells. + If True, return grid node numbers, i.e. index of entry in ``GridIntersect.geoms``. Returns @@ -758,9 +758,9 @@ def intersects( names=["shp_ids", "cellids"], formats=[ int, - "O" - if (return_cellids and self.mfgrid.grid_type == "structured") - else float, + float + if (return_nodenumbers or self.mfgrid.grid_type != "structured") + else "O", ], ) # shp was passed as single geometry @@ -772,8 +772,8 @@ def intersects( rec.shp_ids = qfiltered[0] rec.cellids = qfiltered[1] - if self.mfgrid.grid_type == "structured" and return_cellids: - rec.cellids = self._nodenumber_to_rowcol(rec.cellids) + if self.mfgrid.grid_type == "structured" and not return_nodenumbers: + rec.cellids = self._nodenumbers_to_rowcol(rec.cellids) if dataframe: return DataFrame(rec).set_index("shp_ids") @@ -960,7 +960,8 @@ def plot_linestring(result, ax=None, cmap=None, **kwargs): matplotlib.pyplot.axes returns the axes handle """ - import matplotlib.pyplot as plt + + shapely_plot = import_optional_dependency("shapely.plotting") if ax is None: _, ax = plt.subplots() @@ -988,11 +989,7 @@ def plot_linestring(result, ax=None, cmap=None, **kwargs): c = f"C{i % 10}" else: c = colors[i] - if ishp.geom_type == "MultiLineString": - for part in ishp.geoms: - ax.plot(part.xy[0], part.xy[1], ls="-", c=c, **kwargs) - else: - ax.plot(ishp.xy[0], ishp.xy[1], ls="-", c=c, **kwargs) + shapely_plot.plot_line(ishp, ax=ax, color=c, **kwargs) return ax @@ -1017,24 +1014,20 @@ def plot_point(result, ax=None, **kwargs): matplotlib.pyplot.axes returns the axes handle """ - import matplotlib.pyplot as plt - - shapely_geo = import_optional_dependency("shapely.geometry") + shapely = import_optional_dependency("shapely") + shapely_plot = import_optional_dependency("shapely.plotting") if ax is None: _, ax = plt.subplots() - - x, y = [], [] # allow for result to be geodataframe geoms = ( result.ixshapes if isinstance(result, np.rec.recarray) else result.geometry ) - geo_coll = shapely_geo.GeometryCollection(list(geoms)) - collection = parse_shapely_ix_result([], geo_coll, ["Point"]) - for c in collection: - x.append(c.x) - y.append(c.y) - ax.scatter(x, y, **kwargs) + maskpts = np.isin( + shapely.get_type_id(geoms), + [shapely.GeometryType.POINT, shapely.GeometryType.MULTIPOINT], + ) + shapely_plot.plot_points(geoms[maskpts], ax=ax, **kwargs) return ax @@ -1075,6 +1068,7 @@ def plot_intersection_result(self, result, plot_grid=True, ax=None, **kwargs): [ shapely.GeometryType.LINESTRING, shapely.GeometryType.MULTILINESTRING, + shapely.GeometryType.LINEARRING, ], ).all(): ax = GridIntersect.plot_linestring(result, ax=ax, **kwargs) From 8660978692946076d840db4316714c24ec5bc3cf Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 11:34:50 +0100 Subject: [PATCH 06/16] feat(Gridintersect): support handling z-coordinates for points - add handle_z options: "ignore", "drop" and "return" to intersect() and points_to_cellids() - add method to compute layer position for z coordinates - improve _intersect_point() method by using intersects as starting point. --- flopy/utils/gridintersect.py | 118 ++++++++++++++++++++++++++++++++--- 1 file changed, 109 insertions(+), 9 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index d9acaf2b00..1a9f9dff1b 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -5,9 +5,9 @@ from matplotlib.collections import PatchCollection from matplotlib.patches import PathPatch from matplotlib.path import Path +from numpy.lib import recfunctions as nprecfns from pandas import DataFrame -from .geometry import transform from .geospatial_utils import GeoSpatialUtil from .utl_import import import_optional_dependency @@ -156,6 +156,7 @@ def intersect( return_all_intersections=False, contains_centroid=False, min_area_fraction=None, + handle_z="ignore", geo_dataframe=False, ): """Method to intersect a shape with a model grid. @@ -185,6 +186,13 @@ def intersect( float defining minimum intersection area threshold, if intersection area is smaller than min_frac_area * cell_area, do not store intersection result, only used if shape type is "polygon" + handle_z : str, optional + Method for handling z dimension in intersection results for point + intersections. Default is "ignore" which ignores z-dimension. Other + options are "drop" which only returns results for points within grid + top and bottom, or "return" which returns the computed layer position + for each point. Points above the grid are returned as +np.inf and below + the grid as -np.inf. geo_dataframe : bool, optional if True, return a geopandas GeoDataFrame, default is False @@ -214,6 +222,31 @@ def intersect( sort_by_cellid=sort_by_cellid, return_all_intersections=return_all_intersections, ) + + # handle elevation data for points + # if handle_z is "drop" or "return" + # if shp has z information + # if there are intersection results + if ( + (handle_z != "ignore") + and shapely.has_z(shp).any() + and len(rec.cellids) > 0 + ): + laypos = self.get_layer_from_z(shp, rec.cellids) + if handle_z == "drop": + mask_z = np.isfinite(laypos) + rec = rec[mask_z] + elif handle_z == "return": + # copy data to new array to include layer position + rec = nprecfns.append_fields( + rec, + names="layer", + data=laypos, + dtypes="f8", + usemask=False, + asrecarray=True, + ) + elif shapetype in { "LineString", "MultiLineString", @@ -430,10 +463,8 @@ def _intersect_point( sort_by_cellid=True, return_all_intersections=False, ): - if self.rtree: - qcellids = self.strtree.query(shp, predicate="intersects") - else: - qcellids = self.filter_query_result(self.cellids, shp) + r = self.intersects(shp, return_nodenumbers=True) + qcellids = r.cellids[np.isfinite(r.cellids)].astype(int) if sort_by_cellid: qcellids = np.sort(qcellids) @@ -442,7 +473,10 @@ def _intersect_point( # discard empty intersection results mask_empty = shapely.is_empty(ixresult) # keep only Point and MultiPoint - mask_type = np.isin(shapely.get_type_id(ixresult), [0, 4]) + mask_type = np.isin( + shapely.get_type_id(ixresult), + [shapely.GeometryType.POINT, shapely.GeometryType.MULTIPOINT], + ) ixresult = ixresult[~mask_empty & mask_type] qcellids = qcellids[~mask_empty & mask_type] @@ -801,6 +835,7 @@ def _nodenumbers_to_rowcol(self, nodes): def points_to_cellids( self, pts, + handle_z="ignore", dataframe=False, return_nodenumbers=False, ): @@ -813,6 +848,13 @@ def points_to_cellids( points shape to intersect with the grid dataframe : bool, optional if True, return a pandas.DataFrame, default is False + handle_z : str, optional + Method for handling z dimension in intersection results for point + intersections. Default is "ignore" which ignores z-dimension. Other + options are "drop" which only returns results for points within grid + top and bottom, or "return" which returns the computed layer position + for each point. Points above the grid are returned as +np.inf and below + the grid as -np.inf. return_nodenumbers : bool, optional if False (default), return cellids of intersected grid cells. If True, return grid node numbers, i.e. index of entry in @@ -873,6 +915,22 @@ def points_to_cellids( if self.mfgrid.grid_type == "structured" and not return_nodenumbers: rec.cellids = self._nodenumbers_to_rowcol(rec.cellids) + if handle_z != "ignore": + laypos = self.get_layer_from_z(pts, rec.cellids) + if handle_z == "drop": + mask_z = np.isfinite(laypos) + rec = rec[mask_z] + elif handle_z == "return": + # copy data to new array to include layer position + rec = nprecfns.append_fields( + rec, + names="layer", + data=laypos, + dtypes="f8", + usemask=False, + asrecarray=True, + ) + if dataframe: return DataFrame(rec).set_index("shp_ids") return rec @@ -1080,11 +1138,53 @@ def plot_intersection_result(self, result, plot_grid=True, ax=None, **kwargs): return ax + def get_layer_from_z(self, pts, cellids): + """Method to handle z values for points. -def _polygon_patch(polygon, **kwargs): - from matplotlib.patches import PathPatch - from matplotlib.path import Path + Parameters + ---------- + pts : shapely.geometry + points geometry + cellids : array_like + array of cellids + Returns + ------- + cellids : array_like + array of cellids with z values handled + """ + + def valid_mask(v): + if isinstance(v, tuple): + return True + else: + return not np.isnan(v) + + z_arr = np.atleast_1d(shapely.get_z(pts)) + if self.mfgrid.grid_type == "structured": + mask_valid = list(map(valid_mask, cellids)) + row, col = list(zip(*cellids[mask_valid])) + surface_elevations = self.mfgrid.top_botm[:, row, col] + elif self.mfgrid.grid_type == "vertex": + mask_valid = ~np.isnan(cellids) + surface_elevations = self.mfgrid.top_botm[:, cellids[mask_valid]] + else: + raise NotImplementedError( + "get_layer_from_z() is only implemented for " + "structured and vertex grids." + ) + zb = surface_elevations < z_arr[mask_valid] + mask_above = zb.all(axis=0) + mask_below = (~zb).all(axis=0) + laypos = (np.nanargmax(zb, axis=0) - 1).astype(float) + laypos[mask_above] = np.inf + laypos[mask_below] = -np.inf + laypos_full = np.full_like(z_arr, np.nan, dtype=float) + laypos_full[mask_valid] = laypos + return laypos_full + + +def _polygon_patch(polygon, **kwargs): patch = PathPatch( Path.make_compound_path( Path(np.asarray(polygon.exterior.coords)[:, :2]), From cbcf9c7aba6998fdbe9b68114e207824f4256cb5 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 11:35:07 +0100 Subject: [PATCH 07/16] remove deprecated class from util init --- flopy/utils/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flopy/utils/__init__.py b/flopy/utils/__init__.py index 3481b04e71..bd15a59b86 100644 --- a/flopy/utils/__init__.py +++ b/flopy/utils/__init__.py @@ -29,7 +29,7 @@ from .formattedfile import FormattedHeadFile get_modflow = get_modflow_module.run_main -from .gridintersect import GridIntersect, ModflowGridIndices +from .gridintersect import GridIntersect from .mflistfile import ( Mf6ListBudget, MfListBudget, From bbc5ac31dc0c351f7b551313dbdd403de04fee50 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 11:38:11 +0100 Subject: [PATCH 08/16] test(GridIntersect): update tests - remove tests for deprecated methods - add array tests for intersects() and points_to_cellids() - activate+improve tests for 3D points --- autotest/test_gridintersect.py | 1531 ++++++++++++++++++++------------ 1 file changed, 971 insertions(+), 560 deletions(-) diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py index 15ec8b494b..f094bd76e5 100644 --- a/autotest/test_gridintersect.py +++ b/autotest/test_gridintersect.py @@ -8,9 +8,8 @@ from flopy.utils.gridintersect import GridIntersect from flopy.utils.triangle import Triangle -# TODO: remove all structured tests in v3.10.0, see TODO's in the tests - if has_pkg("shapely", strict=True): + from shapely import linestrings, points, polygons from shapely.geometry import ( LineString, MultiLineString, @@ -117,92 +116,55 @@ def get_rect_vertex_grid(angrot=0.0, xyoffset=0.0): return tgr -# %% test point structured +# %% test point structured shapely @requires_pkg("shapely") def test_rect_grid_3d_point_outside(): botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape((2, 2, 2)) - gr = get_rect_grid(top=np.ones(4).reshape((2, 2)), botm=botm) - ix = GridIntersect(gr, method="vertex") - result = ix.intersect(Point(25.0, 25.0, 0.5)) + gr = get_rect_grid(top=2 * np.ones(4).reshape((2, 2)), botm=botm) + ix = GridIntersect(gr) + result = ix.intersect(Point(25.0, 25.0, 0.5), handle_z="ignore") assert len(result) == 0 - - -# TODO: fix 3D point tests to work when above or below grid -# @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)) -# gr = get_rect_grid(top=np.ones(4).reshape((2, 2)), botm=botm) -# ix = GridIntersect(gr, method="vertex") -# result = ix.intersect(Point(2.0, 2.0, 0.2)) -# assert result.cellids[0] == (1, 0) - - -# @requires_pkg("shapely") -# def test_rect_grid_3d_point_above(): -# botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape((2, 2, 2)) -# gr = get_rect_grid(top=np.ones(4).reshape((2, 2)), botm=botm) -# ix = GridIntersect(gr, method="vertex") -# result = ix.intersect(Point(2.0, 2.0, 2.0)) -# assert len(result) == 0 - - -@requires_pkg("shapely") -def test_rect_grid_point_outside(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - # use GeoSpatialUtil to convert to shapely geometry - result = ix.intersect((25.0, 25.0), shapetype="point") + result = ix.intersect(Point(25.0, 25.0, 0.5), handle_z="drop") + assert len(result) == 0 + result = ix.intersect(Point(25.0, 25.0, 0.5), handle_z="return") assert len(result) == 0 @requires_pkg("shapely") -def test_rect_grid_point_on_outer_boundary(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(Point(20.0, 10.0)) - assert len(result) == 1 - assert np.all(result.cellids[0] == (0, 1)) - - -@requires_pkg("shapely") -def test_rect_grid_point_on_inner_boundary(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(Point(10.0, 10.0)) - assert len(result) == 1 - assert np.all(result.cellids[0] == (0, 0)) - - -@requires_pkg("shapely") -def test_rect_grid_multipoint_in_one_cell(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(2.0, 2.0)])) - assert len(result) == 1 +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)) + gr = get_rect_grid(top=2 * np.ones(4).reshape((2, 2)), botm=botm) + ix = GridIntersect(gr) + result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z="ignore") + assert result.cellids[0] == (1, 0) + result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z="drop") assert result.cellids[0] == (1, 0) + result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z="return") + assert result.cellids[0] == (1, 0) + assert result.layer[0] == 2.0 # returned as float to allow +/-inf @requires_pkg("shapely") -def test_rect_grid_multipoint_in_multiple_cells(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(12.0, 12.0)])) - assert len(result) == 2 +def test_rect_grid_3d_point_above(): + botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape((2, 2, 2)) + gr = get_rect_grid(top=2 * np.ones(4).reshape((2, 2)), botm=botm) + ix = GridIntersect(gr) + result = ix.intersect(Point(2.0, 2.0, 10.0), handle_z="ignore") + assert len(result) == 1 assert result.cellids[0] == (1, 0) - assert result.cellids[1] == (0, 1) + result = ix.intersect(Point(2.0, 2.0, 10.0), handle_z="drop") + assert len(result) == 0 + result = ix.intersect(Point(2.0, 2.0, 10.0), handle_z="return") + assert len(result) == 1 + assert np.isfinite(result.layer[0]) is np.False_ # %% test point shapely @@ -210,9 +172,9 @@ def test_rect_grid_multipoint_in_multiple_cells(): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_point_outside_shapely(rtree): +def test_rect_grid_point_outside(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) # use GeoSpatialUtil to convert to shapely geometry result = ix.intersect((25.0, 25.0), shapetype="point") assert len(result) == 0 @@ -220,9 +182,9 @@ def test_rect_grid_point_outside_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_point_on_outer_boundary_shapely(rtree): +def test_rect_grid_point_on_outer_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(Point(20.0, 10.0)) assert len(result) == 1 assert np.all(result.cellids[0] == (0, 1)) @@ -230,9 +192,9 @@ def test_rect_grid_point_on_outer_boundary_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_point_on_inner_boundary_shapely(rtree): +def test_rect_grid_point_on_inner_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(Point(10.0, 10.0)) assert len(result) == 1 assert np.all(result.cellids[0] == (0, 0)) @@ -240,9 +202,9 @@ def test_rect_grid_point_on_inner_boundary_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_vertex_grid_point_in_one_cell_shapely(rtree): +def test_rect_vertex_grid_point_in_one_cell(rtree): gr = get_rect_vertex_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(Point(4.0, 4.0)) assert len(result) == 1 assert result.cellids[0] == 0 @@ -259,9 +221,9 @@ def test_rect_vertex_grid_point_in_one_cell_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_multipoint_in_one_cell_shapely(rtree): +def test_rect_grid_multipoint_in_one_cell(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(2.0, 2.0)])) assert len(result) == 1 assert result.cellids[0] == (1, 0) @@ -269,9 +231,9 @@ def test_rect_grid_multipoint_in_one_cell_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_multipoint_in_multiple_cells_shapely(rtree): +def test_rect_grid_multipoint_in_multiple_cells(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(12.0, 12.0)])) assert len(result) == 2 assert result.cellids[0] == (0, 1) @@ -341,20 +303,8 @@ def test_tri_grid_multipoint_in_multiple_cells(rtree): @requires_pkg("shapely") @rtree_toggle def test_rect_grid_point_on_all_vertices_return_all_ix(rtree): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured", rtree=rtree) - n_intersections = [1, 2, 1, 2, 4, 2, 1, 2, 1] - for v, n in zip(gr.verts, n_intersections): - r = ix.intersect(Point(*v), return_all_intersections=True) - assert len(r) == n - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_point_on_all_vertices_return_all_ix_shapely(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) n_intersections = [1, 2, 1, 2, 4, 2, 1, 2, 1] for v, n in zip(gr.verts, n_intersections): r = ix.intersect(Point(*v), return_all_intersections=True) @@ -363,157 +313,32 @@ def test_rect_grid_point_on_all_vertices_return_all_ix_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_tri_grid_points_on_all_vertices_return_all_ix_shapely(rtree): +def test_tri_grid_points_on_all_vertices_return_all_ix(rtree): gr = get_tri_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) n_intersections = [2, 2, 2, 2, 8, 2, 2, 2, 2] for v, n in zip(gr.verts, n_intersections): r = ix.intersect(Point(*v), return_all_intersections=True) assert len(r) == n -# %% test linestring structured - - -@requires_pkg("shapely") -def test_rect_grid_linestring_outside(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(LineString([(25.0, 25.0), (21.0, 5.0)])) - assert len(result) == 0 - - -@requires_pkg("shapely") -def test_rect_grid_linestring_in_2cells(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(LineString([(5.0, 5.0), (15.0, 5.0)])) - assert len(result) == 2 - assert result.lengths.sum() == 10.0 - assert result.cellids[0] == (1, 0) - assert result.cellids[1] == (1, 1) - - -@requires_pkg("shapely") -def test_rect_grid_linestring_on_outer_boundary(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(LineString([(15.0, 20.0), (5.0, 20.0)])) - assert len(result) == 2 - assert result.lengths.sum() == 10.0 - assert result.cellids[1] == (0, 0) - assert result.cellids[0] == (0, 1) - - -@requires_pkg("shapely") -def test_rect_grid_linestring_on_inner_boundary(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(LineString([(5.0, 10.0), (15.0, 10.0)])) - assert len(result) == 2 - assert result.lengths.sum() == 10.0 - assert result.cellids[0] == (0, 0) - assert result.cellids[1] == (0, 1) - - -@requires_pkg("shapely") -def test_rect_grid_multilinestring_in_one_cell(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect( - MultiLineString( - [LineString([(1.0, 1), (9.0, 1.0)]), LineString([(1.0, 9.0), (9.0, 9.0)])] - ) - ) - assert len(result) == 1 - assert result.lengths == 16.0 - assert result.cellids[0] == (1, 0) - - -@requires_pkg("shapely") -def test_rect_grid_multilinestring_in_multiple_cells(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect( - MultiLineString( - [ - LineString([(20.0, 0.0), (7.5, 12.0), (2.5, 7.0), (0.0, 4.5)]), - LineString([(5.0, 19.0), (2.5, 7.0)]), - ] - ) - ) - assert len(result) == 3 - assert np.allclose(sum(result.lengths), 40.19197584109293) - - -@requires_pkg("shapely") -def test_rect_grid_linestring_in_and_out_of_cell(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(LineString([(5.0, 9), (15.0, 5.0), (5.0, 1.0)])) - assert len(result) == 2 - assert result.cellids[0] == (1, 0) - assert result.cellids[1] == (1, 1) - assert np.allclose(result.lengths.sum(), 21.540659228538015) - - -@requires_pkg("shapely") -def test_rect_grid_linestring_in_and_out_of_cell2(): - # TODO: remove in 3.10.0 - 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)])) - assert len(result) == 3 - - -@requires_pkg("shapely") -def test_rect_grid_linestrings_on_boundaries_return_all_ix(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - x, y = ix._rect_grid_to_geoms_cellids()[0][0].exterior.xy - n_intersections = [1, 2, 2, 1] - for i in range(4): - ls = LineString([(x[i], y[i]), (x[i + 1], y[i + 1])]) - r = ix.intersect(ls, return_all_intersections=True) - assert len(r) == n_intersections[i] - - -@requires_pkg("shapely") -def test_rect_grid_linestring_starting_on_vertex(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - result = ix.intersect(LineString([(10.0, 10.0), (15.0, 5.0)])) - assert len(result) == 1 - assert np.allclose(result.lengths.sum(), np.sqrt(50)) - assert result.cellids[0] == (1, 1) - - # %% test linestring shapely @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_outside_shapely(rtree): +def test_rect_grid_linestring_outside(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(LineString([(25.0, 25.0), (21.0, 5.0)])) assert len(result) == 0 @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_in_2cells_shapely(rtree): +def test_rect_grid_linestring_in_2cells(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(LineString([(5.0, 5.0), (15.0, 5.0)])) assert len(result) == 2 assert result.lengths.sum() == 10.0 @@ -523,9 +348,9 @@ def test_rect_grid_linestring_in_2cells_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_on_outer_boundary_shapely(rtree): +def test_rect_grid_linestring_on_outer_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(LineString([(15.0, 20.0), (5.0, 20.0)])) assert len(result) == 2 assert result.lengths.sum() == 10.0 @@ -535,9 +360,9 @@ def test_rect_grid_linestring_on_outer_boundary_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_on_inner_boundary_shapely(rtree): +def test_rect_grid_linestring_on_inner_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(LineString([(5.0, 10.0), (15.0, 10.0)])) assert len(result) == 2 assert result.lengths.sum() == 10.0 @@ -547,9 +372,9 @@ def test_rect_grid_linestring_on_inner_boundary_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_multilinestring_in_one_cell_shapely(rtree): +def test_rect_grid_multilinestring_in_one_cell(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect( MultiLineString( [LineString([(1.0, 1), (9.0, 1.0)]), LineString([(1.0, 9.0), (9.0, 9.0)])] @@ -562,9 +387,9 @@ def test_rect_grid_multilinestring_in_one_cell_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_multilinestring_in_multiple_cells_shapely(rtree): +def test_rect_grid_multilinestring_in_multiple_cells(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect( MultiLineString( [ @@ -579,9 +404,9 @@ def test_rect_grid_multilinestring_in_multiple_cells_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_in_and_out_of_cell_shapely(rtree): +def test_rect_grid_linestring_in_and_out_of_cell(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(LineString([(5.0, 9), (15.0, 5.0), (5.0, 1.0)])) assert len(result) == 2 assert result.cellids[0] == (1, 0) @@ -590,17 +415,17 @@ def test_rect_grid_linestring_in_and_out_of_cell_shapely(rtree): @requires_pkg("shapely") -def test_rect_grid_linestring_in_and_out_of_cell2_shapely(): +def test_rect_grid_linestring_in_and_out_of_cell2(): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex") + ix = GridIntersect(gr) result = ix.intersect(LineString([(5, 15), (5.0, 9), (15.0, 5.0), (5.0, 1.0)])) assert len(result) == 3 @requires_pkg("shapely") -def test_rect_grid_linestring_starting_on_vertex_shapely(): +def test_rect_grid_linestring_starting_on_vertex(): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex") + ix = GridIntersect(gr) result = ix.intersect(LineString([(10.0, 10.0), (15.0, 5.0)])) assert len(result) == 1 assert np.allclose(result.lengths.sum(), np.sqrt(50)) @@ -609,9 +434,9 @@ def test_rect_grid_linestring_starting_on_vertex_shapely(): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestrings_on_boundaries_return_all_ix_shapely(rtree): +def test_rect_grid_linestrings_on_boundaries_return_all_ix(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) x, y = ix._rect_grid_to_geoms_cellids()[0][0].exterior.xy n_intersections = [1, 2, 2, 1] for i in range(4): @@ -622,9 +447,9 @@ def test_rect_grid_linestrings_on_boundaries_return_all_ix_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_cell_boundary_shapely(rtree): +def test_rect_grid_linestring_cell_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) ls = LineString(ix._rect_grid_to_geoms_cellids()[0][0].exterior.coords) r = ix.intersect(ls, return_all_intersections=False) assert len(r) == 1 @@ -632,9 +457,9 @@ def test_rect_grid_linestring_cell_boundary_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_rect_grid_linestring_cell_boundary_return_all_ix_shapely(rtree): +def test_rect_grid_linestring_cell_boundary_return_all_ix(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) ls = LineString(ix._rect_grid_to_geoms_cellids()[0][0].exterior.coords) r = ix.intersect(ls, return_all_intersections=True) assert len(r) == 3 @@ -733,7 +558,7 @@ def test_tri_grid_multilinestring_in_multiple_cells(rtree): @rtree_toggle def test_tri_grid_linestrings_on_boundaries_return_all_ix(rtree): tgr = get_tri_grid() - ix = GridIntersect(tgr, method="vertex", rtree=rtree) + ix = GridIntersect(tgr, rtree=rtree) x, y = ix._vtx_grid_to_geoms_cellids()[0][0].exterior.xy n_intersections = [2, 1, 2] for i in range(len(x) - 1): @@ -744,9 +569,9 @@ def test_tri_grid_linestrings_on_boundaries_return_all_ix(rtree): @requires_pkg("shapely") @rtree_toggle -def test_tri_grid_linestring_cell_boundary_shapely(rtree): +def test_tri_grid_linestring_cell_boundary(rtree): tgr = get_tri_grid() - ix = GridIntersect(tgr, method="vertex", rtree=rtree) + ix = GridIntersect(tgr, rtree=rtree) ls = LineString(ix._vtx_grid_to_geoms_cellids()[0][0].exterior.coords) r = ix.intersect(ls, return_all_intersections=False) assert len(r) == 1 @@ -754,9 +579,9 @@ def test_tri_grid_linestring_cell_boundary_shapely(rtree): @requires_pkg("shapely") @rtree_toggle -def test_tri_grid_linestring_cell_boundary_return_all_ix_shapely(rtree): +def test_tri_grid_linestring_cell_boundary_return_all_ix(rtree): tgr = get_tri_grid() - ix = GridIntersect(tgr, method="vertex", rtree=rtree) + ix = GridIntersect(tgr, rtree=rtree) ls = LineString(ix._vtx_grid_to_geoms_cellids()[0][0].exterior.coords) r = ix.intersect(ls, return_all_intersections=True) assert len(r) == 3 @@ -765,7 +590,7 @@ def test_tri_grid_linestring_cell_boundary_return_all_ix_shapely(rtree): @requires_pkg("shapely") def test_rect_vertex_grid_linestring_geomcollection(): gr = get_rect_vertex_grid() - ix = GridIntersect(gr, method="vertex") + ix = GridIntersect(gr) ls = LineString( [ (20.0, 0.0), @@ -782,33 +607,72 @@ def test_rect_vertex_grid_linestring_geomcollection(): assert np.allclose(result.lengths.sum(), ls.length) -# %% test polygon structured +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_polygon_contains_centroid(rtree): + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = Polygon( + [(6.0, 5.0), (4.0, 16.0), (25.0, 14.0), (25.0, -5.0), (6.0, -5.0)], + holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], + ) + result = ix.intersect(p, contains_centroid=True) + assert len(result) == 1 + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_polygon_min_area(rtree): + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = Polygon( + [(5.0, 5.0), (5.0, 15.0), (25.0, 15.0), (25.0, -5.0), (5.0, -5.0)], + holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], + ) + result = ix.intersect(p, min_area_fraction=0.4) + assert len(result) == 2 + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_polygon_centroid_and_min_area(rtree): + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = Polygon( + [(5.0, 5.0), (5.0, 15.0), (25.0, 14.0), (25.0, -5.0), (5.0, -5.0)], + holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], + ) + result = ix.intersect(p, min_area_fraction=0.35, contains_centroid=True) + assert len(result) == 1 + + +# %% test polygon shapely @requires_pkg("shapely") -def test_rect_grid_polygon_outside(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_polygon_outside(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect(Polygon([(21.0, 11.0), (23.0, 17.0), (25.0, 11.0)])) assert len(result) == 0 @requires_pkg("shapely") -def test_rect_grid_polygon_in_2cells(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_polygon_in_2cells(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + ix = GridIntersect(gr, rtree=rtree) 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 @requires_pkg("shapely") -def test_rect_grid_polygon_on_outer_boundary(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_polygon_on_outer_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + ix = GridIntersect(gr, rtree=rtree) result = ix.intersect( Polygon([(20.0, 5.0), (25.0, 5.0), (25.0, 15.0), (20.0, 15.0)]) ) @@ -817,9 +681,8 @@ def test_rect_grid_polygon_on_outer_boundary(): @requires_pkg("shapely") def test_rect_grid_polygon_running_along_boundary(): - # TODO: remove in 3.10.0 gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + ix = GridIntersect(gr) result = ix.intersect( Polygon( [(5.0, 5.0), (5.0, 10.0), (9.0, 10.0), (9.0, 15.0), (1.0, 15.0), (1.0, 5.0)] @@ -828,107 +691,20 @@ def test_rect_grid_polygon_running_along_boundary(): @requires_pkg("shapely") -def test_rect_grid_polygon_on_inner_boundary(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_polygon_on_inner_boundary(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + 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)])) assert len(result) == 2 assert result.areas.sum() == 50.0 - fig, ax = plt.subplots(1, 1, figsize=(8, 8)) - gr.plot(ax=ax) - ix.plot_polygon(result, ax=ax) - # plt.show() - - -@requires_pkg("shapely") -def test_rect_grid_polygon_multiple_polygons(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - p = Polygon( - [ - (0, 0), - (0, 10), - (4, 10), - (4, 0), - (6, 0), - (6, 10), - (9, 10), - (9, -1), - (0, -1), - (0, 0), - ] - ) - - ix = GridIntersect(gr, method="structured") - result = ix.intersect(p) - - fig, ax = plt.subplots(1, 1, figsize=(8, 8)) - gr.plot(ax=ax) - ix.plot_polygon(result, ax=ax) - # plt.show() - - -@requires_pkg("shapely") -def test_rect_grid_multiple_disjoint_polygons_on_inner_boundaries(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - p1 = Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]) - p2 = Polygon([(5.0, 17.5), (15.0, 17.5), (15.0, 12.5), (5.0, 12.5)]) - result = ix.intersect(MultiPolygon([p1, p2])) - - assert len(result) == 4 - assert result.areas.sum() == 100.0 - - fig, ax = plt.subplots(1, 1, figsize=(8, 8)) - gr.plot(ax=ax) - ix.plot_polygon(result, ax=ax) - # plt.show() - - -@requires_pkg("shapely") -@pytest.mark.parametrize("transform", [True, False]) -def test_rect_grid_polygon_reintersects_cell(transform): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - if transform: - gr.set_coord_info(xoff=1, yoff=1, angrot=10.5) - - ix = GridIntersect(gr, method="structured") - p1 = Polygon( - [ - (x, y + 3) - for x, y in [ - (2.5, 2.5), - (2.5, 10.0), - (10.0, 10.0), - (10.0, 2.5), - (7.5, 2.5), - (7.5, 7.5), - (5.0, 7.5), - (5.0, 2.5), - ] - ] - ) - p2 = Polygon([(1, 1), (1, 2), (2, 2), (2, 1)]) - result = ix.intersect(MultiPolygon([p1, p2])) - - assert len(result) == 4 if transform else 2 - assert np.isclose(result.areas.sum(), 44.65733 if transform else 44.75) - - fig, ax = plt.subplots(1, 1, figsize=(8, 8)) - gr.plot(ax=ax) - ix.plot_polygon(result, ax=ax) - # plt.show() - @requires_pkg("shapely") -def test_rect_grid_multipolygon_in_one_cell(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_multipolygon_in_one_cell(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + ix = GridIntersect(gr, rtree=rtree) p1 = Polygon([(1.0, 1.0), (8.0, 1.0), (8.0, 3.0), (1.0, 3.0)]) p2 = Polygon([(1.0, 9.0), (8.0, 9.0), (8.0, 7.0), (1.0, 7.0)]) p = MultiPolygon([p1, p2]) @@ -938,10 +714,10 @@ def test_rect_grid_multipolygon_in_one_cell(): @requires_pkg("shapely") -def test_rect_grid_multipolygon_in_multiple_cells(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_multipolygon_in_multiple_cells(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") + ix = GridIntersect(gr, rtree=rtree) p1 = Polygon([(1.0, 1.0), (19.0, 1.0), (19.0, 3.0), (1.0, 3.0)]) p2 = Polygon([(1.0, 9.0), (19.0, 9.0), (19.0, 7.0), (1.0, 7.0)]) p = MultiPolygon([p1, p2]) @@ -951,145 +727,10 @@ def test_rect_grid_multipolygon_in_multiple_cells(): @requires_pkg("shapely") -def test_rect_grid_polygon_with_hole(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_rect_grid_polygon_with_hole(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="structured") - p = Polygon( - [(5.0, 5.0), (5.0, 15.0), (25.0, 15.0), (25.0, -5.0), (5.0, -5.0)], - holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], - ) - result = ix.intersect(p) - assert len(result) == 3 - assert result.areas.sum() == 104.0 - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_polygon_contains_centroid(rtree): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, rtree=rtree) - p = Polygon( - [(6.0, 5.0), (4.0, 16.0), (25.0, 14.0), (25.0, -5.0), (6.0, -5.0)], - holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], - ) - result = ix.intersect(p, contains_centroid=True) - assert len(result) == 1 - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_polygon_min_area(rtree): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr, rtree=rtree) - p = Polygon( - [(5.0, 5.0), (5.0, 15.0), (25.0, 15.0), (25.0, -5.0), (5.0, -5.0)], - holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], - ) - result = ix.intersect(p, min_area_fraction=0.4) - assert len(result) == 2 - - -@requires_pkg("shapely") -def test_rect_grid_polygon_centroid_and_min_area(): - # TODO: remove in 3.10.0 - gr = get_rect_grid() - ix = GridIntersect(gr) - p = Polygon( - [(5.0, 5.0), (5.0, 15.0), (25.0, 14.0), (25.0, -5.0), (5.0, -5.0)], - holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], - ) - result = ix.intersect(p, min_area_fraction=0.35, contains_centroid=True) - assert len(result) == 1 - - -# %% test polygon shapely - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_polygon_outside_shapely(rtree): - gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) - result = ix.intersect(Polygon([(21.0, 11.0), (23.0, 17.0), (25.0, 11.0)])) - assert len(result) == 0 - - -@requires_pkg("shapely") -@rtree_toggle -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)])) - assert len(result) == 2 - assert result.areas.sum() == 50.0 - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_polygon_on_outer_boundary_shapely(rtree): - gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) - result = ix.intersect( - Polygon([(20.0, 5.0), (25.0, 5.0), (25.0, 15.0), (20.0, 15.0)]) - ) - assert len(result) == 0 - - -@requires_pkg("shapely") -def test_rect_grid_polygon_running_along_boundary_shapely(): - gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex") - result = ix.intersect( - Polygon( - [(5.0, 5.0), (5.0, 10.0), (9.0, 10.0), (9.0, 15.0), (1.0, 15.0), (1.0, 5.0)] - ) - ) - - -@requires_pkg("shapely") -@rtree_toggle -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)])) - assert len(result) == 2 - assert result.areas.sum() == 50.0 - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_multipolygon_in_one_cell_shapely(rtree): - gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) - p1 = Polygon([(1.0, 1.0), (8.0, 1.0), (8.0, 3.0), (1.0, 3.0)]) - p2 = Polygon([(1.0, 9.0), (8.0, 9.0), (8.0, 7.0), (1.0, 7.0)]) - p = MultiPolygon([p1, p2]) - result = ix.intersect(p) - assert len(result) == 1 - assert result.areas.sum() == 28.0 - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_multipolygon_in_multiple_cells_shapely(rtree): - gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) - p1 = Polygon([(1.0, 1.0), (19.0, 1.0), (19.0, 3.0), (1.0, 3.0)]) - p2 = Polygon([(1.0, 9.0), (19.0, 9.0), (19.0, 7.0), (1.0, 7.0)]) - p = MultiPolygon([p1, p2]) - result = ix.intersect(p) - assert len(result) == 2 - assert result.areas.sum() == 72.0 - - -@requires_pkg("shapely") -@rtree_toggle -def test_rect_grid_polygon_with_hole_shapely(rtree): - gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) p = Polygon( [(5.0, 5.0), (5.0, 15.0), (25.0, 15.0), (25.0, -5.0), (5.0, -5.0)], holes=[[(9.0, -1), (9, 11), (21, 11), (21, -1)]], @@ -1103,7 +744,7 @@ def test_rect_grid_polygon_with_hole_shapely(rtree): @rtree_toggle def test_rect_grid_polygon_in_edge_in_cell(rtree): gr = get_rect_grid() - ix = GridIntersect(gr, method="vertex", rtree=rtree) + ix = GridIntersect(gr, rtree=rtree) p = Polygon( [(0.0, 5.0), (3.0, 0.0), (7.0, 0.0), (10.0, 5.0), (10.0, -1.0), (0.0, -1.0)] ) @@ -1240,36 +881,36 @@ def test_tri_grid_polygon_contains_centroid(rtree): @requires_pkg("shapely") -def test_point_offset_rot_structured_grid(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_point_offset_rot_structured_grid(rtree): sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) p = Point(10.0, 10 + np.sqrt(200.0)) - ix = GridIntersect(sgr, method="structured") + ix = GridIntersect(sgr, rtree=rtree) result = ix.intersect(p) assert len(result) == 1 # check empty result when using local model coords - ix = GridIntersect(sgr, method="structured", local=True) + ix = GridIntersect(sgr, rtree=rtree, local=True) result = ix.intersect(p) assert len(result) == 0 @requires_pkg("shapely") -def test_linestring_offset_rot_structured_grid(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_linestring_offset_rot_structured_grid(rtree): sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) - ls = LineString([(5, 25), (15, 25)]) - ix = GridIntersect(sgr, method="structured") + ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))]) + ix = GridIntersect(sgr, rtree=rtree) result = ix.intersect(ls) - assert len(result) == 3 + assert len(result) == 2 # check empty result when using local model coords - ix = GridIntersect(sgr, method="structured", local=True) + ix = GridIntersect(sgr, rtree=rtree, local=True) result = ix.intersect(ls) assert len(result) == 0 @requires_pkg("shapely") -def test_polygon_offset_rot_structured_grid(): - # TODO: remove in 3.10.0 +@rtree_toggle +def test_polygon_offset_rot_structured_grid(rtree): sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) p = Polygon( [ @@ -1279,47 +920,47 @@ def test_polygon_offset_rot_structured_grid(): (5, 10.0 + 1.5 * np.sqrt(200.0)), ] ) - ix = GridIntersect(sgr, method="structured") + ix = GridIntersect(sgr, rtree=rtree) result = ix.intersect(p) assert len(result) == 3 # check empty result when using local model coords - ix = GridIntersect(sgr, method="structured", local=True) + ix = GridIntersect(sgr, rtree=rtree, local=True) result = ix.intersect(p) assert len(result) == 0 @requires_pkg("shapely") @rtree_toggle -def test_point_offset_rot_structured_grid_shapely(rtree): - sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) +def test_point_offset_rot_vertex_grid(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) + ix = GridIntersect(sgr, 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) + ix = GridIntersect(sgr, rtree=rtree, local=True) result = ix.intersect(p) assert len(result) == 0 @requires_pkg("shapely") @rtree_toggle -def test_linestring_offset_rot_structured_grid_shapely(rtree): - sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) +def test_linestring_offset_rot_vertex_grid(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) + ix = GridIntersect(sgr, 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) + ix = GridIntersect(sgr, rtree=rtree, local=True) result = ix.intersect(ls) assert len(result) == 0 @requires_pkg("shapely") @rtree_toggle -def test_polygon_offset_rot_structured_grid_shapely(rtree): - sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) +def test_polygon_offset_rot_vertex_grid(rtree): + sgr = get_rect_vertex_grid(angrot=45.0, xyoffset=10.0) p = Polygon( [ (5, 10.0 + np.sqrt(200.0)), @@ -1328,59 +969,829 @@ def test_polygon_offset_rot_structured_grid_shapely(rtree): (5, 10.0 + 1.5 * np.sqrt(200.0)), ] ) - ix = GridIntersect(sgr, method="vertex", rtree=rtree) + ix = GridIntersect(sgr, 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) + ix = GridIntersect(sgr, rtree=rtree, local=True) result = ix.intersect(p) assert len(result) == 0 +# %% array-based inputs - structured grid points + + @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) +def test_rect_grid_single_point_array_inside(rtree): + """Single point in array inside, returns single intersection.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([1.0], [1.0]) + result = ix.intersects(pts) 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 result.cellids[0] == (1, 0) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_single_point_array_outside(rtree): + """Single point in array outside, returns empty result.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([25.0], [25.0]) + result = ix.intersects(pts) assert len(result) == 0 +@requires_pkg("shapely") +def test_rect_grid_single_point_array_outside_points_to_cellids(): + """Single point in array outside in points_to_cellids, returns single nan result.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + pts = points([25.0], [25.0]) + result = ix.points_to_cellids(pts) + assert len(result) == 1 + assert np.isnan(result.cellids[0]) + + +@requires_pkg("shapely") +def test_rect_grid_single_point_array_on_boundary_points_to_cellids(): + """Single point in array on boundary, returns single result.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + pts = points([10.0], [20.0]) + result = ix.points_to_cellids(pts) + assert len(result) == 1 + assert result.cellids[0] == (0, 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) +def test_rect_grid_single_point_array_on_boundary(rtree): + """Single point in array on boundary, returns multiple results.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([10.0], [20.0]) + result = ix.intersects(pts) 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 result.cellids[0] == (0, 0) + assert result.cellids[1] == (0, 1) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_points_array_in_one_cell(): + """Multiple points in array in one cell, returns results per point.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + pts = points([1.0, 5.0], [2.0, 5.0]) + result = ix.intersects(pts) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 0) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_points_array_in_multiple_cells(): + """Multiple points in array in multiple cells, returns results per point.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + pts = points([1.0, 15.0], [2.0, 15.0]) + result = ix.intersects(pts) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (0, 1) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_points_array_inside_and_outside(): + """Multiple points in array inside and outside, returns one result.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + pts = points([1.0, 25.0], [2.0, 25.0]) + result = ix.intersects(pts) + assert len(result) == 1 + assert result.cellids[0] == (1, 0) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_points_array_inside_and_outside_points_to_cellids(): + """Multiple points in array inside/outside, returns one result and one nan.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + pts = points([1.0, 25.0], [2.0, 25.0]) + result = ix.points_to_cellids(pts) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert np.isnan(result.cellids[1]) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_points_array_with_z_points_to_cellids(): + gr = get_rect_grid( + top=np.ones(4).reshape((2, 2)), botm=np.zeros(4).reshape((1, 2, 2)) + ) + ix = GridIntersect(gr) + pts = points([1.0, 25.0], [2.0, 25.0], [10.0, 0.5]) + result = ix.points_to_cellids(pts, handle_z="ignore") + assert result.cellids[0] == (1, 0) + assert np.isnan(result.cellids[1]) + result = ix.points_to_cellids(pts, handle_z="return") + assert ~np.isfinite(result.layer[0]) + assert np.isnan(result.cellids[1]) + result = ix.points_to_cellids(pts, handle_z="drop") assert len(result) == 0 +# %% array-based input - structured grid linestrings + + @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) +def test_rect_grid_single_linestring_array_in_one_cell(rtree): + """Single linestring in array in 1 cell, returns single intersection.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(5.0, 5.0), (7.5, 5.0)]]) + result = ix.intersects(ls) + assert len(result) == 1 + assert result.cellids[0] == (1, 0) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_single_linestring_array_in_two_cells(rtree): + """Single linestring in array in 2 cells, returns multiple intersections.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(5.0, 5.0), (15.0, 5.0)]]) + result = ix.intersects(ls) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_single_linestring_array_outside(rtree): + """Single linestring in array outside, returns empty result.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(25.0, 5.0), (35.0, 5.0)]]) + result = ix.intersects(ls) assert len(result) == 0 + + +@requires_pkg("shapely") +def test_rect_grid_multiple_linestring_array_in_multiple_cells(): + """Multiple linestrings in array; returns results per linestring, + multiple cellids per linestring.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + ls = linestrings([[(5.0, 5.0), (15.0, 5.0)], [(5.0, 15.0), (15.0, 15.0)]]) + result = ix.intersects(ls) + assert len(result) == 4 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + assert result.cellids[2] == (0, 0) + assert result.cellids[3] == (0, 1) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_linestring_array_inside_outside(): + """Multiple linestrings in array inside/outside; returns results per linestring, + multiple cellids per linestring.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + ls = linestrings([[(5.0, 5.0), (15.0, 5.0)], [(25.0, 15.0), (35.0, 15.0)]]) + result = ix.intersects(ls) + assert len(result) == 2 + assert (result.shp_ids == 0).all() + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + + +# %% array-based input - structured grid polygons + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_single_polygon_array_in_one_cell(rtree): + """Single polygon in array inside, returns single intersection.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons([[(2.5, 5.0), (7.5, 5.0), (7.5, 7.5), (2.5, 7.5)]]) + result = ix.intersects(p) + assert len(result) == 1 + assert result.cellids[0] == (1, 0) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_single_polygon_array_in_two_cells(rtree): + """Single polygon in array in 2 cells, returns multiple intersections.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons([[(2.5, 5.0), (15, 5.0), (15, 7.5), (2.5, 7.5)]]) + result = ix.intersects(p) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_single_polygon_array_outside(rtree): + """Single polygon in array outside, returns empty result.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons([[(25, 5.0), (75, 5.0), (75, 7.5), (25, 7.5)]]) + result = ix.intersects(p) + assert len(result) == 0 + + +@requires_pkg("shapely") +def test_rect_grid_multiple_polygon_array_single_result_per_polygon(): + """Multiple polygons in array; returns results per polygon, + single result per polygon.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + p = polygons( + [ + [(2.5, 5.0), (7.5, 5.0), (7.5, 7.5), (2.5, 7.5)], + [(2.5, 15.0), (7.5, 15.0), (7.5, 17.5), (2.5, 17.5)], + ] + ) + result = ix.intersects(p) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (0, 0) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_polygon_array_multiple_results_per_polygon(): + """Multiple polygons in array; returns results per polygon, + multiple results per polygon.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + p = polygons( + [ + [(2.5, 5.0), (17.5, 5.0), (17.5, 7.5), (2.5, 7.5)], + [(2.5, 15.0), (17.5, 15.0), (17.5, 17.5), (2.5, 17.5)], + ] + ) + result = ix.intersects(p) + assert len(result) == 4 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + assert result.cellids[2] == (0, 0) + assert result.cellids[3] == (0, 1) + + +@requires_pkg("shapely") +def test_rect_grid_multiple_polygon_array_inside_outside(): + """Multiple polygons in array inside/outside; returns results per polygon, + multiple cellids per polygon.""" + gr = get_rect_grid() + ix = GridIntersect(gr) + p = polygons( + [ + [(2.5, 5.0), (7.5, 5.0), (7.5, 7.5), (2.5, 7.5)], + [(25, 15.0), (75, 15.0), (75, 17.5), (25, 17.5)], + ] + ) + result = ix.intersects(p) + assert len(result) == 1 + assert result.cellids[0] == (1, 0) + + +# %% array-based input - structured grid intersection method + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_intersect_single_point_array(rtree): + """Single point in array ok.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([10], [20]) + result = ix.intersect(pts) + assert result.cellids[0] == (0, 0) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_intersect_multiple_points_array(rtree): + """Multiple points in array raises error.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([10, 1], [20, 5]) + with pytest.raises( + ValueError, match="intersect\(\) only accepts arrays containing one" + ): + ix.intersect(pts) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_intersect_single_linestring_array(rtree): + """Single linestring in array ok.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(5.0, 5.0), (15.0, 5.0)]]) + result = ix.intersect(ls) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + assert (result.lengths == 5.0).all() + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_intersect_multiple_linestring_array(rtree): + """Multiple linestrings in array raises error.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings( + [ + [(5.0, 5.0), (15.0, 5.0)], + [(5.0, 15.0), (15.0, 15.0)], + ] + ) + with pytest.raises( + ValueError, match="intersect\(\) only accepts arrays containing one" + ): + ix.intersect(ls) + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_intersect_single_polygon_array(rtree): + """Single polygon in array ok.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons( + [ + [(2.5, 5.0), (17.5, 5.0), (17.5, 7.5), (2.5, 7.5)], + ] + ) + result = ix.intersect(p) + assert len(result) == 2 + assert result.cellids[0] == (1, 0) + assert result.cellids[1] == (1, 1) + assert (result.areas == 18.75).all() + + +@requires_pkg("shapely") +@rtree_toggle +def test_rect_grid_intersect_multiple_polygon_array(rtree): + """Multiple polygons in array input raises error.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons( + [ + [(2.5, 5.0), (17.5, 5.0), (17.5, 7.5), (2.5, 7.5)], + [(2.5, 15.0), (17.5, 15.0), (17.5, 17.5), (2.5, 17.5)], + ] + ) + with pytest.raises( + ValueError, match="intersect\(\) only accepts arrays containing one" + ): + ix.intersect(p) + + +# %% vertex grid points + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_point_array_inside(rtree): + """Single point in array inside returns single intersection.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([9.0], [1.0]) + result = ix.intersects(pts) + assert len(result) == 1 + assert result.cellids[0] == 4 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_point_array_outside(rtree): + """Single point in array outside, returns empty result.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([25.0], [25.0]) + result = ix.intersects(pts) + assert len(result) == 0 + + +@requires_pkg("shapely") +def test_tri_grid_single_point_array_outside_points_to_cellids(): + """Single point in array outside + return_all, returns single nan result.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + pts = points([25.0], [25.0]) + result = ix.points_to_cellids(pts) + assert len(result) == 1 + assert np.isnan(result.cellids[0]) + + +@requires_pkg("shapely") +def test_tri_grid_single_point_array_on_boundary_points_to_cellids(): + """Single point in array on boundary, returns single intersection.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + pts = points([9.0], [1.0]) + result = ix.points_to_cellids(pts) + assert len(result) == 1 + assert result.cellids[0] == 4 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_point_array_on_boundary(rtree): + """Single point in array on boundary, returns multiple intersections.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([5.0], [5.0]) + result = ix.intersects(pts) + assert len(result) == 2 + assert result.cellids[0] == 1 + assert result.cellids[1] == 4 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_points_array_in_one_cell(): + """Multiple points in array in one cell, returns results per point.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + pts = points([9.0, 9.0], [1.0, 8.0]) + result = ix.intersects(pts) + assert len(result) == 2 + assert (result.cellids == 4).all() + + +@requires_pkg("shapely") +def test_tri_grid_multiple_points_array_in_multiple_cells(): + """Multiple points in array in multiple cells, returns results per point.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + pts = points([15.0, 9.0], [3.0, 3.0]) + result = ix.intersects(pts) + assert len(result) == 2 + assert result.cellids[0] == 5 + assert result.cellids[1] == 4 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_points_array_inside_and_outside_points_to_cellids(): + """Multiple points in array inside and outside, returns one result and one nan.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + pts = points([5.0, 25.0], [3.0, 25.0]) + result = ix.points_to_cellids(pts) + assert len(result) == 2 + assert result.cellids[0] == 4 + assert np.isnan(result.cellids[1]) + + +@requires_pkg("shapely") +def test_tri_grid_multiple_points_array_with_z_points_to_cellids(): + gr = get_rect_grid( + top=np.ones(4).reshape((2, 2)), botm=np.zeros(4).reshape((1, 2, 2)) + ) + ix = GridIntersect(gr) + pts = points([1.0, 25.0], [2.0, 25.0], [0.5, 10.0]) + result = ix.points_to_cellids(pts, handle_z="ignore") + assert result.cellids[0] == (1, 0) + assert np.isnan(result.cellids[1]) + result = ix.points_to_cellids(pts, handle_z="return") + assert result.layer[0] == 0.0 + assert np.isnan(result.cellids[1]) + result = ix.points_to_cellids(pts, handle_z="drop") + assert len(result) == 1 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_points_array_inside_and_outside(): + """Multiple points in array inside and outside, returns 1 result.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + pts = points([5.0, 25.0], [3.0, 25.0]) + result = ix.intersects(pts) + assert len(result) == 1 + assert result.cellids[0] == 4 + + +# %% vertex grid linestrings + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_linestring_array_in_one_cell(rtree): + """Single linestring in array in 1 cell, returns single intersection.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(2.0, 1.0), (7.5, 1.0)]]) + result = ix.intersects(ls) + assert len(result) == 1 + assert result.cellids[0] == 4 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_linestring_array_in_two_cells(rtree): + """Single linestring in array in 2 cells, returns multiple intersections.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(2.0, 1.0), (15.0, 1.0)]]) + result = ix.intersects(ls) + assert len(result) == 2 + assert result.cellids[0] == 4 + assert result.cellids[1] == 5 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_linestring_array_outside(rtree): + """Single linestring in array outside, returns empty result.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(25.0, 5.0), (35.0, 5.0)]]) + result = ix.intersects(ls) + assert len(result) == 0 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_linestring_array_in_multiple_cells(): + """Multiple linestrings in array, returns results per linestring, + multiple cellids per linestring.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + ls = linestrings([[(2.0, 1.0), (15.0, 1.0)], [(2.0, 19.0), (15.0, 19.0)]]) + result = ix.intersects(ls) + assert len(result) == 4 + assert result.cellids[0] == 4 + assert result.cellids[1] == 5 + assert result.cellids[2] == 2 + assert result.cellids[3] == 7 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_linestring_array_inside_outside(): + """Multiple linestrings in array inside/outside, returns multiple cellids per + linestring.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + ls = linestrings([[(2.0, 1.0), (15.0, 1.0)], [(25.0, 15.0), (35.0, 15.0)]]) + result = ix.intersects(ls) + assert len(result) == 2 + assert (result.shp_ids == 0).all() + assert result.cellids[0] == 4 + assert result.cellids[1] == 5 + + +# %% vertex grid polygons + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_polygon_array_in_one_cell(rtree): + """Single polygon in array inside, returns single intersection.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons([[(2.0, 1.0), (9.0, 1.0), (9.0, 7.0), (2.0, 1.0)]]) + result = ix.intersects(p) + assert len(result) == 1 + assert result.cellids[0] == 4 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_polygon_array_in_two_cells(rtree): + """Single polygon in array in 2 cells, returns multiple intersections.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons([[(5.0, 1.0), (15.0, 1.0), (15.0, 2.0), (5.0, 2.0)]]) + result = ix.intersects(p) + assert len(result) == 2 + assert result.cellids[0] == 4 + assert result.cellids[1] == 5 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_single_polygon_array_outside(rtree): + """Single polygon in array outside, returns empty result.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons([[(25, 5.0), (75, 5.0), (75, 7.5), (25, 7.5)]]) + result = ix.intersects(p) + assert len(result) == 0 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_polygon_array_single_result_per_polygon(): + """Multiple polygons in array, returns results per polygon, + single result per polygon.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + p = polygons( + [ + [(2.0, 1.0), (9.0, 1.0), (9.0, 7.0), (2.0, 1.0)], + [(2.0, 19.0), (9.0, 19.0), (9.0, 17.0), (2.0, 19.0)], + ] + ) + result = ix.intersects(p) + assert len(result) == 2 + assert result.cellids[0] == 4 + assert result.cellids[1] == 2 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_polygon_array_multiple_results_per_polygon(): + """Multiple polygons in array, returns results per polygon, + multiple results per polygon.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + p = polygons( + [ + [(5.0, 1.0), (15.0, 1.0), (15.0, 2.0), (5.0, 2.0)], + [(5.0, 19.0), (15.0, 19.0), (15.0, 18.0), (5.0, 18.0)], + ] + ) + result = ix.intersects(p) + assert len(result) == 4 + assert result.cellids[0] == 4 + assert result.cellids[1] == 5 + assert result.cellids[2] == 2 + assert result.cellids[3] == 7 + + +@requires_pkg("shapely") +def test_tri_grid_multiple_polygon_array_inside_outside(): + """Multiple polygons in array inside/outside, returns results per polygon, + multiple cellids per polygon.""" + gr = get_tri_grid() + ix = GridIntersect(gr) + p = polygons( + [ + [(2.0, 1.0), (9.0, 1.0), (9.0, 7.0), (2.0, 1.0)], + [(25, 15.0), (75, 15.0), (75, 17.5), (25, 17.5)], + ] + ) + result = ix.intersects(p) + assert len(result) == 1 + assert result.cellids[0] == 4 + + +# %% vertex grid intersection method + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_intersect_single_point_array(rtree): + """Single point in array ok.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([10], [20]) + result = ix.intersect(pts, return_all_intersections=True) + assert len(result.cellids) == 2 + assert result.cellids[0] == 2 + assert result.cellids[1] == 7 + result = ix.intersect(pts, return_all_intersections=False) + assert len(result.cellids) == 1 + assert result.cellids[0] == 2 + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_intersect_multiple_points_array(rtree): + """Multiple points in array raises error.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + pts = points([10, 5], [20, 5]) + with pytest.raises( + ValueError, match="intersect\(\) only accepts arrays containing one" + ): + ix.intersect(pts) + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_intersect_single_linestring_array(rtree): + """Single linestring in array ok.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings([[(5.0, 5.0), (15.0, 5.0)]]) + result = ix.intersect(ls) + assert len(result) == 2 + assert result.cellids[0] == 4 + assert result.cellids[1] == 5 + assert (result.lengths == 5.0).all() + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_intersect_multiple_linestring_array(rtree): + """Multiple linestrings in array raises error.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + ls = linestrings( + [ + [(5.0, 5.0), (15.0, 5.0)], + [(5.0, 15.0), (15.0, 15.0)], + ] + ) + with pytest.raises( + ValueError, match="intersect\(\) only accepts arrays containing one" + ): + ix.intersect(ls) + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_intersect_single_polygon_array(rtree): + """Single polygon in array ok.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons( + [ + [(2.5, 5.0), (17.5, 5.0), (17.5, 7.5), (2.5, 7.5)], + ] + ) + result = ix.intersect(p) + assert len(result) == 4 + assert result.cellids[0] == 1 + assert result.cellids[1] == 4 + assert result.cellids[2] == 5 + assert result.cellids[3] == 6 + assert (result.areas == 9.375).all() + + +@requires_pkg("shapely") +@rtree_toggle +def test_tri_grid_intersect_multiple_polygon_array(rtree): + """Multi-array input raises error.""" + gr = get_tri_grid() + ix = GridIntersect(gr, rtree=rtree) + p = polygons( + [ + [(2.5, 5.0), (17.5, 5.0), (17.5, 7.5), (2.5, 7.5)], + [(2.5, 15.0), (17.5, 15.0), (17.5, 17.5), (2.5, 17.5)], + ] + ) + with pytest.raises( + ValueError, match="intersect\(\) only accepts arrays containing one" + ): + ix.intersect(p) + + +def test_rtree_false_raises_in_points_to_cellids(): + """rtree=False raises error in points_to_cellids.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=False) + pts = points([1.0], [1.0]) + with pytest.raises( + ValueError, + match="points_to_cellids\(\) requires rtree=True when", + ): + ix.points_to_cellids(pts) + + +def test_rtree_false_raises_with_arrays_in_intersects(): + """rtree=False raises error in points_to_cellids.""" + gr = get_rect_grid() + ix = GridIntersect(gr, rtree=False) + pts = points([1.0, 10.0], [1.0, 10.0]) + with pytest.raises( + ValueError, + match="points_to_cellids\(\) requires rtree=True when initializing", + ): + ix.points_to_cellids(pts) + + +# %% +gr = get_rect_grid(top=np.ones(4).reshape((2, 2)), botm=np.zeros(4).reshape((1, 2, 2))) +ix = GridIntersect(gr) +pts = points([1.0, 25.0], [2.0, 25.0], [10.0, 0.5]) +ix.points_to_cellids(pts, handle_z="return", dataframe=True) +# %% +botm = np.concatenate( + [ + np.ones(4), + 0.5 * np.ones(4), + np.zeros(4), + ] +).reshape((3, 2, 2)) +gr = get_rect_grid(top=2 * np.ones(4).reshape((2, 2)), botm=botm) +ix = GridIntersect(gr) +result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z="ignore") +assert result.cellids[0] == (1, 0) +result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z="drop") +assert result.cellids[0] == (1, 0) +result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z="return") +assert result.cellids[0] == (1, 0) +assert result.layer[0] == 2.0 # returned as float to allow +/-inf +# %% From 6cddfc816679121b3f5664b39576cff778bd74ae Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 11:58:47 +0100 Subject: [PATCH 09/16] remove option `handle_z="drop"` from points_to_cellids() --- flopy/utils/gridintersect.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 1a9f9dff1b..80cf067e05 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -851,8 +851,7 @@ def points_to_cellids( handle_z : str, optional Method for handling z dimension in intersection results for point intersections. Default is "ignore" which ignores z-dimension. Other - options are "drop" which only returns results for points within grid - top and bottom, or "return" which returns the computed layer position + options is "return" which returns the computed layer position for each point. Points above the grid are returned as +np.inf and below the grid as -np.inf. return_nodenumbers : bool, optional @@ -917,10 +916,7 @@ def points_to_cellids( if handle_z != "ignore": laypos = self.get_layer_from_z(pts, rec.cellids) - if handle_z == "drop": - mask_z = np.isfinite(laypos) - rec = rec[mask_z] - elif handle_z == "return": + if handle_z == "return": # copy data to new array to include layer position rec = nprecfns.append_fields( rec, From 407378f64aa788e8e1da6f170568bb042a94de39 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 12:58:43 +0100 Subject: [PATCH 10/16] feat(GridIntersect): improve plotting and update example nb - make staticmethods accept shapely geometries as input instead of rec-arrays - update notebook --- .docs/Notebooks/grid_intersection_example.py | 186 ++++++++++++++----- flopy/utils/gridintersect.py | 46 +++-- 2 files changed, 164 insertions(+), 68 deletions(-) diff --git a/.docs/Notebooks/grid_intersection_example.py b/.docs/Notebooks/grid_intersection_example.py index 74218c315a..68da39d794 100644 --- a/.docs/Notebooks/grid_intersection_example.py +++ b/.docs/Notebooks/grid_intersection_example.py @@ -4,17 +4,28 @@ # notebook_metadata_filter: all # text_representation: # extension: .py -# format_name: light -# format_version: '1.5' -# jupytext_version: 1.14.5 +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.18.1 # kernelspec: -# display_name: Python 3 (ipykernel) +# display_name: flopy # language: python # name: python3 +# language_info: +# codemirror_mode: +# name: ipython +# version: 3 +# file_extension: .py +# mimetype: text/x-python +# name: python +# nbconvert_exporter: python +# pygments_lexer: ipython3 +# version: 3.13.7 # metadata: # section: dis # --- +# %% [markdown] # # Intersecting model grids with shapes # # _Note: This feature requires the shapely package (which is an optional FloPy dependency)._ @@ -35,9 +46,10 @@ # - [MultiLineString with triangular grid](#trigrid.2) # - [MultiPoint with triangular grid](#trigrid.3) +# %% [markdown] # Import packages -# + +# %% import sys import matplotlib as mpl @@ -56,8 +68,8 @@ print(f"matplotlib version: {mpl.__version__}") print(f"shapely version: {shapely.__version__}") print(f"flopy version: {flopy.__version__}") -# - +# %% [markdown] # ## [GridIntersect Class](#top) # # The GridIntersect class is constructed by passing a flopy modelgrid object to @@ -67,7 +79,7 @@ # - `rtree`: either `True` (default) or `False`. When True, an STR-tree is built, # which allows for fast spatial queries. Building the STR-tree takes some # time however. Setting the option to False avoids building the STR-tree but requires -# the intersection calculation to loop through all grid cells. It is generally +# the intersection calculation to loop through more candidate grid cells. It is generally # recommended to set this option to True. # - `local`: either `False` (default) or `True`. When True the local model coordinates # are used. When False, real-world coordinates are used. Can be useful if shapes are @@ -76,21 +88,27 @@ # The important methods in the GridIntersect object are: # # - `intersects()`: returns cellids for gridcells that intersect a shape (accepts -# shapely geometry objects, flopy geometry object, shapefile.Shape objects, and -# geojson objects) +# shapely geometry objects, arrays of shapely geometry object, flopy geometry object, +# shapefile.Shape objects, and geojson objects). # - `intersect()`: for intersecting the modelgrid with point, linestrings, and # polygon geometries (accepts shapely geometry objects, flopy geometry object, # shapefile.Shape objects, and geojson objects) +# - `points_to_cellids()`: for quickly locating points in the modelgrid and +# getting a 1:1 mapping of points to cellids. Especially useful when passing in array +# of points. # - `ix.plot_intersection_result()`: for plotting intersection results # # In the following sections examples of intersections are shown for structured # and vertex grids for different types of shapes (Polygon, LineString and Point). +# %% [markdown] # ## [Rectangular regular grid](#top) +# %% delc = 10 * np.ones(10, dtype=float) delr = 10 * np.ones(10, dtype=float) +# %% xoff = 0.0 yoff = 0.0 angrot = 0.0 @@ -98,51 +116,71 @@ delc, delr, top=None, botm=None, xoff=xoff, yoff=yoff, angrot=angrot ) -sgr.plot() +# %% +sgr.plot(); +# %% [markdown] # ### [Polygon with regular grid](#top) # Polygon to intersect with: +# %% p = Polygon( shell=[(15, 15), (20, 50), (35, 80.0), (80, 50), (80, 40), (40, 5), (15, 12)], holes=[[(25, 25), (25, 45), (45, 45), (45, 25)]], ) +# %% # Create the GridIntersect class for our modelgrid. -# TODO: remove method kwarg in v3.9.0 -ix = GridIntersect(sgr, method="vertex") +ix = GridIntersect(sgr) +# %% [markdown] # Do the intersect operation for a polygon +# %% result = ix.intersect(p) +# %% [markdown] # The results are returned as a numpy.recarray containing several fields based on the intersection performed. An explanation of the data in each of the possible fields is given below: # - **cellids**: contains the cell ids of the intersected grid cells -# - **vertices**: contains the vertices of the intersected shape # - **areas**: contains the area of the polygon in that grid cell (only for polygons) # - **lengths**: contains the length of the linestring in that grid cell (only for linestrings) # - **ixshapes**: contains the shapely object representing the intersected shape (useful for plotting the result) # -# Looking at the first few entries of the results of the polygon intersection (convert to pandas.DataFrame for prettier formatting) +# Looking at the first few entries of the results of the polygon intersection. Note that you can convert the result to a GeoDataFrame (if geopandas is installed) with `geo_dataframe=True`. -result[:5] -# pd.DataFrame(result).head() +# %% +ix.intersect(p, geo_dataframe=True).head() +# %% [markdown] # The cellids can be easily obtained +# %% result.cellids +# %% [markdown] # Or the areas +# %% result.areas -# If a user is only interested in which cells the shape intersects (and not the areas or the actual shape of the intersected object) with there is also the `intersects()` method. This method works for all types of shapely geometries. +# %% [markdown] +# If a user is only interested in which cells the shape intersects (and not the areas or +# the actual shape of the intersected object) with there is also the `intersects()` +# method. This method works for all types of shapely geometries including arrays of +# shapely geometries. +# +# This method returns `shp_ids` and `cellids` fields. The `shp_ids` field contains the +# index of the geometry in the original input shape provided by the user. This is useful +# when the input is an array of shapely geometries. In this case we have only one polygon, +# so the `shp_id` is always equal to 0. +# %% ix.intersects(p) +# %% [markdown] # The results of an intersection can be visualized with the `GridIntersect.plot_intersection_result()` method. -# + +# %% # create a figure and plot the grid fig, ax = plt.subplots(1, 1, figsize=(8, 8)) @@ -159,9 +197,9 @@ ) # add legend -ax.legend([h2], [i.get_label() for i in [h2]], loc="best") -# - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); +# %% [markdown] # The `intersect()` method contains several keyword arguments that specifically deal with polygons: # # - `contains_centroid`: only store intersection result if cell centroid is contained within polygon @@ -171,7 +209,7 @@ # # Example with `contains_centroid` set to True, only cells in which centroid is within the intersected polygon are stored. Note the difference with the previous result. -# + +# %% # contains_centroid example result2 = ix.intersect(p, contains_centroid=True) @@ -189,12 +227,12 @@ ) # add legend -ax.legend([h2], [i.get_label() for i in [h2]], loc="best") -# - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); +# %% [markdown] # Example with `min_area_threshold` set to 0.35, the intersection result in a cell should cover 35% or more of the cell area. -# + +# %% # min_area_threshold example result3 = ix.intersect(p, min_area_fraction=0.35) @@ -211,22 +249,25 @@ label="centroids of intersected gridcells", ) # add legend -ax.legend([h2], [i.get_label() for i in [h2]], loc="best") -# - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); +# %% [markdown] # ### [Polyline with regular grid](#top) # MultiLineString to intersect with: +# %% ls1 = LineString([(95, 105), (30, 50)]) ls2 = LineString([(30, 50), (90, 22)]) ls3 = LineString([(90, 22), (0, 0)]) mls = MultiLineString(lines=[ls1, ls2, ls3]) +# %% result = ix.intersect(mls) +# %% [markdown] # Plot the result -# + +# %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) sgr.plot(ax=ax) ix.plot_intersection_result(result, ax=ax, cmap="viridis") @@ -239,27 +280,30 @@ label="centroids of intersected gridcells", ) -ax.legend([h2], [i.get_label() for i in [h2]], loc="best") -# - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); +# %% [markdown] # ### [MultiPoint with regular grid](#top) # # MultiPoint to intersect with +# %% mp = MultiPoint( points=[Point(50.0, 0.0), Point(45.0, 45.0), Point(10.0, 10.0), Point(150.0, 100.0)] ) +# %% [markdown] # For points and linestrings there is a keyword argument `return_all_intersections` which will return multiple intersection results for points or (parts of) linestrings on cell boundaries. As an example, the difference is shown with the MultiPoint intersection. Note the number of red "+" symbols indicating the centroids of intersected cells, in the bottom left case, there are 4 results because the point lies exactly on the intersection between 4 grid cells. +# %% result = ix.intersect(mp) result_all = ix.intersect(mp, return_all_intersections=True) -# + +# %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) sgr.plot(ax=ax) -ix.plot_point(result, ax=ax, s=50, color="C0") -ix.plot_point(result_all, ax=ax, s=50, marker=".", color="C3") +ix.plot_point(result.ixshapes, ax=ax, ms=10, color="C0") +ix.plot_point(result_all.ixshapes, ax=ax, ms=10, marker=".", color="C3") for irow, icol in result.cellids: (h2,) = ax.plot( @@ -279,11 +323,57 @@ label="centroids with `return_all_intersections=True`", ) -ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best") -# - +ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best"); + +# %% [markdown] +# In addition to the `intersect()` and `intersects()` methods there is function to +# quickly get the cellids for many points. This is done with `points_to_cellids()`. +# +# First lets create an array containing shapely Points. + +# %% +n_pts = 10 +rng = np.random.default_rng(seed=42) +x_coords = rng.uniform(0, 100, n_pts) +y_coords = rng.uniform(0, 100, n_pts) +random_points = shapely.points(x_coords, y_coords) + +# %% [markdown] +# Now find the cellid for each point. `points_to_cellids` will only return a single result +# for each point. In case a point is on the boundary between multiple cells, it will +# return the cell with the lowest cellid. +# +# +# In this example we're returning the result as a +# node number to easily select the grid cells for our plot. But by default this method +# returns (row, col) for structured grids. + +# %% +# return cellid as node numbers: so (node,) instead of (row, col) +# this makes it easier to select the grid cells to highlight in the plot +result = ix.points_to_cellids(random_points, return_nodenumbers=True) + +# %% [markdown] +# Plot the result + +# %% +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) +sgr.plot(ax=ax) +ax.plot(x_coords, y_coords, "ro", ms=5, label="random points") +ix.plot_polygon( + ix.geoms[result.cellids.astype(int)], ax=ax, fc="yellow", edgecolor="k", zorder=5 +); +# %% [markdown] +# Note that `points_to_cellids()` returns NaNs for points that lie outside the model grid. + +# %% +ix.points_to_cellids(shapely.points([50, 110], [55, 50]), dataframe=True) + +# %% [markdown] # ## [Vertex Grid](#top) +# %% cell2d = [ [0, 83.33333333333333, 66.66666666666667, 3, 4, 2, 7], [1, 16.666666666666668, 33.333333333333336, 3, 4, 0, 5], @@ -307,17 +397,21 @@ ] tgr = fgrid.VertexGrid(vertices, cell2d) +# %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) pmv = fplot.PlotMapView(modelgrid=tgr) -pmv.plot_grid(ax=ax) +pmv.plot_grid(ax=ax); +# %% [markdown] # ### [Polygon with triangular grid](#top) +# %% ix2 = GridIntersect(tgr) +# %% result = ix2.intersect(p) -# + +# %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) ix.plot_intersection_result(result, ax=ax) @@ -330,14 +424,15 @@ label="centroids of intersected gridcells", ) -ax.legend([h2], [i.get_label() for i in [h2]], loc="best") -# - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); +# %% [markdown] # ### [LineString with triangular grid](#top) +# %% result = ix2.intersect(mls) -# + +# %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) ix2.plot_intersection_result(result, ax=ax, lw=3) @@ -349,17 +444,18 @@ label="centroids of intersected gridcells", ) -ax.legend([h2], [i.get_label() for i in [h2]], loc="best") -# - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); +# %% [markdown] # ### [MultiPoint with triangular grid](#top) +# %% result = ix2.intersect(mp) result_all = ix2.intersect(mp, return_all_intersections=True) -# + +# %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) -ix2.plot_intersection_result(result, ax=ax, color="k", zorder=5, s=80) +ix2.plot_intersection_result(result, ax=ax, color="k", zorder=5, ms=10) for cellid in result.cellids: (h2,) = ax.plot( @@ -378,4 +474,6 @@ label="centroids with return_all_intersections=True", ) -ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best") +ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best"); + +# %% diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 80cf067e05..459a33bb0f 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -932,7 +932,7 @@ def points_to_cellids( return rec @staticmethod - def plot_polygon(result, ax=None, **kwargs): + def plot_polygon(polys, ax=None, **kwargs): """method to plot the polygon intersection results from the resulting numpy.recarray. @@ -972,11 +972,10 @@ def add_poly_patch(poly): ppi = _polygon_patch(poly, facecolor=fc, **kwargs) patches.append(ppi) - # allow for result to be geodataframe - geoms = ( - result.ixshapes if isinstance(result, np.rec.recarray) else result.geometry - ) - for i, ishp in enumerate(geoms): + if isinstance(polys, (shapely.Polygon, shapely.MultiPolygon)): + polys = [polys] + + for i, ishp in enumerate(polys): if hasattr(ishp, "geoms"): for geom in ishp.geoms: add_poly_patch(geom) @@ -992,7 +991,7 @@ def add_poly_patch(poly): return ax @staticmethod - def plot_linestring(result, ax=None, cmap=None, **kwargs): + def plot_linestring(ls, ax=None, cmap=None, **kwargs): """method to plot the linestring intersection results from the resulting numpy.recarray. @@ -1029,15 +1028,14 @@ def plot_linestring(result, ax=None, cmap=None, **kwargs): else: specified_color = False + if isinstance(ls, (shapely.LineString, shapely.MultiLineString)): + ls = [ls] + if cmap is not None: colormap = plt.get_cmap(cmap) - colors = colormap(np.linspace(0, 1, result.shape[0])) + colors = colormap(np.linspace(0, 1, len(ls))) - # allow for result to be geodataframe - geoms = ( - result.ixshapes if isinstance(result, np.rec.recarray) else result.geometry - ) - for i, ishp in enumerate(geoms): + for i, ishp in enumerate(ls): if not specified_color: if cmap is None: c = f"C{i % 10}" @@ -1048,7 +1046,7 @@ def plot_linestring(result, ax=None, cmap=None, **kwargs): return ax @staticmethod - def plot_point(result, ax=None, **kwargs): + def plot_point(pts, ax=None, **kwargs): """method to plot the point intersection results from the resulting numpy.recarray. @@ -1056,8 +1054,8 @@ def plot_point(result, ax=None, **kwargs): Parameters ---------- - result : numpy.recarray or geopandas.GeoDataFrame - record array or GeoDataFrame containing intersection results + pts : array, geopandas.GeoSeries + array or GeoSeries containing geometries ax : matplotlib.pyplot.axes, optional axes to plot onto, if not provided, creates a new figure **kwargs: @@ -1074,14 +1072,14 @@ def plot_point(result, ax=None, **kwargs): if ax is None: _, ax = plt.subplots() # allow for result to be geodataframe - geoms = ( - result.ixshapes if isinstance(result, np.rec.recarray) else result.geometry - ) + if isinstance(pts, (shapely.Point, shapely.MultiPoint)): + pts = [pts] + maskpts = np.isin( - shapely.get_type_id(geoms), + shapely.get_type_id(pts), [shapely.GeometryType.POINT, shapely.GeometryType.MULTIPOINT], ) - shapely_plot.plot_points(geoms[maskpts], ax=ax, **kwargs) + shapely_plot.plot_points(pts[maskpts], ax=ax, **kwargs) return ax @@ -1116,7 +1114,7 @@ def plot_intersection_result(self, result, plot_grid=True, ax=None, **kwargs): shapely.get_type_id(geoms), [shapely.GeometryType.POINT, shapely.GeometryType.MULTIPOINT], ).all(): - ax = GridIntersect.plot_point(result, ax=ax, **kwargs) + ax = GridIntersect.plot_point(geoms, ax=ax, **kwargs) elif np.isin( shapely.get_type_id(geoms), [ @@ -1125,12 +1123,12 @@ def plot_intersection_result(self, result, plot_grid=True, ax=None, **kwargs): shapely.GeometryType.LINEARRING, ], ).all(): - ax = GridIntersect.plot_linestring(result, ax=ax, **kwargs) + ax = GridIntersect.plot_linestring(geoms, ax=ax, **kwargs) elif np.isin( shapely.get_type_id(geoms), [shapely.GeometryType.POLYGON, shapely.GeometryType.MULTIPOLYGON], ).all(): - ax = GridIntersect.plot_polygon(result, ax=ax, **kwargs) + ax = GridIntersect.plot_polygon(geoms, ax=ax, **kwargs) return ax From 3968ce3fed20b5057c5e9e62c4291cc722ef920e Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 13:03:38 +0100 Subject: [PATCH 11/16] ruff --- .docs/Notebooks/grid_intersection_example.py | 40 ++++++++------------ flopy/utils/gridintersect.py | 2 +- 2 files changed, 16 insertions(+), 26 deletions(-) diff --git a/.docs/Notebooks/grid_intersection_example.py b/.docs/Notebooks/grid_intersection_example.py index 68da39d794..a1a8fe99aa 100644 --- a/.docs/Notebooks/grid_intersection_example.py +++ b/.docs/Notebooks/grid_intersection_example.py @@ -88,12 +88,12 @@ # The important methods in the GridIntersect object are: # # - `intersects()`: returns cellids for gridcells that intersect a shape (accepts -# shapely geometry objects, arrays of shapely geometry object, flopy geometry object, +# shapely geometry objects, arrays of shapely geometry object, flopy geometry object, # shapefile.Shape objects, and geojson objects). # - `intersect()`: for intersecting the modelgrid with point, linestrings, and # polygon geometries (accepts shapely geometry objects, flopy geometry object, # shapefile.Shape objects, and geojson objects) -# - `points_to_cellids()`: for quickly locating points in the modelgrid and +# - `points_to_cellids()`: for quickly locating points in the modelgrid and # getting a 1:1 mapping of points to cellids. Especially useful when passing in array # of points. # - `ix.plot_intersection_result()`: for plotting intersection results @@ -117,8 +117,7 @@ ) # %% -sgr.plot(); - +sgr.plot() # %% [markdown] # ### [Polygon with regular grid](#top) # Polygon to intersect with: @@ -197,8 +196,7 @@ ) # add legend -ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best") # %% [markdown] # The `intersect()` method contains several keyword arguments that specifically deal with polygons: # @@ -227,8 +225,7 @@ ) # add legend -ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best") # %% [markdown] # Example with `min_area_threshold` set to 0.35, the intersection result in a cell should cover 35% or more of the cell area. @@ -249,8 +246,7 @@ label="centroids of intersected gridcells", ) # add legend -ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best") # %% [markdown] # ### [Polyline with regular grid](#top) # MultiLineString to intersect with: @@ -280,8 +276,7 @@ label="centroids of intersected gridcells", ) -ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best") # %% [markdown] # ### [MultiPoint with regular grid](#top) # @@ -323,10 +318,9 @@ label="centroids with `return_all_intersections=True`", ) -ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best"); - +ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best") # %% [markdown] -# In addition to the `intersect()` and `intersects()` methods there is function to +# In addition to the `intersect()` and `intersects()` methods there is function to # quickly get the cellids for many points. This is done with `points_to_cellids()`. # # First lets create an array containing shapely Points. @@ -344,7 +338,7 @@ # return the cell with the lowest cellid. # # -# In this example we're returning the result as a +# In this example we're returning the result as a # node number to easily select the grid cells for our plot. But by default this method # returns (row, col) for structured grids. @@ -362,8 +356,7 @@ ax.plot(x_coords, y_coords, "ro", ms=5, label="random points") ix.plot_polygon( ix.geoms[result.cellids.astype(int)], ax=ax, fc="yellow", edgecolor="k", zorder=5 -); - +) # %% [markdown] # Note that `points_to_cellids()` returns NaNs for points that lie outside the model grid. @@ -400,8 +393,7 @@ # %% fig, ax = plt.subplots(1, 1, figsize=(8, 8)) pmv = fplot.PlotMapView(modelgrid=tgr) -pmv.plot_grid(ax=ax); - +pmv.plot_grid(ax=ax) # %% [markdown] # ### [Polygon with triangular grid](#top) @@ -424,8 +416,7 @@ label="centroids of intersected gridcells", ) -ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best") # %% [markdown] # ### [LineString with triangular grid](#top) @@ -444,8 +435,7 @@ label="centroids of intersected gridcells", ) -ax.legend([h2], [i.get_label() for i in [h2]], loc="best"); - +ax.legend([h2], [i.get_label() for i in [h2]], loc="best") # %% [markdown] # ### [MultiPoint with triangular grid](#top) @@ -474,6 +464,6 @@ label="centroids with return_all_intersections=True", ) -ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best"); +ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best") # %% diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 459a33bb0f..2060ab297a 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -1074,7 +1074,7 @@ def plot_point(pts, ax=None, **kwargs): # allow for result to be geodataframe if isinstance(pts, (shapely.Point, shapely.MultiPoint)): pts = [pts] - + maskpts = np.isin( shapely.get_type_id(pts), [shapely.GeometryType.POINT, shapely.GeometryType.MULTIPOINT], From ae92d24a4b84fbab1235f9e91d5a927b0c4ee1c5 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 13:34:49 +0100 Subject: [PATCH 12/16] fix tests --- autotest/test_gridintersect.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py index f094bd76e5..4b6bea841c 100644 --- a/autotest/test_gridintersect.py +++ b/autotest/test_gridintersect.py @@ -1,3 +1,4 @@ +#%% import matplotlib.pyplot as plt import numpy as np import pytest @@ -1099,8 +1100,7 @@ def test_rect_grid_multiple_points_array_with_z_points_to_cellids(): result = ix.points_to_cellids(pts, handle_z="return") assert ~np.isfinite(result.layer[0]) assert np.isnan(result.cellids[1]) - result = ix.points_to_cellids(pts, handle_z="drop") - assert len(result) == 0 + # %% array-based input - structured grid linestrings @@ -1469,8 +1469,6 @@ def test_tri_grid_multiple_points_array_with_z_points_to_cellids(): result = ix.points_to_cellids(pts, handle_z="return") assert result.layer[0] == 0.0 assert np.isnan(result.cellids[1]) - result = ix.points_to_cellids(pts, handle_z="drop") - assert len(result) == 1 @requires_pkg("shapely") From 2a8070ed4fd96f6be20daec59957fea212234674 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 13:54:10 +0100 Subject: [PATCH 13/16] ruff --- autotest/test_gridintersect.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py index 4b6bea841c..dd788dc70f 100644 --- a/autotest/test_gridintersect.py +++ b/autotest/test_gridintersect.py @@ -1,4 +1,4 @@ -#%% +# %% import matplotlib.pyplot as plt import numpy as np import pytest @@ -1102,7 +1102,6 @@ def test_rect_grid_multiple_points_array_with_z_points_to_cellids(): assert np.isnan(result.cellids[1]) - # %% array-based input - structured grid linestrings From 3df8f4b5ae62273971f405691ccaf1aa0808522a Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 14:21:06 +0100 Subject: [PATCH 14/16] test(GridIntersect): remove all references to deprecated method kwarg --- .../groundwater2023_watershed_example.py | 17 ++++++++--------- .../mf6_parallel_model_splitting_example.py | 6 +++--- .docs/md/optional_dependencies.md | 3 +-- autotest/test_model_splitter.py | 2 +- 4 files changed, 13 insertions(+), 15 deletions(-) diff --git a/.docs/Notebooks/groundwater2023_watershed_example.py b/.docs/Notebooks/groundwater2023_watershed_example.py index b6e9e8e49a..cccf526bc0 100644 --- a/.docs/Notebooks/groundwater2023_watershed_example.py +++ b/.docs/Notebooks/groundwater2023_watershed_example.py @@ -94,7 +94,7 @@ def densify_geometry(line, step, keep_internal_nodes=True): # function to set the active and inactive model area def set_idomain(grid, boundary): - ix = GridIntersect(grid, method="vertex", rtree=True) + ix = GridIntersect(grid, rtree=True) result = ix.intersect(Polygon(boundary)) idx = list(result.cellids) idx = np.array(idx, dtype=int) @@ -241,7 +241,7 @@ def set_idomain(grid, boundary): extrapolate_edges=True, ) -ixs = flopy.utils.GridIntersect(struct_grid, method="structured") +ixs = flopy.utils.GridIntersect(struct_grid) cellids = [] for sg in sgs: v = ixs.intersect(LineString(sg), sort_by_cellid=True) @@ -315,7 +315,7 @@ def set_idomain(grid, boundary): extrapolate_edges=True, ) -ixs = flopy.utils.GridIntersect(struct_vrc_grid, method="structured") +ixs = flopy.utils.GridIntersect(struct_vrc_grid) cellids = [] for sg in sgs: v = ixs.intersect(LineString(sg), sort_by_cellid=True) @@ -419,7 +419,7 @@ def set_idomain(grid, boundary): top_nested_grid = [top_ngp, top_ngc] # + -ixs = flopy.utils.GridIntersect(struct_gridp, method="structured") +ixs = flopy.utils.GridIntersect(struct_gridp) cellids = [] for sg in sgs: v = ixs.intersect(LineString(sg), sort_by_cellid=True) @@ -429,7 +429,7 @@ def set_idomain(grid, boundary): intersection_ngp[loc] = 1 intersection_ngp[idomainp[0] == 0] = 0 -ixs = flopy.utils.GridIntersect(struct_gridc, method="structured") +ixs = flopy.utils.GridIntersect(struct_gridc) cellids = [] for sg in sgs: v = ixs.intersect(LineString(sg), sort_by_cellid=True) @@ -515,7 +515,7 @@ def set_idomain(grid, boundary): extrapolate_edges=True, ) -ixs = flopy.utils.GridIntersect(quadtree_grid, method="vertex") +ixs = flopy.utils.GridIntersect(quadtree_grid) cellids = [] for sg in sgs: v = ixs.intersect(LineString(sg), sort_by_cellid=True) @@ -583,7 +583,7 @@ def set_idomain(grid, boundary): extrapolate_edges=True, ) -ixs = flopy.utils.GridIntersect(triangular_grid) # , method="vertex") +ixs = flopy.utils.GridIntersect(triangular_grid) cellids = [] for sg in sgs: v = ixs.intersect( @@ -633,13 +633,12 @@ def set_idomain(grid, boundary): extrapolate_edges=True, ) -ixs = flopy.utils.GridIntersect(voronoi_grid, method="vertex") +ixs = flopy.utils.GridIntersect(voronoi_grid) cellids = [] for sg in sgs: v = ixs.intersect( LineString(sg), return_all_intersections=True, - keepzerolengths=False, sort_by_cellid=True, ) cellids += v["cellids"].tolist() diff --git a/.docs/Notebooks/mf6_parallel_model_splitting_example.py b/.docs/Notebooks/mf6_parallel_model_splitting_example.py index 12595203b1..5ac4944dfa 100644 --- a/.docs/Notebooks/mf6_parallel_model_splitting_example.py +++ b/.docs/Notebooks/mf6_parallel_model_splitting_example.py @@ -354,7 +354,7 @@ def string2geom(geostring, conversion=None): # + # calculate and set idomain -ix = flopy.utils.GridIntersect(modelgrid, method="vertex", rtree=True) +ix = flopy.utils.GridIntersect(modelgrid, rtree=True) result = ix.intersect(Polygon(boundary_polygon)) idxs = tuple(zip(*result.cellids)) idomain = np.zeros((nrow, ncol), dtype=int) @@ -369,7 +369,7 @@ def string2geom(geostring, conversion=None): # Intersect the stream segments with the modelgrid -ixs = flopy.utils.GridIntersect(modelgrid, method="structured") +ixs = flopy.utils.GridIntersect(modelgrid) cellids = [] for seg in segs: v = ixs.intersect(LineString(seg), sort_by_cellid=True) @@ -401,7 +401,7 @@ def string2geom(geostring, conversion=None): # + # intersect stream segs to simulate as drains -ixs = flopy.utils.GridIntersect(modelgrid, method="structured") +ixs = flopy.utils.GridIntersect(modelgrid) drn_cellids = [] drn_lengths = [] for seg in segs: diff --git a/.docs/md/optional_dependencies.md b/.docs/md/optional_dependencies.md index a5a1e17d30..2ce5319c23 100644 --- a/.docs/md/optional_dependencies.md +++ b/.docs/md/optional_dependencies.md @@ -13,9 +13,8 @@ Dependencies for optional features are listed below. These may be installed with | `.resample_to_grid()` in `flopy.utils.rasters` | **scipy.interpolate** | | `.interpolate()` in `flopy.mf6.utils.reference` `StructuredSpatialReference` class | **scipy.interpolate** | | `.get_authority_crs()` in `flopy.utils.crs` | **pyproj** >= 2.2.0 | -| `.generate_classes()` in `flopy.mf6.utils` | [**modflow-devtools**](https://github.com/MODFLOW-ORG/modflow-devtools) | +| `.generate_classes()` in `flopy.mf6.utils` | [**modflow-devtools**](https://github.com/MODFLOW-ORG/modflow-devtools) | | `GridIntersect()` in `flopy.utils.gridintersect` | **shapely** | -| `GridIntersect().plot_polygon()` in `flopy.utils.gridintersect` | **shapely** and **descartes** | | `Raster()` in `flopy.utils.Raster` | **rasterio**, **rasterstats**, **affine**, and **scipy** | | `Raster().sample_polygon()` in `flopy.utils.Raster` | **shapely** | | `Raster().crop()` in `flopy.utils.Raster` | **shapely** | diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py index 895a87aca9..31560e1566 100644 --- a/autotest/test_model_splitter.py +++ b/autotest/test_model_splitter.py @@ -1008,7 +1008,7 @@ def string2geom(geostring, conversion=None): botm=np.full((1, nrow, ncol), -100.0), ) - ixs = flopy.utils.GridIntersect(modelgrid, method="vertex", rtree=True) + ixs = flopy.utils.GridIntersect(modelgrid, rtree=True) result = ixs.intersect( [ boundary, From a2bc6990cf853739098f0ce821d70dd5be4defbd Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 16:17:10 +0100 Subject: [PATCH 15/16] forgot one --- .docs/Notebooks/groundwater2023_watershed_example.py | 1 - 1 file changed, 1 deletion(-) diff --git a/.docs/Notebooks/groundwater2023_watershed_example.py b/.docs/Notebooks/groundwater2023_watershed_example.py index cccf526bc0..619df91f92 100644 --- a/.docs/Notebooks/groundwater2023_watershed_example.py +++ b/.docs/Notebooks/groundwater2023_watershed_example.py @@ -589,7 +589,6 @@ def set_idomain(grid, boundary): v = ixs.intersect( LineString(sg), return_all_intersections=True, - keepzerolengths=False, sort_by_cellid=True, ) cellids += v["cellids"].tolist() From 6953bfb81c7de3aa4d437ff75c18f715a0e3513b Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Thu, 27 Nov 2025 17:26:43 +0100 Subject: [PATCH 16/16] refactor(GridIntersect): more cleanup - replace 3 instances where new nodenumber to rowcol method should be used --- flopy/utils/gridintersect.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 2060ab297a..69109ae63f 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -508,9 +508,7 @@ def _intersect_point( # if structured calculate (i, j) cell address if self.mfgrid.grid_type == "structured": - rec.cellids = list( - zip(*self.mfgrid.get_lrc([self.cellids[keep_cid]])[0][1:]) - ) + rec.cellids = self._nodenumbers_to_rowcol(self.cellids[keep_cid]) else: rec.cellids = self.cellids[keep_cid] rec.ixshapes = keep_pts @@ -631,9 +629,7 @@ def parse_linestrings_in_geom_collection(gc): rec = np.recarray(len(ixresult), names=names, formats=formats) # if structured grid calculate (i, j) cell address if self.mfgrid.grid_type == "structured": - rec.cellids = list( - zip(*self.mfgrid.get_lrc([self.cellids[qcellids]])[0][1:]) - ) + rec.cellids = self._nodenumbers_to_rowcol(self.cellids[qcellids]) else: rec.cellids = self.cellids[qcellids] rec.ixshapes = ixresult @@ -718,9 +714,7 @@ def parse_polygons_in_geom_collection(gc): rec = np.recarray(len(ixresult), names=names, formats=formats) # if structured calculate (i, j) cell address if self.mfgrid.grid_type == "structured": - rec.cellids = list( - zip(*self.mfgrid.get_lrc([self.cellids[qcellids]])[0][1:]) - ) + rec.cellids = self._nodenumbers_to_rowcol(self.cellids[qcellids]) else: rec.cellids = self.cellids[qcellids] rec.ixshapes = ixresult