如何用Python库（如Shapely）计算多个相交多边形的平均值？

针对多个存在大量交集的多边形，以下几种基于Python库的方案可以实现符合直观预期的“平均”效果：

1. 顶点坐标平均法（适合形状相似的多边形） #

直接对所有多边形的对应顶点坐标取均值，构建新多边形。适合顶点数量、顺序相近的多边形集合。

from shapely.geometry import Polygon
import numpy as np

def average_polygons_by_vertices(polygons):
    # 提取每个多边形的外轮廓顶点（忽略最后一个重复的闭合点）
    vertices_list = [np.array(poly.exterior.coords)[:-1] for poly in polygons]
    
    # 若顶点数量不一致，通过角度排序对齐（以质心为原点）
    def align_vertices(vertices):
        centroid = np.mean(vertices, axis=0)
        # 计算每个顶点相对于质心的角度
        angles = np.arctan2(vertices[:,1]-centroid[1], vertices[:,0]-centroid[0])
        # 按角度排序顶点
        return vertices[np.argsort(angles)]
    
    aligned_vertices = [align_vertices(v) for v in vertices_list]
    # 统一顶点数量（通过线性插值补点，示例补到最大数量）
    max_len = max(len(v) for v in aligned_vertices)
    resampled_vertices = []
    for v in aligned_vertices:
        if len(v) == max_len:
            resampled_vertices.append(v)
        else:
            # 线性插值补点
            indices = np.linspace(0, len(v)-1, max_len)
            resampled = np.array([np.interp(indices, np.arange(len(v)), v[:,0]),
                                 np.interp(indices, np.arange(len(v)), v[:,1])]).T
            resampled_vertices.append(resampled)
    
    # 计算顶点均值
    avg_vertices = np.mean(resampled_vertices, axis=0)
    return Polygon(avg_vertices)

2. 缓冲区交集迭代法（贴合重叠区域的直观平均） #

从多边形交集出发，逐步扩展缓冲区，保留与多数多边形重叠的区域，最终得到介于交集和并集之间的“平均”形状。

from shapely.ops import unary_union, intersection

def average_polygons_by_buffer(polygons, step=0.1, overlap_threshold=0.7):
    total_polys = len(polygons)
    # 初始形状为所有多边形的交集
    current_shape = intersection(*polygons)
    
    # 迭代扩展缓冲区，直到重叠比例低于阈值
    while True:
        buffered_shape = current_shape.buffer(step)
        # 统计与buffered_shape相交的多边形数量
        overlap_count = sum(1 for poly in polygons if buffered_shape.intersects(poly))
        if overlap_count / total_polys < overlap_threshold:
            break
        current_shape = buffered_shape
    
    # 最后用所有多边形的并集裁剪，避免过度扩展
    final_shape = intersection(current_shape, unary_union(polygons))
    return final_shape

• 参数说明：step控制每次扩展的步长，overlap_threshold控制保留的重叠比例（如0.7表示保留至少70%多边形覆盖的区域）。

3. 栅格化矢量化法（鲁棒性强，适合复杂形状） #

将多边形转换为栅格，计算每个栅格点的覆盖比例，提取符合比例阈值的区域后再矢量化为多边形，适合处理顶点不规则的复杂多边形。

import rasterio
from rasterio.features import rasterize, shapes
import numpy as np
from shapely.geometry import shape
from shapely.ops import unary_union

def average_polygons_by_raster(polygons, resolution=100, coverage_threshold=0.5):
    # 获取所有多边形的边界范围
    union_shape = unary_union(polygons)
    bounds = union_shape.bounds
    # 计算栅格尺寸
    width = int((bounds[2] - bounds[0]) * resolution)
    height = int((bounds[3] - bounds[1]) * resolution)
    transform = rasterio.transform.from_bounds(*bounds, width, height)
    
    # 初始化栅格并累加每个多边形的覆盖
    raster = np.zeros((height, width), dtype=np.int32)
    for poly in polygons:
        raster += rasterize([poly], out_shape=(height, width), transform=transform)
    
    # 提取覆盖比例达标区域
    coverage_mask = (raster / len(polygons)) >= coverage_threshold
    # 矢量化回多边形
    shape_gen = shapes(coverage_mask.astype(np.uint8), transform=transform)
    avg_polygons = [shape(geom) for geom, val in shape_gen if val == 1]
    
    # 合并多个多边形（如果有）
    return unary_union(avg_polygons) if avg_polygons else None

• 参数说明：resolution控制栅格精度（数值越高精度越高），coverage_threshold控制覆盖比例（如0.5表示保留至少被一半多边形覆盖的区域）。

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