Optimal way of aggregating geographic points with Python/Shapely

Question

I would like to convert a long list of lat/long coordinates to the US state (or county) they belong to. One possible solution, given that I have the state geometries, is to check each points against all the states.

for point in points:
    for state in states:
        if point.within(state['shape']):
            print state.name

Is there a more optimized way to do this, possibly in O(1)?

Answer 1

Use Rtree as a spatial index to very quickly identify the point in a bounding box of zero or more polygons, then use Shapely to determine the polygon that the point is within.

Similar to this example https://stackoverflow.com/a/14804366/327026

from shapely.geometry import Polygon, Point
from rtree import index

# List of non-overlapping polygons
polygons = [
    Polygon([(0, 0), (0, 1), (1, 1), (0, 0)]),
    Polygon([(0, 0), (1, 0), (1, 1), (0, 0)]),
]

# Populate R-tree index with bounds of polygons
idx = index.Index()
for pos, poly in enumerate(polygons):
    idx.insert(pos, poly.bounds)

# Query a point to see which polygon it is in
# using first Rtree index, then Shapely geometry's within
point = Point(0.5, 0.2)
poly_idx = [i for i in idx.intersection((point.coords[0]))
            if point.within(polygons[i])]
for num, idx in enumerate(poly_idx, 1):
    print("%d:%d:%s" % (num, idx, polygons[idx]))

If you dissect the list comprehension, you will see that list(idx.intersection((point.coords[0]))) actually matches the bounding boxes of two polygons. Also, take note that points on the boundary like Point(0.5, 0.5) will not match anything with within, but will match both for intersects. So be prepared to match 0, 1 or more polygons.