从多边形列表中减去内部环路

Question

我有一个包含大量多边形的列表(>10^6)，其中大多数是非交叉的，但其中一些多边形是另一个多边形的孔洞(约10^3个案例)。下面是一张图片来解释这个问题，小多边形是大多边形中的一个孔洞，但两者都是多边形列表中的独立多边形。 现在，我希望能够高效地确定哪些多边形是孔洞并减去这些孔洞，即减去完全位于另一个多边形内部的小多边形，并返回一个“清理过”的多边形列表。 孔洞和父多边形对应关系应该转换成如下形式(也就是将孔洞从父多边形中减去)： 虽然stackoverflow和gis.stackexchange.com上有很多类似的问题，但我没有找到一个实际解决这个问题的问题。以下是一些相关的问题: 1. https://gis.stackexchange.com/questions/5405/using-shapely-translating-between-polygons-and-multipolygons 2. https://gis.stackexchange.com/questions/319546/converting-list-of-polygons-to-multipolygon-using-shapely 下面是一段示例代码。

from shapely.geometry import Point
from shapely.geometry import MultiPolygon
from shapely.ops import unary_union
import numpy as np

#Generate a list of polygons, where some are holes in others; 
def generateRandomPolygons(polygonCount = 100, areaDimension = 1000, holeProbability = 0.5):
    pl = []
    radiusLarge = 2 #In the real dataset the size of polygons can vary
    radiusSmall = 1 #Size of holes can also vary

    for i in range(polygonCount):
        x, y = np.random.randint(0,areaDimension,(2))
        rn1 = np.random.random(1)
        pl.append(Point(x, y).buffer(radiusLarge))
        if rn1 < holeProbability: #With a holeProbability add a hole in the large polygon that was just added to the list
            pl.append(Point(x, y).buffer(radiusSmall))
    return pl

polygons = generateRandomPolygons()
print(len(pl))

输出结果如下：

现在我该如何创建一个新的多边形列表，将孔洞去除。Shapely提供了从一个多边形中减去另一个多边形（difference）的函数，但是是否有类似于多个多边形列表的函数（可能类似于unary_union），可以去除重叠的部分？另外，如何高效地确定哪些是孔洞，并从大的多边形中删除它们？

- KarateKid

---

这些是内环还是孔洞？没有孔洞，内环就没有意义。把孔洞想象成你土地表面上的湖泊，内环就像是湖泊上的岛屿。你不能从土地中减去岛屿，但可以减去湖泊。它仍然是一个单一的多边形，除非你再添加一个岛屿，这将使它成为多重多边形。 - Michael Entin

正如@MichaelEntin所说，您似乎混淆了多边形（例如带有离岸岛屿的区域）和具有洞的多边形（例如公园中的湖泊）。如何处理它们取决于您的意思是哪一个。 - Ian Turton

它们是漏洞。我已经修正了我的问题。 - KarateKid

Answer 1

你的问题是你不知道哪些是“空洞”，对吧？要“有效地确定哪些多边形是空洞”，你可以使用r树来加速相交检查：

from rtree.index import Index

# create an rtree for efficient spatial queries
rtree = Index((i, p.bounds, None) for i, p in enumerate(polygons))
donuts = []

for i, this_poly in enumerate(polygons):
    # loop over indices of approximately intersecting polygons
    for j in rtree.intersection(this_poly.bounds):
        # ignore the intersection of this polygon with itself
        if i == j:
            continue
        other_poly = polygons[j]
        # ensure the polygon fully contains our match
        if this_poly.contains(other_poly):
            donut = this_poly.difference(other_poly)
            donuts.append(donut)
            break  # quit searching

print(len(donuts))

- berardig