LDT*_*LDT 8 math voronoi r ggplot2 tidyverse
我有一个看起来像这样的数据框。它包含了每个国家的葵花籽生产力。我想在这个数据旁边添加多边形数据,这样我就可以用 ggplot2 绘制它。
我被告知使用这个网站: https: //observablehq.com/@ladataviz/wip-voronoi-data-generator,我想了解如何创建多边形并绘制圆形 voronoi 图。
我过去曾创建过类似的帖子,但我在这里的问题非常不同。我想找到一种方法来创建多边形数据
df <- data.frame(country = c("Ukraine", "Russia", "Argentina", "China", "Romania", "Other"),
prod = c(11.0, 10.6, 3.1, 2.4, 2.1, 15.3))
df
#> country prod
#> 1 Ukraine 11.0
#> 2 Russia 10.6
#> 3 Argentina 3.1
#> 4 China 2.4
#> 5 Romania 2.1
#> 6 Other 15.3
Run Code Online (Sandbox Code Playgroud)
创建于 2023-01-20,使用reprex v2.0.2
如果将多边形添加到我的数据中应该如下所示:
x y path split group value
1 472.0117 220.08122253 0 Ukraine Ukraine 11
2 471.8336 217.18476868 1 Ukraine Ukraine 11
3 471.6556 214.28833008 2 Ukraine Ukraine 11
4 471.4776 211.39187622 3 Ukraine Ukraine 11
5 471.2996 208.49542236 4 Ukraine Ukraine 11
6 471.1216 205.59896851 5 Ukraine Ukraine 11
Run Code Online (Sandbox Code Playgroud)
我希望我的数据看起来像这样。
可能有一个聪明的算法可以做到这一点,但这里是你如何通过暴力来制作这样的图表。
您的数据
df <- data.frame(country = c("Ukraine", "Russia", "Argentina", "China", "Romania", "Other"),
prod = c(11.0, 10.6, 3.1, 2.4, 2.1, 15.3))
Run Code Online (Sandbox Code Playgroud)
通过优化找到解决方案的函数
library(terra)
vtreeMap <- function(d, p=NULL) {
if (is.null(p)) {
p <- vect(cbind(0,0), crs="+proj=utm +zone=1") |> buffer(1)
}
A <- expanse(p) * d / sum(d)
f <- function(xy) {
if (any(xy > 1) || any(xy < -1)) return(Inf)
xy <- vect(matrix(xy, ncol=2), crs=crs(p))
e <- extract(p, xy)
if (any(is.na(e[,2]))) return(Inf)
v <- crop(voronoi(xy, bnd=p), p)
mean( (A - expanse(v))^2 )
}
xy <- spatSample(p, length(A)) |> crds() |> as.vector()
opt <- optim(xy, f)
print(paste("MSE:", round(opt$value, 5)))
vp <- vect(matrix(opt$par, ncol=2), crs=crs(p))
crop(voronoi(vp, bnd=p), p)
}
Run Code Online (Sandbox Code Playgroud)
调用函数
set.seed(3)
vtm <- vtreeMap(df$prod)
[1] "MSE: 0.01187"
Run Code Online (Sandbox Code Playgroud)
并绘制它
library(RColorBrewer)
vtm$country <- df$country
plot(vtm, col=brewer.pal(6, "Set2"), axes=FALSE, lwd=4, border="white", mar=rep(0.1, 4))
text(vtm, "country", halo=TRUE)
Run Code Online (Sandbox Code Playgroud)
您可能需要稍微调整优化过程(不同的算法、附加选项)以获得最佳结果(低 MSE)。
例如,您可以使用
opt <- optim(xy, f, method="BFGS", control=list(abstol=0.001, maxit=500))
Run Code Online (Sandbox Code Playgroud)
如果您不喜欢这个特定的解决方案,请更改种子并重试,直到找到令您满意的解决方案。
如果你想使用 ggplot2 你可以这样做
library(tidyterra)
library(ggplot2)
ggplot(vp) + geom_spatvector(aes(fill = country)) + theme_void()
Run Code Online (Sandbox Code Playgroud)
通过这些示例数据
set.seed(1)
regions <- data.frame(country=rep(df$country,each=3), region=c("A","B","C"),
value=sample(10,nrow(df), replace=TRUE))
Run Code Online (Sandbox Code Playgroud)
您可以按国家内的区域细分多边形
out <- vector(nrow(df), mode="list")
for (i in 1:length(df$country)) {
vv <- vtm[vtm$country==df$country[i], ]
regv <- regions$value[regions$country==df$country[i]]
out[[i]] <- vtreeMap(regv, vv)
}
vtm2 <- vect(out)
values(vtm2) <- regions
Run Code Online (Sandbox Code Playgroud)
在 R 中制作 Voronoi 穗状结构相对容易,但制作 Voronoi树形图比较困难。链接的问答是通过使用voronoiTreemap包来实现的,该包本质上只是 JavaScript 库的包装器。据我所知,这是唯一发布的生成 Voronoi 树形图的 R 包。
我们的两个选择是从头开始计算多边形,或者以某种方式从 的 SVG 输出中提取多边形voronoiTreemap。
对于第一个选项,这不是一个小问题。要了解它有多复杂,并在 R 中获得完整的解决方案,您可以查看Paul Murrell 撰写的这篇精彩文章。该代码运行了好几页,并且已经有十多年的历史了,所以我不确定所有依赖项是否仍然有效。令人失望的是,没有人将所有这些都整合到 CRAN 上的一个包中,但也许它有点小众。
如果您对 Paul Murrell 的方法感到困惑,那么您只能尝试从 的输出中获取多边形voronoiTreemap。虽然这个包工作得很好,但输出并不适合收集多边形,而且我们无法访问允许我们在 R 中自己生成多边形的中间计算。这并非不可能,并且有几种方法来解决这个问题,但它们都相当复杂。
以下方法首先使用 正常绘制树状图voronoiTreemap,但不带标签:
library(voronoiTreemap)
library(terra)
library(tidyverse)
df <- data.frame(country = c("Ukraine", "Russia", "Argentina",
"China", "Romania", "Other"),
prod = c(11.0, 10.6, 3.1, 2.4, 2.1, 15.3))
vor <- data.frame(h1 = 'World',
h2 = c('Europe', 'Europe', 'Americas', 'Asia',
'Europe', 'Other'),
h3 = df$country,
color = hcl.colors(nrow(df), palette = 'TealRose'),
weight = df$prod,
codes = "")
vt <- vt_input_from_df(vor)
v <- vt_d3(vt_export_json(vt))
v
Run Code Online (Sandbox Code Playgroud)
现在单击Export -> Save as image并将您的图另存为Rplot.png
现在我们可以做
polygons <- rast('Rplot02.png')[[2]] %>%
app(fun = function(x) ifelse(x > 220, 255, 0)) %>%
as.polygons() %>%
sf::st_as_sf() %>%
filter(lyr.1 == 0) %>%
sf::st_buffer(dist = -0.002) %>%
sf::st_coordinates() %>%
as.data.frame() %>%
mutate(country = df$country[L2], prod = df$prod[L2]) %>%
select(-(L1:L3))
Run Code Online (Sandbox Code Playgroud)
使用我们的多边形生成以下数据框:
head(polygons)
#> X Y country prod
#> 1 0.6460000 0.3970068 Ukraine 11
#> 2 0.6460000 0.4054322 Ukraine 11
#> 3 0.6460501 0.4054499 Ukraine 11
#> 4 0.6461468 0.4054900 Ukraine 11
#> 5 0.6462413 0.4055351 Ukraine 11
Run Code Online (Sandbox Code Playgroud)
通过执行以下操作,我们可以看到这是 Voronoi 树形图的多边形数据框:
ggplot(polygons, aes(X, Y, fill = country)) +
geom_polygon() +
coord_fixed(0.52) +
theme_void()
Run Code Online (Sandbox Code Playgroud)
被这个问题迷住了,我写了一个小包来解决它,你可以通过安装
devtools::install_github("AllanCameron/VoronoiPlus")
Run Code Online (Sandbox Code Playgroud)
它可以通过调用 来处理 Voronoi 映射(如本页问题中所示)voronoi_map,该调用采用权重和组标签。它还可以采用任意形状作为图块的边界,但如果缺少,则默认为单位圆。
library(VoronoiPlus)
res <- voronoi_map(values = df$prod, groups = df$country)
plot(res)
Run Code Online (Sandbox Code Playgroud)
您可以使用以下命令从此对象中提取多边形作为数据框:
polys <- get_polygons(res)
head(polys)
#> geom x y group value
#> 1 1 -0.6436006 -0.7649495 Ukraine 11
#> 2 1 -0.6691306 -0.7431448 Ukraine 11
#> 3 1 -0.7071068 -0.7071068 Ukraine 11
#> 4 1 -0.7431448 -0.6691306 Ukraine 11
#> 5 1 -0.7771460 -0.6293204 Ukraine 11
#> 6 1 -0.8090170 -0.5877853 Ukraine 11
Run Code Online (Sandbox Code Playgroud)
该包还可以处理任意嵌套的组,以通过该voronoi_treemap函数生成真正的树形图,该函数采用公式接口(左侧的权重和右侧的分组变量)
df$region <- c("Europe", "Europe", "Other", "Other", "Europe", "Other")
dat <- voronoi_treemap(prod ~ region + country, data = df)
head(dat)
#> x y group value parent level
#> 1 0.6657716 -0.7460137 Europe 23.7 root 1
#> 2 0.6293204 -0.7771460 Europe 23.7 root 1
#> 3 0.5877853 -0.8090170 Europe 23.7 root 1
#> 4 0.5446390 -0.8386706 Europe 23.7 root 1
#> 5 0.5000000 -0.8660254 Europe 23.7 root 1
#> 6 0.4539905 -0.8910065 Europe 23.7 root 1
Run Code Online (Sandbox Code Playgroud)
这允许嵌套树形图,如下所示:
library(tidyverse)
ggplot(dat[dat$level == 2,], aes(x, y, label = group)) +
geom_polygon(aes(fill = parent)) +
geom_polygon(fill = "white", aes(group = group, alpha = group),
color = "black") +
geom_text(data = . %>% group_by(group) %>%
summarize(x = mean(x), y = mean(y))) +
scale_alpha_discrete(guide = "none") +
coord_equal() +
theme_void()
Run Code Online (Sandbox Code Playgroud)
目前使用的算法也是一种暴力方法,尽管与 Robert Hijmans 演示的算法略有不同。我正在研究一种更有针对性的方法来缩短收敛时间。
一个主要的警告是,该软件包还处于初级阶段,在撰写本文时尚未经过适当的测试或记录。
| 归档时间: |
|
| 查看次数: |
621 次 |
| 最近记录: |