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

Question

我想使用世界的球形性质(不是它的投影)制作带有 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 ?

Answer 1

(这个答案没有告诉你怎么做，但确实告诉你出了什么问题。)
当我运行这段代码时，我得到了

Warning message:In st_voronoi.sfc(sf::st_union(points)) :st_voronoi does not correctly triangulate longitude/latitude data

从深入研究代码来看，这是一个已知的限制。查看 CPL_geos_voronoi 的 C++ 代码，看起来它直接调用了一个 GEOS 方法来构建 Voronoi 图。可能值得开一个 sf issue表明这是一个您会重视的功能(如果没有人告诉开发人员特定功能将是有用的，它们不会被优先考虑......) GEOS 不会自动进行计算并不让我感到惊讶考虑球面几何。尽管 S2 代码库 mentions Voronoi diagrams in a variety of places ，看起来好像没有 GEOS 算法的替代品......球面 Voronoi 图有其他语言的各种实现(例如 Python)，但有人可能不得不将它们移植到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(抑制警告)。我不知道为什么 - 从干净的 session 中运行确实给了我警告。也许管道中的某些东西(例如 rnaturalearth 告诉我我需要安装 rnaturalearthdata 包)不小心设置了该选项？