多角形とぴったりのやり方

Question

多角形（〜100000）があり、通常の格子セルとの交差領域を計算するスマートな方法を見つけようとします。

現在、ポリゴンとグリッドセルを（コーナー座標に基づいて）滑らかに作成しています。次に、単純なfor-loopを使って、各ポリゴンを通り、近くのグリッドセルと比較します。

ポリゴン/グリッドセルを示す小さな例です。

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 

（BTW：グリッドセルは、寸法0.25x0.25とmaxのポリゴンの1x1のを持っている）

実はこれは周り0.003秒で、個々のポリゴン/グリッドセルのコンボのために非常に高速です。しかし、このコードを膨大な量のポリゴン（各グリッドセルに交差する可能性があります）で実行すると、マシン上で約15 +分（交差グリッドセルの数に応じて最大30+分）がかかります。残念ながら、ポリゴンの交差のコードを記述してオーバーラップ領域を取得することは、どのように可能であるかわかりません。何かヒントはありますか？洗練されたものに代わるものはありますか？

Answer 1

Rtreeを使用して、ポリゴンが交差する可能性のあるグリッドセルを特定することを検討してください。こうすることで、おそらく遅い部分であるlat/lonsの配列で使用される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 

Answer 2

むしろ、期限が切れている利用可能 The Shapely User Manual ように見えるが、2014分の2013以来、格好の良いクラスSTRtreeと strtree.py がありました。私はそれを使いました、そして、それはうまくいくようです。ここ

は、ドキュメント文字列からの抜粋である：

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

>>> from shapely.geometry import Polygon 
>>> 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