使用shapely进行多边形交集的更快方法
使用shapely进行多边形交集的更快方法
我有大量的多边形(约100,000个),并试图找到一种智能的方法来计算它们与正规网格单元的相交面积。\n目前,我使用shapely创建多边形和网格单元(基于它们的角坐标)。然后,使用简单的for循环,我遍历每个多边形并将其与附近的网格单元进行比较。\n这里是一个小例子,用来说明多边形/网格单元的情况。\n
from shapely.geometry import box, Polygon # 例子多边形 xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]] polygon_shape = Polygon(xy) # 例子网格单元 gridcell_shape = box(129.5, -27.0, 129.75, 27.25) # 相交面积 polygon_shape.intersection(gridcell_shape).area
\n(顺便说一下:网格单元的尺寸为0.25x0.25,多边形的最大尺寸为1x1)\n对于一个单独的多边形/网格单元组合,这实际上是相当快的,大约需要0.003秒。然而,如果在大量的多边形上运行此代码(每个多边形可能与几十个网格单元相交),则需要大约15分钟以上(取决于相交网格单元的数量),这在我的机器上是不可接受的。不幸的是,我不知道如何编写一个计算多边形相交并得到重叠面积的代码。你有什么建议吗?是否有shapely的替代方案?
快速计算多边形交集的问题出现的原因是需要在速度较快的情况下计算多边形的交集。为了解决这个问题,可以使用shapely库来合并相交的多边形,并且只有当交并比大于0.5时才进行合并。使用STRTree可以加速计算过程。
解决方法是首先将输入的多边形列表转换成shapely多边形对象。然后创建一个STRTree来管理这些多边形。接下来,遍历每个多边形,找到与其相交的其他多边形。对于每个相交的多边形,计算其与当前多边形的交集并计算交并比。如果交并比大于0.5,则合并这两个多边形。最后,将合并后的多边形转换回列表形式并返回。
以下是具体的代码实现:
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
通过以上方法,可以快速计算多边形的交集,并且只合并交并比大于0.5的多边形。使用STRTree可以提高计算速度。
Shapely自2013/2014年以来就有了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
这非常有帮助。你知道STRtree是否可以使用pickle或marshall库进行序列化,以便以后使用?
不,我不熟悉STRtree的序列化能力。我相信它完全依赖于`shapely.geos.GEOSSTRtree_create(max(2, len(geoms)))`返回的_tree_handle的序列化。
使用Rtree可以帮助识别多边形可能相交的网格单元格。这样,您可以删除使用纬度/经度数组的for循环,这可能是速度较慢的部分。
代码结构如下:
from shapely.ops import cascaded_union from rtree import index idx = index.Index() # 用网格单元的边界填充R树索引 for pos, cell in enumerate(grid_cells): # 假设cell是一个shapely对象 idx.insert(pos, cell.bounds) # 遍历每个Shapely多边形 for poly in polygons: # 合并具有重叠边界框的单元格 merged_cells = cascaded_union([grid_cells[pos] for pos in idx.intersection(poly.bounds)]) # 现在进行实际的相交操作 print(poly.intersection(merged_cells).area)
这仍然是一个非常有用的答案-它应该被接受。我遇到了类似的问题,Rtree使算法运行速度提高了约5000倍。
请注意,Rtree只能用于矩形(4个点),而不是复杂的多边形。
对于“真实”的多边形,只需为每个边界相交添加适当的几何形状相交检查即可。Rtree可以帮助您减少搜索空间,速度非常快。
虽然这是一个老问题的答案,但我发现这是一个提醒我使用rtree的结果。如果您有一个复杂的多边形,只需将其包裹起来,然后先与之进行检查,然后再进行复杂的检查。这样会再次加快速度。对于简单的示例,我给出了+1。
顺便说一句,shapely有一个Rtree实现,如下面的答案所示:shapely.readthedocs.io/en/stable/manual.html#str-packed-r-tree