使用Shapely按线分割多边形却返回原多边形的问题

使用Shapely 2.0版本，尝试用直线L将多边形P分割为两部分：

1. 先获取P与L的交集，得到LineString类型的L2；
2. 调用split函数传入P和L2，预期返回两个子多边形，但实际仅返回原多边形P。

处理数百个同类多边形时，仅约10%的案例分割成功，其余均出现上述问题。怀疑是浮点精度误差导致Shapely判定L2的边界未与P的边界精确相交。

相关数据 #

多边形P的WKT：

POLYGON ((292718.0381447676 6638193.414029885, 292694.50537013356 6637994.803334004, 292718.0381447676 6638193.414029885, 292718.9708331647 6638193.303518486, 292722.0992155936 6638192.97038651, 292725.48975053953 6638192.648911643, 292729.16966360115 6638192.34135049, 292733.16283984896 6638192.050457474, 292737.48817684915 6638191.779469689, 292742.1600771896 6638191.532042792, 292747.1894678763 6638191.312155023, 292752.58358723175 6638191.124185786, 292758.34681013395 6638190.973003162, 292753.9574101781 6637991.021175833, 292758.34681013395 6638190.973003162, 292753.9574101781 6637991.021175833, 292694.50537013356 6637994.803334004, 292718.0381447676 6638193.414029885))

交集L2的WKT：

LINESTRING (292756.2221414001 6638094.18724655, 292743.0611 6638095.284, 292709.9881 6638095.284, 292706.4410368586 6638095.537357549)

代码示例 #

from shapely.ops import split
from shapely import wkt

box = wkt.loads(my_polygon_as_str)
line = wkt.loads(my_line_as_str)
result = split(box, line)  # 预期返回2个多边形，实际返回原多边形

针对浮点精度导致的分割失败问题，可采用以下几种方法处理：

1. 对分割线应用极小缓冲 #

通过给分割线L2添加极小的缓冲（比如1e-8量级，根据数据精度调整），让它确保与多边形边界产生有效交集，从而触发分割。

from shapely.ops import split
from shapely import wkt

box = wkt.loads(my_polygon_as_str)
line = wkt.loads(my_line_as_str)

# 添加极小缓冲，将线转化为极窄的多边形
buffered_line = line.buffer(1e-8)
result = split(box, buffered_line)

2. 使用snap函数对齐端点 #

利用Shapely的snap函数，将分割线的端点对齐到多边形的边界上，消除浮点误差导致的端点偏移。

from shapely.ops import split, snap
from shapely import wkt

box = wkt.loads(my_polygon_as_str)
line = wkt.loads(my_line_as_str)

# 将分割线端点对齐到多边形边界，精度阈值设为1e-8
snapped_line = snap(line, box.boundary, 1e-8)
result = split(box, snapped_line)

3. 验证并修正分割线的端点 #

检查分割线的端点是否真的在多边形边界上，若不在则替换为边界上的最近点：

from shapely.ops import split, nearest_points
from shapely import wkt

box = wkt.loads(my_polygon_as_str)
line = wkt.loads(my_line_as_str)
boundary = box.boundary

# 获取分割线的端点
coords = list(line.coords)
new_coords = []
for coord in coords:
    # 找到边界上离当前端点最近的点
    nearest_pt = nearest_points(coord, boundary)[1]
    new_coords.append((nearest_pt.x, nearest_pt.y))

# 重新创建修正后的分割线
corrected_line = wkt.loads(f"LINESTRING ({' '.join([f'{x} {y}' for x,y in new_coords])})")
result = split(box, corrected_line)

4. 调整Shapely的精度阈值 #

Shapely底层依赖GEOS库，可通过设置环境变量调整精度：

import os
# 设置GEOS的精度阈值，单位与数据坐标一致
os.environ['GEOS_PRECISION'] = '1e-8'

from shapely.ops import split
from shapely import wkt

box = wkt.loads(my_polygon_as_str)
line = wkt.loads(my_line_as_str)
result = split(box, line)

内容的提问来源于stack exchange，提问作者Beinje