使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）

之前也给大家分享过类似的课程，不过感觉方法不够范式，最近终于开发出了一套范式的方法，使用该方法可以非常轻松的绘制历年的中国各省市区县的地图。附件中我提供了 1949～2021 年每年的省市区县地图数据，包含 mini 和 long 两套：

今天的课程中我们将一起学习这两套数据的使用。

数据介绍

该套绘图数据是经过预先编辑的（添加了秦岭淮河线、胡焕庸线等元素）。在之前的课程「使用 R 语言操作地理矢量数据、使用 Stata 绘制精美地图」（可以从平台上搜索下）中我向大家介绍过这种使用 R 语言进行地理矢量数据编辑的方法。

附件中提供了三个文件夹：

• provmapdata
• citymapdata
• countymapdata

每个文件夹中都包含了 longshp 和 minishp 两个文件夹，这两个文件夹下面又分别包含了 1949～2021 年每年的省市区县地图数据，long 表示长版，mini 表示小地图版本。

每份数据又包含两种 shp 数据，例如 1992 年 long 版本的：

• chinaprov1992long_line：线条元素，例如省界、国界线、九段线、海岸线、秦岭淮河线、胡焕庸线等。
• chinaprov1992long：多边形元素，也就是各省的区域范围。

小地图版本

1. 在地图上添加散点图

首先加载所需的 R 包：

library(tidyverse)
library(sf)

中国地图通常使用这样的坐标系：

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("provmapdata/minishp/chinaprov2019mini/chinaprov2019mini.shp") %>%
  filter(!is.na(省代码)) -> provmap

provmap

#> Simple feature collection with 35 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -2625586 ymin: 1876787 xmax: 2962719 ymax: 5921583
#> 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: 35 × 5
#>    省代码 省           class objid                                      geometry
#>  *  <dbl> <chr>        <chr> <dbl>                            <MULTIPOLYGON [m]>
#>  1 110000 北京市       省份      1 (((964331.3 4474014, 964405 4474018, 964455.…
#>  2 120000 天津市       省份      2 (((1038620 4393376, 1039708 4392960, 1040590…
#>  3 130000 河北省       省份      3 (((1173559 4280994, 1173178 4280926, 1172773…
#>  4 140000 山西省       省份      4 (((758535.4 4417694, 758671.5 4417536, 75873…
#>  5 150000 内蒙古自治区 省份      5 (((1134469 5872436, 1134764 5871947, 1134790…
#>  6 150000 中朝共有     省份      6 (((1616450 4451680, 1616434 4451670, 1616458…
#>  7 210000 辽宁省       省份      7 (((1361070 4282431, 1361146 4282390, 1361289…
#>  8 220000 吉林省       省份      8 (((2052712 4841890, 2052677 4841862, 2052498…
#>  9 230000 黑龙江省     省份      9 (((1259864 5917473, 1260015 5917498, 1260301…
#> 10 310000 上海市       省份     10 (((1529226 3496447, 1529303 3496452, 1529371…
#> # ℹ 25 more rows

# 线条
read_sf("provmapdata/minishp/chinaprov2019mini/chinaprov2019mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> provlinemap

provlinemap

#> Simple feature collection with 3 features and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -2625102 ymin: 1825139 xmax: 2983002 ymax: 5921059
#> 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: 3 × 2
#>   class                                                                 geometry
#>   <chr>                                                    <MULTILINESTRING [m]>
#> 1 九段线     ((-2608129 4545302, -2625102 4558824, -2622683 4591937, -2615168 4…
#> 2 海岸线     ((953456.7 2016301, 951363.1 2015404, 949199.8 2015293, 946810 201…
#> 3 小地图框格 ((2160000 1825139, 2983002 1825139, 2983002 2944032, 2160000 29440…

在地图上添加散点图，这里以瞪羚企业分布为例：

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

然后就可以绘图了：

代码里面用到的 cnfont 是预先配置的字体，可以学习平台上的系列课程「R 语言数据科学」第一课时了解如何配置。

library(ggspatial)
library(ggnewscale)
ggplot(provmap) +
  geom_sf(fill = "white", color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(label = 省), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 3) +
  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)) +
  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)

2. 绘制填充地图

下面我们再计算各城市的瞪羚企业数量，虽然绘制市级地图，但是我们仍然使用省级的线条（突出省份区域）：

# 读取 2021 年城市地图
read_sf("citymapdata/minishp/chinacity2021mini/chinacity2021mini.shp") %>%
  filter(!is.na(省代码)) -> citymap

# 在市级地图上强调省级边界
read_sf("provmapdata/minishp/chinaprov2021mini/chinaprov2021mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格", "省份")) %>%
  select(class) -> provlinemap

统计每个城市的公司数量：

pointdf %>%
  count(市, 市代码) -> citydf

citydf

和地图数据合并：

citymap %>%
  left_join(citydf) %>%
  mutate(n = if_else(is.na(n), 0, n)) -> citymap2

然后就可以绘图了：

citymap2 %>%
  mutate(n = n + 1) %>%
  ggplot() +
  geom_sf(aes(fill = n), color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black",
               "省份" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3,
               "省份" = 0.3)
  ) +
  scico::scale_fill_scico(palette = "acton",
                          name = "公司数量",
                          trans = "log10") +
  new_scale_color() +
  geom_sf(data = pointdfsfall, aes(color = 省),
          size = 0.1, show.legend = F) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(label = 省), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 3) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  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("pic2.png", width = 10, height = 8.5, device = png)

使用连续变量绘制填充地图的时候，默认情况下会使用渐变色，我们也可以先把连续变量切分成分类变量（这里是有序因子变量）再绘图：

# 使用分段填色
citymap2$n %>%
  quantile(probs = 1:10/10, digits = 1) %>%
  as.integer() %>%
  unique() -> cutlist

cutlist

citymap2 %>%
  mutate(group = cut(n, breaks = unique(cutlist),
                     include.lowest = T,
                     labels = c("<= 1", "1～3", "3～9", "9～18",
                                "19～39", "39～92", ">= 92"))) -> citymap3

citymap3 %>%
  ggplot() +
  geom_sf(aes(fill = group), color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black",
               "省份" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3,
               "省份" = 0.3)
  ) +
  scico::scale_fill_scico_d(palette = "acton",
                            name = "公司数量",
                            end = 0.95) +
  new_scale_color() +
  geom_sf(data = pointdfsfall, aes(color = 省),
          size = 0.1, show.legend = F) +
  scale_color_manual(values = paletteer::paletteer_d("ggsci::default_igv", n = 32)) +
  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
    )
  ) +
  theme(legend.position = c(0.12, 0.2)) +
  guides(fill = guide_legend(ncol = 2)) -> p3

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

3. 使用栅格数据绘制地图

使用栅格数据绘制地图的场景也很常见，首先我们看一下高分辨率的数据如何绘图。

高分辨率的栅格数据可以先转换成密集的散点数据：

library(raster)
raster("cn2022.tif") -> cn2022
cn2022 %>%
  aggregate(fact = 10, fun = mean) %>%
  rasterToPoints(spatial = TRUE) %>%
  st_as_sf() -> cn2022points

然后绘图方法就和散点图一样了：

cn2022points %>%
  st_transform(mycrs) -> cn2022points

# 提取小地图内的点
cn2022points %>%
  st_intersection(small_bbox) -> cn2022points_small

# 移动小地图内的点到恰当位置：
cn2022points_small %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs) -> cn2022points_small

# 合并两个部分
bind_rows(cn2022points, cn2022points_small) -> cn2022points_all

# 绘图
read_sf("provmapdata/minishp/chinaprov2019mini/chinaprov2019mini_line.shp") %>%
  filter(class %in% c("九段线", "海岸线", "小地图框格")) %>%
  select(class) -> provlinemap

ggplot(provmap) +
  geom_sf(data = cn2022points_all, aes(color = cn2022),
          size = 0.01, show.legend = F) +
  scico::scale_color_scico(
    palette = "lajolla",
    direction = -1,
    trans = "log10"
  ) +
  new_scale_color() +
  geom_sf(fill = NA, color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(label = 省), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 3) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "2022 年中国夜间灯光地图",
       subtitle = "数据处理&绘制：微信公众号 RStata",
       caption = "数据来源：An extended time-series (2000-2020) of global NPP-VIIRS-like nighttime light data - Harvard Dataverse\nhttps://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/YGIVCD") +
  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
    )
  ) -> p4

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

如果是低分辨率的栅格数据，使用上面的方法就不太好看了，可以考虑转换成 polygon 绘图：

# 低分辨率
raster("NH3_em_anthro_2015_sector_ENE.tif") %>%
  rasterToPolygons() %>%
  st_as_sf() -> cn2022polygons

cn2022polygons %>%
  st_transform(mycrs) -> cn2022polygons

# 提取小地图内的点
cn2022polygons %>%
  st_intersection(small_bbox) -> cn2022polygons_small

# 移动小地图内的点到恰当位置：
cn2022polygons_small %>%
  mutate(geometry = geometry * 0.5 + c(2100000, 1665139)) %>%
  sf::st_set_crs(mycrs) -> cn2022polygons_small

# 合并两个部分
bind_rows(cn2022polygons, cn2022polygons_small) -> cn2022polygons_all

# 提取中国范围的
cn2022polygons_all %>%
  st_intersection(provmap) -> cn2022polygons_all

cn2022polygons_all$NH3_em_anthro_2015_sector_ENE -> numlist
range(log10(numlist), na.rm = T, finite = T) -> rlist
10^(mean(rlist)) -> median
10^(rlist[1] + 0.02 * (rlist[2] - rlist[1])) -> min
10^(rlist[2] - 0.02 * (rlist[2] - rlist[1])) -> max

min;median;max

ggplot(provmap) +
  geom_sf(data = cn2022polygons_all,
          aes(fill = NH3_em_anthro_2015_sector_ENE),
          linewidth = 0.001) +
  scico::scale_fill_scico(
    palette = "imola", trans = "log10",
    direction = -1,
    breaks = c(min, median, max),
    labels = c("低", latex2exp::TeX("$\\leftarrow$ NH3 排放($kg/m^2/yr$) $\\rightarrow$"), "高"),
    name = "2015 年各地能源开采部门",
    guide = guide_colorbar(
      direction = "horizontal",
      barheight = unit(3, units = "mm"),
      barwidth = unit(60, units = "mm"),
      draw.ulim = FALSE,
      ticks.colour = "transparent",
      title.position = 'top',
      title.hjust = 0.5,
      label.hjust = 0.5
    )
  ) +
  new_scale_color() +
  geom_sf(fill = NA, color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(label = 省), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 3) +
  scale_color_manual(
    values = c("九段线" = "#A29AC4",
               "海岸线" = "#0055AA",
               "小地图框格" = "black")
  ) +
  scale_linewidth_manual(
    values = c("九段线" = 0.6,
               "海岸线" = 0.3,
               "小地图框格" = 0.3)
  ) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = latex2exp::TeX("各地区 $NH_{3}$ 排放分布(单位：$kg/m^2/yr$)"),
       subtitle = "数据处理&绘制：微信公众号 RStata",
       caption = "数据来源：中国历史空气污染物排放数据集（1990-2015 逐年）\n<https://zenodo.org/record/4741285>") +
  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
    )
  ) +
  theme(legend.position = c(0.15, 0.12)) -> p5

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

长版地图

长版地图数据的使用更简单，这里仅仅演示如何在地图上添加散点图：

# 读取小地图版本的中国省级地图
read_sf("provmapdata/longshp/chinaprov2019long/chinaprov2019long.shp") %>%
  filter(!is.na(省代码)) -> provmap

# 线条
read_sf("provmapdata/longshp/chinaprov2019long/chinaprov2019long_line.shp") %>%
  filter(class %in% c("九段线", "海岸线")) %>%
  select(class) -> provlinemap

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

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

# 范围
st_bbox(provmap)

# 绘图
library(ggspatial)
library(ggnewscale)
ggplot(provmap) +
  geom_sf(fill = "white", color = "gray", linewidth = 0.01) +
  geom_sf(data = provlinemap,
          aes(color = class, linewidth = class),
          show.legend = F) +
  stat_sf_coordinates(data = provmap,
                      geom = "text", color = "gray",
                      aes(label = 省), family = cnfont,
                      fun.geometry = st_point_on_surface,
                      size = 3) +
  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)) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  coord_sf(xlim = c(-2625585.8, 2206964.7) * 1.2) +
  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
    )
  ) -> p6

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

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

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言绘制历年中国省市区县地图（小地图版本+长版）