48

多数のポリゴン (〜 100000) があり、通常のグリッド セルを使用して交差領域を計算するスマートな方法を見つけようとしています。

現在、(コーナー座標に基づいて) shapely を使用してポリゴンとグリッド セルを作成しています。次に、単純な for ループを使用して各ポリゴンを調べ、近くのグリッド セルと比較します。

ポリゴン/グリッド セルを説明するためのほんの一例です。

from shapely.geometry import box, Polygon
# Example polygon 
xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]]
polygon_shape = Polygon(xy)
# Example grid cell
gridcell_shape = box(129.5, -27.0, 129.75, 27.25)
# The intersection
polygon_shape.intersection(gridcell_shape).area

(ところで: グリッド セルの寸法は 0.25x0.25 で、ポリゴンは最大で 1x1 です)

実際、これは個々のポリゴン/グリッド セルの組み合わせで約 0.003 秒と非常に高速です。ただし、このコードを大量のポリゴン (それぞれが数十のグリッド セルと交差する可能性がある) で実行すると、私のマシンでは約 15 分以上 (交差するグリッド セルの数によっては最大 30 分以上) かかります。残念ながら、ポリゴン交差のコードを記述して重複領域を取得する方法がわかりません。ヒントはありますか？格好良いに代わるものはありますか？

4

2 に答える 2

74

Rtreeを使用して、ポリゴンが交差する可能性のあるグリッド セルを特定することを検討してください。このようにして、緯度/経度の配列で使用される for ループを削除できます。これはおそらく遅い部分です。

コードを次のように構成します。

from shapely.ops import cascaded_union
from rtree import index
idx = index.Index()

# Populate R-tree index with bounds of grid cells
for pos, cell in enumerate(grid_cells):
    # assuming cell is a shapely object
    idx.insert(pos, cell.bounds)

# Loop through each Shapely polygon
for poly in polygons:
    # Merge cells that have overlapping bounding boxes
    merged_cells = cascaded_union([grid_cells[pos] for pos in idx.intersection(poly.bounds)])
    # Now do actual intersection
    print(poly.intersection(merged_cells).area)
于 2013-02-11T00:37:55.253 に答える
27

2013/2014 から Shapely には STRtreeがあります。私はそれを使用しましたが、うまく機能しているようです。

docstring のスニペットを次に示します。

STRtree は、Sort-Tile-Recursive アルゴリズムを使用して作成される R ツリーです。STRtree は、一連のジオメトリ オブジェクトを初期化パラメーターとして受け取ります。初期化後、query メソッドを使用して、これらのオブジェクトに対して空間クエリを作成できます。

>>> from shapely.geometry import Polygon
>>> from shapely.strtree import STRtree
>>> polys = [Polygon(((0, 0), (1, 0), (1, 1))), Polygon(((0, 1), (0, 0), (1, 0))), Polygon(((100, 100), (101, 100), (101, 101)))]
>>> s = STRtree(polys)
>>> query_geom = Polygon([(-1, -1), (2, 0), (2, 2), (-1, 2)])
>>> result = s.query(query_geom)
>>> polys[0] in result
True
于 2017-03-29T22:54:47.413 に答える