numpy 用shapely实现多边形求交的快速方法

yqlxgs2m  于 2023-04-21  发布在  其他
关注(0)|答案(3)|浏览(343)

我有大量的多边形(~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

(BTW:网格单元格的尺寸为0.25x0.25,多边形的尺寸最大为1x 1)
实际上,这对于一个单独的多边形/网格单元组合来说是相当快的,大约需要0.003秒。(每一个都可以与几十个网格单元相交)需要大约15分钟以上(长达30+分钟,取决于交叉网格单元的数量)在我的机器上,这是不能接受的。不幸的是,我不知道如何写一个多边形相交的代码来得到重叠的面积。你有什么提示吗?有没有替代shapely的方法?

xdyibdwo

xdyibdwo1#

考虑使用Rtree来帮助识别多边形可能与哪些网格单元相交。这样,您可以删除与lat/隆恩数组一起使用的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)
uujelgoq

uujelgoq2#

从2013/2014年开始Shapely就有STRtree了,我用过,看起来效果不错。
下面是文档字符串的一个片段:
STRtree是一个R树,它是使用Sort-Tile-Recursive算法创建的。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
woobm2wo

woobm2wo3#

下面是另一个版本的答案,但对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

相关问题