从多边形列表中减去内环答案

【问题标题】：Substracting inner rings from a list of polygons从多边形列表中减去内环
【发布时间】：2020-03-25 03:59:42
【问题描述】：

我有很多多边形 (>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，但在哪里消除了重叠）？替代如何有效地确定哪些是孔，然后从较大的多边形中减去它们？

【问题讨论】：

这些真的是内圈还是孔？没有孔的内圈没有意义。把一个洞想象成你地表内的一个湖，把内环想象成这个湖上的岛屿。您不能从陆地中减去岛屿，但可以减去湖泊。除非您随后添加一个岛，否则它仍然是一个多边形，这将使其成为多多边形。
正如@MichaelEntin 所说，您似乎将多面体（一组多边形，如离岸岛屿的区域）和带孔的多边形（一个有湖的公园）混为一谈。如何处理它们取决于您的意思
它们是漏洞。我已经更正了我的问题。

标签： python gis polygon shapely

【解决方案1】：

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

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))

【讨论】：