6

我想使用世界的球形性质（不是它的投影）制作一个带有 voronoi 镶嵌的世界地图，类似于使用 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 ？

4

1 回答 1

3

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

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

警告消息：在 st_voronoi.sfc(sf::st_union(points)) 中：st_voronoi 未正确对经度/纬度数据进行三角测量

从深入研究代码来看，这是一个已知的限制。查看CPL_geos_voronoi的 C++ 代码，看起来它直接调用了 GEOS 方法来构建 Voronoi 图。可能值得打开一个sf issue来表明这是一个你会重视的功能（如果没有人告诉开发人员特定功能会有用，它们不会被优先考虑......）这并不让我感到惊讶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软件包）意外设置了选项？

于 2021-07-11T17:50:55.457 回答