bounty还有7天到期。回答此问题可获得+100声望奖励。Arthur Welle正在寻找这个问题的最新答案：这个问题是从2021年开始的，因为也许有一种更新的方法来解决这个问题（我仍然对答案感兴趣，如果有的话，用R完成）。
我想用voronoi镶嵌来制作一个世界Map，它使用了世界的球形特性（而不是它的投影），类似于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 *？
1条答案
按热度按时间hjzp0vay1#
(This答案不会告诉你如何去做，但会告诉你哪里出了问题。）
当我运行这个代码时
警告消息：在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图，在其它语言中存在多种实现（例如，Python），但可能得有人把它们移植到R（或C）...
如果我真的需要这样做，我可能会尝试找出如何从R中调用Python代码（将数据从
sf导出为Python需要的任何格式，然后将结果重新导入为适当的sf...）打印
sf:::st_voronoi.sfc的代码：换句话说，如果GEOS版本小于3.5.0，则操作完全失败。如果它>= 3.5.0（
sf:::CPL_geos_version()报告我的版本是3.8.1），并且正在使用long-lat数据，则假定将发出警告（但无论如何都会完成计算）。我第一次运行时没有收到警告;我检查并将
options("warn")设置为-1（抑制警告）。我不知道为什么-从一个干净的会议运行确实给予了我警告。可能是管道中的东西（例如rnaturalearth告诉我我需要安装rnaturalearthdata包）不小心设置了选项？