如何使用 sf 真正计算球形 voronoi 图？

问题描述

我想使用世界的球形性质（不是它的投影）制作一个带有 voronoi 镶嵌的世界地图，类似于 this using D3.js，但使用 R。

据我所知（“Goodbye flat Earth,welcome S2 spherical geometry”），sf 包现在完全基于 s2 包并且应该按照我的需要执行。但我不认为我得到了预期的结果。一个可重现的例子：

library(tidyverse)
library(sf)
library(rnaturalearth)
library(tidygeocoder)

# just to be sure
sf::sf_use_s2(TRUE)

# download map 
world_map <- rnaturalearth::ne_countries(
                scale = 'small',type = 'map_units',returnclass = 'sf')

# addresses that you want to find lat long and to become centroids of the voronoi tessellation 
addresses <- tribble(
~addr,"Juneau,Alaska","Saint Petersburg,Russia","Melbourne,Australia" 
)

# retrive lat long using tidygeocoder
points <- addresses %>% 
          tidygeocoder::geocode(addr,method = 'osm')

# Transform lat long in a single geometry point and join with sf-base of the world
points <- points %>% 
          dplyr::rowwise() %>% 
          dplyr::mutate(point = list(sf::st_point(c(long,lat)))) %>% 
          sf::st_as_sf() %>% 
          sf::st_set_crs(4326)

# voronoi tessellation
voronoi <- sf::st_voronoi(sf::st_union( points ) ) %>% 
     sf::st_as_sf() %>% 
     sf::st_set_crs(4326)

# plot
ggplot2::ggplot() +
    geom_sf(data = world_map,mapping = aes(geometry = geometry),fill = "gray95") +
    geom_sf(data = points,mapping = aes(geometry = point),colour = "red") +
    geom_sf(data = voronoi,mapping = aes(geometry = x),colour = "red",alpha = 0.5)

整个南极洲应该比其他两点更靠近墨尔本。我在这里缺少什么？如何使用sf计算球体上的voronoi？

解决方法

（这个答案不会告诉你怎么做，但会告诉你出了什么问题。）

当我运行这段代码时，我得到了

警告信息：在 st_voronoi.sfc(sf::st_union(points)) 中： st_voronoi 不能正确三角测量经度/纬度数据

从深入研究代码来看，这是一个已知的限制。查看 CPL_geos_voronoi 的 C++ 代码，它看起来像是直接调用 GEOS 方法来构建 Voronoi 图。可能值得打开一个 sf issue 以表明这是一个您会重视的功能（如果没有人告诉开发人员特定功能会很有用，它们不会被优先考虑......）它没有令我惊讶的是 GEOS 不会自动进行考虑球面几何的计算。尽管 S2 代码库 mentions Voronoi diagrams in a variety of places，但它看起来并没有替代 GEOS 算法......对于球形 Voronoi 图（例如 {{3} })，但有人可能不得不将它们移植到 R（或 C++）...

如果我真的需要这样做，我可能会尝试弄清楚如何从 R 中调用 Python 代码（将数据从 sf 格式导出到 Python 需要的任何内容，然后将结果重新导入适当的 sf 格式...)

打印 sf:::st_voronoi.sfc 的代码：

function (x,envelope = st_polygon(),dTolerance = 0,bOnlyEdges = FALSE) 
{
    if (compareVersion(CPL_geos_version(),"3.5.0") > -1) {
        if (isTRUE(st_is_longlat(x))) 
            warning("st_voronoi does not correctly triangulate longitude/latitude data")
        st_sfc(CPL_geos_voronoi(x,st_sfc(envelope),dTolerance = dTolerance,bOnlyEdges = as.integer(bOnlyEdges)))
    }
    else stop("for voronoi,GEOS version 3.5.0 or higher is required")
}

也就是说，如果 GEOS 版本低于 3.5.0，则操作完全失败。如果 >= 3.5.0（sf:::CPL_geos_version() 报告我有 3.8.1 版），并且正在使用长纬度数据，则应该会发出警告（但无论如何计算都完成了）。

我第一次运行时没有收到警告；我检查了一下，options("warn") 被设置为 -1（抑制警告）。我不知道为什么 - 从干净的会话中运行确实给了我警告。也许管道中的某些东西（例如 rnaturalearth 告诉我我需要安装 rnaturalearthdata 包）不小心设置了该选项？

geospatial ggplot2 ggspatial r r sf