使用Shapely实现多边形相交的更快方法

Question

使用Shapely实现多边形相交的更快方法

53

我有大量的多边形(~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分钟或更长时间（取决于相交格网单元的数量），这是不可接受的。不幸的是，我不知道如何编写用于多边形相交以获取重叠区域面积的代码。你有什么技巧吗？是否有比shapely更好的替代方案？

- HyperCube

1

我很好奇你是如何循环和相交你的多边形的。你能展示更多关于这个过程的代码吗？这样会更容易找出如何进行优化。 - tdedecko

我基本上会取一个纬度/经度角落值的数组，并在for循环中将它们转换为多边形。然后，我再次在for循环中将每个多边形与特定的网格单元进行比较。请参见此链接：https://dev59.com/I2Yr5IYBdhLWcg3wH23K#13956110 - HyperCube

3个回答

32

自2013/2014年以来，Shapely已经包含了STRtree。我使用过它，似乎运行良好。

以下是文档字符串的摘录：

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

- Phil

2

这非常有帮助。您知道STRtree是否可以使用pickle或marshall库进行序列化，以便稍后使用吗？ - eguaio

1

不，我不熟悉STRtree的序列化能力。我认为它完全依赖于shapely.geos.GEOSSTRtree_create(max(2, len(geoms)))返回的_tree_handle的序列化。 - Phil

0

这里是另一个答案的版本，但对 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

```

- Achyut Sarma

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mike T · Accepted Answer

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