Python中基于重叠阈值构建Shapely多边形关联索引字典及合并

一、构建重叠索引字典 #

直接对12000个多边形两两计算重叠率会产生O(n²)的时间复杂度，完全无法实用。这里用Shapely的STRtree空间索引快速筛选潜在重叠的多边形，大幅减少计算量。

实现步骤 #

1. 自定义重叠阈值（示例设为0.5，可按需调整）
2. 用STRtree创建空间索引，快速获取每个多边形的候选重叠对象
3. 对候选对象计算重叠率（交集面积÷较小多边形面积），筛选符合阈值的索引并构建字典

代码示例 #

from shapely.strtree import STRtree
from shapely.geometry import Polygon

# 假设你的多边形列表为polygons，包含12000个Shapely Polygon实例
polygons = [Polygon(...), ...]
threshold = 0.5  # 自定义重叠阈值

# 优化：绑定索引与多边形，避免后续O(n)查找
indexed_polys = list(enumerate(polygons))
# 创建空间索引
tree = STRtree([poly for _, poly in indexed_polys])
# 初始化结果字典
overlap_dict = {idx: [] for idx in range(len(polygons))}

for idx, poly in indexed_polys:
    # 获取所有潜在重叠的候选多边形
    candidate_polys = tree.query(poly)
    for cand_poly in candidate_polys:
        # 找到候选多边形对应的原始索引
        cand_idx = next(i for i, p in indexed_polys if p is cand_poly)
        # 跳过自身
        if cand_idx == idx:
            continue
        # 计算重叠率
        intersection_area = poly.intersection(cand_poly).area
        min_area = min(poly.area, cand_poly.area)
        # 跳过面积为0的无效多边形
        if min_area == 0:
            continue
        overlap_ratio = intersection_area / min_area
        # 符合阈值则双向添加（若仅需单向关联，可移除其中一行）
        if overlap_ratio >= threshold:
            if cand_idx not in overlap_dict[idx]:
                overlap_dict[idx].append(cand_idx)
            if idx not in overlap_dict[cand_idx]:
                overlap_dict[cand_idx].append(idx)

二、获取合并后的多边形列表 #

合并重叠多边形有两种常用方案，可根据场景选择：

方案1：直接全局合并 #

利用unary_union直接合并所有多边形，再拆分结果为单个多边形列表。优点是代码简洁，适合大多数场景。

from shapely.ops import unary_union

# 全局合并所有多边形
merged_union = unary_union(polygons)
# 拆分结果为多边形列表
merged_polygons = []
if merged_union.geom_type == 'MultiPolygon':
    merged_polygons = list(merged_union.geoms)
elif merged_union.geom_type == 'Polygon':
    merged_polygons = [merged_union]

方案2：按重叠分组合并 #

基于之前的overlap_dict，用并查集（Union-Find）找出所有互相重叠的多边形组，再对每组单独合并。适合多边形数量极大、重叠区域集中的场景，能减少内存占用。

from shapely.ops import unary_union

# 实现并查集类
class UnionFind:
    def __init__(self, size):
        self.parent = list(range(size))
    
    def find(self, x):
        if self.parent[x] != x:
            self.parent[x] = self.find(self.parent[x])
        return self.parent[x]
    
    def union(self, x, y):
        xr = self.find(x)
        yr = self.find(y)
        if xr != yr:
            self.parent[yr] = xr

# 构建连通分量（互相重叠的多边形组）
uf = UnionFind(len(polygons))
for idx, neighbors in overlap_dict.items():
    for neighbor in neighbors:
        uf.union(idx, neighbor)

# 按连通分量分组
poly_groups = {}
for idx in range(len(polygons)):
    root = uf.find(idx)
    if root not in poly_groups:
        poly_groups[root] = []
    poly_groups[root].append(polygons[idx])

# 每组合并并收集结果
merged_polygons = []
for group in poly_groups.values():
    group_union = unary_union(group)
    if group_union.geom_type == 'MultiPolygon':
        merged_polygons.extend(list(group_union.geoms))
    elif group_union.geom_type == 'Polygon':
        merged_polygons.append(group_union)

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