如何用sf真实地计算球面voronoi图？

Question

我想用世界的球形(不是它的投影)绘制一张世界地图，类似于使用D3.js，但是用R。

据我所知("再见地球，欢迎S2球面几何学")，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

(这个答案没有告诉你怎么做，而是告诉你出了什么问题。)

当我运行这段代码时

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

从深入研究代码看，这似乎是已知的限制。查看C++代码沃罗内，它看起来像是直接调用一个GEOS方法来构建Voronoi图。打开一个sf发行可能是值得的，以表明这是一个您将重视的特性(如果没有人告诉开发人员特定的特性是有用的，那么它们不会得到优先排序.)我一点也不惊讶，GEOS不会自动进行计算，从而解释球面几何学。虽然S2代码基在不同的地方提到Voronoi图，但似乎没有一个插入替代的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版)，并且正在使用long-lat数据，则警告是要发出假定的(但计算还是完成了)。

第一次运行这个命令时，我没有收到警告；我检查并将options("warn")设置为-1 (抑制警告)。我不知道为什么-逃避干净的治疗确实给了我警告。也许正在进行中的一些事情(例如，rnaturalearth告诉我我需要安装rnaturalearthdata包)意外地设置了选项？