与多个 MultiPolygons 的形状相交，其中一定比例的 Multipolygons 相交？

Question

我将 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...

我可以在之前的图上画出这个交点，我们看到：

这看起来不错。但是，此方法仅 returns 所有 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])

在这四个多边形上，交集方法将returns没有交集，因为没有一个点被所有四个多边形占据。然而，用黄色标记突出的蓝色正方形被蓝色、橙色和绿色的多边形占据。因此，75% 的多面体在此位置重叠。

有没有办法（最好使用 Shapely）获取 75% 的多边形重叠的位置？

接受的答案似乎在某些情况下崩溃了。我发现了问题：包含形状的边界有时会与内部多边形重叠。我可以用 if geom.contains(polygon.buffer(-1)): 替换 if geom.contains(polygon): 并且它完全按预期工作。

Answer 1

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

首先，我们需要解决一个问题。您放在一起的示例有几个无效的几何图形，这将导致 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]'

基本上，当它们应该是 multi-polygons:

时，像这样的形状被表示为多边形并不令人高兴

因此，为了使您的示例有效，您的一些 multi-polygons 将具有 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 实现，但目前还没有已发布），否则下面描述的方法将不起作用。

除此之外，第一步是从形状列表中获取 non-intersecting 区域的平面列表。为此，我们从原始的相交形状列表开始（注意多个形状重叠处的深色阴影）：

使用 linemerge (in combination with unaryunion):

将它们转换为行

然后polygonize结果：

从图片上看可能不太清楚，但想法是 none 这些几何形状相交 - 这些多边形中的一些有孔（在一个形状以前包含另一个形状的情况下）。所以这代表了我一开始提到的“XY-plane 上 non-intersecting 区域的平面列表”。

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

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)

与多个 MultiPolygons 的形状相交，其中一定比例的 Multipolygons 相交？

Shapely intersection with multiple MultiPolygons where a percentage of Multipoloygons intersect?

python

polygon

shapely