更快速的多边形交叉方式

Question

我有大量的多边形(~100000),并尝试找到一种智能的方法来计算它们与常规网格单元的交叉区域.

目前,我正在使用形状(基于它们的角坐标)创建多边形和网格单元.然后,使用一个简单的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分钟,具体取决于交叉网格单元的数量),这是不可接受的.不幸的是,我不知道如何编写多边形交集的代码来获得重叠区域.你有什么建议吗？有一种形状的替代品吗？

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用户手册 已经过时了,但自2013/2014以来,Shapely已经使用了 STRtree 类的strtree.py.我用过它似乎运作良好.

这是docstring的一个片段:

STRtree是使用Sort-Tile-Recursive算法创建的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