附件下载(点击文末的阅读原文即可跳转)：https://rstata./#/brief/course/0e80355d6df245e3aaa51e5860ed9381
为了让大家更好的理解本文的代码，欢迎各位培训班会员参加明晚 8 点的直播课：「如何使用 R 语言绘制双变量填充中国地图？」
在过去的教程中，我一直给大家分享的都是单变量填充地图的绘制方法，但是有时候我们想在一张地图上展示两个变量，这个时候该怎么绘制地图呢？我们可以选择双变量填充地图。如果我们在 Observable 上搜索“Bivariate Choropleth” 就可以看到一大堆使用 D3 绘制的双变量填充地图，非常炫酷，例如：https:///@d3/bivariate-choropleth
如何使用 R 绘制这种双变量填充地图，Timo Grossenbacher 做了一个非常炫酷的示范：
该图的绘制方法可以参考：https:///2019/04/bivariate-maps-with-ggplot2-and-sf/ 学习。不过代码过于复杂，所以我这里对里面的代码进行了简化，并且介绍了如何绘制更多分阶的双变量填充地图（该作者只给出了三阶的画法）。
我暂时也没有时候绘制这种地图的数据，所以下面的演示中我是用随机生成的变量。
准备工作
数据
- 2019行政区划：中国省市区县级矢量地图（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
#> Simple feature collection with 2900 features and 7 fields (with 8 geometries empty)
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: -2577078 ymin: 2093713 xmax: 2095448 ymax: 6388032
#> CRS: +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
#> # A tibble: 2,900 × 8
#> 县代码 县 省代码 省 市代码 市 类型 geometry
#> * <dbl> <chr> <dbl> <chr> <dbl> <chr> <chr> <GEOMETRY [m]>
#> 1 110101 东城区 110000 北京市 110000 北京市 市辖区 POLYGON ((939697.1 486078…
#> 2 110102 西城区 110000 北京市 110000 北京市 市辖区 POLYGON ((937736.1 486023…
#> 3 110105 朝阳区 110000 北京市 110000 北京市 市辖区 MULTIPOLYGON (((943956.6 …
#> 4 110106 丰台区 110000 北京市 110000 北京市 市辖区 POLYGON ((933688.3 485266…
#> 5 110107 石景山区 110000 北京市 110000 北京市 市辖区 POLYGON ((917889.2 486087…
#> 6 110108 海淀区 110000 北京市 110000 北京市 市辖区 POLYGON ((920458.5 487719…
#> 7 110109 门头沟区 110000 北京市 110000 北京市 市辖区 POLYGON ((887385 4873900,…
#> 8 110111 房山区 110000 北京市 110000 北京市 市辖区 POLYGON ((912222.9 484403…
#> 9 110112 通州区 110000 北京市 110000 北京市 市辖区 POLYGON ((956386.6 486952…
#> 10 110113 顺义区 110000 北京市 110000 北京市 市辖区 POLYGON ((972163.9 490020…
#> # … with 2,890 more rows
# plot(county_simple[2])
# 生成示例数据
set.seed(1234)
county_simple %>%
st_drop_geometry() %>%
select(code = 县代码, name = 县) %>%
mutate(v1 = runif(2900, 1, 100),
v2 = runif(2900, 1, 100)) -> df
df
#> # A tibble: 2,900 × 4
#> code name v1 v2
#> <dbl> <chr> <dbl> <dbl>
#> 1 110101 东城区 12.3 80.7
#> 2 110102 西城区 62.6 82.4
#> 3 110105 朝阳区 61.3 34.3
#> 4 110106 丰台区 62.7 4.61
#> 5 110107 石景山区 86.2 38.6
#> 6 110108 海淀区 64.4 63.3
#> 7 110109 门头沟区 1.94 48.2
#> 8 110111 房山区 24.0 72.7
#> 9 110112 通州区 66.9 3.24
#> 10 110113 顺义区 51.9 50.9
#> # … with 2,890 more rows
数据非常简单，就是每个区县对于两个变量的两个值，v1 和 v2 就是我们需要的 “bivariate”。
上面我展示的图片中使用的是三阶，我自然得比他的更复杂点，我用 4 阶！下面我们将变量 v1 和 v2 使用分位数各分成 4 段：
# 对 df 的两个变量进行分组
no_classes <- 4
# 计算 v1 分位数
df %>%
pull(v1) %>%
quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
as.vector() -> quantiles1
quantiles1
#> [1] 1.033839 26.453382 51.048444 74.977486 99.931000
# 计算 v2 分位数
df %>%
pull(v2) %>%
quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
as.vector() -> quantiles2
quantiles2
#> [1] 1.089039 25.761603 49.748025 75.341730 99.959791
对 df 里的 v1 v2 进行分组：
# 对 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
#> # A tibble: 2,900 × 5
#> code name v1 v2 class
#> <dbl> <chr> <dbl> <dbl> <chr>
#> 1 110101 东城区 12.3 80.7 1 - 4
#> 2 110102 西城区 62.6 82.4 3 - 4
#> 3 110105 朝阳区 61.3 34.3 3 - 2
#> 4 110106 丰台区 62.7 4.61 3 - 1
#> 5 110107 石景山区 86.2 38.6 4 - 2
#> 6 110108 海淀区 64.4 63.3 3 - 3
#> 7 110109 门头沟区 1.94 48.2 1 - 2
#> 8 110111 房山区 24.0 72.7 1 - 3
#> 9 110112 通州区 66.9 3.24 3 - 1
#> 10 110113 顺义区 51.9 50.9 3 - 3
#> # … with 2,890 more rows
这样我们就把数据都分好组了，下面最重要的事情来了？两个变量，那每组应该是什么颜色？不用担心，前人已经把工具开发好了：Bivariate Choropleth Color Generator(https:///@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 和地图数据集合并起来：
# 合并地图数据和 df
county_simple %>%
left_join(df, by = c("县代码" = "code", "县" = "name")) -> mapdf
mapdf
#> Simple feature collection with 2900 features and 11 fields (with 8 geometries empty)
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: -2577078 ymin: 2093713 xmax: 2095448 ymax: 6388032
#> CRS: +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
#> # A tibble: 2,900 × 12
#> 县代码 县 省代码 省 市代码 市 类型 geometry v1
#> <dbl> <chr> <dbl> <chr> <dbl> <chr> <chr> <GEOMETRY [m]> <dbl>
#> 1 110101 东城区 110000 北京… 110000 北京… 市辖… POLYGON ((939697.1 48607… 12.3
#> 2 110102 西城区 110000 北京… 110000 北京… 市辖… POLYGON ((937736.1 48602… 62.6
#> 3 110105 朝阳区 110000 北京… 110000 北京… 市辖… MULTIPOLYGON (((943956.6… 61.3
#> 4 110106 丰台区 110000 北京… 110000 北京… 市辖… POLYGON ((933688.3 48526… 62.7
#> 5 110107 石景… 110000 北京… 110000 北京… 市辖… POLYGON ((917889.2 48608… 86.2
#> 6 110108 海淀区 110000 北京… 110000 北京… 市辖… POLYGON ((920458.5 48771… 64.4
#> 7 110109 门头… 110000 北京… 110000 北京… 市辖… POLYGON ((887385 4873900… 1.94
#> 8 110111 房山区 110000 北京… 110000 北京… 市辖… POLYGON ((912222.9 48440… 24.0
#> 9 110112 通州区 110000 北京… 110000 北京… 市辖… POLYGON ((956386.6 48695… 66.9
#> 10 110113 顺义区 110000 北京… 110000 北京… 市辖… POLYGON ((972163.9 49002… 51.9
#> # … with 2,890 more rows, and 3 more variables: v2 <dbl>, class <chr>,
#> # fill <chr>
再读入省级地图数据，这个是覆盖在县级地图上面让地图更好看的：
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
#> # A tibble: 4 × 7
#> x y xend yend xla yla label
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 -1826165. 5463075. -1926165. 5863075. -1926165. 6163075. "low v1\nlow v2"
#> 2 -134148. 2930179. 65852. 2130179. 65852. 1930179. "low v1\nhigh v2"
#> 3 1389428. 3988377. 1889428. 4088377. 2289428. 4088377. "high v1\nlow v2"
#> 4 1491219. 5171978. 1991219. 5171978. 2391219. 5171978. "high v1\nhigh v2"
这里为了满足“四角”的要求，我从新疆、云南、江苏和辽宁四个省选取了四个极端组的县进行标注。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,
# bar_cols = c("#be64ac", "#5ac8c8"),
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(
# line_col = "#be64ac",
# text_col = "#5ac8c8",
# fill = c("#be64ac", "#5ac8c8"),
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
#> # A tibble: 16 × 3
#> v1 v2 fill
#> <int> <int> <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
绘制图例：
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)
使用类似的方法还可以绘制城市的：
