如何使用 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, it looks like it directly calls a GEOS method for building Voronoi diagrams. It might be worth opening an sf issue to indicate that this is a feature you would value (if no-one tells the developer that particular features would be useful, they don't get prioritized ...) It doesn't surprise me that GEOS doesn't automatically do computations that account for spherical geometry. Although the S2 code base mentions Voronoi diagrams in a variety of places, it doesn't look like there is a drop-in replacement for the GEOS algorithm ... there are a variety of implementations in other languages for spherical Voronoi diagrams (e.g. Python) 的 C++ 代码，但有人可能不得不将它们移植到 R（或 C++）...

如果我真的需要这样做，我可能会尝试弄清楚如何从 R 中调用 Python 代码（从 [=11 导出数据） =] 格式为任何 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 数据，则 supposed 发出警告（但无论如何计算都完成了）。

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

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

How to truly calculate a spherical voronoi diagram using sf?

r

ggplot2

geospatial

sf

ggspatial