使用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]))) 实际上是在匹配两个多边形的边界框。另外要注意，像 Point(0.5, 0.5) 这样的边界上的点，使用 within 是不会匹配到任何东西的，但使用 intersects 时会匹配到两个多边形。所以要做好准备，可能会匹配到0个、1个或多个多边形。