如何使用 sf 真正计算球形 voronoi 图？
Art*_*hur 6 r geospatial ggplot2 r-sf ggspatial
我想使用世界的球形性质（不是它的投影）制作带有 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？
这是一个基于 St\xc3\xa9phane Laurent\ 方法构建的方法，但输出sf对象。
让我们获取sf世界所有首都的对象：
library(sf)\n#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE\n\ncapitals <- do.call(rbind,\n subset(maps::world.cities, capital == 1, select = c("long", "lat")) |>\n as.matrix() |>\n asplit(1) |>\n lapply(st_point) |>\n lapply(st_sfc) |>\n lapply(st_sf, crs = \'WGS84\')) |>\n `st_geometry<-`(\'geometry\') |>\n cbind(city = subset(maps::world.cities, capital == 1, select = c("name")))\n\ncapitals\n#> Simple feature collection with 230 features and 1 field\n#> Geometry type: POINT\n#> Dimension: XY\n#> Bounding box: xmin: -176.13 ymin: -51.7 xmax: 179.2 ymax: 78.21\n#> Geodetic CRS: WGS 84\n#> First 10 features:\n#> name geometry\n#> 1 \'Amman POINT (35.93 31.95)\n#> 2 Abu Dhabi POINT (54.37 24.48)\n#> 3 Abuja POINT (7.17 9.18)\n#> 4 Accra POINT (-0.2 5.56)\n#> 5 Adamstown POINT (-130.1 -25.05)\n#> 6 Addis Abeba POINT (38.74 9.03)\n#> 7 Agana POINT (144.75 13.47)\n#> 8 Algiers POINT (3.04 36.77)\n#> 9 Alofi POINT (-169.92 -19.05)\n#> 10 Amsterdam POINT (4.89 52.37)\n还有我们的世界地图：
\nworld_map <- rnaturalearth::ne_countries(\n scale = \'small\', \n type = \'map_units\',\n returnclass = \'sf\')\n现在我们使用 St\xc3\xa9phane Laurent\ 的方法对球体进行细分，然后将投影反转回球坐标。这允许转换回sf，但我们必须小心分割“环绕”180/-180 经度线的任何对象：
voronoi <- capitals %>%\n st_coordinates() %>%\n `*`(pi/180) %>%\n cbind(1) %>%\n pracma::sph2cart() %>%\n sphereTessellation::VoronoiOnSphere() %>%\n lapply(\\(x) rbind(t(x$cell), t(x$cell)[1,])) %>%\n lapply(\\(x) {\n n <- nrow(x) - 1\n lapply(seq(n), function(i) {\n a <- approx(x[i + 0:1, 1], x[i + 0:1, 2], n = 1000)\n b <- approx(x[i + 0:1, 1], x[i + 0:1, 3], n = 1000)\n d <- cbind(a$x, a$y, b$y) |> pracma::cart2sph() \n d <- d[,1:2] * 180/pi\n if(max(abs(diff(d[,1]))) > 180) {\n s <- which.max(abs(diff(d[,1])))\n d <- list(d[1:s, ], d[(s+1):nrow(d),])\n }\n d })}) |>\n lapply(\\(x) {\n st_geometrycollection(lapply(x, \\(y) {\n if(class(y)[1] == "list") {\n st_multilinestring(y)\n } else {\n st_linestring(y)\n }}))}) %>%\n lapply(st_sfc) %>%\n lapply(st_sf, crs = \'WGS84\') %>%\n {do.call(rbind, .)} %>%\n `st_geometry<-`(\'geometry\')\n现在我们将 Voronoi 网格作为sf对象，因此我们可以使用以下方法绘制它ggplot：
library(ggplot2)\n\nggplot() +\n geom_sf(data = world_map, fill = "cornsilk", color = NA) +\n geom_sf(data = voronoi, color = "gray40") +\n geom_sf(data = capitals, color = "black", size = 0.2) + \n coord_sf(crs = "ESRI:53011") +\n theme(panel.background = element_rect(fill = "lightblue"))\n
\n
附录
\n虽然上述解决方案适用于在整个地球上绘制曲面细分，但如果我们只想获得陆地区域的多边形，我们可以按如下方式执行：
\n首先，我们将世界地图上的所有陆地合并起来
\nwm <- st_make_valid(world_map) |> st_union()\n现在我们得到了 Voronoi 瓦片顶点的坐标：
\npieces <- capitals %>%\n st_coordinates() %>%\n `*`(pi/180) %>%\n cbind(1) %>%\n pracma::sph2cart() %>%\n sphereTessellation::VoronoiOnSphere() %>%\n lapply(\\(x) rbind(t(x$cell), t(x$cell)[1,])) %>%\n lapply(pracma::cart2sph) %>%\n lapply(\\(x) x[,1:2] * 180/pi)\n现在我们需要找到跨越 -180 / 180 线的图块：
\ncomplete <- pieces %>% sapply(\\(x) abs(diff(c(min(x[,1]), max(x[,1])))) < 180)\n现在我们将它们分开并将它们变成多边形，找到它们与世界地图的交集：
\norphans <- pieces[!complete] %>%\n lapply(\\(x) {x[,1] + 180 -> x[,1]; x}) %>%\n lapply(\\(x) st_polygon(list(x)) |> st_sfc(crs = "WGS84")) %>%\n lapply(\\(x) {\n west <- st_intersection(x, matrix(c(-180, -0.001, -0.001, -180, -180, \n -89, -89, 89, 89, -89), ncol = 2) |>\n list() |> st_polygon() |> st_sfc(crs = "WGS84"))\n east <- st_intersection(x, matrix(c(0, 180, 180, 0, 0, \n -89, -89, 89, 89, -89), ncol = 2) |>\n list() |> st_polygon() |> st_sfc(crs = "WGS84"))\n west <- st_coordinates(west)[,1:2]\n east <- st_coordinates(east)[,1:2]\n west[,1] <- west[,1] + 180\n east[,1] <- east[,1] - 180\n w <- st_polygon(list(west)) |> st_sfc(crs = "WGS84") |> st_intersection(wm)\n e <- st_polygon(list(east)) |> st_sfc(crs = "WGS84") |> st_intersection(wm)\n st_combine(st_union(e, w))\n }) %>%\n lapply(st_sf) %>%\n lapply(\\(x) { if(nrow(x) > 0) { st_segmentize(x, 100000) } else {\n st_point(matrix(c(0, 0), ncol = 2)) |> \n st_sfc(crs = "WGS84") |> st_sf()\n }\n }) %>% \n lapply(\\(x) `st_geometry<-`(x, \'geometry\')) %>%\n {do.call(rbind, .)} %>%\n cbind(city = capitals$name[!complete])\n我们可以像这样对非环绕图块进行交叉：
\nnon_orphans <- pieces %>%\n subset(complete) %>%\n lapply(list) %>%\n lapply(st_polygon) %>%\n lapply(st_sfc, crs = "WGS84") %>%\n lapply(st_intersection, y = wm) %>%\n lapply(st_sf) %>%\n lapply(\\(x) { if(nrow(x) > 0) { st_segmentize(x, 100000) } else {\n st_point(matrix(c(0, 0), ncol = 2)) |> \n st_sfc(crs = "WGS84") |> st_sf()\n }\n }) %>%\n lapply(\\(x) `st_geometry<-`(x, \'geometry\')) %>%\n {do.call(rbind, .)} %>%\n cbind(city = capitals$name[complete])\n最后，我们将所有这些绑定到一个sf对象中：
voronoi <- rbind(orphans, non_orphans)\nvoronoi <- voronoi[!st_is_empty(voronoi),]\nvoronoi <- voronoi[sapply(voronoi$geometry, \\(x) class(x)[2] != "POINT"),]\n现在我们准备好绘制了。让我们定义一个调色板函数，它给出的结果与链接的示例类似：
\nf <- colorRampPalette(c("#dae7b4", "#c5b597", "#f3dca8", "#b4b6e7", "#d6a3a4"))\n我们还将创建一个背景“地球”和一个平滑的网格来绘制我们的地图，如示例所示：
\ngrid <- lapply(seq(-170, 170, 10), \\(x) rbind(c(x, -89), c(x, 0), c(x, 89))) |>\n lapply(st_linestring) |>\n lapply(\\(x) st_sfc(x, crs = "WGS84")) |>\n lapply(\\(x) st_segmentize(x, dfMaxLength = 100000)) |>\n c(\n lapply(seq(-80, 80, 10), \\(x) rbind(c(-179, x), c(0, x), c(179, x))) |>\n lapply(st_linestring) |>\n lapply(\\(x) st_sfc(x, crs = "WGS84"))\n ) |>\n lapply(st_sf) |>\n lapply(\\(x) `st_geometry<-`(x, \'geometry\')) %>%\n {do.call(rbind, .)}\n\nglobe <- st_polygon(list(cbind(c(-179, 179, 179, -179, -179), \n c(-89, -89, 89, 89, -89)))) |> \n st_sfc(crs = "WGS84") |> \n st_segmentize(100000)\n最终结果是sf链接示例的忠实版本：
ggplot() +\n geom_sf(data = globe, fill = "#4682b4", color = "black") +\n geom_sf(data = voronoi, color = "black", aes(fill = city)) +\n geom_sf(data = capitals, color = "black", size = 1) + \n geom_sf(data = grid, color = "black", linewidth = 0.2) +\n coord_sf(crs = "ESRI:53011") +\n scale_fill_manual(values = f(nrow(voronoi))) +\n theme(panel.background = element_blank(),\n legend.position = "none",\n panel.grid = element_blank())\n
创建于 2023-06-24，使用reprex v2.0.2
\n- 这是惊人的！stackoverflow R 英雄 Stéphane Laurent 和 Allan Cameron 编写代码的力量！谢谢！我希望我能把赏金分给你们两个。 (3认同)
对于sf我不知道，但你可以使用sphereTessellation包代替。
\nlibrary(pracma) # for the sph2cart function\nlibrary(maps)\ndata(world.cities)\ndata(worldMapEnv)\n\nworld <- map("world", plot = FALSE)\ncountries <- sph2cart(cbind(world$x*pi/180, world$y*pi/180, 1))\n\ncapitals_ll <- as.matrix(\n subset(world.cities, capital == 1, select = c("long", "lat"))\n) * pi / 180\ncapitals <- sph2cart(cbind(capitals_ll, 1))\n\n\nlibrary(rgl)\nopen3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.7)\nsphereMesh <- Rvcg::vcgSphere()\nshade3d(sphereMesh, color = "cyan", polygon_offset = 1)\nlines3d(countries)\npoints3d(capitals, size = 4, color = "red")\nsnapshot3d("world.png", webshot = FALSE)\n
library(sphereTessellation)\nvor <- VoronoiOnSphere(capitals)\nopen3d(windowRect = 50 + c(0, 0, 512, 512), zoom = 0.7)\nplotVoronoiOnSphere(\n vor, colors = "white", edges = TRUE, ecolor = "green",\n sites = TRUE, scolor = "red", polygon_offset = 1\n)\nlines3d(countries)\nsnapshot3d("world_voronoi.png", webshot = FALSE)\n
但我不知道如何将 Vorono\xc3\xaf 单元格剪辑到国家/地区。
\n（这个答案并没有告诉你如何做，但确实告诉你出了什么问题。）
\n当我运行这段代码时，我得到了
\n\n\n警告消息：\n在 st_voronoi.sfc(sf::st_union(points)) 中：\nst_voronoi 未正确对经度/纬度数据进行三角测量
\n
从深入研究代码来看，这似乎是一个已知的限制。查看CPL_geos_voronoi的 C++ 代码，它看起来像是直接调用 GEOS 方法来构建 Voronoi 图。可能值得打开一个sf 问题来表明这是一个您会重视的功能（如果没有人告诉开发人员特定功能会很有用，那么它们就不会被优先考虑......）这并不奇怪我认为 GEOS 不会自动进行球面几何的计算。尽管 S2 代码库在多个地方提到了 Voronoi 图，但看起来并没有 GEOS 算法的直接替代品...其他语言有多种球形 Voronoi 图的实现（例如Python），但有人可能不得不将它们移植到 R （或 C++）......
\n如果我真的需要这样做，我可能会尝试弄清楚如何从 R 中调用 Python 代码（将数据从sf格式导出到 Python 需要的任何格式，然后将结果重新导入到适当的sf格式......）
打印代码sf:::st_voronoi.sfc：
function (x, envelope = st_polygon(), dTolerance = 0, bOnlyEdges = FALSE) \n{\n if (compareVersion(CPL_geos_version(), "3.5.0") > -1) {\n if (isTRUE(st_is_longlat(x))) \n warning("st_voronoi does not correctly triangulate longitude/latitude data")\n st_sfc(CPL_geos_voronoi(x, st_sfc(envelope), dTolerance = dTolerance, \n bOnlyEdges = as.integer(bOnlyEdges)))\n }\n else stop("for voronoi, GEOS version 3.5.0 or higher is required")\n}\n也就是说，如果GEOS版本低于3.5.0，则操作完全失败。如果它>= 3.5.0（sf:::CPL_geos_version()报告我有版本3.8.1），并且正在使用长纬度数据，则应该发出警告（但无论如何计算都会完成）。
我第一次运行这个时，我没有收到警告；我检查并options("warn")设置为-1（抑制警告）。我不确定为什么从干净会话运行的 \xe2\x80\x94 确实给了我警告。也许管道中的某些东西（例如rnaturalearth告诉我我需要安装该rnaturalearthdata包）意外设置了该选项？
| 归档时间： |
|
| 查看次数： |
76 次 |
| 最近记录： |