perf(Gridintersect): optimize intersection methods for shapely 2.0 (#1666)

Author: dbrakenhoff

CITATION.cff

@@ -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

autotest/test_gridintersect.py

@@ -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
 
 
@@ -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):
@@ -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
@@ -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")
@@ -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")
@@ -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")
@@ -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")
@@ -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")
@@ -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
@@ -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()