如何在Python中生成带最大半径限制的Voronoi单元格？

我的需求和R+ggplot中的一个类似问题一致，但需要Python实现：如何生成单元格带有最大扩展/生长限制的二维Voronoi图？生成的图中部分单元格会带有圆弧边界，同时存在未被任何单元格填充的空间区域。

从图像处理角度理解，所需输出相当于对每个Voronoi单元格，与以该单元格中心为圆心、指定半径为半径的圆形执行交集操作。

注：我正在处理地理定位数据，目前可以接受多数Voronoi库将经纬度视为笛卡尔坐标带来的偏差，这是另一个独立问题。

注意：本问题并非关于裁剪或处理“无限”Voronoi单元格，因此以下类型的解决方案不适用：

• 裁剪Voronoi图（Python）
• 在Python中限制无限Voronoi单元格

相关主题参考：

• 改变Voronoi单元格的生长速率
• R语言中的加权Voronoi多边形
• 创建加权泰森多边形
• PyQGIS实现加权Voronoi

R+ggplot解决方案示例效果： #

核心思路是利用Shapely库的几何交集操作，将每个Voronoi单元格与对应种子点的圆形区域做交集，从而限制单元格的最大扩展范围。以下是修改后的完整代码：

from typing import Tuple, List, Union
from pathlib import Path

import geovoronoi
import matplotlib.colors as clr
import matplotlib.pyplot as plt
import numpy as np
from geographiclib.geodesic import Geodesic
from shapely.geometry import Polygon, Point, MultiPolygon
from simplekml import LineStyle, Color, PolyStyle, Container, Kml

WGS84_Tool = Geodesic.WGS84

T_Color_List = List[Union[clr.Colormap, clr.LinearSegmentedColormap, clr.ListedColormap]]


def generate_n_colors(n, cmap_name='tab20') -> T_Color_List:
    """
    获取matplotlib色图`cmap_name`中的`n`种颜色。如果`n`超过色图的颜色数量，会循环复用颜色。
    :param n: 需要生成的颜色数量
    :param cmap_name: matplotlib色图名称
    :return: 包含`n`种颜色的列表
    """
    pt_region_colormap = plt.get_cmap(cmap_name)
    max_i = len(pt_region_colormap.colors)
    return [pt_region_colormap(i % max_i) for i in range(n)]


def _plot_cell_and_seed(kml: Container,
                        name: str,
                        region_polygon: Union[Polygon, MultiPolygon],
                        seed_coords: Tuple[float, float],
                        _kml_color: str):
    # 处理MultiPolygon情况，拆分多个多边形绘制
    if isinstance(region_polygon, MultiPolygon):
        for idx, poly in enumerate(region_polygon.geoms):
            _p = kml.newpolygon(name=f"{name} zone_part_{idx}",
                                outerboundaryis=[(lon, lat)
                                                for lat, lon
                                                in poly.exterior.coords], )
            _p.style.linestyle = LineStyle(color=Color.darkgrey, width=1.)
            _p.style.polystyle = PolyStyle(color=_kml_color)
    else:
        _p = kml.newpolygon(name=f"{name} zone",
                            outerboundaryis=[(lon, lat)
                                            for lat, lon
                                            in region_polygon.exterior.coords], )
        _p.style.linestyle = LineStyle(color=Color.darkgrey, width=1.)
        _p.style.polystyle = PolyStyle(color=_kml_color)
    
    p = kml.newpoint(coords=[seed_coords],
                     name=name)
    p.style.iconstyle.icon.href = "http://maps.google.com/mapfiles/kml/shapes/placemark_circle.png"
    p.style.iconstyle.scale = 0.5
    p.style.labelstyle.scale = 0.5


def plot_regions(region_polys,
                 region_pts,
                 seeds_coords_list: np.ndarray,
                 seeds_names: List[str],
                 kml: Container,
                 colors: T_Color_List,
                 ):
    assert (len(seeds_names) == len(seeds_coords_list))
    index = 0
    for region_id, region_polygon in region_polys.items():
        _cell_point_indexes = region_pts[region_id]
        _cell_seed_coords = seeds_coords_list[_cell_point_indexes][0]
        name = seeds_names[_cell_point_indexes[0]]
        _kml_airport_coords = (_cell_seed_coords[-1], _cell_seed_coords[0])
        _mpl_color = colors[index]
        _mpl_hexa_color = clr.to_hex(_mpl_color, keep_alpha=True)
        _hexa_color_no_sharp = _mpl_hexa_color.split("#")[-1]
        _kml_color = Color.hexa(_hexa_color_no_sharp)
        _kml_color = Color.changealphaint(alpha=7 * 255 // 10, gehex=_kml_color)
        _plot_cell_and_seed(kml=kml,
                            name=name,
                            region_polygon=region_polygon,
                            seed_coords=_kml_airport_coords,
                            _kml_color=_kml_color)
        index += 1


# 地理边界框
geo_boundaries = {"min": {"lat": +30, "lon": -12},
                  "max": {"lat": +75, "lon": +35}, }
# 生成随机种子点
n = 150
seeds_coords_list = np.dstack(
    [np.random.uniform(low=geo_boundaries["min"]["lat"], high=geo_boundaries["max"]["lat"], size=n),
     np.random.uniform(low=geo_boundaries["min"]["lon"], high=geo_boundaries["max"]["lon"], size=n), ]).reshape((n, 2))
seeds_names = [f"{lat:+_.2f};{lon:+_.2f}" for lat, lon in seeds_coords_list]

boundary_points = geovoronoi.coords_to_points([[geo_boundaries["min"]["lat"], geo_boundaries["min"]["lon"]],
                                               [geo_boundaries["min"]["lat"], geo_boundaries["max"]["lon"]],
                                               [geo_boundaries["max"]["lat"], geo_boundaries["max"]["lon"]],
                                               [geo_boundaries["max"]["lat"], geo_boundaries["min"]["lon"]],
                                               [geo_boundaries["min"]["lat"], geo_boundaries["min"]["lon"]]])
boundary_polygon = Polygon(boundary_points)

# 生成原始Voronoi区域
region_polys, region_pts = geovoronoi.voronoi_regions_from_coords(seeds_coords_list, boundary_polygon)

# --------------------------
# 核心处理：限制单元格最大扩展
# --------------------------
MAX_RADIUS = 5.0  # 设置最大半径（单位：度，因经纬度被视为笛卡尔坐标）
restricted_region_polys = {}

for region_id, poly in region_polys.items():
    # 获取当前单元格对应的种子点
    seed_idx = region_pts[region_id][0]
    seed_lat, seed_lon = seeds_coords_list[seed_idx]
    # 创建以种子点为圆心、指定半径的圆形
    seed_point = Point(seed_lon, seed_lat)  # Shapely的Point顺序为(lon, lat)
    seed_circle = seed_point.buffer(MAX_RADIUS)
    # 计算Voronoi单元格与圆形的交集
    restricted_poly = poly.intersection(seed_circle)
    # 保存处理后的区域
    restricted_region_polys[region_id] = restricted_poly

# --------------------------

# 生成KML文件
kdoc = Kml()
p_kml = Path.cwd() / "voronoi_restricted.kml"
colors: T_Color_List = generate_n_colors(len(restricted_region_polys))

# 绘制处理后的区域
plot_regions(region_polys=restricted_region_polys,
             region_pts=region_pts,
             seeds_coords_list=seeds_coords_list,
             seeds_names=seeds_names,
             kml=kdoc.newfolder(name="restricted regions"),
             colors=colors)

# 可选：绘制原始区域用于对比
plot_regions(region_polys=region_polys,
             region_pts=region_pts,
             seeds_coords_list=seeds_coords_list,
             seeds_names=seeds_names,
             kml=kdoc.newfolder(name="raw regions"),
             colors=colors)

print("保存KML文件")
kdoc.save(p_kml.as_posix())
print(p_kml.as_uri())

关键说明： #

1. 核心操作：使用shapely.geometry.Point.buffer()创建圆形区域，再通过intersection()方法将Voronoi单元格与圆形做交集，得到受限的单元格区域。
2. 半径单位：由于将经纬度视为笛卡尔坐标，这里的半径单位是度；如果需要实际距离（如公里），需先将经纬度转换为投影坐标系（如UTM）后再操作。
3. MultiPolygon处理：交集操作可能生成MultiPolygon（比如单元格被圆形切割成多个部分），绘制时需拆分处理。
4. 对比查看：代码中保留了原始区域的绘制，方便对比受限前后的效果。

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