Optimal way of aggregating geographic points with Python/Shapely
python, shapely
Solution
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.
Problem
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)?