如何在R中为Voronoi树图添加多边形?

apeeds0o  于 2023-01-22  发布在  其他
关注(0)|答案(2)|浏览(105)

我有一个像这样的数据框,它包含了每个国家的向日葵籽产量,我想在这个数据旁边添加多边形数据,这样我就可以用ggplot2来绘制它。
我被告知使用这个网站:https://observablehq.com/@ladataviz/wip-voronoi-data-generator,我想了解如何创建多边形并绘制圆形Voronoi图。
我以前也创建过类似的post,但这里的问题非常不同,我想找到创建面数据的方法

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

创建于2023年1月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

我希望我的数据看起来像这样。

rhfm7lfc

rhfm7lfc1#

在R中创建Voronoi流苏相对容易,但创建Voronoi * 树图 * 就比较难了。链接的Q&A使用voronoiTreemap包来完成,它本质上只是一个JavaScript库的 Package 器。据我所知,这是唯一一个发布的生成Voronoi树图的R包。
我们有两种选择,一种是从头开始计算多边形,另一种是从voronoiTreemap的SVG输出中提取多边形。
关于第一个选项,这不是一个小问题。要了解它有多复杂,也要在R中得到一个完整的解决方案,你可以查看this fantastic article by 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

现在单击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))

生成以下多边形数据框:

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

我们可以看到这是一个Voronoi树图多边形的数据框,通过执行以下操作:

ggplot(polygons, aes(X, Y, fill = country)) + 
  geom_polygon() +
  coord_fixed(0.52) +
  theme_void()

2vuwiymt

2vuwiymt2#

也许有一个聪明的算法,但这里是你如何可以作出这样一个图表的蛮力。
您的数据

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))

通过优化找到解的函数

library(terra)

vtreeMap <- function(d) {

    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))
    crop(voronoi(vp, bnd=p), p)
}

调用函数

set.seed(3)
vp <- vtreeMap(df$prod)
[1] "MSE: 0.01187"

并绘制它

library(RColorBrewer)
vp$country <- df$country
plot(vp, col=brewer.pal(6, "Set2"), axes=FALSE, lwd=4, border="white", mar=rep(0.1, 4))
text(vp, "country", halo=TRUE)

您可能需要稍微调整优化过程(不同的算法、附加选项)以获得最佳结果(低MSE)。
例如,您可以使用

opt <- optim(xy, f, method="BFGS", control=list(abstol=0.001, maxit=500))

如果你不喜欢这个特定的解决方案,改变种子,然后再试一次,直到你找到一个让你满意的。
如果您想使用ggplot2,可以执行以下操作

library(tidyterra)
library(ggplot2)
ggplot(vp) + geom_spatvector(aes(fill = country)) + theme_void()

相关问题