在过去的教程中，我一直给大家分享的都是单变量填充地图的绘制方法，但是有时候我们想在一张地图上展示两个变量，这个时候该怎么绘制地图呢？我们可以选择双变量填充地图。如果我们在 Observable 上搜索“Bivariate Choropleth” 就可以看到一大堆使用 D3 绘制的双变量填充地图，非常炫酷，例如：https://observablehq.com/@d3/bivariate-choropleth
如何使用 R 绘制这种双变量填充地图，Timo Grossenbacher 做了一个非常炫酷的示范：
该图的绘制方法可以参考：https://timogrossenbacher.ch/2019/04/bivariate-maps-with-ggplot2-and-sf/ 学习。不过代码过于复杂，所以我这里对里面的代码进行了简化，并且介绍了如何绘制更多分阶的双变量填充地图（该作者只给出了三阶的画法）。
我暂时也没有时候绘制这种地图的数据，所以下面的演示中我是用随机生成的变量。
准备工作
数据
- 2019行政区划：中国省市区县级矢量地图（shp 文件）；
- 海岸线：海岸线矢量数据（shp 文件）；
- 九段线.geojson：九段线和陆地国界线矢量数据（GEOJSON 文件）。
加载 R 包：
library(sf)
library(ggspatial)
library(readr)
library(stringr)
library(tidyr)
library(purrr)
library(hrbrthemes)
library(tidyverse)
library(showtext)
showtext_auto(enable = TRUE)
font_add("songti", regular = "song.otf")
|
这里使用的是 showtext 包设置字体，“song.otf” 文件可以在附件中找到。
准备数据
我随机生成了两个变量作为绘制双变量填充地图的“双变量”：
mycrs <- "+proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"
read_sf("2019行政区划/县.shp") %>%
st_transform(mycrs) %>%
rename(县 = NAME, 县代码 = PAC) %>%
st_simplify(dTolerance = 2000) -> county_simple
county_simple
set.seed(1234)
county_simple %>%
st_drop_geometry() %>%
select(code = 县代码, name = 县) %>%
mutate(v1 = runif(2900, 1, 100),
v2 = runif(2900, 1, 100)) -> df
df
|
数据非常简单，就是每个区县对于两个变量的两个值，v1 和 v2 就是我们需要的 “bivariate”。
上面我展示的图片中使用的是三阶，我自然得比他的更复杂点，我用 4 阶！下面我们将变量 v1 和 v2 使用分位数各分成 4 段：
no_classes <- 4
df %>%
pull(v1) %>%
quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
as.vector() -> quantiles1
quantiles1
df %>%
pull(v2) %>%
quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
as.vector() -> quantiles2
quantiles2
|
对 df 里的 v1 v2 进行分组：
df %>%
mutate(vgroup1 = cut(v1, breaks = quantiles1,
labels = c(1, 2, 3, 4)),
vgroup2 = cut(v2, breaks = quantiles2,
labels = c(1, 2, 3, 4)),
vgroup1 = as.character(vgroup1),
vgroup2 = as.character(vgroup2),
class = paste(vgroup1, "-", vgroup2)) %>%
select(-starts_with("vgroup")) -> df
df
|
这样我们就把数据都分好组了，下面最重要的事情来了？两个变量，那每组应该是什么颜色？不用担心，前人已经把工具开发好了：Bivariate Choropleth Color Generator(https://observablehq.com/@czxa/bivariate-choropleth-color-generator)
这个工具可以让你自由的进行颜色选择、色阶选择，完美！这里我们选择 4 阶的：对应的颜色是：
["#e8e8e8", "#bddede", "#8ed4d4", "#5ac8c8", "#dabdd4", "#bdbdd4", "#8ebdd4", "#5abdc8", "#cc92c1", "#bd92c1", "#8e92c1", "#5a92c1", "#be64ac", "#bd64ac", "#8e64ac", "#5a64ac"]
|
可以用 R 查看这些颜色：
c("#e8e8e8", "#bddede", "#8ed4d4", "#5ac8c8", "#dabdd4", "#bdbdd4", "#8ebdd4", "#5abdc8", "#cc92c1", "#bd92c1", "#8e92c1", "#5a92c1", "#be64ac", "#bd64ac", "#8e64ac", "#5a64ac") %>% scales::show_col()
|
我们使用这些颜色构建一个数据框，这个数据框是用于接下来绘制图例的：
tibble(
class = c("1 - 1","1 - 2","1 - 3","1 - 4",
"2 - 1","2 - 2","2 - 3","2 - 4",
"3 - 1","3 - 2","3 - 3","3 - 4",
"4 - 1","4 - 2","4 - 3","4 - 4"),
fill = c("#e8e8e8", "#bddede", "#8ed4d4", "#5ac8c8",
"#dabdd4", "#bdbdd4", "#8ebdd4", "#5abdc8",
"#cc92c1", "#bd92c1", "#8e92c1", "#5a92c1",
"#be64ac", "#bd64ac", "#8e64ac", "#5a64ac")
) -> bivariate_color_scale
bivariate_color_scale
#> # A tibble: 16 × 2
#> class fill
#> <chr> <chr>
#> 1 1 - 1 #e8e8e8
#> 2 1 - 2 #bddede
#> 3 1 - 3 #8ed4d4
#> 4 1 - 4 #5ac8c8
#> 5 2 - 1 #dabdd4
#> 6 2 - 2 #bdbdd4
#> 7 2 - 3 #8ebdd4
#> 8 2 - 4 #5abdc8
#> 9 3 - 1 #cc92c1
#> 10 3 - 2 #bd92c1
#> 11 3 - 3 #8e92c1
#> 12 3 - 4 #5a92c1
#> 13 4 - 1 #be64ac
#> 14 4 - 2 #bd64ac
#> 15 4 - 3 #8e64ac
#> 16 4 - 4 #5a64ac
|
我们把这个数据框和 df 合并起来：
df %>%
left_join(bivariate_color_scale) -> df
|
再把 df 和地图数据集合并起来：
county_simple %>%
left_join(df, by = c("县代码" = "code", "县" = "name")) -> mapdf
mapdf
|
再读入省级地图数据，这个是覆盖在县级地图上面让地图更好看的：
read_sf("2019行政区划/省.shp") %>%
st_transform(crs = mycrs) %>%
st_simplify(dTolerance = 2000) -> prov
|
生成一个标签数据框，这个数据框是用于给地图添加注解的，因为我想像上面的瑞典地图一样在地图的四角选择四个县通过箭头添加注解文本：
rbind(
slice(subset(mapdf, class == "1 - 1" & str_detect(省, "新疆")), 3),
slice(subset(mapdf, class == "1 - 4" & str_detect(省, "云南")), 3),
slice(subset(mapdf, class == "4 - 1" & str_detect(省, "江苏")), 3),
slice(subset(mapdf, class == "4 - 4" & str_detect(省, "辽宁")), 3)
) %>%
st_centroid() %>%
st_coordinates() %>%
as_tibble() %>%
set_names(c("x", "y")) %>%
mutate(xend = x + c(-100000, 200000, 500000, 500000),
yend = y + c(400000, -800000, 100000, 0),
xla = x + c(-100000, 200000, 900000, 900000),
yla = y + c(700000, -1000000, 100000, 0)) %>%
mutate(label = c("low v1\nlow v2",
"low v1\nhigh v2",
"high v1\nlow v2",
"high v1\nhigh v2")) -> labeldf
labeldf
|
这里为了满足“四角”的要求，我从新疆、云南、江苏和辽宁四个省选取了四个极端组的县进行标注。labeldf 数据框中的 x 和 y 就是选取的四个极端县的质心，同时也是绘图时箭头的箭尾坐标，xend、yend 是箭头坐标，是我试出来的。xla 和 yla 是标签放置的坐标，也是我试出来的。label 则是用于标注的文本，\n 表示换行。
然后我们就可以绘制地图了：
read_sf("海岸线/海岸线.shp") %>%
st_transform(mycrs) -> hax
read_sf("九段线.geojson") %>%
st_transform(mycrs) -> jdx
mapdf %>%
ggplot() +
geom_sf(aes(fill = I(fill)), color = "black", size = 0.01) +
geom_sf(data = prov, fill = NA, color = "white", size = 0.2) +
geom_sf(data = hax, color = "#0055AA", size = 0.5) +
geom_sf(data = jdx, color = "black", size = 0.5) +
theme_ipsum(base_family = "songti") +
geom_curve(data = labeldf,
aes(x = x, y = y, xend = xend, yend = yend),
curvature = 0.1,
arrow = arrow(length = unit(0.01, "npc"))) +
geom_text(data = labeldf,
aes(x = xla, y = yla, label = label),
family = "songti") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
coord_sf(crs = mycrs, xlim = c(-2577078-500000, 2095963+500000)) +
annotation_scale(location = "bl", width_hint = 0.3,
text_family = "songti") +
annotation_north_arrow(location = "tr", which_north = "false",
pad_x = unit(0.5, "cm"),
pad_y = unit(0.5, "cm"),
style = north_arrow_fancy_orienteering(
text_family = "songti"
)) +
labs(title = "中国地图双变量填充") -> map
|
绘制图例需要这样的数据：
bivariate_color_scale %>%
separate(class, into = c("v1", "v2"), sep = " - ") %>%
mutate(v1 = as.integer(v1),
v2 = as.integer(v2)) -> bivariate_color_scale
bivariate_color_scale
|
绘制图例：
library(latex2exp)
ggplot() +
geom_tile(
data = bivariate_color_scale,
mapping = aes(
x = v1,
y = v2,
fill = fill)
) +
scale_fill_identity() +
labs(x = TeX("Higher v1 \\rightarrow"),
y = TeX("Higher v2 \\rightarrow")) +
coord_fixed() +
theme_minimal(base_family = "songti") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title = element_text(size = 10)) -> legend
legend
|
最后就是把这两个图组合在一起了：
library(cowplot)
ggdraw() +
draw_plot(map, 0, 0, 1, 1) +
draw_plot(legend, 0.2, 0.15, 0.22, 0.22)
|
使用类似的方法还可以绘制城市的：
点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制双变量填充中国地图（拼图方法）
评论