Skip to content

Commit

Permalink
perf(Gridintersect): optimize intersection methods for shapely 2.0 (#…
Browse files Browse the repository at this point in the history 
…1666)
  • Loading branch information
@dbrakenhoff
dbrakenhoff authored Dec 23, 2022
1 parent 19b3daa commit 1114c93
Show file tree
Hide file tree
Showing 4 changed files with 526 additions and 299 deletions.
3 changes: 2 additions & 1 deletion CITATION.cff
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,8 @@ authors:
given-names: Jon Jeffrey
orcid: https://orcid.org/0000-0001-5909-0010
- family-names: Brakenhoff
given-names: Davíd
given-names: Davíd A.
orcid: https://orcid.org/0000-0002-2993-2202
- family-names: Bonelli
given-names: Wesley P.
orcid: https://orcid.org/0000-0002-2665-5078
Expand Down
154 changes: 47 additions & 107 deletions autotest/test_gridintersect.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,42 +120,6 @@ def get_rect_vertex_grid(angrot=0.0, xyoffset=0.0):
return tgr


def plot_structured_grid(sgr):
_, ax = plt.subplots(1, 1, figsize=(8, 8))
sgr.plot(ax=ax)
return ax


def plot_vertex_grid(tgr):
_, ax = plt.subplots(1, 1, figsize=(8, 8))
pmv = fplot.PlotMapView(modelgrid=tgr)
pmv.plot_grid(ax=ax)
return ax


def plot_ix_polygon_result(rec, ax):
from descartes import PolygonPatch

for i, ishp in enumerate(rec.ixshapes):
ppi = PolygonPatch(ishp, facecolor=f"C{i % 10}")
ax.add_patch(ppi)


def plot_ix_linestring_result(rec, ax):
for i, ishp in enumerate(rec.ixshapes):
if ishp.type == "MultiLineString":
for part in ishp:
ax.plot(part.xy[0], part.xy[1], ls="-", c=f"C{i % 10}")
else:
ax.plot(ishp.xy[0], ishp.xy[1], ls="-", c=f"C{i % 10}")


def plot_ix_point_result(rec, ax):
x = [ip.x for ip in rec.ixshapes]
y = [ip.y for ip in rec.ixshapes]
ax.scatter(x, y)


# %% test point structured


Expand Down Expand Up @@ -597,6 +561,26 @@ def test_rect_grid_linestrings_on_boundaries_return_all_ix_shapely(rtree):
assert len(r) == n_intersections[i]


@requires_pkg("shapely")
@rtree_toggle
def test_rect_grid_linestring_cell_boundary_shapely(rtree):
gr = get_rect_grid()
ix = GridIntersect(gr, method="vertex", rtree=rtree)
ls = LineString(ix._rect_grid_to_geoms_cellids()[0][0].exterior.coords)
r = ix.intersect(ls, return_all_intersections=False)
assert len(r) == 1


@requires_pkg("shapely")
@rtree_toggle
def test_rect_grid_linestring_cell_boundary_return_all_ix_shapely(rtree):
gr = get_rect_grid()
ix = GridIntersect(gr, method="vertex", rtree=rtree)
ls = LineString(ix._rect_grid_to_geoms_cellids()[0][0].exterior.coords)
r = ix.intersect(ls, return_all_intersections=True)
assert len(r) == 3


@requires_pkg("shapely")
@rtree_toggle
def test_tri_grid_linestring_outside(rtree):
Expand Down Expand Up @@ -681,7 +665,26 @@ def test_tri_grid_linestrings_on_boundaries_return_all_ix(rtree):
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]
return


@requires_pkg("shapely")
@rtree_toggle
def test_tri_grid_linestring_cell_boundary_shapely(rtree):
tgr = get_tri_grid()
ix = GridIntersect(tgr, method="vertex", rtree=rtree)
ls = LineString(ix._vtx_grid_to_geoms_cellids()[0][0].exterior.coords)
r = ix.intersect(ls, return_all_intersections=False)
assert len(r) == 1


@requires_pkg("shapely")
@rtree_toggle
def test_tri_grid_linestring_cell_boundary_return_all_ix_shapely(rtree):
tgr = get_tri_grid()
ix = GridIntersect(tgr, method="vertex", rtree=rtree)
ls = LineString(ix._vtx_grid_to_geoms_cellids()[0][0].exterior.coords)
r = ix.intersect(ls, return_all_intersections=True)
assert len(r) == 3


# %% test polygon structured
Expand Down Expand Up @@ -1064,7 +1067,7 @@ def test_point_offset_rot_structured_grid():
p = Point(10.0, 10 + np.sqrt(200.0))
ix = GridIntersect(sgr, method="structured")
result = ix.intersect(p)
# assert len(result) == 1.
assert len(result) == 1


@requires_pkg("shapely")
Expand All @@ -1073,7 +1076,8 @@ def test_linestring_offset_rot_structured_grid():
ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))])
ix = GridIntersect(sgr, method="structured")
result = ix.intersect(ls)
# assert len(result) == 2.
# NOTE: in shapely 2.0, this returns a Linestring with length 10^-15 in cell (0, 1)
assert len(result) == 2 or len(result) == 3


@requires_pkg("shapely")
Expand All @@ -1089,7 +1093,7 @@ def test_polygon_offset_rot_structured_grid():
)
ix = GridIntersect(sgr, method="structured")
result = ix.intersect(p)
# assert len(result) == 3.
assert len(result) == 3


@requires_pkg("shapely")
Expand All @@ -1099,7 +1103,7 @@ def test_point_offset_rot_structured_grid_shapely(rtree):
p = Point(10.0, 10 + np.sqrt(200.0))
ix = GridIntersect(sgr, method="vertex", rtree=rtree)
result = ix.intersect(p)
# assert len(result) == 1.
assert len(result) == 1


@requires_pkg("shapely")
Expand All @@ -1109,7 +1113,7 @@ def test_linestring_offset_rot_structured_grid_shapely(rtree):
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.
assert len(result) == 2


@requires_pkg("shapely")
Expand All @@ -1126,66 +1130,7 @@ def test_polygon_offset_rot_structured_grid_shapely(rtree):
)
ix = GridIntersect(sgr, method="vertex", rtree=rtree)
result = ix.intersect(p)
# assert len(result) == 3.


# %% test non strtree shapely intersect


@requires_pkg("shapely")
def test_all_intersections_shapely_no_strtree():
"""avoid adding separate tests for rtree=False"""
# Points
# regular grid
test_rect_grid_point_on_inner_boundary_shapely(rtree=False)
test_rect_grid_point_on_outer_boundary_shapely(rtree=False)
test_rect_grid_point_outside_shapely(rtree=False)
test_rect_grid_multipoint_in_one_cell_shapely(rtree=False)
test_rect_grid_multipoint_in_multiple_cells_shapely(rtree=False)
# vertex grid
test_tri_grid_point_on_inner_boundary(rtree=False)
test_tri_grid_point_on_outer_boundary(rtree=False)
test_tri_grid_point_outside(rtree=False)
test_tri_grid_multipoint_in_multiple_cells(rtree=False)
test_tri_grid_multipoint_in_one_cell(rtree=False)

# LineStrings
# regular grid
test_rect_grid_linestring_on_inner_boundary_shapely(rtree=False)
test_rect_grid_linestring_on_outer_boundary_shapely(rtree=False)
test_rect_grid_linestring_outside_shapely(rtree=False)
test_rect_grid_linestring_in_2cells_shapely(rtree=False)
test_rect_grid_linestring_in_and_out_of_cell_shapely(rtree=False)
test_rect_grid_multilinestring_in_one_cell_shapely(rtree=False)
# vertex grid
test_tri_grid_linestring_on_inner_boundary(rtree=False)
test_tri_grid_linestring_on_outer_boundary(rtree=False)
test_tri_grid_linestring_outside(rtree=False)
test_tri_grid_linestring_in_2cells(rtree=False)
test_tri_grid_multilinestring_in_one_cell(rtree=False)

# Polygons
# regular grid
test_rect_grid_polygon_on_inner_boundary_shapely(rtree=False)
test_rect_grid_polygon_on_outer_boundary_shapely(rtree=False)
test_rect_grid_polygon_outside_shapely(rtree=False)
test_rect_grid_polygon_in_2cells_shapely(rtree=False)
test_rect_grid_polygon_with_hole_shapely(rtree=False)
test_rect_grid_multipolygon_in_one_cell_shapely(rtree=False)
test_rect_grid_multipolygon_in_multiple_cells_shapely(rtree=False)
# vertex grid
test_tri_grid_polygon_on_inner_boundary(rtree=False)
test_tri_grid_polygon_on_outer_boundary(rtree=False)
test_tri_grid_polygon_outside(rtree=False)
test_tri_grid_polygon_in_2cells(rtree=False)
test_tri_grid_polygon_with_hole(rtree=False)
test_tri_grid_multipolygon_in_multiple_cells(rtree=False)
test_tri_grid_multipolygon_in_one_cell(rtree=False)

# offset and rotated grids
test_point_offset_rot_structured_grid_shapely(rtree=False)
test_linestring_offset_rot_structured_grid_shapely(rtree=False)
test_polygon_offset_rot_structured_grid_shapely(rtree=False)
assert len(result) == 3


# %% test rasters
Expand Down Expand Up @@ -1298,8 +1243,3 @@ def test_raster_sampling_methods(example_data_path):
raise AssertionError(
f"{method} resampling returning incorrect values"
)


if __name__ == "__main__":

test_all_intersections_shapely_no_strtree()
Loading

0 comments on commit 1114c93

Please sign in to comment.