如何解决如何使用 sf 真正计算球形 voronoi 图?
我想使用世界的球形性质(不是它的投影)制作一个带有 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?
解决方法
(这个答案不会告诉你怎么做,但会告诉你出了什么问题。)
当我运行这段代码时,我得到了
警告信息: 在 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 图(例如 {{3} }),但有人可能不得不将它们移植到 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(抑制警告)。我不知道为什么 - 从干净的会话中运行确实给了我警告。也许管道中的某些东西(例如 rnaturalearth
告诉我我需要安装 rnaturalearthdata
包)不小心设置了该选项?
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。