这是一种基于Stéphane Laurent方法的方法，但输出sf对象。
让我们获取所有世界首都的sf对象：
library(sf)
capitals <- do.call(rbind,
subset(maps::world.cities, capital == 1, select = c("long", "lat")) |>
as.matrix() |>
asplit(1) |>
lapply(st_point) |>
lapply(st_sfc) |>
lapply(st_sf, crs = 'WGS84')) |>
`st_geometry<-`('geometry') |>
cbind(city = subset(maps::world.cities, capital == 1, select = c("name")))
capitals
我们的世界地图：
world_map <- rnaturalearth::ne_countries(
scale = 'small',
type = 'map_units',
returnclass = 'sf')
现在我们使用 Stéphane Laurent 的方法来将球体网格化，然后将投影反转回球面坐标。这允许回到
sf，但我们必须小心地拆分任何跨越 180/-180 经度线的对象。
voronoi <- capitals %>%
st_coordinates() %>%
`*`(pi/180) %>%
cbind(1) %>%
pracma::sph2cart() %>%
sphereTessellation::VoronoiOnSphere() %>%
lapply(\(x) rbind(t(x$cell), t(x$cell)[1,])) %>%
lapply(\(x) {
n <- nrow(x) - 1
lapply(seq(n), function(i) {
a <- approx(x[i + 0:1, 1], x[i + 0:1, 2], n = 1000)
b <- approx(x[i + 0:1, 1], x[i + 0:1, 3], n = 1000)
d <- cbind(a$x, a$y, b$y) |> pracma::cart2sph()
d <- d[,1:2] * 180/pi
if(max(abs(diff(d[,1]))) > 180) {
s <- which.max(abs(diff(d[,1])))
d <- list(d[1:s, ], d[(s+1):nrow(d),])
}
d })}) |>
lapply(\(x) {
st_geometrycollection(lapply(x, \(y) {
if(class(y)[1] == "list") {
st_multilinestring(y)
} else {
st_linestring(y)
}}))}) %>%
lapply(st_sfc) %>%
lapply(st_sf, crs = 'WGS84') %>%
{do.call(rbind, .)} %>%
`st_geometry<-`('geometry')
现在我们有了Voronoi网格作为一个
sf对象，所以我们可以使用
ggplot来绘制它。
library(ggplot2)
ggplot() +
geom_sf(data = world_map, fill = "cornsilk", color = NA) +
geom_sf(data = voronoi, color = "gray40") +
geom_sf(data = capitals, color = "black", size = 0.2) +
coord_sf(crs = "ESRI:53011") +
theme(panel.background = element_rect(fill = "lightblue"))
附录
虽然上述解决方案适用于在整个地球上绘制镶嵌图案，但如果我们只想得到陆地区域的多边形，可以按照以下步骤进行：
首先，我们从世界地图中合并所有陆地区域。
wm <- st_make_valid(world_map) |> st_union()
现在我们得到了Voronoi图块顶点的坐标：
pieces <- capitals %>%
st_coordinates() %>%
`*`(pi/180) %>%
cbind(1) %>%
pracma::sph2cart() %>%
sphereTessellation::VoronoiOnSphere() %>%
lapply(\(x) rbind(t(x$cell), t(x$cell)[1,])) %>%
lapply(pracma::cart2sph) %>%
lapply(\(x) x[,1:2] * 180/pi)
现在我们需要找到跨越-180 / 180线的瓦片。
complete <- pieces %>% sapply(\(x) abs(diff(c(min(x[,1]), max(x[,1])))) < 180)
现在我们将它们拆分并转换为多边形，找到它们与世界地图的交集：
orphans <- pieces[!complete] %>%
lapply(\(x) {x[,1] + 180 -> x[,1]; x}) %>%
lapply(\(x) st_polygon(list(x)) |> st_sfc(crs = "WGS84")) %>%
lapply(\(x) {
west <- st_intersection(x, matrix(c(-180, -0.001, -0.001, -180, -180,
-89, -89, 89, 89, -89), ncol = 2) |>
list() |> st_polygon() |> st_sfc(crs = "WGS84"))
east <- st_intersection(x, matrix(c(0, 180, 180, 0, 0,
-89, -89, 89, 89, -89), ncol = 2) |>
list() |> st_polygon() |> st_sfc(crs = "WGS84"))
west <- st_coordinates(west)[,1:2]
east <- st_coordinates(east)[,1:2]
west[,1] <- west[,1] + 180
east[,1] <- east[,1] - 180
w <- st_polygon(list(west)) |> st_sfc(crs = "WGS84") |> st_intersection(wm)
e <- st_polygon(list(east)) |> st_sfc(crs = "WGS84") |> st_intersection(wm)
st_combine(st_union(e, w))
}) %>%
lapply(st_sf) %>%
lapply(\(x) { if(nrow(x) > 0) { st_segmentize(x, 100000) } else {
st_point(matrix(c(0, 0), ncol = 2)) |>
st_sfc(crs = "WGS84") |> st_sf()
}
}) %>%
lapply(\(x) `st_geometry<-`(x, 'geometry')) %>%
{do.call(rbind, .)} %>%
cbind(city = capitals$name[!complete])
我们可以这样对非环绕瓦片进行交叉处理：
non_orphans <- pieces %>%
subset(complete) %>%
lapply(list) %>%
lapply(st_polygon) %>%
lapply(st_sfc, crs = "WGS84") %>%
lapply(st_intersection, y = wm) %>%
lapply(st_sf) %>%
lapply(\(x) { if(nrow(x) > 0) { st_segmentize(x, 100000) } else {
st_point(matrix(c(0, 0), ncol = 2)) |>
st_sfc(crs = "WGS84") |> st_sf()
}
}) %>%
lapply(\(x) `st_geometry<-`(x, 'geometry')) %>%
{do.call(rbind, .)} %>%
cbind(city = capitals$name[complete])
最后，我们将所有这些绑定到一个单独的
sf 对象中：
voronoi <- rbind(orphans, non_orphans)
voronoi <- voronoi[!st_is_empty(voronoi),]
voronoi <- voronoi[sapply(voronoi$geometry, \(x) class(x)[2] != "POINT"),]
现在我们准备好绘图了。让我们定义一个调色板函数，其结果类似于你提供的示例链接:
f <- colorRampPalette(c("#dae7b4", "#c5b597", "#f3dca8", "#b4b6e7", "#d6a3a4"))
我们还将创建一个背景“地球仪”和一个平滑的网格，以在我们的地图上绘制，就像示例中一样。
grid <- lapply(seq(-170, 170, 10), \(x) rbind(c(x, -89), c(x, 0), c(x, 89))) |>
lapply(st_linestring) |>
lapply(\(x) st_sfc(x, crs = "WGS84")) |>
lapply(\(x) st_segmentize(x, dfMaxLength = 100000)) |>
c(
lapply(seq(-80, 80, 10), \(x) rbind(c(-179, x), c(0, x), c(179, x))) |>
lapply(st_linestring) |>
lapply(\(x) st_sfc(x, crs = "WGS84"))
) |>
lapply(st_sf) |>
lapply(\(x) `st_geometry<-`(x, 'geometry')) %>%
{do.call(rbind, .)}
globe <- st_polygon(list(cbind(c(-179, 179, 179, -179, -179),
c(-89, -89, 89, 89, -89)))) |>
st_sfc(crs = "WGS84") |>
st_segmentize(100000)
最终结果是一个忠实的链接示例的
sf版本。
ggplot() +
geom_sf(data = globe, fill = "#4682b4", color = "black") +
geom_sf(data = voronoi, color = "black", aes(fill = city)) +
geom_sf(data = capitals, color = "black", size = 1) +
geom_sf(data = grid, color = "black", linewidth = 0.2) +
coord_sf(crs = "ESRI:53011") +
scale_fill_manual(values = f(nrow(voronoi))) +
theme(panel.background = element_blank(),
legend.position = "none",
panel.grid = element_blank())
2023-06-24创建，使用reprex v2.0.2
