使用shapely的多边形交集更快的方法

Question

我有很多多边形（大约10万个），想找到一个聪明的方法来计算它们与规则网格单元的交集面积。

目前，我是用shapely这个工具来创建多边形和网格单元，主要是根据它们的角落坐标来做的。然后，我用一个简单的循环，逐个检查每个多边形和附近的网格单元。

这里有个小例子来说明多边形和网格单元的关系。

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分钟以上），这个速度是不可接受的。可惜的是，我不知道怎么写代码来计算多边形的交集，以得到重叠的面积。你有什么建议吗？有没有比shapely更好的选择？

Answer 1

这是另一个版本的答案，不过对IoU的控制更多一些。

def merge_intersecting_polygons(list_of_polygons, image_width, image_height):
    """Merge intersecting polygons with shapely library.

    Merge only if Intersection over Union is greater than 0.5 
    speed up with STRTree.
    """
    # create shapely polygons
    shapely_polygons = []
    for polygon in list_of_polygons:
        shapely_polygons.append(Polygon(polygon))
    # create STRTree
    tree = STRtree(shapely_polygons)
    # merge polygons
    merged_polygons = []
    for i, polygon in enumerate(shapely_polygons):
        # find intersecting polygons
        intersecting_polygons = tree.query(polygon)
        # merge intersecting polygons
        for intersecting_polygon in intersecting_polygons:
            if polygon != intersecting_polygon:
                # compute intersection over union
                intersection = polygon.intersection(intersecting_polygon).area
                union = polygon.union(intersecting_polygon).area
                iou = intersection/union
                if iou > 0.5:
                    # merge polygons
                    polygon = polygon.union(intersecting_polygon)
        # add merged polygon to list
        merged_polygons.append(polygon)
    # remove duplicates
    merged_polygons = list(set(merged_polygons))
    # convert shapely polygons to list of polygons
    list_of_polygons = []
    for polygon in merged_polygons:
        list_of_polygons.append(np.array(polygon.exterior.coords).tolist())
    return list_of_polygons

Answer 2

自从2013/2014年，Shapely就有了一个叫做STRtree的功能。我用过这个功能，感觉效果不错。

下面是文档中的一段说明：

STRtree是一种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 3

可以考虑使用 Rtree 来帮助找出一个多边形可能会与哪些网格单元相交。这样，你就可以去掉用经纬度数组进行的循环，这部分可能是比较慢的。

你可以把代码结构安排成下面这样：

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)