使用 R 语言绘制双变量填充中国地图（插件方法）

在之前「使用 R 语言绘制历年中国各省市区县地图（小地图版本 + 长版）」课程的基础上，我们今天再来学习下如何绘制双变量填充地图，也同样包含两种版本：

小地图版本

在本课程中我们会讲解如何生成双变量色块图例的数据，但是实际上大家之后绘图的时候并不用每次都生成，而是直接读取附件中生成好的 shp 文件即可。

首先加载所需 R 包，读取地图数据和瞪羚数据：

library(tidyverse)
library(sf)
library(raster)
library(ggspatial)
library(ggnewscale)
library(ggtext)

# 中国地图通常使用这样的坐标系
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 读取小地图版本的中国城市地图
read_sf("chinacity2019mini/chinacity2019mini.shp") %>%
  filter(!is.na(省代码)) -> citymap

# 线条
read_sf("chinacity2019mini/chinacity2019mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> citylinemap

# 在地图上添加散点图，这里以瞪羚企业分布为例：
readxl::read_xlsx("瞪羚、独角兽、创新型企业经纬度数据.xlsx") %>%
  filter(!is.na(经度)) -> pointdf

# 转换成 sf 对象
pointdf %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326) %>%
  st_transform(mycrs) -> pointdfsf

# 小地图上的点
# 小地图的范围
small_bbox <- st_bbox(c(xmin = 120000,
                        xmax = 1766004.1,
                        ymax = 2557786.0,
                        ymin = 320000),
                      crs = st_crs(mycrs)) %>%
  st_as_sfc()

# 提取这个范围的点
pointdfsf %>%
  st_intersection(small_bbox) -> pointdfsf_small

# 把这些点移动到小地图的位置
pointdfsf_small %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs) -> pointdfsf_small

# 合并两部分
bind_rows(pointdfsf, pointdfsf_small) -> pointdfsfall

pointdfsfall

#> Simple feature collection with 40820 features and 9 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -1489994 ymin: 2072412 xmax: 2720657 ymax: 5261103
#> Projected CRS: +proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#> # A tibble: 40,820 × 10
#>    公司编号             公司名称 省    省代码 市    市代码 县    县代码 认定分类
#>  * <chr>                <chr>    <chr>  <dbl> <chr>  <dbl> <chr>  <dbl> <chr>
#>  1 00073fd4d5204e1db2c… 北京领…  北京… 110000 北京… 110000 东城… 110101 ,瞪羚企…
#>  2 0007fc40037c446c8c7… 北京华…  北京… 110000 北京… 110000 密云… 110118 ,瞪羚企…
#>  3 000e2dda3c0f45c8b99… 广州穗…  广东… 440000 广州… 440100 黄埔… 440112 ,潜在瞪…
#>  4 0013abf442d84d9caa4… 武汉Ｔ…  湖北… 420000 武汉… 420100 江夏… 420115 ,瞪羚企…
#>  5 001434ab5e154531abd… 北京启…  北京… 110000 北京… 110000 昌平… 110114 ,瞪羚企…
#>  6 00155fb6ee9544c7af9… 中通铁…  北京… 110000 北京… 110000 西城… 110102 ,瞪羚企…
#>  7 0015cb8bb70e4b60b77… 广东国…  广东… 440000 广州… 440100 天河… 440106 ,瞪羚企…
#>  8 0017ac5d15584cfca9a… 思科系…  北京… 110000 北京… 110000 海淀… 110108 ,瞪羚企…
#>  9 00190b6fdb964662b72… 合肥英…  安徽… 340000 合肥… 340100 蜀山… 340104 ,瞪羚企…
#> 10 001f5df0e4ad43e5b47… 北京好…  北京… 110000 北京… 110000 昌平… 110114 ,瞪羚企…
#> # ℹ 40,810 more rows
#> # ℹ 1 more variable: geometry <POINT [m]>

citymap 的坐标范围是这样的：

citymap %>%
  st_bbox()

#>     xmin     ymin     xmax     ymax
#> -2625586  1876787  2962719  5921583

后面的图例位置就是根据这个范围尝试出来的。

这里我们以 4x4 的双变量填充地图为例，使用下面的代码即可得到所有的 16 种类别标签：

crossing(1:4, 1:4) %>%
  set_names("v1", "v2") %>%
  unite("v", c("v1", "v2"), sep = " - ", remove = T) -> groupclassdf

groupclassdf$v

#>  [1] "1 - 1" "1 - 2" "1 - 3" "1 - 4" "2 - 1" "2 - 2" "2 - 3" "2 - 4" "3 - 1"
#> [10] "3 - 2" "3 - 3" "3 - 4" "4 - 1" "4 - 2" "4 - 3" "4 - 4"

然后就可以构建数据框了：

tibble(
  groupclass = groupclassdf$v,
  fill = c("#e8e8e8", "#bddede", "#8ed4d4", "#5ac8c8",
           "#dabdd4", "#bdbdd4", "#8ebdd4", "#5abdc8",
           "#cc92c1", "#bd92c1", "#8e92c1", "#5a92c1",
           "#be64ac", "#bd64ac", "#8e64ac", "#5a64ac")
) -> bivariate_color_scale

这里面的 15 种渐变色是使用这个网页应用生成的：https://observablehq.com/@czxa/bivariate-choropleth-color-generator

该应用可以生成 2～20 阶的双变量填充颜色，并且颜色可以自定义。

然后我们需要生成一个 4x4 的 multipolygon 类的 sf 对象，这里我是先创建的栅格数据，然后再把栅格数据转换成 sf 对象。

首先看一下直接转换得到的结果：

x0 <- -2525586
y0 <- 2076787

raster(resolution = c(187500, 187500), nrows = 4, ncols = 4,
       crs = mycrs, xmn = x0, xmx = x0 + 750000,
       ymn = y0, ymx = y0 + 750000) %>%
  rasterToPolygons() %>%
  st_as_sf() %>%
  select(-layer) %>%
  mutate(ID = row_number()) %>%
  bind_cols(bivariate_color_scale) %>%
  ggplot() +
  geom_sf(aes(fill = I(fill))) +
  geom_sf_text(aes(label = paste0(groupclass, "\n", "(", ID, ")")))

这样的结果并不是我们想要的，应该把这个图例逆时针旋转 90 度才对，所以我们需要改变下 polygon 的顺序，注意到现在色块现在的编号顺序是 1-16，如果逆时针旋转 90 度，对应的编号应该是 4、8、12、16、3、7、11、15、2、6、10、14、1、5、9、13。我们可以通过编程生成这样的编号顺序的：

matrix(1:16, nrow = 4, byrow = T) %>%
  raster() %>%
  t() %>%
  flip() %>%
  .[] -> sortid

sortid

#>  [1]  4  8 12 16  3  7 11 15  2  6 10 14  1  5  9 13

raster(resolution = c(187500, 187500), nrows = 4, ncols = 4,
       crs = mycrs, xmn = x0, xmx = x0 + 750000,
       ymn = y0, ymx = y0 + 750000) %>%
  rasterToPolygons() %>%
  st_as_sf() %>%
  select(-layer) %>%
  mutate(ID = sortid) %>%
  arrange(ID) %>%
  bind_cols(bivariate_color_scale) -> bivariate_color_scale_sf

bivariate_color_scale_sf %>%
  ggplot() +
  geom_sf(aes(fill = I(fill))) +
  geom_sf_text(aes(label = groupclass))

这样的结果就正确了。

另外我们还需要给图例的标签生成位置数据：

# 标签的位置
x0 + 2 * 187500
y0 + 2 * 187500

library(latex2exp)
tribble(
  ~x, ~y, ~label, ~angle,
  x0 + 2 * 187500, y0 - 90000, "v1 →", 0,
  x0 - 90000, y0 + 2 * 187500, "v2 →", 90,
) -> labeldf

为了绘制双变量填充地图，我们需要准备两个连续变量，这里我们统计每个城市的瞪羚企业和其他企业的数量：

# 统计每个城市的瞪羚企业数量和独角兽企业数量
pointdf %>%
  mutate(认定分类 = case_when(
    str_detect(认定分类, "瞪羚") ~ "瞪羚企业",
    T ~ "独角兽和创新企业"
  )) %>%
  count(市, 市代码, 认定分类) %>%
  spread(认定分类, n, fill = 0) -> countdf

countdf
#> # A tibble: 232 × 4
#>    市     市代码 瞪羚企业 独角兽和创新企业
#>    <chr>   <dbl>    <dbl>            <dbl>
#>  1 安康市 610900       16                0
#>  2 安顺市 520400        1                0
#>  3 安阳市 410500       10                0
#>  4 鞍山市 210300       35                1
#>  5 白银市 620400        5                0
#>  6 百色市 451000       10                0
#>  7 蚌埠市 340300       59                0
#>  8 包头市 150200       27                1
#>  9 宝鸡市 610300       56                0
#> 10 保定市 130600       15                2
#> # ℹ 222 more rows

然后再把每个变量分成四组，这里使用分位数进行分组：

# 计算两个变量的分位数
no_classes <- 4
# 计算 v1 分位数
countdf %>%
  pull(瞪羚企业) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
  as.vector() -> quantiles1
quantiles1

#> [1]     0.00     4.00    15.50    51.25 16433.00

# 计算 v2 分位数
countdf %>%
  filter(独角兽和创新企业 != 0) %>%
  pull(独角兽和创新企业) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
  as.vector() -> quantiles2
quantiles2

#> [1]   1.00   1.00   2.00  10.25 341.00

quantiles2[1] <- 0
quantiles2

#> [1]   0.00   1.00   2.00  10.25 341.00

# 对 countdf 里的两个进行分组
countdf %>%
  mutate(vgroup1 = cut(瞪羚企业, breaks = quantiles1,
                       labels = c(1, 2, 3, 4),
                       include.lowest = T),
         vgroup2 = cut(独角兽和创新企业, breaks = quantiles2,
                       labels = c(1, 2, 3, 4),
                       include.lowest = T),
         vgroup1 = as.character(vgroup1),
         vgroup2 = as.character(vgroup2),
         groupclass = paste(vgroup1, "-", vgroup2)) %>%
  select(-starts_with("vgroup")) -> countdf

countdf

#> # A tibble: 232 × 5
#>    市     市代码 瞪羚企业 独角兽和创新企业 groupclass
#>    <chr>   <dbl>    <dbl>            <dbl> <chr>
#>  1 安康市 610900       16                0 3 - 1
#>  2 安顺市 520400        1                0 1 - 1
#>  3 安阳市 410500       10                0 2 - 1
#>  4 鞍山市 210300       35                1 3 - 1
#>  5 白银市 620400        5                0 2 - 1
#>  6 百色市 451000       10                0 2 - 1
#>  7 蚌埠市 340300       59                0 4 - 1
#>  8 包头市 150200       27                1 3 - 1
#>  9 宝鸡市 610300       56                0 4 - 1
#> 10 保定市 130600       15                2 2 - 2
#> # ℹ 222 more rows

然后再和地图数据以及 bi—scale 数据匹配：

citymap %>%
  left_join(countdf) %>%
  mutate(groupclass = if_else(is.na(groupclass), "1 - 1", groupclass)) %>%
  left_join(bivariate_color_scale) -> citymap2

label 数据也需要根据自己的需要设置下：

labeldf %>%
  mutate(label = c("瞪羚企业更多 →",
                   "独角兽、创新企业更多 →")) -> labeldf2

然后就可以绘图了：

ggplot(citymap2) +
  geom_sf(aes(fill = I(fill)), color = "gray", linewidth = 0.01) +
  geom_sf(data = citylinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  new_scale_color() +
  geom_sf(data = pointdfsfall, aes(color = 省),
          size = 0.5, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  geom_sf(data = bivariate_color_scale_sf, aes(fill = I(fill)),
          linewidth = 0.01) +
  geom_richtext(data = labeldf2, aes(x, y, label = label, angle = angle),
                fill = NA, label.color = NA, family = cnfont, size = 3) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "瞪羚、独角兽、创新型企业城市分布",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：瞪羚云网站") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p1

ggsave("pic1.png", width = 10, height = 8.5, device = png)

最后为了方便使用，保存为 shp 文件，然后 label 数据也保存成了 xlsx 文件，这样大家之后直接读取就可以使用了：

# 保存结果
dir.create("biscale4x4")
bivariate_color_scale_sf %>%
  st_make_valid() %>%
  st_collection_extract("POLYGON") %>%
  st_write(paste0("biscale4x4/biscale4x4.shp"), layer_options = "ENCODING=UTF-8", delete_layer = TRUE, layer = "MULTIPOLYGON")

labeldf %>%
  writexl::write_xlsx("biscale4x4/biscale4x4_label.xlsx")

附件中还提供了 2x2、3x3、5x5、8x8 的生成方法和结果。

长版地图

方法类似：

library(tidyverse)
library(sf)
library(raster)
library(ggspatial)
library(ggnewscale)
library(ggtext)

# 中国地图通常使用这样的坐标系
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# 读取长版的中国城市地图
read_sf("chinacity2019long/chinacity2019long.shp") %>%
  filter(!is.na(省代码)) -> citymap

# 线条
read_sf("chinacity2019long/chinacity2019long_line.shp") %>%
  filter(class %in% c("九段线", "海岸线")) %>%
  select(class) -> citylinemap

# 在地图上添加散点图，这里以瞪羚企业分布为例：
readxl::read_xlsx("瞪羚、独角兽、创新型企业经纬度数据.xlsx") %>%
  filter(!is.na(经度)) -> pointdf

# 转换成 sf 对象
pointdf %>%
  st_as_sf(coords = c("经度", "纬度"), crs = 4326) %>%
  st_transform(mycrs) -> pointdfsf

# 城市地图的范围
citymap %>%
  st_bbox()

#>       xmin       ymin       xmax       ymax
#> -2625585.8   667046.3  2206964.7  5921583.0

crossing(1:4, 1:4) %>%
  set_names("v1", "v2") %>%
  unite("v", c("v1", "v2"), sep = " - ", remove = T) -> groupclassdf

groupclassdf$v

#>  [1] "1 - 1" "1 - 2" "1 - 3" "1 - 4" "2 - 1" "2 - 2" "2 - 3" "2 - 4" "3 - 1"
#> [10] "3 - 2" "3 - 3" "3 - 4" "4 - 1" "4 - 2" "4 - 3" "4 - 4"

# 4X4
tibble(
  groupclass = groupclassdf$v,
  fill = c("#e8e8e8", "#bddede", "#8ed4d4", "#5ac8c8",
           "#dabdd4", "#bdbdd4", "#8ebdd4", "#5abdc8",
           "#cc92c1", "#bd92c1", "#8e92c1", "#5a92c1",
           "#be64ac", "#bd64ac", "#8e64ac", "#5a64ac")
) -> bivariate_color_scale

x0 <- -2425586
y0 <-   767047

matrix(1:16, nrow = 4, byrow = T) %>%
  raster() %>%
  t() %>%
  flip() %>%
  .[] -> sortid

raster(resolution = c(300000, 300000), nrows = 4, ncols = 4,
       crs = mycrs, xmn = x0, xmx = x0 + 1200000,
       ymn = y0, ymx = y0 + 1200000) %>%
  rasterToPolygons() %>%
  st_as_sf() %>%
  select(-layer) %>%
  mutate(ID = sortid) %>%
  arrange(ID) %>%
  bind_cols(bivariate_color_scale) -> bivariate_color_scale_sf

bivariate_color_scale_sf %>%
  ggplot() +
  geom_sf(aes(fill = I(fill))) +
  geom_sf_text(aes(label = groupclass))

# 标签的位置
x0 + 2 * 300000
y0 + 2 * 300000

library(latex2exp)
tribble(
  ~x, ~y, ~label, ~angle,
  x0 + 2 * 300000, y0 - 120000, "v1 →", 0,
  x0 - 120000, y0 + 2 * 300000, "v2 →", 90,
) -> labeldf

# 统计每个城市的瞪羚企业数量和独角兽企业数量
pointdf %>%
  mutate(认定分类 = case_when(
    str_detect(认定分类, "瞪羚") ~ "瞪羚企业",
    T ~ "独角兽和创新企业"
  )) %>%
  count(市, 市代码, 认定分类) %>%
  spread(认定分类, n, fill = 0) -> countdf

countdf

# 计算两个变量的分位数
no_classes <- 4
# 计算 v1 分位数
countdf %>%
  pull(瞪羚企业) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
  as.vector() -> quantiles1
quantiles1

# 计算 v2 分位数
countdf %>%
  filter(独角兽和创新企业 != 0) %>%
  pull(独角兽和创新企业) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1)) %>%
  as.vector() -> quantiles2
quantiles2

quantiles2[1] <- 0
quantiles2

# 对 countdf 里的两个进行分组
countdf %>%
  mutate(vgroup1 = cut(瞪羚企业, breaks = quantiles1,
                       labels = c(1, 2, 3, 4),
                       include.lowest = T),
         vgroup2 = cut(独角兽和创新企业, breaks = quantiles2,
                       labels = c(1, 2, 3, 4),
                       include.lowest = T),
         vgroup1 = as.character(vgroup1),
         vgroup2 = as.character(vgroup2),
         groupclass = paste(vgroup1, "-", vgroup2)) %>%
  select(-starts_with("vgroup")) -> countdf

countdf

# 和 citymap 数据匹配
citymap %>%
  left_join(countdf) %>%
  mutate(groupclass = if_else(is.na(groupclass), "1 - 1", groupclass)) %>%
  left_join(bivariate_color_scale) -> citymap2

labeldf %>%
  mutate(label = c("瞪羚企业更多 &rarr;",
                   "独角兽、创新企业更多 &rarr;")) -> labeldf2

ggplot(citymap2) +
  geom_sf(aes(fill = I(fill)), color = "gray", linewidth = 0.01) +
  geom_sf(data = citylinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3)
  ) +
  new_scale_color() +
  geom_sf(data = pointdfsf, aes(color = 省),
          size = 0.5, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  geom_sf(data = bivariate_color_scale_sf, aes(fill = I(fill)),
          linewidth = 0.01) +
  geom_richtext(data = labeldf2, aes(x, y, label = label, angle = angle),
                fill = NA, label.color = NA, family = cnfont, size = 3) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "瞪羚、独角兽、创新型企业城市分布",
       subtitle = "数据爬取&绘制：微信公众号 RStata",
       caption = "数据来源：瞪羚云网站") +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p2

ggsave("pic4x4_long.png", width = 9, height = 9, device = png)

# 保存结果
dir.create("biscale4x4_long")
bivariate_color_scale_sf %>%
  st_make_valid() %>%
  st_collection_extract("POLYGON") %>%
  st_write(paste0("biscale4x4_long/biscale4x4_long.shp"), layer_options = "ENCODING=UTF-8", delete_layer = TRUE, layer = "MULTIPOLYGON")

labeldf %>%
  writexl::write_xlsx("biscale4x4_long/biscale4x4_long_label.xlsx")

同样，附件中也提供了其他阶数的结果。

更多关于地图绘制的内容可以学习下平台上的系列课程「使用 R 语言进行地理计算」。

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制双变量填充中国地图（插件方法）