python - 使用 Python/Shapely 聚合地理点的最佳方式

Question

我想将一长串经纬度坐标转换为它们所属的美国州（或县）。鉴于我有状态几何，一种可能的解决方案是针对所有状态检查每个点。

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

有没有更优化的方法来做到这一点，可能在 O(1) 中？

Answer 1

使用Rtree作为空间索引可以非常快速地识别零个或多个多边形的边界框中的点，然后使用 Shapely 确定该点所在的多边形。

与此示例类似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]))

如果你剖析列表推导，你会发现它list(idx.intersection((point.coords[0])))实际上匹配了两个多边形的边界框。另外，请注意，边界上的点 likePoint(0.5, 0.5)不会与 匹配within，但会与 匹配intersects。因此，请准备好匹配 0、1 个或更多多边形。