与多个MultiPolygons的形状相交，其中有一定百分比的Multipoloygons相交？

问题描述

我将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,1),(0,1)])
P2 = polygon([(2.5,2),(3,3),(2.5,3)])
multi_a = Multipolygon([P1,P2])

P1 = polygon([(-1,-1),(-1,(2,-1)])
P2 = polygon([(3,(4,4),4)])    
multi_b = Multipolygon([P1,P2])  

P1 = polygon([(-2,(2.2,-2),(-2,-2)])
P2 = polygon([(-1.5,(-1.1,(-1.5,-1)])    
multi_c = Multipolygon([P1,P2])

P1 = polygon([(2.5,(3.2,-2)])
P2 = polygon([(3,1)])
multi_d = Multipolygon([P1,P2])

在这四个多边形上，交集方法将不返回交集，因为所有四个多边形都没有一个点。但是，带有黄色标记的蓝色正方形被蓝色，橙色和绿色多边形占据。因此，有75％的多边形在此位置重叠。

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

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

解决方法

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

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

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

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

>>> P2.is_valid
False

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

基本上，对于像这样的形状表示为多边形的情况感到不满意：

为使您的示例有效，您的某些多面多边形将包含3个而非2个组成多边形：

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

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

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

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

希望您的真实源数据具有有效的几何形状（或者您可以通过某种方式将其转换为有效的几何形状-顺便说一句，这是Shapely 1.8中即将推出的功能，通过make_valid实现，但尚未发布），因为否则将无法使用下面描述的方法。

首先，第一步是从形状列表中获得不相交区域的平面列表。为此，我们从相交形状的原始列表开始（请注意多个形状重叠的深色阴影）：

使用linemerge（结合unaryunion）将它们转换为行：

然后polygonize结果：

从图片上可能不清楚，但想法是这些几何形状都不相交-这些多边形中的一些在其中具有孔（在一种形状先前包含另一种形状的情况下）。因此，这代表了我一开始所指的“ 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)

polygon python shapely