【Python】在地球球面、椭球面上构建Voronoi多边形

介绍

1. 在局部区域研究中，我们将地球椭球面投影到平面上，在二维平面上构建Voronoi图（或泰森多边形），已经有成熟的算法及各种编程语言程序库的API。
   以Python为例，有scipy.spatial.Voronoi库。

2. 但是在涉及全球范围内的研究上，投影到平面后变形很大。在精度要求不高时可以将地球表面近似为球面，构建球面Voronoi。
   在SO问答上有相关讨论。
   以Python为例，有scipy.spatial.SphericalVoronoi库，以及一个开源库：py_sphere_Voronoi。

3. 在精度要求较高时，则直接构建地球椭球面上的Voronoi（将地球近似为椭球面）。
   该算法还没有可以直接调用的库。只能参考论文学习。

地球球面上构建Voronoi图

虽然Python有scipy.spatial.SphericalVoronoi库可以用，但是该库支持的输出数据为笛卡尔坐标系的$(x,y,z)$值。

所以，如果输出数据为经纬度数据，则需要参考下文论文[1]中的转换公式：

• 经度不变，先将纬度从椭球面纬度转为球面纬度：
  

• 再将球面经纬度转为笛卡尔坐标：
  

• 然后就可以使用scipy.spatial.SphericalVoronoi库了。

地球椭球面上构建Voronoi图

基于矢量的方法

参考论文[1]：Voronoi tessellation on the ellipsoidal earth for vector data

该论文先以球面的角度生成球面Voronoi的近似点，再通过椭球面上的三点问题，将近似点精确化。

关于“椭球面上的三点问题”以及相关知识的补充与扩展，可以参考以下论文：

[2]The Three-point Problem of The Median Line Turning Point- on the Solutions for the Sphere and Ellipsoid
[3]New Solutions to Classical Geodetic Problems on the Ellipsoid
[4]Intersections on the sphere and ellipsoid
[5]Precise Determination of the Clairaut Constant in Ellipsoidal Geodesy
[6]A Method for the Delimitation of an Equidistant Boundary Between Coastal States on the Surface of a Geodetic Ellipsoid

基于栅格的方法

参考论文[7]：Voronoi diagram generation on the ellipsoidal earth

该论文提出了一种在椭球体上生成误差为$\sqrt{2}\ast r$，其中 $r$ 是网格的空间分辨率。
在他们的方法中，点由网格单元近似，Voronoi 顶点通过计算单元中心之间的距离生成。
用于确定 Voronoi 顶点的单元到单元距离的精度取决于选定的单元大小，无法达到使用矢量方法在椭球体上计算的距离的精度。

矢栅混合方法

参考论文[8]：Determination of circumcenter of triangle on ellipsoidal surface based on map algebra

该方法提高了顶点位置的准确性。然而，对于从密集点数据集生成 Voronoi图，例如，由中分辨率和高分辨率海岸线组成的点，必须使用过小的单元大小，以防止每个单元中包含多个生成点。

尽管基于矢栅混合的方法有其无可置疑的优势，但如论文[7]指出的：“精确的 Voronoi 图只能通过矢量方法获得，在矢量方法中，空间是无限可辨别的，并且每个点都由其精确的坐标来表示”。