Python计算泰森多边形?Voronoi咋画?

wowwwai GIS研习社 · 工具流程与项目排障

很多 GIS 同学搜索“Python计算泰森多边形?Voronoi咋画?”，其实是在解决同一个问题：手里有一批点数据，想用 Python 自动生成每个点的服务范围、影响范围或最近邻分区。本文用 GeoPandas、SciPy 和 Shapely 演示一套可复现的 Python 计算泰森多边形流程，并说明 Voronoi 结果为什么必须裁剪到研究区边界。

引言：Python计算泰森多边形适合哪些 GIS 场景

泰森多边形，也常叫 Voronoi 多边形，是一种根据点之间距离划分空间的方法。每个多边形内部的任意位置，都离该多边形对应的点最近。

在 GIS 工作中，Python计算泰森多边形常用于以下场景：

• 根据医院、学校、门店、基站等点位划分理论服务范围。
• 对采样点、气象站、水文站建立最近邻控制区。
• 为点数据生成空间分区，再统计每个分区内的人口、道路、土地利用等指标。
• 批量处理多组点数据，避免在桌面 GIS 中重复点击工具。

需要注意，泰森多边形不是行政区、实际服务区或交通可达范围。它只表达“欧氏距离最近”这个规则。如果你的分析依赖道路时间、通行成本或行政边界，就不能直接把 Voronoi 结果当作真实服务范围。

背景：为什么直接画 Voronoi 经常得到奇怪结果

很多人在 Python 里第一次画 Voronoi，会遇到三个常见问题：

• 边缘点的多边形无限延伸，无法直接保存为面数据。
• 坐标是经纬度，结果面积和距离判断不可靠。
• 生成的 Voronoi 超出研究区范围，看起来不像最终 GIS 成果。

这些问题不是代码写错，而是 Voronoi 算法本身和 GIS 数据环境共同造成的。SciPy 的 Voronoi 计算结果主要是几何拓扑结构，其中边界外侧区域可能是无穷区域；而 GIS 制图和空间统计通常需要有限的面，因此必须结合研究区边界进行裁剪。

原理：泰森多边形和 Voronoi 的关系

泰森多边形与 Voronoi 图本质上是同一个概念在 GIS 和计算几何中的不同叫法。给定一组点，空间会被分割成若干区域，每个区域对应一个点，区域内所有位置到该点的距离小于到其他点的距离。

在 Python GIS 工作流里，可以这样理解：

• 输入点：门店、采样站、设施点等点图层。
• Voronoi 计算：根据点坐标生成分割线和区域。
• 构造面：把计算结果转换为 Shapely polygon。
• 边界裁剪：用研究区范围把无限或过大的多边形裁成有限结果。
• 输出结果：保存为 GeoPackage、Shapefile 或 GeoJSON。

GIS 分析中最重要的一点是：Python计算泰森多边形前，应先使用适合距离计算的投影坐标系，而不是直接使用经纬度坐标。

步骤：用 Python 计算泰森多边形并裁剪到研究区

1. 准备 Python 环境

建议使用 Conda 创建环境，避免 GeoPandas、Shapely、Fiona、PyProj 等库之间版本冲突。

conda create -n gis-voronoi python=3.11
conda activate gis-voronoi
conda install -c conda-forge geopandas scipy shapely pyproj matplotlib

如果你已经在使用 Jupyter Notebook 或 VS Code，只要确认下面这些库可以正常导入即可。

import geopandas as gpd
import numpy as np
from scipy.spatial import Voronoi
from shapely.geometry import Polygon
import matplotlib.pyplot as plt

2. 准备点数据和研究区边界

假设我们有两个文件：

• points.gpkg：点图层，例如门店、站点或采样点。
• boundary.gpkg：研究区边界，例如城市范围或项目区范围。

读取数据：

points = gpd.read_file("points.gpkg")
boundary = gpd.read_file("boundary.gpkg")

检查几何类型：

print(points.geom_type.value_counts())
print(boundary.geom_type.value_counts())
print(points.crs)
print(boundary.crs)

点图层必须是真正的 Point 几何。如果你的数据是 CSV，需要先用经纬度字段构造点图层，再设置正确坐标系。

3. 投影到适合距离计算的坐标系

Voronoi 的距离计算依赖平面坐标。如果数据是 WGS84 经纬度，单位是度，不适合直接计算泰森多边形面积和距离。

如果研究区在中国某个城市，可以优先使用当地 CGCS2000 高斯投影、UTM 分带投影，或项目指定的米制投影。下面示例使用 EPSG:3857 只是为了演示，正式项目不建议无脑使用。

target_crs = "EPSG:3857"

points_proj = points.to_crs(target_crs)
boundary_proj = boundary.to_crs(target_crs)

如果你不确定用哪个坐标系，可以先问自己三个问题：

• 研究区范围是否很大？跨省、跨国家时不能随便用一个局部投影。
• 结果是否要计算面积？如果要，优先选择适合面积或本地测绘标准的投影。
• 是否需要和项目中的其他图层叠加？如果需要，应保持同一项目坐标系。

4. 提取点坐标并生成 Voronoi

把点几何转换为 NumPy 坐标数组：

coords = np.array([[geom.x, geom.y] for geom in points_proj.geometry])

vor = Voronoi(coords)

这一步完成的是 Voronoi 拓扑计算，但还没有得到可以直接写入 GIS 文件的面图层。接下来需要把无限区域处理成有限多边形。

5. 将无限 Voronoi 区域转成有限多边形

SciPy 官方示例中常用一个辅助函数把无限 Voronoi 区域延伸到指定半径内。下面函数可直接复制使用。

def voronoi_finite_polygons_2d(vor, radius=None):
    if vor.points.shape[1] != 2:
        raise ValueError("只支持二维点数据")

    new_regions = []
    new_vertices = vor.vertices.tolist()

    center = vor.points.mean(axis=0)

    if radius is None:
        radius = np.ptp(vor.points, axis=0).max() * 2

    all_ridges = {}

    for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
        all_ridges.setdefault(p1, []).append((p2, v1, v2))
        all_ridges.setdefault(p2, []).append((p1, v1, v2))

    for p1, region_index in enumerate(vor.point_region):
        vertices = vor.regions[region_index]

        if all(v >= 0 for v in vertices):
            new_regions.append(vertices)
            continue

        ridges = all_ridges[p1]
        new_region = [v for v in vertices if v >= 0]

        for p2, v1, v2 in ridges:
            if v2 < 0:
                v1, v2 = v2, v1

            if v1 >= 0:
                continue

            tangent = vor.points[p2] - vor.points[p1]
            tangent = tangent / np.linalg.norm(tangent)

            normal = np.array([-tangent[1], tangent[0]])

            midpoint = vor.points[[p1, p2]].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, normal)) * normal

            far_point = vor.vertices[v2] + direction * radius

            new_vertices.append(far_point.tolist())
            new_region.append(len(new_vertices) - 1)

        vs = np.asarray([new_vertices[v] for v in new_region])
        centroid = vs.mean(axis=0)
        angles = np.arctan2(vs[:, 1] - centroid[1], vs[:, 0] - centroid[0])
        new_region = np.array(new_region)[np.argsort(angles)]

        new_regions.append(new_region.tolist())

    return new_regions, np.asarray(new_vertices)

使用该函数生成多边形：

regions, vertices = voronoi_finite_polygons_2d(vor)

polygons = []

for region in regions:
    polygon = Polygon(vertices[region])
    polygons.append(polygon)

voronoi_gdf = gpd.GeoDataFrame(points_proj.drop(columns="geometry"), geometry=polygons, crs=points_proj.crs)

这里保留了原始点的属性字段，因此输出的每个泰森多边形可以对应回原始点。例如门店编号、站点名称、采样点 ID 都可以继续使用。

6. 按研究区边界裁剪泰森多边形

直接生成的 Voronoi 多边形通常比研究区大很多。实际 GIS 成果一般要裁剪到项目区边界。

boundary_union = boundary_proj.dissolve()

voronoi_clip = gpd.overlay(voronoi_gdf, boundary_union, how="intersection")

如果边界图层有多个面，先 dissolve 成一个整体，可以避免裁剪结果被边界内部属性拆成过多碎片。

7. 检查并修复几何

裁剪后建议检查几何有效性。无效几何可能导致后续面积统计、空间叠加或文件导出失败。

voronoi_clip["is_valid"] = voronoi_clip.geometry.is_valid
print(voronoi_clip["is_valid"].value_counts())

voronoi_clip["geometry"] = voronoi_clip.geometry.make_valid()

如果你的 Shapely 版本不支持 make_valid，可以尝试：

voronoi_clip["geometry"] = voronoi_clip.geometry.buffer(0)

不过 buffer(0) 是经验性修复方法，不保证所有复杂几何都能正确修好。正式项目中建议记录修复前后的面数量、面积变化和异常对象。

8. 计算面积并导出结果

如果坐标系单位是米，可以直接计算平方米面积，也可以换算为平方公里。

voronoi_clip["area_m2"] = voronoi_clip.geometry.area
voronoi_clip["area_km2"] = voronoi_clip["area_m2"] / 1000000

voronoi_clip.to_file("voronoi_result.gpkg", layer="voronoi", driver="GPKG")

GeoPackage 比 Shapefile 更推荐，原因是字段名限制少、中文兼容更好，也更适合保存现代 GIS 数据。

9. 快速可视化结果

可以用 GeoPandas 简单看一下结果是否合理：

ax = boundary_proj.boundary.plot(figsize=(8, 8), color="black", linewidth=1)
voronoi_clip.plot(ax=ax, column="area_km2", edgecolor="white", alpha=0.7, legend=True)
points_proj.plot(ax=ax, color="red", markersize=10)
plt.show()

检查重点不是颜色是否好看，而是每个点是否位于对应的泰森多边形内部或边界附近，以及边缘区域是否已经被研究区边界正确裁剪。

常见坑：Python画 Voronoi 时最容易出错的地方

坑 1：直接用经纬度计算泰森多边形

经纬度坐标单位是度，不是米。直接使用经纬度做 Voronoi，结果在小范围内可能看起来“差不多”，但面积、距离和形状都不严谨。Python计算泰森多边形前，应先投影到适合分析区域的平面坐标系。

坑 2：点数量太少

Voronoi 至少需要多个非共线点。如果点数量太少，或者所有点几乎在一条直线上，SciPy 可能报错。

if len(points_proj) < 3:
    raise ValueError("生成二维 Voronoi 至少需要 3 个点")

如果点共线，需要检查数据来源，或者考虑该任务是否适合用泰森多边形表达。

坑 3：重复点没有处理

重复点会让 Voronoi 计算不稳定。计算前建议删除完全重合的点，或按业务规则聚合。

points_proj["x"] = points_proj.geometry.x
points_proj["y"] = points_proj.geometry.y

points_proj = points_proj.drop_duplicates(subset=["x", "y"])

如果重复点代表多个设施位于同一位置，应先决定属性如何合并，而不是简单删除。

坑 4：裁剪边界存在无效几何

研究区边界如果有自相交、缝隙或重复面，overlay 裁剪可能失败。建议先检查并修复边界。

boundary_proj["geometry"] = boundary_proj.geometry.make_valid()
boundary_proj = boundary_proj[~boundary_proj.geometry.is_empty]

坑 5：把泰森多边形误当真实服务范围

Voronoi 只考虑直线距离，不考虑道路、河流、山体、交通拥堵、行政边界和用户选择行为。如果你要分析门店实际覆盖范围，可能更适合使用网络分析、等时圈或基于订单数据的服务区模型。

方法比较：Python、QGIS、ArcGIS Pro 该怎么选

方法	适合场景	优点	限制
Python + GeoPandas + SciPy	批量生成、自动化流程、可重复分析	可写脚本、便于集成数据清洗和统计	需要处理无限区域、投影和几何修复
QGIS Voronoi polygons 工具	单次处理、教学演示、快速出图	界面友好，上手快	批量自动化和复杂属性处理不如 Python 灵活
ArcGIS Pro Create Thiessen Polygons	企业 GIS 项目、ArcGIS 工作流	工具成熟，和地理处理模型集成好	依赖授权环境，自动化通常需要 ArcPy
PostGIS ST_VoronoiPolygons	数据库内空间处理、大数据流程	适合和空间查询、索引、服务端流程结合	SQL 写法和几何集合处理有一定门槛

如果你只是偶尔画一次 Voronoi，QGIS 或 ArcGIS Pro 更快。如果你要对很多城市、很多批次点位反复生成结果，Python计算泰森多边形更容易形成稳定流程。

检查清单：生成结果前后务必确认这些项

• 点数据是否为 Point 几何：不要把线、面或空几何误传给 Voronoi 计算。
• 是否存在重复点：重复点要按业务规则删除或合并。
• 是否已经投影：不要直接用 WGS84 经纬度计算面积和距离。
• 研究区边界是否有效：先修复无效几何，再做 overlay 裁剪。
• 输出结果是否完整覆盖研究区：检查是否有空洞、缺口或异常碎片。
• 每个多边形是否能对应原始点：保留点 ID、名称或编号字段。
• 面积单位是否正确：确认投影坐标单位，再解释 area_m2 或 area_km2。
• 是否符合业务含义：泰森多边形只是最近邻分区，不等于真实交通服务区。

FAQ：Python计算泰森多边形常见问题

1. Python 画 Voronoi 必须用 SciPy 吗？

不一定。SciPy 是常见选择，适合计算 Voronoi 的几何结构。GIS 场景中还需要 GeoPandas 和 Shapely 负责坐标系、属性表、空间裁剪和文件输出。也可以使用 Shapely 的相关能力或 PostGIS 的 ST_VoronoiPolygons，但工作流会有所不同。

2. 为什么我的 Voronoi 多边形没有覆盖整个研究区？

通常是因为生成有限区域时半径设置太小，或者研究区边界比点分布范围大很多。可以增大 radius，或先用研究区外扩范围作为辅助包络。最终仍然需要用研究区边界裁剪。

3. Python计算泰森多边形后为什么面积不对？

最常见原因是使用了经纬度坐标系。经纬度单位是度，不能直接解释为平方米。请先把点和边界投影到米制坐标系，再计算 geometry.area。

4. Voronoi咋画才能在 QGIS 里继续编辑？

建议导出为 GeoPackage。示例代码中的 voronoi_result.gpkg 可以直接拖入 QGIS。相比 Shapefile，GeoPackage 对中文字段、长字段名和多图层管理更友好。

5. 点在边界外会影响泰森多边形吗？

会。边界外的点也会参与最近邻分割，从而影响边界内的 Voronoi 结果。是否保留边界外点取决于业务含义。例如分析城市内门店服务范围时，城市边界附近的外部门店可能确实会影响居民最近选择。

6. 泰森多边形可以用于人口分配吗？

可以作为一种简化方法，但要谨慎。它只按最近点分配空间，不考虑人口实际流动、道路阻隔和服务能力。如果用于人口、客户或订单分配，最好结合真实业务数据验证。

结论：Python 画 Voronoi 的关键不是一行代码，而是完整 GIS 流程

Python计算泰森多边形的核心步骤并不复杂：读取点数据，投影到合适坐标系，用 SciPy 生成 Voronoi，将结果转成面，再按研究区边界裁剪并导出。

真正容易出错的地方在 GIS 细节：坐标系是否正确、重复点是否处理、边界是否有效、结果是否被正确裁剪、面积单位是否可信。只要把这些环节检查清楚，Voronoi 就能从一个计算几何结果，变成可以用于 GIS 分析和制图的泰森多边形图层。

如果你的任务是批量生成多个区域的最近邻分区，建议把本文代码整理成函数，并固定输入字段、坐标系、输出格式和质量检查步骤。这样比每次在软件界面中手动操作更稳定，也更容易复现。