更快速的多邊形相交的方式與形狀

我有大量的多邊形（〜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

（BTW：網格單元具有尺寸0.25x0.25和最大的多邊形的1x1）

其實這是相當快的周圍0.003秒一個單獨的多邊形/網格單元組合。但是，在我的機器上運行此代碼的大量多邊形（每個可能會交叉幾十個網格單元）在我的機器上花費大約15分鐘以上（取決於交叉網格單元的數量，最多爲30+分鐘），這是不可接受的。不幸的是，我不知道如何編寫多邊形交集的代碼來獲得重疊區域。你有什麼建議嗎？有沒有一個適合的選擇？

來源

2013-02-04 HyperCube

我很好奇你是如何循環和相交你的多邊形。你能在這個過程中展示更多的代碼嗎？弄清楚這可以如何優化會更容易。 – tdedecko

我基本上是取一個緯度/經度角值的數組，並將它們轉換爲for循環到多邊形。然後，我將每個多邊形與特定的網格單元進行比較，這是在for循環中完成的。看到這個：http://stackoverflow.com/a/13956110/1740928 – HyperCube

考慮使用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

來源

2013-02-11 00:37:55

這仍然是一個非常有用的答案 - 它應該被接受。我有一個類似的問題，'Rtree'使算法運行速度提高了5000倍。 – Gabriel

請注意，'Rtree'只能用於框（4個點），而不能用於複雜的多邊形。 –

看起來可用 The Shapely User Manual 是相當過時，但由於2013/2014，身材勻稱已 strtree.py 與類STRtree。我已經使用它，它似乎運作良好。

下面是從文檔字符串片段：

STRtree是使用排序瓦片遞歸算法創建一個R-tree。 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

來源

2017-03-29 22:54:47 Phil

這非常有用。你知道STRtree是否可以用pickle或marshall庫序列化以保存以備後用？ – eguaio

不，我不熟悉STRtree的序列化功能。我相信它完全依賴於由''shapely.geos.GEOSSTRtree_create（max（2，len（geoms））' – Phil

'返回的_tree_handle的序列化。請注意，更新的Shapely手冊位於此處：http： //shapely.readthedocs.io/en/stable/ –

更快速的多邊形相交的方式與形狀

回答

相關問題