Intersecting model grids with shapes
Note: This feature requires the shapely package (which is an optional FloPy dependency).
This notebook shows the grid intersection functionality in flopy. The intersection methods are available through the GridIntersect object. A flopy modelgrid is passed to instantiate the object. Then the modelgrid can be intersected with Points, LineStrings and Polygons and their Multi variants.
Table of Contents
GridIntersect Class
Rectangular regular grid
Polygon with regular grid
MultiLineString with regular grid
MultiPoint with regular grid
Vertex grid
Polygon with triangular grid
MultiLineString with triangular grid
MultiPoint with triangular grid
Import packages
[1]:
import sys
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import shapely
from shapely.geometry import LineString, MultiLineString, MultiPoint, Point, Polygon
import flopy
import flopy.discretization as fgrid
import flopy.plot as fplot
from flopy.utils import GridIntersect
print(sys.version)
print(f"numpy version: {np.__version__}")
print(f"matplotlib version: {mpl.__version__}")
print(f"shapely version: {shapely.__version__}")
print(f"flopy version: {flopy.__version__}")
3.12.8 | packaged by conda-forge | (main, Dec 5 2024, 14:24:40) [GCC 13.3.0]
numpy version: 2.2.1
matplotlib version: 3.10.0
shapely version: 2.0.6
flopy version: 3.10.0.dev1
GridIntersect Class
The GridIntersect class is constructed by passing a flopy modelgrid object to the constructor. There are options users can select to change how the intersection is calculated.
rtree: eitherTrue(default) orFalse. 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 recommended to set this option to True.local: eitherFalse(default) orTrue. When True the local model coordinates are used. When False, real-world coordinates are used. Can be useful if shapes are defined in local coordinates.
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)intersect(): for intersecting the modelgrid with point, linestrings, and polygon geometries (accepts shapely geometry objects, flopy geometry object, shapefile.Shape objects, and geojson objects)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).
Rectangular regular grid
[2]:
delc = 10 * np.ones(10, dtype=float)
delr = 10 * np.ones(10, dtype=float)
[3]:
xoff = 0.0
yoff = 0.0
angrot = 0.0
sgr = fgrid.StructuredGrid(
delc, delr, top=None, botm=None, xoff=xoff, yoff=yoff, angrot=angrot
)
[4]:
sgr.plot()
[4]:
<matplotlib.collections.LineCollection at 0x7f1325b433e0>
Polygon with regular grid
Polygon to intersect with:
[5]:
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)]],
)
[6]:
# Create the GridIntersect class for our modelgrid.
# TODO: remove method kwarg in v3.9.0
ix = GridIntersect(sgr, method="vertex")
Do the intersect operation for a polygon
[7]:
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)
Looking at the first few entries of the results of the polygon intersection (convert to pandas.DataFrame for prettier formatting)
[8]:
result[:5]
# pd.DataFrame(result).head()
[8]:
rec.array([((np.int64(2), np.int64(3)), <POLYGON ((35 80, 40 76.667, 40 70, 30 70, 35 80))>, 66.66666667),
((np.int64(2), np.int64(4)), <POLYGON ((50 70, 40 70, 40 76.667, 50 70))>, 33.33333333),
((np.int64(3), np.int64(2)), <POLYGON ((30 70, 30 60, 25 60, 30 70))>, 25. ),
((np.int64(3), np.int64(3)), <POLYGON ((40 70, 40 60, 30 60, 30 70, 40 70))>, 100. ),
((np.int64(3), np.int64(4)), <POLYGON ((50 70, 50 60, 40 60, 40 70, 50 70))>, 100. )],
dtype=[('cellids', 'O'), ('ixshapes', 'O'), ('areas', '<f8')])
The cellids can be easily obtained
[9]:
result.cellids
[9]:
array([(np.int64(2), np.int64(3)), (np.int64(2), np.int64(4)),
(np.int64(3), np.int64(2)), (np.int64(3), np.int64(3)),
(np.int64(3), np.int64(4)), (np.int64(3), np.int64(5)),
(np.int64(3), np.int64(6)), (np.int64(4), np.int64(2)),
(np.int64(4), np.int64(3)), (np.int64(4), np.int64(4)),
(np.int64(4), np.int64(5)), (np.int64(4), np.int64(6)),
(np.int64(4), np.int64(7)), (np.int64(5), np.int64(1)),
(np.int64(5), np.int64(2)), (np.int64(5), np.int64(3)),
(np.int64(5), np.int64(4)), (np.int64(5), np.int64(5)),
(np.int64(5), np.int64(6)), (np.int64(5), np.int64(7)),
(np.int64(6), np.int64(1)), (np.int64(6), np.int64(2)),
(np.int64(6), np.int64(4)), (np.int64(6), np.int64(5)),
(np.int64(6), np.int64(6)), (np.int64(6), np.int64(7)),
(np.int64(7), np.int64(1)), (np.int64(7), np.int64(2)),
(np.int64(7), np.int64(3)), (np.int64(7), np.int64(4)),
(np.int64(7), np.int64(5)), (np.int64(7), np.int64(6)),
(np.int64(8), np.int64(1)), (np.int64(8), np.int64(2)),
(np.int64(8), np.int64(3)), (np.int64(8), np.int64(4)),
(np.int64(8), np.int64(5)), (np.int64(9), np.int64(2)),
(np.int64(9), np.int64(3)), (np.int64(9), np.int64(4))],
dtype=object)
Or the areas
[10]:
result.areas
[10]:
array([ 66.66666667, 33.33333333, 25. , 100. ,
100. , 66.66666667, 8.33333333, 75. ,
100. , 100. , 100. , 91.66666667,
33.33333333, 7.14285714, 75. , 50. ,
75. , 100. , 100. , 100. ,
21.42857143, 50. , 50. , 100. ,
99.10714286, 43.75 , 35.71428571, 75. ,
50. , 75. , 96.42857143, 32.14285714,
41.71428571, 99.35714286, 100. , 91.96428571,
22.32142857, 8.64285714, 36. , 14.28571429])
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.
[11]:
ix.intersects(p)
[11]:
rec.array([((np.int64(8), np.int64(1)),), ((np.int64(7), np.int64(2)),),
((np.int64(7), np.int64(1)),), ((np.int64(6), np.int64(1)),),
((np.int64(6), np.int64(2)),), ((np.int64(9), np.int64(4)),),
((np.int64(9), np.int64(2)),), ((np.int64(9), np.int64(3)),),
((np.int64(8), np.int64(2)),), ((np.int64(8), np.int64(3)),),
((np.int64(8), np.int64(4)),), ((np.int64(7), np.int64(3)),),
((np.int64(7), np.int64(4)),), ((np.int64(6), np.int64(4)),),
((np.int64(5), np.int64(1)),), ((np.int64(4), np.int64(1)),),
((np.int64(4), np.int64(2)),), ((np.int64(3), np.int64(2)),),
((np.int64(5), np.int64(2)),), ((np.int64(5), np.int64(3)),),
((np.int64(5), np.int64(4)),), ((np.int64(4), np.int64(3)),),
((np.int64(4), np.int64(4)),), ((np.int64(3), np.int64(4)),),
((np.int64(3), np.int64(3)),), ((np.int64(2), np.int64(4)),),
((np.int64(2), np.int64(3)),), ((np.int64(2), np.int64(2)),),
((np.int64(1), np.int64(3)),), ((np.int64(8), np.int64(5)),),
((np.int64(7), np.int64(5)),), ((np.int64(7), np.int64(6)),),
((np.int64(6), np.int64(5)),), ((np.int64(6), np.int64(6)),),
((np.int64(6), np.int64(7)),), ((np.int64(6), np.int64(8)),),
((np.int64(5), np.int64(5)),), ((np.int64(5), np.int64(7)),),
((np.int64(5), np.int64(6)),), ((np.int64(4), np.int64(6)),),
((np.int64(4), np.int64(5)),), ((np.int64(3), np.int64(6)),),
((np.int64(3), np.int64(5)),), ((np.int64(2), np.int64(5)),),
((np.int64(5), np.int64(8)),), ((np.int64(4), np.int64(8)),),
((np.int64(4), np.int64(7)),)],
dtype=[('cellids', 'O')])
The results of an intersection can be visualized with the GridIntersect.plot_intersection_result() method.
[12]:
# 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)
# add black x at cell centers
for irow, icol in result.cellids:
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
"kx",
label="centroids of intersected gridcells",
)
# add legend
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
[12]:
<matplotlib.legend.Legend at 0x7f12fec1b800>
The intersect() method contains several keyword arguments that specifically deal with polygons:
contains_centroid: only store intersection result if cell centroid is contained within polygonmin_area_fraction: minimal intersecting cell area (expressed as a fraction of the total cell area) to include cells in intersection result
Two examples showing the usage of these keyword arguments are shown below.
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.
[13]:
# contains_centroid example
result2 = ix.intersect(p, contains_centroid=True)
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
ix.plot_intersection_result(result2, ax=ax)
# add black x at cell centers
for irow, icol in result2.cellids:
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
"kx",
label="centroids of intersected gridcells",
)
# add legend
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
[13]:
<matplotlib.legend.Legend at 0x7f12fea1eff0>
Example with min_area_threshold set to 0.35, the intersection result in a cell should cover 35% or more of the cell area.
[14]:
# min_area_threshold example
result3 = ix.intersect(p, min_area_fraction=0.35)
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
ix.plot_intersection_result(result3, ax=ax)
# add black x at cell centers
for irow, icol in result3.cellids:
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
"kx",
label="centroids of intersected gridcells",
)
# add legend
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
[14]:
<matplotlib.legend.Legend at 0x7f12feab9c70>
Polyline with regular grid
MultiLineString to intersect with:
[15]:
ls1 = LineString([(95, 105), (30, 50)])
ls2 = LineString([(30, 50), (90, 22)])
ls3 = LineString([(90, 22), (0, 0)])
mls = MultiLineString(lines=[ls1, ls2, ls3])
[16]:
result = ix.intersect(mls)
Plot the result
[17]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
sgr.plot(ax=ax)
ix.plot_intersection_result(result, ax=ax, cmap="viridis")
for irow, icol in result.cellids:
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
"kx",
label="centroids of intersected gridcells",
)
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
[17]:
<matplotlib.legend.Legend at 0x7f12f67d36b0>
MultiPoint with regular grid
MultiPoint to intersect with
[18]:
mp = MultiPoint(
points=[Point(50.0, 0.0), Point(45.0, 45.0), Point(10.0, 10.0), Point(150.0, 100.0)]
)
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.
[19]:
result = ix.intersect(mp)
result_all = ix.intersect(mp, return_all_intersections=True)
[20]:
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")
for irow, icol in result.cellids:
(h2,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
"kx",
ms=15,
label="centroids of intersected cells",
)
for irow, icol in result_all.cellids:
(h3,) = ax.plot(
sgr.xcellcenters[0, icol],
sgr.ycellcenters[irow, 0],
"C3+",
ms=15,
label="centroids with `return_all_intersections=True`",
)
ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best")
[20]:
<matplotlib.legend.Legend at 0x7f12f66fd760>
Vertex Grid
[21]:
cell2d = [
[0, 83.33333333333333, 66.66666666666667, 3, 4, 2, 7],
[1, 16.666666666666668, 33.333333333333336, 3, 4, 0, 5],
[2, 33.333333333333336, 83.33333333333333, 3, 1, 8, 4],
[3, 16.666666666666668, 66.66666666666667, 3, 5, 1, 4],
[4, 33.333333333333336, 16.666666666666668, 3, 6, 0, 4],
[5, 66.66666666666667, 16.666666666666668, 3, 4, 3, 6],
[6, 83.33333333333333, 33.333333333333336, 3, 7, 3, 4],
[7, 66.66666666666667, 83.33333333333333, 3, 8, 2, 4],
]
vertices = [
[0, 0.0, 0.0],
[1, 0.0, 100.0],
[2, 100.0, 100.0],
[3, 100.0, 0.0],
[4, 50.0, 50.0],
[5, 0.0, 50.0],
[6, 50.0, 0.0],
[7, 100.0, 50.0],
[8, 50.0, 100.0],
]
tgr = fgrid.VertexGrid(vertices, cell2d)
[22]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
pmv = fplot.PlotMapView(modelgrid=tgr)
pmv.plot_grid(ax=ax)
[22]:
<matplotlib.collections.LineCollection at 0x7f12f6699190>
Polygon with triangular grid
[23]:
ix2 = GridIntersect(tgr)
[24]:
result = ix2.intersect(p)
[25]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
ix.plot_intersection_result(result, ax=ax)
# only cells that intersect with shape
for cellid in result.cellids:
(h2,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
"kx",
label="centroids of intersected gridcells",
)
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
[25]:
<matplotlib.legend.Legend at 0x7f12f673fbc0>
LineString with triangular grid
[26]:
result = ix2.intersect(mls)
[27]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
ix2.plot_intersection_result(result, ax=ax, lw=3)
for cellid in result.cellids:
(h2,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
"kx",
label="centroids of intersected gridcells",
)
ax.legend([h2], [i.get_label() for i in [h2]], loc="best")
[27]:
<matplotlib.legend.Legend at 0x7f12f65d64b0>
MultiPoint with triangular grid
[28]:
result = ix2.intersect(mp)
result_all = ix2.intersect(mp, return_all_intersections=True)
[29]:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
ix2.plot_intersection_result(result, ax=ax, color="k", zorder=5, s=80)
for cellid in result.cellids:
(h2,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
"kx",
ms=15,
label="centroids of intersected cells",
)
for cellid in result_all.cellids:
(h3,) = ax.plot(
tgr.xcellcenters[cellid],
tgr.ycellcenters[cellid],
"r+",
ms=15,
label="centroids with return_all_intersections=True",
)
ax.legend([h2, h3], [i.get_label() for i in [h2, h3]], loc="best")
[29]:
<matplotlib.legend.Legend at 0x7f13306a45c0>
