Efficiently check if polygon contains any point from a list

Question

I have a large, static list of Shapely Points that I want to check against a large list of Shapely Polygons, filtering out any polygons that contain one or more of the points.

The polygons do not overlap one another.

The same list of points is checked against all polygons.

A simplified example of the overall shape of the data is pictured below. The real data is magnified to a much larger extent.

Currently, I am using the Python any and the Shapely contains functions:

from shapely.geometry.point import Point
from shapely.geometry.polygon import Polygon

# A large list of thousands of points
points: list[Point] = []

# A large list of thousands of polygons
polygons: list[Polygon] = []
...

# Want to speed up this loop, which can take many seconds to complete.
# Retain polygons that do not contain any points.
retained_polygons: list[Polygon] = []
for polygon in polygons:
    if not any(polygon.contains(point) for point in points):
        retained_polygons.append(polygon)

return retained_polygons

For performance reasons, we were hoping there might be a native method within Shapely that could do this more efficiently. We have thousands of polygons to check against thousands of points, and have profiled slow down to this part of the code.

Answer 1

This is a typical case to use a spatial index.

Shapely has support for spatial indexes via shapely.STRtree.query

I did a quick test with 1000 polygons and 750 points derived of them... and even with those relatively small numbers using the rtree index speeds it up by a factor 1000. The larger the number... the bigger the gain.

• Retained 499 polygons with tree, took 0:00:00.006998
• Retained 499 polygons with loop, took 0:00:08.774927

Sample script:

import time
import geopandas as gpd
import numpy as np
import shapely
from shapely import Polygon


# Load some polygons as test data
poly_path = r"C:\temp\file.gpkg"
polygons_gdf = gpd.read_file(poly_path, columns=[], rows=1000, engine="pyogrio")
polygons = polygons_gdf.geometry.tolist()

# For the first half of the polygons, use the centroid as points
points = (polygons_gdf.loc[: len(polygons) / 2].geometry.centroid).tolist()
# For the first quarter of the polygons, add the representative_point as points as well,
# so we have multiple points for those polygons
points2 = polygons_gdf.loc[: len(polygons) / 2].geometry.representative_point().tolist()
points.extend(points2)

# Test the performance using spatial index
start = time.perf_counter()
polygons_tree = shapely.STRtree(polygons)
result = polygons_tree.query(points, predicate="intersects")
mask = np.zeros(len(polygons), dtype=bool)
mask[np.unique(result[1])] = True
retained_polygons_tree = polygons_tree.geometries[~mask]
elapsed = time.perf_counter() - start
print(f"Retained {len(retained_polygons_tree)} polygons with tree, took {elapsed:.6f}")

# Test the performance using a simple loop
start = time.perf_counter()
retained_polygons: list[Polygon] = []
for polygon in polygons:
    if not any(polygon.contains(point) for point in points):
        retained_polygons.append(polygon)
elapsed = time.perf_counter() - start
print(f"Retained {len(retained_polygons)} polygons with loop, took {elapsed:.6f}")

# Make sure the results are the same
assert all(retained_polygons_tree == retained_polygons)