多个MultiPolygons的Shapely相交，其中一个百分比的Multipoloygons相交？

Question

多个MultiPolygons的Shapely相交，其中一个百分比的Multipoloygons相交？

5

我正在使用Shapely的多边形来处理人类生成的数据。多个人被要求在图像中的某些特定特征周围绘制多边形。因此，对于每个图像，我们都有n个MultiPolygon，其中n等于每个图像的参与者数量。

我可以绘制这些Multipolygon中的每一个。

fig, ax = plt.subplots()

for ii, multi_poly in enumerate(multi_polys):
    for poly in multi_poly.geoms:
        x,y = poly.exterior.xy
        plt.plot(x,y, c=colors[ii])

我们可以看到，在某些位置，Multipolygon 重叠，而在其他位置则没有重叠。

我希望得到这些多边形的重叠部分或交集。这应该很简单，我可以做类似以下的操作：

intersection = multi_a.intersection(multi_b) \
            .intersection(multi_c) \
            .intersection(multi_d) \
            .inters...

我可以在之前的图表上绘制这个交点，我们可以看到：

这看起来很不错。然而，此方法仅返回所有 Multipoloygon 重叠的区域。是否有一种方法可以获得多边形重叠 75% 的交集？或者重叠 50% 的交集？

代码示例：以下虚拟数据给出了此图形：

P1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
P2 = Polygon([(2.5, 2), (3, 2), (3, 3), (2.5, 3)])
multi_a = MultiPolygon([P1, P2])

P1 = Polygon([(-1, -1), (-1, 2),(2, 2), (2, -1)])
P2 = Polygon([(3, 3), (4, 2), (4, 4), (3, 4)])    
multi_b = MultiPolygon([P1,P2])  

P1 = Polygon([(-2, 4), (2.2, 4),(2.2, -2), (-2, -2)])
P2 = Polygon([(-1.5, 3), (-1.1, 3), (-1.5, -1), (-1.1, -1)])    
multi_c = MultiPolygon([P1,P2])

P1 = Polygon([(2.5, -1), (3.2, -1),(3.2, -2), (-2, -2)])
P2 = Polygon([(3, 0), (4, 0), (3, 1), (4, 1)])
multi_d = MultiPolygon([P1,P2])

在这四个多边形中，由于没有一个点被所有四个多边形占据，因此交集方法将返回无交集。然而，用黄色标记强调的蓝色正方形被蓝色、橙色和绿色多边形占据。因此，在该位置上，75%的多边形重叠。

有没有一种方法（最好使用Shapely）来获取多边形重叠75％的位置？

接受的答案在某些情况下似乎会失效。我找到了问题所在：包含形状的边界有时会重叠在内部多边形上。我可以将if geom.contains(polygon):替换为if geom.contains(polygon.buffer(-1)):，它会按预期正常工作。

- Mitchell van Zuylen

1个回答

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

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

- Lukas S. · Accepted Answer

一种方法是将所有几何图形拆分，以获得在XY平面上的非交叉区域的平面列表，然后查看原始几何图形中包含每个区域的数量。任何被至少一定阈值数量的原始几何图形所包含的区域都可以添加到结果中。这更容易通过代码和插图的组合来解释。

首先，我们需要解决一个问题。你提供的示例有一些无效的几何图形，这将导致Shapely在尝试查询空间关系时抛出错误（例如调用contains或intersects）。您可以使用is_valid属性检查，通过调用explain_validity获取更详细的信息：

from shapely.geometry import Polygon
from shapely.validation import explain_validity

P2 = Polygon([(-1.5, 3), (-1.1, 3), (-1.5, -1), (-1.1, -1)])

>>> P2.is_valid
False

>>> explain_validity(P2)
'Self-intersection[-1.3 1]'

基本上，它不喜欢像这样的形状被表示为多边形，而应该是多多边形：

因此，为了使您的示例有效，您的一些多边形将具有3个而不是2个组成多边形：

P1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
P2 = Polygon([(2.5, 2), (3, 2), (3, 3), (2.5, 3)])
multi_a = MultiPolygon([P1, P2])

P1 = Polygon([(-1, -1), (-1, 2),(2, 2), (2, -1)])
P2 = Polygon([(3, 3), (4, 2), (4, 4), (3, 4)])    
multi_b = MultiPolygon([P1,P2])  

P1 = Polygon([(-2, 4), (2.2, 4),(2.2, -2), (-2, -2)])
P2 = Polygon([(-1.5, 3), (-1.1, 3), (-1.3, 1)])
P3 = Polygon([(-1.5, -1), (-1.3, 1), (-1.1, -1)])
multi_c = MultiPolygon([P1,P2,P3])

P1 = Polygon([(2.5, -1), (3.2, -1),(3.2, -2), (-2, -2)])
P2 = Polygon([(3, 0), (4, 0), (3.5, 0.5)])
P3 = Polygon([(3.5, 0.5), (3, 1), (4, 1)])
multi_d = MultiPolygon([P1,P2,P3])

希望您的真实源数据具有有效的几何形状（或者您有一些方法将它们转换为有效的几何形状——顺便说一下，这是 Shapely 1.8 中即将推出的一个功能，通过 make_valid 实现，但目前尚未发布），否则下面描述的方法将无法工作。

解决了这个问题之后，第一步是从您的形状列表中获取一个不相交区域的平坦列表。为此，我们从原始的相交形状列表开始（请注意多个形状重叠的较深阴影）。

使用linemerge（与unaryunion结合使用）将它们转换为线条：

然后将结果进行多边形化：

从图片上可能不太清楚，但是这个想法是这些几何图形中没有一个相交 - 其中一些多边形有洞（在一个形状先前包含另一个形状的情况下）。因此，这代表了我一开始提到的“XY平面上非相交区域的平面列表”。

到目前为止，该过程的代码如下：

from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import linemerge, unary_union, polygonize

# Original shape definitions here (snipped)...

shapes = [multi_a, multi_b, multi_c, multi_d]

lines = unary_union(linemerge([geom.exterior for shape in shapes for geom in shape.geoms]))
polygons = list(polygonize(lines))

接下来，我们检查polygons列表中的每个结果区域，并检查它与原始列表中多少形状相交。如果它超过了阈值（在这里定义为0.75 * len(shapes)），那么我们将其添加到结果中：

threshold = 0.75 * len(shapes)

def overlaps(polygon, shape):
    for geom in shape.geoms:
        if geom.contains(polygon):
            return True
    return False

result = []

for polygon in polygons:
    containing_shapes = []
    for shape in shapes:
        if overlaps(polygon, shape):
            containing_shapes.append(shape)
    if len(containing_shapes) >= threshold:
        result.append(polygon)

如果你正在处理一个大数据集，像这样在嵌套循环中检查交集可能会非常慢（O(N^2)），因此你可以使用 STRtree 来加速它：

from shapely.strtree import STRtree

# (Previous code here to get the flattened list of polygons...)

tree = STRtree([geom for shape in shapes for geom in shape.geoms])
result = []

for polygon in polygons:
    matches = [o.wkt for o in tree.query(polygon) if o.contains(polygon)]
    if len(matches) >= threshold:
        result.append(polygon)