地球不是平的，你的Voronoi图也不应该是

关于精确性的故事，揭示使用Python生成球面地理Voronoi图的强大力量

你可能对Voronoi图及其在地理空间分析中的应用熟悉。如果不熟悉，这里是快速总结：它将平面分割为区域，这些区域由平面上距离给定种子点更近而不是其他任何点的所有点组成。它以数学家Georgy Voronoy命名。你可以在维基百科上了解更多相关信息。

它如何应用在地理空间领域？使用Voronoi图，你可以在更大的尺度上快速找到给定城市居民最近的公共交通站点，比分别为每栋建筑物计算要快。或者你还可以将其用于不同品牌之间的市场份额分析。

在本文中，我想展示使用平面坐标计算的典型Voronoi图与球面Voronoi图之间的差异，并希望展示球面Voronoi图的优势。

尺寸和投影－为什么它很重要？

如果我们想在地图上显示数据，我们必须使用投影。为了在二维平面上显示东西，我们必须将坐标从地球的三维坐标投影到二维平面上。我们都知道并使用最流行的投影是墨卡托投影（Web墨卡托投影或WGS84墨卡托投影，因为大多数地图供应商使用它），最流行的坐标系统是世界大地测量系统1984 – WGS84（或EPSG:4326）。该系统基于角度，经度范围从-180°到180°（从西到东），纬度范围从-90°到90°（从南到北）。

每种二维平面投影都会有一些失真。墨卡托投影是一种保角度投影，这意味着地球上物体之间的角度应该保持不变。纬度越高（或越低），区域和距离的失真就越大。由于Voronoi图严重依赖种子之间的距离，生成图时会传递相同的失真错误。

地球是一个不规则的椭球体，但为了我们的目的，我们可以近似为球体形状。通过在球体上生成Voronoi图，我们可以根据球体表面的弧线准确计算距离。然后，我们可以将生成的球面多边形映射到投影的二维坐标上，我们可以确保分隔相邻Voronoi单元格的线段与定义这些单元格的两个种子点之间的连线垂直。

下面你可以看到我上面描述的角度和距离问题。尽管线在同一点交叉，但Voronoi单元格的形状和角度不同。

另一个问题是，如果使用2D Voronoi图，无法比较世界上不同地区（即不在同一纬度上）的区域，因为这些区域会严重失真。

完整的Jupyter笔记本，包含下面示例中使用的代码，可以在GitHub上找到。这里为了简洁起见跳过了一些函数。

先决条件

安装所需的库

pip install -q srai[voronoi,osm,plotting] geodatasets

导入所需的模块和函数

import geodatasetsimport geopandas as gpdimport matplotlib.pyplot as pltimport plotly.express as pxfrom shapely.geometry import MultiPoint, Pointfrom shapely.ops import voronoi_diagramfrom srai.regionalizers import VoronoiRegionalizer, geocode_to_region_gdf

第一个示例

让我们在地球上定义六个点：北极和南极，以及赤道上的四个点。

earth_points_gdf = gpd.GeoDataFrame(    geometry=[        Point(0, 0),        Point(90, 0),        Point(180, 0),        Point(-90, 0),        Point(0, 90),        Point(0, -90),    ],    index=[1, 2, 3, 4, 5, 6],    crs="EPSG:4326",)

使用Shapely库中的voronoi_diagram生成Voronoi图

def generate_flat_voronoi_diagram_regions(    seeds_gdf: gpd.GeoDataFrame,) -> gpd.GeoDataFrame:    points = MultiPoint(seeds_gdf.geometry.values)    # Generate 2D diagram    regions = voronoi_diagram(points)    # Map geometries to GeoDataFrame    flat_voronoi_regions = gpd.GeoDataFrame(        geometry=list(regions.geoms),        crs="EPSG:4326",    )    # Apply indexes from the seeds dataframe    flat_voronoi_regions.index = gpd.pd.Index(        flat_voronoi_regions.sjoin(seeds_gdf)["index_right"],        name="region_id",    )    # Clip to Earth boundaries    flat_voronoi_regions.geometry = flat_voronoi_regions.geometry.clip_by_rect(        xmin=-180, ymin=-90, xmax=180, ymax=90    )    return flat_voronoi_regions

earth_poles_flat_voronoi_regions = generate_flat_voronoi_diagram_regions(    earth_points_gdf)

使用srai库中的VoronoiRegionalizer生成Voronoi图。它使用scipy库中的SphericalVoronoi实现，并且能够正确将WGS84坐标转换为球坐标系统。

earth_points_spherical_voronoi_regions = VoronoiRegionalizer(    seeds=earth_points_gdf).transform()

让我们在图表上看一下这两者之间的区别。

首先可以看到的是，2D维诺图不在地球上环绕一圈，因为它在一个扁平的笛卡尔平面上工作。球形维诺图正确覆盖了地球，并且不会在反子午线（经度从180°变为-180°的位置）断裂。

为了数值化量化差异，我们可以计算IoU（交并比）指标（或Jaccard指数）来衡量多边形形状之间的差异。该指标的值在0和1之间，其中0表示没有重叠，1表示完全重叠。

def calculate_iou(    flat_regions: gpd.GeoDataFrame, spherical_regions: gpd.GeoDataFrame) -> float:    total_intersections_area = 0    total_unions_area = 0    # 遍历所有区域    for index in spherical_regions.index:        # 找到相匹配的球形和扁平维诺伊区域        spherical_region_geometry = spherical_regions.loc[index].geometry        flat_region_geometry = flat_regions.loc[index].geometry        # 计算它们的交集面积        intersections_area = spherical_region_geometry.intersection(            flat_region_geometry        ).area        # 计算它们的并集面积        # 或者代码：        # spherical_region_geometry.union(flat_region_geometry).area        unions_area = (            spherical_region_geometry.area            + flat_region_geometry.area            - intersections_area        )        # 添加到总和中        total_intersections_area += intersections_area        total_unions_area += unions_area    # 将交集面积除以并集面积    return round(total_intersections_area / total_unions_area, 3)

calculate_iou(    earth_points_flat_voronoi_regions, earth_points_spherical_voronoi_regions)

计算出的值为0.423，这是比较低的，在大尺度上，这些多边形彼此不同，可以在上面的图中很容易地看到。

真实数据示例：使用自动体外除颤器（AEDs）位置分割地球

此示例中使用的数据来自OpenAEDMap，基于OpenStreetMap数据。预处理的文件具有经过滤的位置（确切地说为80694个），其中没有重复定义在彼此之上的节点。

# 将AED位置加载到GeoDataFrameaed_world_gdf = gpd.read_file(    "https://raw.githubusercontent.com/RaczeQ/medium-articles/main/articles/spherical-geovoronoi/aed_world.geojson")

为AED生成维诺伊图

aed_flat_voronoi_regions = generate_flat_voronoi_diagram_regions(aed_world_gdf)aed_spherical_voronoi_regions = VoronoiRegionalizer(    seeds=aed_world_gdf, max_meters_between_points=1_000).transform()

让我们比较这些维诺伊图。

在图表上看这些差异非常明显。 2D版本的边境都是直的，而球面版本的边境在WGS84坐标系中看起来非常弯曲。 您还可以清楚地看到，在平面版本中，许多区域在极点（正交投影聚焦在南极）收敛，而球面版本则没有。 另一个明显的差异是经度线周围的连续性，这在第一个示例中提到过。 从新西兰出现的区域在平面版本上被突然切断。

让我们来看一下IoU值：

calculate_iou(aed_flat_voronoi_regions, aed_spherical_voronoi_regions)

计算出的值是0.511，比第一个示例稍好，但多边形仍大致匹配50％。

缩放至城市规模

让我们在较小的尺度上看看差异。 我们可以选择所有位于伦敦的AED并绘制出来。

greater_london_area = geocode_to_region_gdf("伦敦市")aeds_in_london = aed_world_gdf.sjoin(greater_london_area)

calculate_iou(    aed_flat_voronoi_regions.loc[aeds_in_london.index],    aed_spherical_voronoi_regions.loc[aeds_in_london.index],)

该值为0.675。 这是一个明显的差异。 由于AED分布更密集，形状和距离变小，因此在投影的二维平面和球体上计算的Voronoi图之间的差异变小。

让我们看看一些重叠在一起的个别示例。

这些多边形的区域大部分是匹配的，但您可以看到角度和形状上的差异。 这些差异可能在空间分析中很重要，并可能改变其结果。感兴趣的区域越大，则差异越大。

总结

我希望现在您能明白为什么球面Voronoi图比平面Voronoi图更适合地理空间领域的使用。

目前，该领域中的大多数分析都是使用在投影的二维平面上计算的Voronoi图，这可能导致错误结果。

在Python中，长时间以来没有供地理空间数据科学家和分析师使用的球面Voronoi图的简单解决方案。 现在只需安装一个库就可以轻松实现了。虽然它要比平面解决方案计算的时间略长，因为它必须将点投影到球面坐标系中并从中投射，同时正确剪裁与经度线交叉的多边形，但在您的分析中保持精度应该并不重要。对于JavaScript用户，已经有一个可用的球面Voronoi图的D3.js实现。

免责声明

我是srai库的维护者之一。