如何获取使用shapely.simplify简化后多边形的原始顶点索引？

先看示例代码：

from shapely import simplify, points, contains, Point

circle = Point(0, 0).buffer(1.0, quad_segs=8).exterior
simple = simplify(circle, 0.1)

将原始多边形circle（红色）通过simplify方法简化后得到多边形simple（蓝色），简化后的多边形顶点是原始多边形顶点的子集：

我们需要得到简化后顶点对应的原始顶点索引列表，比如示例中的[0, 4, 8, 12, 16, 20, 24, 28, 32]。当前的实现方式是遍历原始顶点，逐个判断是否被简化后的多边形包含：

iCircle = []

for i, p in enumerate(points(circle.coords)):
  if contains(simple, p):
    iCircle.append(i) 

这种方式效率极低，尤其是面对顶点数量多的多边形时。如何避免这种低效查找，高效计算该索引列表？

注：圆形仅为示例，该问题适用于任意多边形。

高效解决方案 #

方法1：坐标精确匹配（无浮点误差场景） #

如果简化后的顶点是原始顶点的精确副本，可直接通过坐标映射快速查找索引：

# 提取原始与简化后的坐标列表
original_coords = list(circle.coords)
simplified_coords = list(simple.coords)

# 构建坐标到原始索引的映射（若有重复坐标需额外处理，此处假设原始坐标唯一）
coord_to_idx = {tuple(coord): idx for idx, coord in enumerate(original_coords)}

# 批量获取索引
iCircle = [coord_to_idx[tuple(coord)] for coord in simplified_coords]

该方法时间复杂度为O(n)+O(m)，远优于原方法的O(n*m)。若存在浮点精度问题，可先对坐标保留固定小数位后再构建映射。

方法2：用snap对齐坐标（处理浮点误差） #

若浮点运算导致坐标存在微小差异，先用snap将简化后的顶点对齐到原始顶点，再进行匹配：

from shapely.ops import snap

# 将简化后的多边形顶点对齐到原始多边形顶点，tolerance根据误差调整
snapped_simple = snap(simple, circle, tolerance=1e-9)
simplified_coords = list(snapped_simple.coords)

original_coords = list(circle.coords)
coord_to_idx = {tuple(coord): idx for idx, coord in enumerate(original_coords)}

iCircle = [coord_to_idx[tuple(coord)] for coord in simplified_coords]

方法3：自定义简化算法跟踪索引 #

若要完全规避坐标匹配问题，可自行实现Douglas-Peucker算法（Shapely默认简化算法），在简化过程中直接记录保留的顶点索引：

def douglas_peucker_with_indices(points, epsilon):
    def perpendicular_distance(point, line_start, line_end):
        if line_start == line_end:
            return ((point[0]-line_start[0])**2 + (point[1]-line_start[1])**2)**0.5
        numerator = abs((line_end[0]-line_start[0])*(line_start[1]-point[1]) - (line_start[0]-point[0])*(line_end[1]-line_start[1]))
        denominator = ((line_end[0]-line_start[0])**2 + (line_end[1]-line_start[1])**2)**0.5
        return numerator / denominator

    max_dist = 0.0
    index = 0
    for i in range(1, len(points)-1):
        dist = perpendicular_distance(points[i], points[0], points[-1])
        if dist > max_dist:
            max_dist = dist
            index = i

    indices = []
    if max_dist > epsilon:
        left_indices = douglas_peucker_with_indices(points[:index+1], epsilon)
        right_indices = douglas_peucker_with_indices(points[index:], epsilon)
        indices = left_indices[:-1] + right_indices
    else:
        indices = [0, len(points)-1]

    return indices

# 应用自定义算法获取索引
original_coords = list(circle.coords)
iCircle = douglas_peucker_with_indices(original_coords, 0.1)

该方法直接在简化流程中记录索引，无需后续匹配，精度与效率都有保障。

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