多边形相交的更快方式

Question

我有大量的多边形（~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分钟），这是不可接受的。不幸的是，我不知道如何编写多边形相交的代码来获得重叠面积。你有什么建议吗？有没有一个替代品来塑造身材？

Answer 1

自2013/2014年以来，Shapely拥有STRtree。我已经用过了，而且看起来效果很好。

以下是docstring中的一段代码：

STRTree是一个使用Sort-Tile-Recursion算法创建的R树。STRTree将一系列几何对象作为初始化参数。初始化后，查询方法可用于对这些对象进行空间查询。

>>> 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

Answer 2

考虑使用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)