如何获取GeoDataFrame的Voronoi面积？关联初始点与多边形方法

问题说明 #

用Shapely生成Voronoi图时，输出的多边形顺序和初始点顺序不对应，直接按顺序给GeoDataFrame赋值面积会出错。现有代码直接将voronoi_diagram输出的多边形面积按列表顺序赋值，但数据量达数百万级，逐点判断"点是否在多边形内"的方式效率极低，需要更高效的方案，同时减少不必要的格式转换。

当前问题代码：

# gdf = is a GeoDataFrame
minx, miny, maxx, maxy = gdf.total_bounds
bound = Polygon([(minx, miny),
                    (maxx, miny),
                    (maxx, maxy),
                    (minx, maxy)])
points = MultiPoint(gdf.geometry.to_list())
parcels = voronoi_diagram(points , envelope=bound)
areas = [p.area for p in parcels]
gdf['area'] = areas

解决方案 #

方案一：基于Scipy Voronoi的索引映射（百万级数据首选） #

Shapely的voronoi_diagram底层依赖Scipy的Voronoi实现，可直接通过索引映射关联初始点与对应多边形，时间复杂度O(n)，无需逐点判断：

from scipy.spatial import Voronoi
import numpy as np
from shapely.geometry import Polygon

# 提取与gdf行顺序一致的坐标数组
coords = np.array([(p.x, p.y) for p in gdf.geometry])
vor = Voronoi(coords)

# 构建边界多边形
minx, miny, maxx, maxy = gdf.total_bounds
bound_poly = Polygon([(minx, miny), (maxx, miny), (maxx, maxy), (minx, maxy)])

areas = []
for i in range(len(coords)):
    # 获取当前点对应的Voronoi区域索引
    region_idx = vor.point_region[i]
    region = vor.regions[region_idx]
    
    # 处理无效区域（边界外或无顶点的情况）
    if not region or -1 in region:
        poly_coords = [vor.vertices[j] for j in region if j != -1]
        if poly_coords:
            vor_poly = Polygon(poly_coords)
            # 裁剪到边界内计算面积
            clipped_poly = vor_poly.intersection(bound_poly)
            areas.append(clipped_poly.area)
        else:
            areas.append(0.0)
    else:
        # 构建有效Voronoi多边形并计算面积
        vor_poly = Polygon([vor.vertices[j] for j in region])
        areas.append(vor_poly.area)

# 赋值回GeoDataFrame
gdf['area'] = areas

方案二：基于Shapely空间索引的批量匹配 #

如果不想引入Scipy，可利用Shapely的STRtree构建空间索引，将查询复杂度降至O(n log n)：

from shapely.strtree import STRtree

points = MultiPoint(gdf.geometry.to_list())
parcels = list(voronoi_diagram(points, envelope=bound))

# 构建多边形空间索引
tree = STRtree(parcels)
# 批量查询每个点对应的包含它的多边形
matches = tree.query(gdf.geometry, predicate='contains')

# 整理点与多边形的对应关系（Voronoi特性保证一个点对应唯一多边形）
point_to_poly = {match[1]: match[0] for match in matches}
# 按原始点顺序提取面积
areas = [parcels[point_to_poly[i]].area for i in range(len(gdf))]

gdf['area'] = areas

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