4

緯度/経度座標の長いリストを、それらが属する米国の州 (または郡) に変換したいと思います。状態のジオメトリがある場合、考えられる解決策の 1 つは、すべての状態に対して各ポイントをチェックすることです。

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

おそらくO（1）で、これを行うためのより最適化された方法はありますか？

4

1 に答える 1

8

Rtreeを空間インデックスとして使用して、0 個以上のポリゴンのバウンディング ボックス内のポイントを非常に迅速に識別し、次に 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])))、2 つのポリゴンのバウンディング ボックスと実際に一致することがわかります。Point(0.5, 0.5)また、 のような境界上のポイントはとは何にも一致しませんwithinが、 の場合は両方に一致することに注意してくださいintersects。したがって、0、1、またはそれ以上のポリゴンに一致するように準備してください。

于 2014-06-09T06:10:06.027 に答える