最快实现批量点缓冲区几何合并(union)的方法有哪些？

你当前使用的三种方案性能基线如下：

• cascaded_union是shapely已废弃的旧接口，性能落后于unary_union，不建议继续使用
• 原生unary_union是shapely官方推荐的矢量合并接口，geopandas的GeoSeries.unary_union底层直接调用该接口，二者性能几乎一致

以下是速度更快的实现方案，按优化效果排序： #

1. 替换循环创建缓冲逻辑，使用向量化操作 #

你当前用列表循环创建Point再调用buffer的逻辑是最大的性能瓶颈，点数量超过1万时这部分耗时会占总耗时的60%以上。用shapely的向量化接口可以直接批量处理坐标：

from shapely import multipoints, buffer, unary_union

# 直接基于坐标数组批量生成缓冲，完全避免循环
points_buffered = buffer(multipoints(coords), 0.01)
mpoly = unary_union(points_buffered)

如果是shapely 2.0+版本，这个优化可以把缓冲生成的速度提升10倍以上。 #

2. 空间分块合并 #

如果你的点量在10万级以上，且点分布比较分散、很多缓冲不存在重叠，可以先按空间网格对坐标分组，每个分组内单独做合并，最后再合并所有分组的结果：

import numpy as np
from shapely import unary_union, buffer, multipoints

# 按坐标整数倍分块，网格大小可以根据你的缓冲半径调整
grid_size = 0.5
grid_keys = np.floor(coords / grid_size).astype(int)
unique_keys = np.unique(grid_keys, axis=0)

part_unions = []
for key in unique_keys:
    # 筛选当前网格内的坐标
    mask = (grid_keys == key).all(axis=1)
    part_coords = coords[mask]
    part_unions.append(unary_union(buffer(multipoints(part_coords), 0.01)))

# 最后合并所有分块的结果
mpoly = unary_union(part_unions)

这个方案可以减少单次合并的几何数量，性能提升幅度在2~5倍左右，具体取决于点的分布密集程度。 #

3. 栅格化近似方案 #

如果对合并结果的精度要求不高，可以用栅格化转轮廓的方案，点量越大性能优势越明显，100万级点的场景下速度比纯矢量合并快10倍以上：

import numpy as np
from shapely.geometry import Polygon
from skimage.measure import find_contours

# 计算坐标范围
x_min, y_min = coords.min(axis=0)
x_max, y_max = coords.max(axis=0)
# 栅格分辨率，根据精度要求调整
res = 0.001

# 生成栅格掩码
x_bins = np.arange(x_min - 0.01, x_max + 0.01, res)
y_bins = np.arange(y_min - 0.01, y_max + 0.01, res)
mask = np.zeros((len(y_bins), len(x_bins)), dtype=bool)

# 批量计算每个点对应的栅格位置，标记缓冲区域
radius = int(0.01 / res)
x_idx = np.clip(np.digitize(coords[:, 0], x_bins), radius, len(x_bins) - radius - 1)
y_idx = np.clip(np.digitize(coords[:, 1], y_bins), radius, len(y_bins) - radius - 1)

for x, y in zip(x_idx, y_idx):
    mask[y-radius:y+radius+1, x-radius:x+radius+1] = True

# 提取轮廓转多边形
contours = find_contours(mask, 0.5)
polygons = []
for cnt in contours:
    cnt_coords = np.column_stack((x_bins[cnt[:, 1].astype(int)], y_bins[cnt[:, 0].astype(int)]))
    if len(cnt_coords) >= 3:
        polygons.append(Polygon(cnt_coords))

mpoly = unary_union(polygons)

小提示：如果还在使用shapely 1.x版本，升级到shapely 2.0+可以直接获得30%以上的性能提升，不需要修改业务代码。

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