diff --git a/.docs/Notebooks/grid_intersection_example.py b/.docs/Notebooks/grid_intersection_example.py
index 74218c315a..1b947faa2d 100644
--- a/.docs/Notebooks/grid_intersection_example.py
+++ b/.docs/Notebooks/grid_intersection_example.py
@@ -1,20 +1,32 @@
+# %%
# ---
# jupyter:
# jupytext:
# 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 +47,10 @@
# - [MultiLineString with triangular grid](#trigrid.2)
# - [MultiPoint with triangular grid](#trigrid.3)
+# %% [markdown]
# Import packages
-# +
+# %%
import sys
import matplotlib as mpl
@@ -56,8 +69,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 +80,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,81 +89,117 @@
# 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 objects, 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
sgr = fgrid.StructuredGrid(
- delc, delr, top=None, botm=None, xoff=xoff, yoff=yoff, angrot=angrot
+ delc,
+ delr,
+ top=np.ones(100).reshape((10, 10)),
+ botm=np.zeros(100).reshape((1, 10, 10)),
+ xoff=xoff,
+ yoff=yoff,
+ angrot=angrot,
)
+# %%
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")
+gi = GridIntersect(sgr)
+# %% [markdown]
# Do the intersect operation for a polygon
-result = ix.intersect(p)
-
-# 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)
+# %%
+result = gi.intersect(p, geo_dataframe=False)
+
+# %% [markdown]
+# The results are returned as a geopandas.GeoDataFrame or numpy.recarray containing several fields based on the intersection performed. An explanation of the data in each of the possible columns is given below:
+# - `cellid`: contains the cell ids of the intersected grid cells
+# - `row`: contains the row index of the intersected grid cells (only for structured grids)
+# - `col`: contains the column index of the intersected grid cells (only for structured grids)
+# - `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`/`geometry`: contains the shapely object representing the intersected shape (useful for plotting the result)
+#
+# __Note__: The `cellids` column is deprecated since flopy 3.11 but still included in the result for backward compatibility. It contains (row, column) tuples for structured grids and cellids for vertex grids.
#
-# 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()
+# %%
+gi.intersect(p, geo_dataframe=True).head()
-# The cellids can be easily obtained
+# %% [markdown]
+# The rows and columns can be easily obtained.
-result.cellids
+# %%
+result.row, result.col
+# %% [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_id` and `cellid` columns. The `shp_id` column 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. For structured grids the row and column indices
+# are returned in the `row` and `col` columns.
-ix.intersects(p)
+# %%
+gi.intersects(p, dataframe=True)
+# %% [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))
# the intersection object contains some helpful plotting commands
-ix.plot_intersection_result(result, ax=ax)
+gi.plot_intersection_result(result, ax=ax)
# add black x at cell centers
-for irow, icol in result.cellids:
+for irow, icol in zip(result.row, result.col):
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
@@ -160,8 +209,7 @@
# add legend
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,16 +219,16 @@
#
# 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)
+result2 = gi.intersect(p, contains_centroid=True, geo_dataframe=True)
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
-ix.plot_intersection_result(result2, ax=ax)
+gi.plot_intersection_result(result2, ax=ax)
# add black x at cell centers
-for irow, icol in result2.cellids:
+for irow, icol in zip(result2.row, result2.col):
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
@@ -190,20 +238,19 @@
# add legend
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)
+result3 = gi.intersect(p, min_area_fraction=0.35)
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
-ix.plot_intersection_result(result3, ax=ax)
+gi.plot_intersection_result(result3, ax=ax)
# add black x at cell centers
-for irow, icol in result3.cellids:
+for irow, icol in zip(result3.row, result3.col):
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
@@ -212,26 +259,28 @@
)
# add legend
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)
+# %%
+result = gi.intersect(mls, geo_dataframe=True)
+# %% [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")
+gi.plot_intersection_result(result, ax=ax, cmap="tab20")
-for irow, icol in result.cellids:
+for irow, icol in zip(result.row, result.col):
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
@@ -240,28 +289,30 @@
)
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)
+# %%
+result = gi.intersect(mp, geo_dataframe=True)
+result_all = gi.intersect(mp, return_all_intersections=True, geo_dataframe=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")
+gi.plot_point(result.geometry, ax=ax, ms=10, color="C0")
+gi.plot_point(result_all.geometry, ax=ax, ms=10, marker=".", color="C3")
-for irow, icol in result.cellids:
+for irow, icol in zip(result.row, result.col):
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
@@ -270,7 +321,7 @@
label="centroids of intersected cells",
)
-for irow, icol in result_all.cellids:
+for irow, icol in zip(result_all.row, result_all.col):
(h3,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
@@ -280,10 +331,47 @@
)
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.
+
+# %%
+result = gi.points_to_cellids(random_points)
+result
+
+# %% [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")
+gi.plot_polygon(gi.geoms[result.cellid], ax=ax, fc="yellow", edgecolor="k", zorder=5)
+# %% [markdown]
+# Note that `points_to_cellids()` returns NA for points that lie outside the model grid.
+
+# %%
+gi.points_to_cellids(shapely.points([50, 110], [55, 50]))
+
+# %% [markdown]
# ## [Vertex Grid](#top)
+# %%
cell2d = [
[0, 83.33333333333333, 66.66666666666667, 3, 4, 2, 7],
[1, 16.666666666666668, 33.333333333333336, 3, 4, 0, 5],
@@ -307,22 +395,25 @@
]
tgr = fgrid.VertexGrid(vertices, cell2d)
+# %%
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
pmv = fplot.PlotMapView(modelgrid=tgr)
pmv.plot_grid(ax=ax)
-
+# %% [markdown]
# ### [Polygon with triangular grid](#top)
-ix2 = GridIntersect(tgr)
+# %%
+gi2 = GridIntersect(tgr)
-result = ix2.intersect(p)
+# %%
+result = gi2.intersect(p, geo_dataframe=True)
-# +
+# %%
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
-ix.plot_intersection_result(result, ax=ax)
+gi2.plot_intersection_result(result, ax=ax)
# only cells that intersect with shape
-for cellid in result.cellids:
+for cellid in result.cellid:
(h2,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
@@ -331,17 +422,17 @@
)
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
-# -
-
+# %% [markdown]
# ### [LineString with triangular grid](#top)
-result = ix2.intersect(mls)
+# %%
+result = gi2.intersect(mls, geo_dataframe=True)
-# +
+# %%
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
-ix2.plot_intersection_result(result, ax=ax, lw=3)
+gi2.plot_intersection_result(result, ax=ax, lw=3)
-for cellid in result.cellids:
+for cellid in result.cellid:
(h2,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
@@ -350,18 +441,18 @@
)
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)
+# %%
+result = gi2.intersect(mp, geo_dataframe=True)
+result_all = gi2.intersect(mp, return_all_intersections=True, geo_dataframe=True)
-# +
+# %%
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
-ix2.plot_intersection_result(result, ax=ax, color="k", zorder=5, s=80)
+gi2.plot_intersection_result(result, ax=ax, color="k", zorder=5, ms=10)
-for cellid in result.cellids:
+for cellid in result.cellid:
(h2,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
@@ -369,7 +460,7 @@
ms=15,
label="centroids of intersected cells",
)
-for cellid in result_all.cellids:
+for cellid in result_all.cellid:
(h3,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
@@ -379,3 +470,5 @@
)
ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best")
+
+# %%
diff --git a/.docs/Notebooks/groundwater2023_watershed_example.py b/.docs/Notebooks/groundwater2023_watershed_example.py
index b6e9e8e49a..4d64baeca4 100644
--- a/.docs/Notebooks/groundwater2023_watershed_example.py
+++ b/.docs/Notebooks/groundwater2023_watershed_example.py
@@ -94,8 +94,8 @@ 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)
- result = ix.intersect(Polygon(boundary))
+ ix = GridIntersect(grid, rtree=True)
+ result = ix.intersect(Polygon(boundary), geo_dataframe=False)
idx = list(result.cellids)
idx = np.array(idx, dtype=int)
nr = idx.shape[0]
@@ -241,10 +241,10 @@ 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)
+ v = ixs.intersect(LineString(sg), sort_by_cellid=True, geo_dataframe=False)
cellids += v["cellids"].tolist()
intersection_sg = np.zeros(struct_grid.shape[1:])
for loc in cellids:
@@ -315,10 +315,10 @@ 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)
+ v = ixs.intersect(LineString(sg), sort_by_cellid=True, geo_dataframe=False)
cellids += v["cellids"].tolist()
intersection_sg_vrc = np.zeros(struct_vrc_grid.shape[1:])
for loc in cellids:
@@ -419,20 +419,20 @@ 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)
+ v = ixs.intersect(LineString(sg), sort_by_cellid=True, geo_dataframe=False)
cellids += v["cellids"].tolist()
intersection_ngp = np.zeros(struct_gridp.shape[1:])
for loc in cellids:
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)
+ v = ixs.intersect(LineString(sg), sort_by_cellid=True, geo_dataframe=False)
cellids += v["cellids"].tolist()
intersection_ngc = np.zeros(struct_gridc.shape[1:])
for loc in cellids:
@@ -515,10 +515,10 @@ 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)
+ v = ixs.intersect(LineString(sg), sort_by_cellid=True, geo_dataframe=False)
cellids += v["cellids"].tolist()
intersection_qg = np.zeros(quadtree_grid.shape[1:])
for loc in cellids:
@@ -583,14 +583,14 @@ 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(
LineString(sg),
return_all_intersections=True,
- keepzerolengths=False,
sort_by_cellid=True,
+ geo_dataframe=False,
)
cellids += v["cellids"].tolist()
intersection_tg = np.zeros(triangular_grid.shape[1:])
@@ -633,14 +633,14 @@ 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,
+ geo_dataframe=False,
)
cellids += v["cellids"].tolist()
intersection_vg = np.zeros(voronoi_grid.shape[1:])
diff --git a/.docs/Notebooks/mf6_parallel_model_splitting_example.py b/.docs/Notebooks/mf6_parallel_model_splitting_example.py
index 12595203b1..16be74a4be 100644
--- a/.docs/Notebooks/mf6_parallel_model_splitting_example.py
+++ b/.docs/Notebooks/mf6_parallel_model_splitting_example.py
@@ -354,8 +354,8 @@ def string2geom(geostring, conversion=None):
# +
# calculate and set idomain
-ix = flopy.utils.GridIntersect(modelgrid, method="vertex", rtree=True)
-result = ix.intersect(Polygon(boundary_polygon))
+ix = flopy.utils.GridIntersect(modelgrid, rtree=True)
+result = ix.intersect(Polygon(boundary_polygon), geo_dataframe=False)
idxs = tuple(zip(*result.cellids))
idomain = np.zeros((nrow, ncol), dtype=int)
idomain[idxs] = 1
@@ -369,10 +369,10 @@ 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)
+ v = ixs.intersect(LineString(seg), sort_by_cellid=True, geo_dataframe=False)
cellids += v["cellids"].tolist()
intersection_rg = np.zeros(modelgrid.shape[1:])
for loc in cellids:
@@ -401,11 +401,11 @@ 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:
- v = ixs.intersect(LineString(seg), sort_by_cellid=True)
+ v = ixs.intersect(LineString(seg), sort_by_cellid=True, geo_dataframe=False)
drn_cellids += v["cellids"].tolist()
drn_lengths += v["lengths"].tolist()
# -
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_gridintersect.py b/autotest/test_gridintersect.py
index 15ec8b494b..2a1a876fd2 100644
--- a/autotest/test_gridintersect.py
+++ b/autotest/test_gridintersect.py
@@ -1,5 +1,7 @@
+# %%
import matplotlib.pyplot as plt
import numpy as np
+import pandas as pd
import pytest
from modflow_devtools.markers import requires_pkg
from modflow_devtools.misc import has_pkg
@@ -8,9 +10,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,
@@ -21,6 +22,7 @@
)
rtree_toggle = pytest.mark.parametrize("rtree", [True, False])
+df_toggle = False # set False to silence warnings, remove when dataframe is default
def get_tri_grid(angrot=0.0, xyoffset=0.0, triangle_exe=None):
@@ -117,92 +119,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=False, geo_dataframe=df_toggle
+ )
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=True, geo_dataframe=df_toggle
+ )
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))
+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=False, geo_dataframe=df_toggle)
+ assert result["cellids"][0] == (1, 0)
+ result = ix.intersect(Point(2.0, 2.0, 0.2), handle_z=True, geo_dataframe=df_toggle)
+ assert result["cellids"][0] == (1, 0)
+ assert result["layer"][0] == 2
@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)]))
+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=False, geo_dataframe=df_toggle
+ )
assert len(result) == 1
assert result.cellids[0] == (1, 0)
-
-
-@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
- assert result.cellids[0] == (1, 0)
- assert result.cellids[1] == (0, 1)
+ result = ix.intersect(Point(2.0, 2.0, 10.0), handle_z=True, geo_dataframe=df_toggle)
+ assert len(result) == 1
+ assert np.ma.is_masked(result["layer"][0])
# %% test point shapely
@@ -210,69 +175,73 @@ 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")
+ result = ix.intersect((25.0, 25.0), shapetype="point", geo_dataframe=df_toggle)
assert len(result) == 0
@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)
- result = ix.intersect(Point(20.0, 10.0))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(Point(20.0, 10.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert np.all(result.cellids[0] == (0, 1))
@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)
- result = ix.intersect(Point(10.0, 10.0))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(Point(10.0, 10.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert np.all(result.cellids[0] == (0, 0))
@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)
- result = ix.intersect(Point(4.0, 4.0))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(Point(4.0, 4.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert result.cellids[0] == 0
- result = ix.intersect(Point(4.0, 6.0))
+ result = ix.intersect(Point(4.0, 6.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert result.cellids[0] == 0
- result = ix.intersect(Point(6.0, 6.0))
+ result = ix.intersect(Point(6.0, 6.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert result.cellids[0] == 0
- result = ix.intersect(Point(6.0, 4.0))
+ result = ix.intersect(Point(6.0, 4.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert result.cellids[0] == 0
@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)
- result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(2.0, 2.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ MultiPoint([Point(1.0, 1.0), Point(2.0, 2.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 1
assert result.cellids[0] == (1, 0)
@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)
- result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(12.0, 12.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ MultiPoint([Point(1.0, 1.0), Point(12.0, 12.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.cellids[0] == (0, 1)
assert result.cellids[1] == (1, 0)
@@ -285,7 +254,7 @@ def test_tri_grid_point_outside(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(Point(25.0, 25.0))
+ result = ix.intersect(Point(25.0, 25.0), geo_dataframe=df_toggle)
assert len(result) == 0
@@ -296,7 +265,7 @@ def test_tri_grid_point_on_outer_boundary(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(Point(20.0, 10.0))
+ result = ix.intersect(Point(20.0, 10.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert np.all(result.cellids[0] == 0)
@@ -308,7 +277,7 @@ def test_tri_grid_point_on_inner_boundary(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(Point(10.0, 10.0))
+ result = ix.intersect(Point(10.0, 10.0), geo_dataframe=df_toggle)
assert len(result) == 1
assert np.all(result.cellids[0] == 0)
@@ -320,7 +289,9 @@ def test_tri_grid_multipoint_in_one_cell(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(2.0, 2.0)]))
+ result = ix.intersect(
+ MultiPoint([Point(1.0, 1.0), Point(2.0, 2.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 1
assert result.cellids[0] == 1
@@ -332,7 +303,9 @@ def test_tri_grid_multipoint_in_multiple_cells(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(MultiPoint([Point(1.0, 1.0), Point(12.0, 12.0)]))
+ result = ix.intersect(
+ MultiPoint([Point(1.0, 1.0), Point(12.0, 12.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.cellids[0] == 0
assert result.cellids[1] == 1
@@ -341,160 +314,27 @@ 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)
+ r = ix.intersect(
+ Point(*v), return_all_intersections=True, geo_dataframe=df_toggle
+ )
assert len(r) == n
@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)]),
- ]
+ r = ix.intersect(
+ Point(*v), return_all_intersections=True, geo_dataframe=df_toggle
)
- )
- 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)
+ assert len(r) == n
# %% test linestring shapely
@@ -502,19 +342,23 @@ def test_rect_grid_linestring_starting_on_vertex():
@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)
- result = ix.intersect(LineString([(25.0, 25.0), (21.0, 5.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ LineString([(25.0, 25.0), (21.0, 5.0)]), geo_dataframe=df_toggle
+ )
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)
- result = ix.intersect(LineString([(5.0, 5.0), (15.0, 5.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ LineString([(5.0, 5.0), (15.0, 5.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.lengths.sum() == 10.0
assert result.cellids[0] == (1, 0)
@@ -523,10 +367,12 @@ 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)
- result = ix.intersect(LineString([(15.0, 20.0), (5.0, 20.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ LineString([(15.0, 20.0), (5.0, 20.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.lengths.sum() == 10.0
assert result.cellids[0] == (0, 0)
@@ -535,10 +381,12 @@ 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)
- result = ix.intersect(LineString([(5.0, 10.0), (15.0, 10.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ LineString([(5.0, 10.0), (15.0, 10.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.lengths.sum() == 10.0
assert result.cellids[0] == (0, 0)
@@ -547,13 +395,14 @@ 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)])]
- )
+ ),
+ geo_dataframe=df_toggle,
)
assert len(result) == 1
assert result.lengths == 16.0
@@ -562,16 +411,17 @@ 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(
[
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)]),
]
- )
+ ),
+ geo_dataframe=df_toggle,
)
assert len(result) == 3
assert np.allclose(sum(result.lengths), 40.19197584109293)
@@ -579,10 +429,12 @@ 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)
- result = ix.intersect(LineString([(5.0, 9), (15.0, 5.0), (5.0, 1.0)]))
+ ix = GridIntersect(gr, rtree=rtree)
+ result = ix.intersect(
+ LineString([(5.0, 9), (15.0, 5.0), (5.0, 1.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.cellids[0] == (1, 0)
assert result.cellids[1] == (1, 1)
@@ -590,18 +442,23 @@ 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")
- result = ix.intersect(LineString([(5, 15), (5.0, 9), (15.0, 5.0), (5.0, 1.0)]))
+ ix = GridIntersect(gr)
+ result = ix.intersect(
+ LineString([(5, 15), (5.0, 9), (15.0, 5.0), (5.0, 1.0)]),
+ geo_dataframe=df_toggle,
+ )
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")
- result = ix.intersect(LineString([(10.0, 10.0), (15.0, 5.0)]))
+ ix = GridIntersect(gr)
+ result = ix.intersect(
+ LineString([(10.0, 10.0), (15.0, 5.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 1
assert np.allclose(result.lengths.sum(), np.sqrt(50))
assert result.cellids[0] == (1, 1)
@@ -609,34 +466,34 @@ 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):
ls = LineString([(x[i], y[i]), (x[i + 1], y[i + 1])])
- r = ix.intersect(ls, return_all_intersections=True)
+ r = ix.intersect(ls, return_all_intersections=True, geo_dataframe=df_toggle)
assert len(r) == n_intersections[i]
@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)
+ r = ix.intersect(ls, return_all_intersections=False, geo_dataframe=df_toggle)
assert len(r) == 1
@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)
+ r = ix.intersect(ls, return_all_intersections=True, geo_dataframe=df_toggle)
assert len(r) == 3
@@ -647,7 +504,9 @@ def test_tri_grid_linestring_outside(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(LineString([(25.0, 25.0), (21.0, 5.0)]))
+ result = ix.intersect(
+ LineString([(25.0, 25.0), (21.0, 5.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 0
@@ -658,7 +517,9 @@ def test_tri_grid_linestring_in_2cells(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(LineString([(5.0, 5.0), (5.0, 15.0)]))
+ result = ix.intersect(
+ LineString([(5.0, 5.0), (5.0, 15.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.lengths.sum() == 10.0
assert result.cellids[0] == 1
@@ -672,7 +533,9 @@ def test_tri_grid_linestring_on_outer_boundary(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(LineString([(15.0, 20.0), (5.0, 20.0)]))
+ result = ix.intersect(
+ LineString([(15.0, 20.0), (5.0, 20.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.lengths.sum() == 10.0
assert result.cellids[0] == 2
@@ -686,7 +549,9 @@ def test_tri_grid_linestring_on_inner_boundary(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(LineString([(5.0, 10.0), (15.0, 10.0)]))
+ result = ix.intersect(
+ LineString([(5.0, 10.0), (15.0, 10.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 2
assert result.lengths.sum() == 10.0
assert result.cellids[0] == 0
@@ -703,7 +568,8 @@ def test_tri_grid_multilinestring_in_one_cell(rtree):
result = ix.intersect(
MultiLineString(
[LineString([(1.0, 1), (9.0, 1.0)]), LineString([(2.0, 2.0), (9.0, 2.0)])]
- )
+ ),
+ geo_dataframe=df_toggle,
)
assert len(result) == 1
assert result.lengths == 15.0
@@ -723,7 +589,8 @@ def test_tri_grid_multilinestring_in_multiple_cells(rtree):
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)]),
]
- )
+ ),
+ geo_dataframe=df_toggle,
)
assert len(result) == 5
assert np.allclose(sum(result.lengths), 40.19197584109293)
@@ -733,39 +600,39 @@ 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):
ls = LineString([(x[i], y[i]), (x[i + 1], y[i + 1])])
- r = ix.intersect(ls, return_all_intersections=True)
+ r = ix.intersect(ls, return_all_intersections=True, geo_dataframe=df_toggle)
assert len(r) == n_intersections[i]
@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)
+ r = ix.intersect(ls, return_all_intersections=False, geo_dataframe=df_toggle)
assert len(r) == 1
@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)
+ r = ix.intersect(ls, return_all_intersections=True, geo_dataframe=df_toggle)
assert len(r) == 3
@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),
@@ -777,324 +644,152 @@ def test_rect_vertex_grid_linestring_geomcollection():
(10.0, 20.0),
]
)
- result = ix.intersect(ls)
+ result = ix.intersect(ls, geo_dataframe=df_toggle)
assert len(result) == 3
assert np.allclose(result.lengths.sum(), ls.length)
-# %% test polygon structured
-
-
@requires_pkg("shapely")
-def test_rect_grid_polygon_outside():
- # TODO: remove in 3.10.0
+@rtree_toggle
+def test_rect_grid_polygon_contains_centroid(rtree):
gr = get_rect_grid()
- ix = GridIntersect(gr, method="structured")
- result = ix.intersect(Polygon([(21.0, 11.0), (23.0, 17.0), (25.0, 11.0)]))
- assert len(result) == 0
+ 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, geo_dataframe=df_toggle)
+ assert len(result) == 1
@requires_pkg("shapely")
-def test_rect_grid_polygon_in_2cells():
- # TODO: remove in 3.10.0
+@rtree_toggle
+def test_rect_grid_polygon_min_area(rtree):
gr = get_rect_grid()
- ix = GridIntersect(gr, method="structured")
- result = ix.intersect(Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.0), (2.5, 15.0)]))
+ 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, geo_dataframe=df_toggle)
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_centroid_and_min_area(rtree):
gr = get_rect_grid()
- ix = GridIntersect(gr, method="structured")
+ 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(
- Polygon([(20.0, 5.0), (25.0, 5.0), (25.0, 15.0), (20.0, 15.0)])
+ p, min_area_fraction=0.35, contains_centroid=True, geo_dataframe=df_toggle
)
- assert len(result) == 0
+ assert len(result) == 1
+
+
+# %% test polygon shapely
@requires_pkg("shapely")
-def test_rect_grid_polygon_running_along_boundary():
- # 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(
- [(5.0, 5.0), (5.0, 10.0), (9.0, 10.0), (9.0, 15.0), (1.0, 15.0), (1.0, 5.0)]
- )
+ Polygon([(21.0, 11.0), (23.0, 17.0), (25.0, 11.0)]), geo_dataframe=df_toggle
)
+ assert len(result) == 0
@requires_pkg("shapely")
-def test_rect_grid_polygon_on_inner_boundary():
- # 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")
- result = ix.intersect(Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]))
+ 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)]),
+ geo_dataframe=df_toggle,
+ )
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
+@rtree_toggle
+def test_rect_grid_polygon_on_outer_boundary(rtree):
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, rtree=rtree)
+ result = ix.intersect(
+ Polygon([(20.0, 5.0), (25.0, 5.0), (25.0, 15.0), (20.0, 15.0)]),
+ geo_dataframe=df_toggle,
)
-
- 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()
+ assert len(result) == 0
@requires_pkg("shapely")
-def test_rect_grid_multiple_disjoint_polygons_on_inner_boundaries():
- # TODO: remove in 3.10.0
+def test_rect_grid_polygon_running_along_boundary():
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()
+ 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)]
+ ),
+ geo_dataframe=df_toggle,
+ )
@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
- gr = get_rect_grid()
- ix = GridIntersect(gr, method="structured")
- 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")
-def test_rect_grid_multipolygon_in_multiple_cells():
- # TODO: remove in 3.10.0
- gr = get_rect_grid()
- ix = GridIntersect(gr, method="structured")
- 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")
-def test_rect_grid_polygon_with_hole():
- # TODO: remove in 3.10.0
- 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
+@rtree_toggle
+def test_rect_grid_polygon_on_inner_boundary(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")
-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)]
- )
+ Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]),
+ geo_dataframe=df_toggle,
)
-
-
-@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):
+def test_rect_grid_multipolygon_in_one_cell(rtree):
gr = get_rect_grid()
- ix = GridIntersect(gr, method="vertex", rtree=rtree)
+ 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])
- result = ix.intersect(p)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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):
+def test_rect_grid_multipolygon_in_multiple_cells(rtree):
gr = get_rect_grid()
- ix = GridIntersect(gr, method="vertex", rtree=rtree)
+ 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])
- result = ix.intersect(p)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 2
assert result.areas.sum() == 72.0
@requires_pkg("shapely")
@rtree_toggle
-def test_rect_grid_polygon_with_hole_shapely(rtree):
+def test_rect_grid_polygon_with_hole(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)]],
)
- result = ix.intersect(p)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 3
assert result.areas.sum() == 104.0
@@ -1103,11 +798,11 @@ 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)]
)
- result = ix.intersect(p)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 1
assert result.areas.sum() == 15.0
@@ -1119,7 +814,9 @@ def test_tri_grid_polygon_outside(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(Polygon([(21.0, 11.0), (23.0, 17.0), (25.0, 11.0)]))
+ result = ix.intersect(
+ Polygon([(21.0, 11.0), (23.0, 17.0), (25.0, 11.0)]), geo_dataframe=df_toggle
+ )
assert len(result) == 0
@@ -1130,7 +827,10 @@ def test_tri_grid_polygon_in_2cells(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(Polygon([(2.5, 5.0), (5.0, 5.0), (5.0, 15.0), (2.5, 15.0)]))
+ result = ix.intersect(
+ Polygon([(2.5, 5.0), (5.0, 5.0), (5.0, 15.0), (2.5, 15.0)]),
+ geo_dataframe=df_toggle,
+ )
assert len(result) == 2
assert result.areas.sum() == 25.0
@@ -1143,7 +843,8 @@ def test_tri_grid_polygon_on_outer_boundary(rtree):
return
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)])
+ Polygon([(20.0, 5.0), (25.0, 5.0), (25.0, 15.0), (20.0, 15.0)]),
+ geo_dataframe=df_toggle,
)
assert len(result) == 0
@@ -1155,7 +856,10 @@ def test_tri_grid_polygon_on_inner_boundary(rtree):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
- result = ix.intersect(Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]))
+ result = ix.intersect(
+ Polygon([(5.0, 10.0), (15.0, 10.0), (15.0, 5.0), (5.0, 5.0)]),
+ geo_dataframe=df_toggle,
+ )
assert len(result) == 4
assert result.areas.sum() == 50.0
@@ -1170,7 +874,7 @@ def test_tri_grid_multipolygon_in_one_cell(rtree):
p1 = Polygon([(1.0, 1.0), (8.0, 1.0), (8.0, 3.0), (3.0, 3.0)])
p2 = Polygon([(5.0, 5.0), (8.0, 5.0), (8.0, 8.0)])
p = MultiPolygon([p1, p2])
- result = ix.intersect(p)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 1
assert result.areas.sum() == 16.5
@@ -1185,7 +889,7 @@ def test_tri_grid_multipolygon_in_multiple_cells(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)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 4
assert result.areas.sum() == 72.0
@@ -1201,7 +905,7 @@ def test_tri_grid_polygon_with_hole(rtree):
[(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)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 6
assert result.areas.sum() == 104.0
@@ -1217,7 +921,7 @@ def test_tri_grid_polygon_min_area(rtree):
[(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.5)
+ result = ix.intersect(p, min_area_fraction=0.5, geo_dataframe=df_toggle)
assert len(result) == 2
@@ -1232,7 +936,7 @@ def test_tri_grid_polygon_contains_centroid(rtree):
[(5.0, 5.0), (6.0, 14.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, contains_centroid=True)
+ result = ix.intersect(p, contains_centroid=True, geo_dataframe=df_toggle)
assert len(result) == 2
@@ -1240,36 +944,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")
- result = ix.intersect(p)
+ ix = GridIntersect(sgr, rtree=rtree)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 1
# check empty result when using local model coords
- ix = GridIntersect(sgr, method="structured", local=True)
- result = ix.intersect(p)
+ ix = GridIntersect(sgr, rtree=rtree, local=True)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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")
- result = ix.intersect(ls)
- assert len(result) == 3
+ 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, geo_dataframe=df_toggle)
+ assert len(result) == 2
# check empty result when using local model coords
- ix = GridIntersect(sgr, method="structured", local=True)
- result = ix.intersect(ls)
+ ix = GridIntersect(sgr, rtree=rtree, local=True)
+ result = ix.intersect(ls, geo_dataframe=df_toggle)
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 +983,47 @@ def test_polygon_offset_rot_structured_grid():
(5, 10.0 + 1.5 * np.sqrt(200.0)),
]
)
- ix = GridIntersect(sgr, method="structured")
- result = ix.intersect(p)
+ ix = GridIntersect(sgr, rtree=rtree)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
assert len(result) == 3
# check empty result when using local model coords
- ix = GridIntersect(sgr, method="structured", local=True)
- result = ix.intersect(p)
+ ix = GridIntersect(sgr, rtree=rtree, local=True)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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)
- result = ix.intersect(p)
+ ix = GridIntersect(sgr, rtree=rtree)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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)
+ ix = GridIntersect(sgr, rtree=rtree, local=True)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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)
- result = ix.intersect(ls)
+ ix = GridIntersect(sgr, rtree=rtree)
+ result = ix.intersect(ls, geo_dataframe=df_toggle)
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)
+ ix = GridIntersect(sgr, rtree=rtree, local=True)
+ result = ix.intersect(ls, geo_dataframe=df_toggle)
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 +1032,797 @@ 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)
- result = ix.intersect(p)
+ ix = GridIntersect(sgr, rtree=rtree)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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)
+ ix = GridIntersect(sgr, rtree=rtree, local=True)
+ result = ix.intersect(p, geo_dataframe=df_toggle)
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, dataframe=df_toggle)
assert len(result) == 1
- # check empty result when using local model coords
- ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
- result = ix.intersect(p)
- assert len(result) == 0
+ assert result.cellids[0] == (1, 0)
@requires_pkg("shapely")
@rtree_toggle
-def test_linestring_offset_rot_vertex_grid_shapely(rtree):
- sgr = get_rect_vertex_grid(angrot=45.0, xyoffset=10.0)
- ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))])
- ix = GridIntersect(sgr, method="vertex", rtree=rtree)
- result = ix.intersect(ls)
- assert len(result) == 2
- # check empty result when using local model coords
- ix = GridIntersect(sgr, method="vertex", rtree=rtree, local=True)
- result = ix.intersect(ls)
+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, dataframe=df_toggle)
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 pd.isna(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_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)),
- ]
+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, dataframe=df_toggle)
+ assert len(result) == 2
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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 pd.isna(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(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)
+ 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=False)
+ assert result.cellids[0] == (1, 0)
+ assert pd.isna(result.cellids[1])
+ result = ix.points_to_cellids(pts, handle_z=True)
+ assert pd.isna(result.layer[0])
+ assert pd.isna(result.cellids[1])
+
+
+# %% array-based input - structured grid linestrings
+
+
+@requires_pkg("shapely")
+@rtree_toggle
+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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ assert len(result) == 2
+ assert (result.shp_id == 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+
+
+@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, geo_dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+
+
+@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, geo_dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+
+
+# %% 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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 pd.isna(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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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 pd.isna(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=True)
+ assert result.layer[0] == 0.0
+ assert pd.isna(result.cellids[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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ assert len(result) == 2
+ assert (result.shp_id == 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+ assert len(result.cellids) == 2
+ assert result.cellids[0] == 2
+ assert result.cellids[1] == 7
+ result = ix.intersect(pts, return_all_intersections=False, geo_dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+
+
+@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, geo_dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+
+
+@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, geo_dataframe=df_toggle)
+ 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, geo_dataframe=df_toggle)
+
+
+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)
diff --git a/autotest/test_model_splitter.py b/autotest/test_model_splitter.py
index 895a87aca9..624c01b7ee 100644
--- a/autotest/test_model_splitter.py
+++ b/autotest/test_model_splitter.py
@@ -1008,12 +1008,13 @@ 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,
],
shapetype="Polygon",
+ geo_dataframe=False,
)
r, c = list(zip(*list(result.cellids)))
idomain = np.zeros(modelgrid.shape, dtype=int)
@@ -1033,7 +1034,12 @@ def string2geom(geostring, conversion=None):
lengths = []
for sg in stream_segs:
sg = string2geom(sg)
- v = ixs.intersect(sg, shapetype="LineString", sort_by_cellid=True)
+ v = ixs.intersect(
+ sg,
+ shapetype="LineString",
+ sort_by_cellid=True,
+ geo_dataframe=False,
+ )
cellids += v["cellids"].tolist()
lengths += v["lengths"].tolist()
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,
diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py
index 9effe75dc1..f8c9f46c8a 100644
--- a/flopy/utils/gridintersect.py
+++ b/flopy/utils/gridintersect.py
@@ -1,27 +1,18 @@
import warnings
+import matplotlib.pyplot as plt
import numpy as np
-from pandas import DataFrame
+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 NA, DataFrame
-from .geometry import transform
from .geospatial_utils import GeoSpatialUtil
from .utl_import import import_optional_dependency
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
@@ -72,35 +63,33 @@ class GridIntersect:
Notes
-----
+ - Supports structured and vertex grids. No support for unstructured grids.
+ If grid is layered-unstructured, creating a single layer vertex-grid can be
+ used as a workaround.
+ - For linestrings and polygons only 2D intersection operations are supported.
+ - Point intersections can optionally return layer position based on the
+ z-coordinate.
+ - The STR-tree can be disabled, but this is generally not recommended as some
+ functions will not work without it.
- The STR-tree query is based on the bounding box of the shape or
collection, if the bounding box of the shape covers nearly the entire
grid, the query won't be able to limit the search space much, resulting
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,69 +99,75 @@ 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
+ self.rtree = rtree
+
+ # build arrays of geoms and cellids
+ 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()
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,
+ raise NotImplementedError(
+ f"Grid type {self.mfgrid.grid_type} not supported"
)
- 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()
+ # 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)
- if self.method == "vertex":
- # build arrays of geoms and cellids
- self.geoms, self.cellids = self._get_gridshapes()
+ def _parse_input_shape(self, shp, shapetype=None):
+ """Internal method to parse input shape.
- # 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)
+ Allows numpy arrays containing shapely geometries; otherwise delegates to
+ ``GeoSpatialUtil`` for parsing.
- 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
+ Parameters
+ ----------
+ shp : shapely geometry, geojson object, shapefile.Shape, np.ndarray,
+ or a FloPy geometry object
+ Shape to intersect with the grid. If an ``np.ndarray`` is provided, all
+ elements must be shapely geometries of the same type.
+ shapetype : str, optional
+ Type of shape ("point", "linestring", "polygon" or their
+ multi-variants). Used by ``GeoSpatialUtil`` if shp is passed as a list
+ of vertices. Default is None.
- else:
- raise ValueError(
- f"Method '{self.method}' not recognized or not supported "
- f"for grid_type '{self.mfgrid.grid_type}'!"
+ Returns
+ -------
+ tuple
+ (geom, gtype) where geom is a shapely geometry or an array of
+ shapely geometries, and gtype is the corresponding shapely
+ GeometryType (or a matching string value).
+ """
+ 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,
shapetype=None,
sort_by_cellid=True,
- keepzerolengths=False,
return_all_intersections=False,
contains_centroid=False,
min_area_fraction=None,
- geo_dataframe=False,
- shapely2=None,
+ handle_z=False,
+ geo_dataframe=None,
):
- """Method to intersect a shape with a model grid.
+ """Intersect a shape with a model grid.
Parameters
----------
@@ -185,9 +180,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
@@ -202,71 +194,117 @@ 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 : bool, optional
+ Method for handling z-dimension in intersection results for point
+ intersections. Default is False which ignores z-dimension. If True
+ returns the layer index for each point. Points below/above the grid
+ are returned as masked values or pd.NA if returned as dataframe.
geo_dataframe : bool, optional
- if True, return a geopandas GeoDataFrame, default is False
+ If True, return a geopandas ``GeoDataFrame``, otherwise return a
+ numpy recarray. Default will be set to True in the future;
+ currently, None triggers a deprecation warning and returns a
+ recarray (legacy behavior).
Returns
-------
- numpy.recarray or gepandas.GeoDataFrame
- a record array containing information about the intersection or
- a geopandas.GeoDataFrame if geo_dataframe=True
+ numpy.recarray or geopandas.GeoDataFrame
+ Intersection results. Result contains the following columns:
+ - cellid: cellid of intersected grid cell
+ - cellids: legacy column containing (row, col) tuples for structured grids,
+ or cellids for vertex grids (deprecated, use "cellid" column instead)
+ - ixshapes or geometry: shapely geometry of the intersection result
+ And optionally:
+ - layer: layer index of for points with z-dimension when handle_z is True
+ - row: row index of intersected grid cell, for structured grids
+ - col: column index of intersected grid cell, for structured grids
+ - lengths: length of intersection results for linestrings
+ - areas: areas of intersection results for polygons
+
"""
- if shapely2 is not None:
- warnings.warn(
- "The shapely2 keyword argument is deprecated. "
- "Shapely<2 support was dropped in flopy version 3.9.0."
+ shp, shapetype = self._parse_input_shape(shp, shapetype=shapetype)
+
+ # 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."
)
- 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,
- )
- 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,
- )
- 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,
+ 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,
+ )
+
+ # handle elevation data for points
+ # if handle_z is True
+ # if shp has z information
+ # if there are intersection results
+ if handle_z and shapely.has_z(shp).any() and len(rec.cellid) > 0:
+ laypos = self.get_layer_from_z(shp, rec.cellid)
+ # copy data to new array to include layer position
+ rec = nprecfns.append_fields(
+ rec,
+ names="layer",
+ data=laypos,
+ dtypes=int,
+ usemask=False,
+ asrecarray=True,
)
+ rec = np.ma.masked_where(rec.layer < 0, rec)
+
+ 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 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 is None:
+ warnings.warn(
+ "In the future this function will return a GeoDataFrame by default. "
+ "Set geo_dataframe=True to adopt future behavior and silence this "
+ "warning. Set geo_dataframe=False to silence this warning and maintain "
+ "old behavior",
+ DeprecationWarning,
+ )
if geo_dataframe:
gpd = import_optional_dependency("geopandas")
+ if gpd is None:
+ raise ModuleNotFoundError(
+ "GeoDataFrame output requires geopandas to be installed."
+ )
gdf = (
gpd.GeoDataFrame(rec)
.rename(columns={"ixshapes": "geometry"})
.set_geometry("geometry")
+ .replace(999999, NA)
)
if self.mfgrid.crs is not None:
gdf = gdf.set_crs(self.mfgrid.crs)
@@ -274,17 +312,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.
@@ -335,7 +362,6 @@ def _rect_grid_to_geoms_cellids(self):
indices=np.repeat(cellids, 4),
)
)
-
return geoms, cellids
def _usg_grid_to_geoms_cellids(self):
@@ -383,63 +409,112 @@ def _vtx_grid_to_geoms_cellids(self):
]
return np.array(geoms), np.arange(self.mfgrid.ncpl)
- def query_grid(self, shp):
- """Perform spatial query on grid with shapely geometry. If no spatial
- query is possible returns all grid cells.
+ def query_grid(self, shp, predicate=None):
+ """Perform spatial query on the grid using a shapely geometry.
+
+ If no spatial query is possible (``rtree=False``), returns all grid cells.
Parameters
----------
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
-------
- array_like
- array containing cellids of grid cells in query result
+ numpy.ndarray
+ For a single geometry, a 1-D array of cellids. For an array of geometries,
+ a 2xN array where the first row is the input geometry index (``shp_id``)
+ and the second row contains the matching cellids.
"""
if self.rtree:
- result = self.strtree.query(shp)
+ result = self.strtree.query(shp, predicate=predicate)
else:
# no spatial query
result = self.cellids
return result
- def filter_query_result(self, cellids, shp):
- """Filter array of geometries to obtain grid cells that intersect with
- shape.
+ def filter_query_result(self, shp, cellids):
+ """Filter a query result to cells that truly intersect a shape.
- Used to (further) reduce query result to cells that intersect with
- shape.
+ Used to (further) reduce the query result to cells that intersect with
+ shp. This is primarily used when ``rtree=False``.
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
+ numpy.ndarray
+ Array of cellids that intersect shp.
"""
+ # 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
- def _intersect_point_shapely(
+ def _intersect_point_shapely(self, *args, **kwargs):
+ """Deprecated method, use _intersect_point instead."""
+ 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,
return_all_intersections=False,
):
- if self.rtree:
- qcellids = self.strtree.query(shp, predicate="intersects")
- else:
- qcellids = self.filter_query_result(self.cellids, shp)
+ """Intersect a Point or MultiPoint with the grid.
+
+ Parameters
+ ----------
+ shp : shapely Point or MultiPoint (or array with a single point)
+ Geometry to intersect.
+ sort_by_cellid : bool, optional
+ If True (default), sort candidate cells by cellid
+ return_all_intersections : bool, optional
+ If True, return multiple intersection results for points on cell boundaries.
+ Default is False.
+
+ Returns
+ -------
+ numpy.recarray
+ Intersection results. Result contains the following columns:
+ - cellid: cellid of intersected grid cell
+ - cellids: legacy column containing (row, col) tuples for structured grids,
+ or cellids for vertex grids (deprecated, use "cellid" column instead)
+ - ixshapes or geometry: shapely geometry of the intersection result
+ And optionally:
+ - layer: layer index of for points with z-dimension when handle_z is True
+ - row: row index of intersected grid cell, for structured grids
+ - col: column index of intersected grid cell, for structured grids
+ """
+ r = self.intersects(shp, dataframe=False)
+ qcellids = r.cellid[r.cellid >= 0]
if sort_by_cellid:
qcellids = np.sort(qcellids)
@@ -448,7 +523,10 @@ def _intersect_point_shapely(
# 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]
@@ -474,39 +552,71 @@ def _intersect_point_shapely(
keep_pts = ixresult
keep_cid = qcellids
- names = ["cellids", "ixshapes"]
- formats = ["O", "O"]
+ if self.mfgrid.grid_type == "structured":
+ names = ["cellids", "cellid", "row", "col", "ixshapes"]
+ formats = ["O", int, int, int, "O"]
+ else:
+ names = ["cellids", "cellid", "ixshapes"]
+ formats = [int, int, "O"]
rec = np.recarray(len(keep_pts), names=names, formats=formats)
+ rec.cellid = self.cellids[keep_cid]
+ rec.ixshapes = keep_pts
# 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, rec.row, rec.col = self._cellid_to_rowcol(
+ self.cellids[keep_cid]
)
else:
- rec.cellids = self.cellids[keep_cid]
- rec.ixshapes = keep_pts
-
+ rec.cellids = rec.cellid # NOTE: legacy support for cellids column
return rec
- def _intersect_linestring_shapely(
+ def _intersect_linestring_shapely(self, *args, **kwargs):
+ """Deprecated method, use _intersect_linestring instead."""
+
+ warnings.warn(
+ "_intersect_linestring_shapely is deprecated, "
+ "use _intersect_linestring instead.",
+ DeprecationWarning,
+ )
+ return self._intersect_linestring(*args, **kwargs)
+
+ def _intersect_linestring(
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,
- )
+ """Intersect a LineString or MultiLineString with the grid.
+
+ Parameters
+ ----------
+ shp : shapely LineString or MultiLineString
+ Geometry to intersect.
+ sort_by_cellid : bool, optional
+ If True (default), sort candidate cells by cellid before
+ intersecting.
+ return_all_intersections : bool, optional
+ If True, keep overlapping boundary segments as separate results.
+ Default is False.
+ Returns
+ -------
+ numpy.recarray
+ Intersection results. Result contains the following columns:
+ - cellid: cellid of intersected grid cell
+ - cellids: legacy column containing (row, col) tuples for structured grids,
+ or cellids for vertex grids (deprecated, use "cellid" column instead)
+ - ixshapes or geometry: shapely geometry of the intersection result
+ And optionally:
+ - row: row index of intersected grid cell, for structured grids
+ - col: column index of intersected grid cell, for structured grids
+ - lengths: length of intersection results for linestrings
+ """
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)
@@ -595,33 +705,80 @@ def parse_linestrings_in_geom_collection(gc):
ixresult = ixresult[mask_keep]
qcellids = qcellids[mask_keep]
- names = ["cellids", "ixshapes", "lengths"]
- formats = ["O", "O", "f8"]
+ if self.mfgrid.grid_type == "structured":
+ names = ["cellids", "cellid", "row", "col", "ixshapes", "lengths"]
+ formats = ["O", int, int, int, "O", float]
+ else:
+ names = ["cellids", "cellid", "ixshapes", "lengths"]
+ formats = ["O", int, "O", float]
rec = np.recarray(len(ixresult), names=names, formats=formats)
+ rec.cellid = self.cellids[qcellids]
+ rec.ixshapes = ixresult
+ rec.lengths = shapely.length(ixresult)
+
# 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, rec.row, rec.col = self._cellid_to_rowcol(
+ self.cellids[qcellids]
)
else:
- rec.cellids = self.cellids[qcellids]
- rec.ixshapes = ixresult
- rec.lengths = shapely.length(ixresult)
+ rec.cellids = rec.cellid # NOTE: legacy support for cellids column
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,
contains_centroid=False,
min_area_fraction=None,
):
+ """Intersect a Polygon or MultiPolygon with the grid.
+
+ Parameters
+ ----------
+ shp : shapely Polygon or MultiPolygon
+ Geometry to intersect.
+ sort_by_cellid : bool, optional
+ If True (default), sort candidate cells by cellid before
+ intersecting.
+ contains_centroid : bool, optional
+ If True, only keep results where the cell centroid is contained in
+ (or touches) the intersection. Default is False.
+ min_area_fraction : float, optional
+ Minimum fractional area threshold relative to the cell area. Cells with
+ an intersection area smaller than ``min_area_fraction * cell_area`` are
+ discarded. Default is None (no threshold).
+
+ Returns
+ -------
+ numpy.recarray
+ Intersection results. Result contains the following columns:
+ - cellid: cellid of intersected grid cell
+ - cellids: legacy column containing (row, col) tuples for structured grids,
+ or cellids for vertex grids (deprecated, use "cellid" column instead)
+ - ixshapes or geometry: shapely geometry of the intersection result
+ And optionally:
+ - layer: layer index of for points with z-dimension when handle_z is True
+ - row: row index of intersected grid cell, for structured grids
+ - col: column index of intersected grid cell, for structured grids
+ - areas: areas of intersection results for polygons
+ """
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)
@@ -673,23 +830,34 @@ def parse_polygons_in_geom_collection(gc):
qcellids = qcellids[mask_area_frac]
# fill rec array
- names = ["cellids", "ixshapes", "areas"]
- formats = ["O", "O", "f8"]
+ if self.mfgrid.grid_type == "structured":
+ names = ["cellids", "cellid", "row", "col", "ixshapes", "areas"]
+ formats = ["O", int, int, int, "O", float]
+ else:
+ names = ["cellids", "cellid", "ixshapes", "areas"]
+ formats = ["O", int, "O", float]
rec = np.recarray(len(ixresult), names=names, formats=formats)
+ rec.cellid = self.cellids[qcellids]
+ rec.ixshapes = ixresult
+ rec.areas = shapely.area(ixresult)
+
# 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, rec.row, rec.col = self._cellid_to_rowcol(
+ self.cellids[qcellids]
)
else:
- rec.cellids = self.cellids[qcellids]
- rec.ixshapes = ixresult
- rec.areas = shapely.area(ixresult)
+ rec.cellids = rec.cellid # NOTE: legacy support for cellids column
return rec
- def intersects(self, shp, shapetype=None, dataframe=False):
- """Return cellids for grid cells that intersect with shape.
+ def intersects(
+ self,
+ shp,
+ shapetype=None,
+ dataframe=None,
+ ):
+ """Return candidate grid cellids that intersect with shape(s).
Parameters
----------
@@ -701,952 +869,216 @@ def intersects(self, shp, shapetype=None, dataframe=False):
their multi-variants), used by GeoSpatialUtil if shp is
passed as a list of vertices, default is None
dataframe : bool, optional
- if True, return a pandas.DataFrame, default is False
+ If True, return a ``pandas.DataFrame``; otherwise return a numpy
+ recarray. Default will be set to True in the future; currently,
+ None triggers a deprecation warning and returns a recarray (legacy
+ behavior).
Returns
-------
numpy.recarray or pandas.DataFrame
- a record array or pandas.DataFrame containing cell IDs of the gridcells
- the shape intersects with.
+ Grid cell candidates for intersections. Result contains the following
+ columns:
+ - cellid: cellid of candidate grid cell
+ - cellids: legacy column containing (row, col) tuples for structured grids,
+ or cellids for vertex grids (deprecated, use "cellid" column instead)
+ And optionally:
+ - row: row index of intersected grid cell, for structured grids
+ - col: column index of intersected grid cell, for structured grids
"""
- shp = GeoSpatialUtil(shp, shapetype=shapetype).shapely
- qfiltered = self.strtree.query(shp, predicate="intersects")
-
- # 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:]))
+ 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:
- rec.cellids = qfiltered
-
- if dataframe:
- 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
+ qfiltered = qcellids
- if isinstance(shp, shapely_geo.Point):
- shp = [shp]
- elif isinstance(shp, shapely_geo.MultiPoint):
- shp = list(shp.geoms)
+ # sort cellids
+ if qfiltered.ndim == 1:
+ qfiltered = np.sort(qfiltered)
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
+ qfiltered = qfiltered[:, np.lexsort((qfiltered[1], qfiltered[0]))]
- 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.
+ # determine size of output array
+ nr = len(qfiltered) if qfiltered.ndim == 1 else qfiltered.shape[1]
- 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])
+ # build rec-array
+ # use float dtype to allow nans in row/col/cellid
+ if self.mfgrid.grid_type == "structured":
+ names = ["shp_id", "cellids", "cellid", "row", "col"]
+ formats = [int, "O", int, int, int]
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,
- )
+ names = ["shp_id", "cellids", "cellid"]
+ formats = [int, int, int]
+ rec = np.recarray(nr, names=names, formats=formats)
+
+ # shp was passed as single geometry
+ if qfiltered.ndim == 1:
+ rec.shp_id[:] = 0
+ rec.cellid = qfiltered
+ # shape passed as array of geometries
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)
+ rec.shp_id = qfiltered[0]
+ rec.cellid = qfiltered[1]
+
+ if self.mfgrid.grid_type == "structured":
+ rec.cellids, rec.row, rec.col = self._cellid_to_rowcol(rec.cellid)
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
- )
+ rec.cellids = rec.cellid # NOTE: legacy support for cellids column
+ if dataframe is None:
+ warnings.warn(
+ "In the future this function will return a dataframe by default. "
+ "Set dataframe=True to adopt future behavior and silence this warning. "
+ "Set dataframe=False to silence this warning and maintain old behavior",
+ DeprecationWarning,
)
+ if dataframe:
+ return DataFrame(rec).set_index("shp_id")
+ return rec
- 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.
+ def _cellid_to_rowcol(self, cellids):
+ """Convert cellid (node number) to row, col.
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
+ cellids : array_like
+ array of cellids to convert
Returns
-------
- node, length, verts: lists
- lists containing nodes, lengths and vertices of
- intersections with adjacent cells relative to the
- current cell (i, j)
+ tuple of array_like
+ array of (row, col) tuples, array of rows, array of columns
"""
- 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))
+ idx = np.nonzero(cellids >= 0)
+ row = np.full_like(cellids, -1, dtype=int) # -1 is invalid
+ col = np.full_like(cellids, -1, dtype=int) # -1 is invalid
+ row[idx], col[idx] = self.mfgrid.get_lrc([cellids[idx]])[0][1:]
+ # NOTE: build tuple for backward compatibility
+ rctuple = np.full_like(cellids, -1, dtype=object)
+ rctuple[idx] = list(zip(row[idx], col[idx]))
+ return rctuple, row, col
+
+ def points_to_cellids(
+ self,
+ pts,
+ handle_z=False,
+ dataframe=True,
+ ):
+ """Return cellids for points intersecting the grid.
Parameters
----------
- rectangle : list of tuples
- list of lower-left coordinate and upper-right
- coordinate: [(xmin, ymin), (xmax, ymax)]
+ pts : shapely geometry, geojson geometry, shapefile.Shape, np.ndarray,
+ or FloPy geometry object
+ Point(s) to intersect with the grid. array inputs must contain
+ shapely point (or multipoint) geometries.
+ handle_z : bool, optional
+ Handle z-dimension for points. If True, returns a "layer" column with
+ the computed layer index for each point (NA is returned where the
+ points lie above/below the model grid). Default is False.
+ dataframe : bool, optional
+ If True, return a ``pandas.DataFrame``; otherwise return a numpy
+ recarray. Default is True.
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()`.
-
+ pandas.DataFrame or numpy.recarray
+ Grid cell indices per point. Result contains the following
+ columns:
+ - cellid: cellid of grid cell containing point
+ - cellids: legacy column containing (row, col) tuples for structured grids,
+ or cellids for vertex grids (deprecated, use "cellid" column instead)
+ And optionally:
+ - row: row index of intersected grid cell, for structured grids
+ - col: column index of intersected grid cell, for structured grids
+ - layer: layer index of for points with z-dimension when handle_z is True
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
+ Requires ``rtree=True`` when initializing ``GridIntersect``.
"""
- 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 not self.rtree:
+ raise ValueError(
+ "points_to_cellids() requires rtree=True when"
+ " initializing GridIntersect"
)
- 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]),
- ]
+ if not isinstance(pts, np.ndarray):
+ if shapely.get_type_id(pts) == shapely.GeometryType.MULTIPOINT:
+ pts = np.array(shapely.get_parts(pts))
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
+ pts = np.array([pts])
- return rec
+ # query grid or strtree
+ qfiltered = self.query_grid(pts, predicate="intersects")
- def _transform_geo_interface_polygon(self, polygon):
- """Internal method, helper function to transform geometry
- __geo_interface__.
+ # sort cellids
+ if qfiltered.ndim == 1:
+ qfiltered = np.sort(qfiltered)
+ else:
+ qfiltered = qfiltered[:, np.lexsort((qfiltered[1], qfiltered[0]))]
- .. deprecated:: 3.9.0
- method="structured" is deprecated. Only used by
- `_intersect_polygon_structured()`
+ # 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
- Used for translating intersection result coordinates back into
- real-world coordinates.
+ # build rec-array
+ if self.mfgrid.grid_type == "structured":
+ names = ["shp_id", "cellids", "cellid", "row", "col"]
+ formats = [int, "O", int, int, int] # float to allow nans in row/col
+ else:
+ names = ["shp_id", "cellids", "cellid"]
+ formats = [int, int, int] # float to allow nans
+ rec = np.recarray(nr, names=names, formats=formats)
- Parameters
- ----------
- polygon : shapely.geometry.Polygon
- polygon to transform coordinates for
+ rec.cellid = -1 # invalid value by default
+ # return only 1 gr id cell intersection result for each point
+ uniques, idx = np.unique(qfiltered[0], return_index=True)
+ rec.shp_id = np.arange(nr)
+ rec.cellid[uniques] = qfiltered[1, idx]
- Returns
- -------
- geom_list : list
- list containing transformed coordinates in same structure as
- the original __geo_interface__.
- """
+ if self.mfgrid.grid_type == "structured":
+ rec.cellids, rec.row, rec.col = self._cellid_to_rowcol(rec.cellid)
+ else:
+ rec.cellids = rec.cellid # NOTE: legacy support for cellids column
+
+ if handle_z and shapely.has_z(pts).any() and len(rec.cellid) > 0:
+ laypos = self.get_layer_from_z(pts, rec.cellid)
+ # copy data to new array to include layer position
+ rec = nprecfns.append_fields(
+ rec,
+ names="layer",
+ data=laypos,
+ dtypes=int,
+ usemask=False,
+ asrecarray=True,
+ )
- 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
+ if dataframe:
+ # replace invalid indices with NA, replace invalid layer with NA
+ return (
+ DataFrame(rec).replace(-1, NA).replace(999_999, NA).set_index("shp_id")
+ )
+ return np.ma.masked_where(rec.cellid < 0, rec)
@staticmethod
- def plot_polygon(result, ax=None, **kwargs):
- """method to plot the polygon intersection results from the resulting
- numpy.recarray.
-
- Note: only works when recarray has 'ixshapes' column!
+ def plot_polygon(polys, ax=None, **kwargs):
+ """Plot polygons.
Parameters
----------
- result : numpy.recarray or geopandas.GeoDataFrame
- record array or GeoDataFrame containing intersection results
+ polys : collection of polygons
+ list, array or GeoSeries containing polygons
ax : matplotlib.pyplot.axes, optional
axes to plot onto, if not provided, creates a new figure
**kwargs:
@@ -1657,10 +1089,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")
@@ -1681,11 +1109,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)
@@ -1701,16 +1128,13 @@ def add_poly_patch(poly):
return ax
@staticmethod
- def plot_linestring(result, ax=None, cmap=None, **kwargs):
- """method to plot the linestring intersection results from the
- resulting numpy.recarray.
-
- Note: only works when recarray has 'ixshapes' column!
+ def plot_linestring(ls, ax=None, cmap=None, **kwargs):
+ """Plot linestrings.
Parameters
----------
- result : numpy.recarray or geopandas.GeoDataFrame
- record array or GeoDataFrame containing intersection results
+ ls : collection of linestrings
+ list, array or GeoSeries containing linestrings
ax : matplotlib.pyplot.axes, optional
axes to plot onto, if not provided, creates a new figure
cmap : str
@@ -1723,7 +1147,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()
@@ -1737,39 +1162,31 @@ 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}"
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
@staticmethod
- def plot_point(result, ax=None, **kwargs):
- """method to plot the point intersection results from the resulting
- numpy.recarray.
-
- Note: only works when recarray has 'ixshapes' column!
+ def plot_point(pts, ax=None, **kwargs):
+ """Plot points.
Parameters
----------
- result : numpy.recarray or geopandas.GeoDataFrame
- record array or GeoDataFrame containing intersection results
+ pts : collection of points
+ array or GeoSeries containing point geometries
ax : matplotlib.pyplot.axes, optional
axes to plot onto, if not provided, creates a new figure
**kwargs:
@@ -1780,24 +1197,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
+ if isinstance(pts, (shapely.Point, shapely.MultiPoint)):
+ pts = [pts]
+
+ maskpts = np.isin(
+ shapely.get_type_id(pts),
+ [shapely.GeometryType.POINT, shapely.GeometryType.MULTIPOINT],
)
- 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)
+ shapely_plot.plot_points(pts[maskpts], ax=ax, **kwargs)
return ax
@@ -1832,168 +1245,64 @@ 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),
[
shapely.GeometryType.LINESTRING,
shapely.GeometryType.MULTILINESTRING,
+ 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
-
-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.
+ def get_layer_from_z(self, pts, cellids):
+ """Compute layer indices from point z-values.
Parameters
----------
- arr : A one dimensional array (such as Xe) that contains
- coordinates for the left cell edge.
+ pts : shapely geometry
+ Points geometry (single or array-like) with z-values.
+ cellids : array_like
+ Array of candidate cellids for the points.
- x : float
- The x position to find in arr.
+ Returns
+ -------
+ numpy.ndarray
+ layer index for each point. Points below/above the grid are returned
+ with index -1.
"""
- 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]
+ z_arr = np.atleast_1d(shapely.get_z(pts))
+ mask_valid = cellids >= 0
+ if self.mfgrid.grid_type == "structured":
+ row, col = self.mfgrid.get_lrc([cellids[mask_valid].astype(int)])[0][1:]
+ surface_elevations = self.mfgrid.top_botm[:, row, col]
+ elif self.mfgrid.grid_type == "vertex":
+ surface_elevations = self.mfgrid.top_botm[:, cellids[mask_valid]]
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)
+ 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(int)
+ laypos[mask_above] = -1
+ laypos[mask_below] = -1
+ laypos_full = np.full_like(z_arr, -1, dtype=int)
+ laypos_full[mask_valid] = laypos
+ return laypos_full
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]),