在Shapely中，拆分自相交的多边形只返回一个多边形

Question

我在Windows 7 64位上使用Python 3.5 64位，shapely版本为1.5.13。

以下是我的代码，返回了一个自相交的多边形：

import numpy as np
from shapely.geometry import Polygon, MultiPolygon
import matplotlib.pyplot as plt

x = np.array([ 0.38517325,  0.40859912,  0.43296919,  0.4583215 ,  0.4583215 ,
               0.43296919,  0.40859912,  0.38517325,  0.36265506,  0.34100929])
y = np.array([ 62.5       ,  56.17977528,  39.39698492,   0.        ,
               0.        ,  17.34605377,  39.13341671,  60.4180932 ,
               76.02574417,  85.47008547])
polygon = Polygon(np.c_[x, y])
plt.plot(*polygon.exterior.xy)

这是正确的。然后我尝试使用buffer(0)方法来获取两个单独的多边形：

split_polygon = polygon.buffer(0)
plt.plot(*polygon.exterior.xy)
print(type(split_polygon))
plt.fill(*split_polygon.exterior.xy)

很遗憾，它只返回了两个多边形中的一个：

请问有人可以帮忙吗？谢谢！

- Yuxiang Wang

Answer 1

第一步是关闭LineString以创建LinearRing，这是构成多边形的基础。

from shapely.geometry import LineString, MultiPolygon
from shapely.ops import polygonize, unary_union

# original data
ls = LineString(np.c_[x, y])
# closed, non-simple
lr = LineString(ls.coords[:] + ls.coords[0:1])
lr.is_simple  # False

然而，请注意它不是简单的，因为线条交叉形成了一个蝴蝶结。 （根据我的经验，通常使用的缓冲区（0）技巧无法修复蝴蝶结）。这对于LinearRing不合适，因此需要进一步处理。使用unary_union使其变得简单并成为MultiLineString：

mls = unary_union(lr)
mls.geom_type  # MultiLineString'

然后使用polygonize从线条中找到多边形：

for polygon in polygonize(mls):
    print(polygon)

或者，如果你想要一个MultiPolygon几何图形：

mp = MultiPolygon(list(polygonize(mls)))

- Mike T

---

谢谢，这似乎在大多数情况下都有效。实际上，我必须在之后进行缓冲（0），才能完全清理我的要素，但至少它没有丢失大块。小细节：在制作多边形时不需要调用“list()”，MultiPolygon对象将接受生成器。 - berardig

@Mike T：非常感谢！我正在大学的论文中使用它，我很高兴我找到了它。我还需要做一些关于实现方面的笔记，但我仍然没有完全理解你的代码，你介意再解释一次吗？ - Mauritius

使用 shapely 1.8.0，您可以使用 shapely.validation.make_valid(shapely.geometry.Polygon(np.c_[x, y])) 一行代码来获取相同的多边形。 - RunOrVeith

Answer 2

我在2020年仍然为此苦苦挣扎，最终编写了一种清理自相交的方法。

这需要使用Shapely v1.2.1的explain_validity()方法才能实现。

def clean_bowtie_geom(base_linearring):
    base_polygon = Polygon(base_linearring)

    invalidity = explain_validity(base_polygon)
    invalid_regex = re.compile('^(Self-intersection)[[](.+)\s(.+)[]]$')
    match = invalid_regex.match(invalidity)
    if match:
        groups = match.groups()
        intersect_point = (float(groups[1]), float(groups[2]))
        new_linring_coords1 = []
        new_linring_coords2 = []
        pop_new_linring = False

        for i in range(0, len(base_linearring.coords)):
            if i == len(base_linearring.coords) - 1:
                end_point = base_linearring.coords[0]
            else:
                end_point = base_linearring.coords[i + 1]
            start_point = base_linearring.coords[i]

            if not pop_new_linring:
                if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
                    new_linring_coords2.append(intersect_point)
                    new_linring_coords1.append(intersect_point)
                    pop_new_linring = True
                else:
                    new_linring_coords1.append(start_point)

            else:
                new_linring_coords2.append(start_point)
                if is_point_on_line_and_between(start=start_point, end=end_point, pt=intersect_point):
                    new_linring_coords2.append(intersect_point)
                    pop_new_linring = False

        corrected_linear_ring1 = LinearRing(coordinates=new_linring_coords1)
        corrected_linear_ring2 = LinearRing(coordinates=new_linring_coords2)

        polygon1 = Polygon(corrected_linear_ring1)
        polygon2 = Polygon(corrected_linear_ring2)
        
def is_point_on_line_and_between(start, end, pt, tol=0.0005):
    """
    Checks to see if pt is directly in line and between start and end coords
    :param start: list or tuple of x, y coordinates of start point of line
    :param end: list or tuple of x, y coordinates of end point of line
    :param pt: list or tuple of x, y coordinates of point to check if it is on the line
    :param tol: Tolerance for checking if point on line
    :return: True if on the line, False if not on the line
    """
    v1 = (end[0] - start[0], end[1] - start[1])
    v2 = (pt[0] - start[0], pt[1] - start[1])
    cross = cross_product(v1, v2)
    if cross <= tol:
        # The point lays on the line, but need to check if in between
        if ((start[0] <= pt[0] <= end[0]) or (start[0] >= pt[0] >= end[0])) and ((start[1] <= pt[1] <= end[1]) or (start[1] >= pt[1] >= end[1])):
            return True
    return False

这段代码不是最干净的，但对我来说完成了任务。

输入是一个具有自相交几何形状（is_simple=False）的LinearRing，输出可以是2个LinearRings或两个多边形，任选其一（或者有条件选择其中一个），真的很随意。

---

编辑

---

在Shapely 1.8.0中，添加了新功能。 shapely.validation.make_valid()将获取一个自相交的多边形并返回由拆分自相交点创建的每个多边形的MultiPolygon。

- Chris Reed

---

我想您是指点积，而不是叉积？ - Shariq