与多个多重多边形的形状相交，其中一定百分比的多重多边形相交？

Question

我正在使用 Shapely 的多边形来获取人类生成的数据。要求多人围绕图像中的某些特征绘制多边形。因此，对于每个图像，我们有 nMultiPolygon个，其中 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)):并且它完全按照预期工作。

Answer 1

一种方法是分割所有几何图形，以获得 XY 平面上不相交区域的平面列表，然后查看有多少原始几何图形包含每个区域。至少包含某个阈值数量的原始几何图形的任何区域都可以添加到结果中。使用代码和插图的组合更容易解释这一点。

首先，我们需要解决一个问题。contains您放在一起的示例有几个无效的几何图形，这将导致 Shapely 在尝试查询空间关系（例如，调用或）时抛出错误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 平面上不相交区域的平面列表”。

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

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

接下来，我们检查列表中的每个结果区域polygons，并检查它与原始列表中的多少个形状相交。如果它高于阈值（此处定义为0.75 * len(shapes)，则我们将其添加到结果中：

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

如果您正在处理大型数据集，检查像这样的嵌套循环中的交集可能会非常慢（O(N^2)），因此您可以使用以下方法加快速度STRtree：

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)