R 语言：如何为区县名称添加行政区划代码

在之前的课程「R 语言：如何整理 2022 年县域统计年鉴：caj 文件转 pdf、文本识别与数据清洗」中我们讲解了如何从 caj 文件中提取表格数据的方法，今天我们再来学习下如何根据区县名称匹配行政区划代码，另外在该过程中还可以检查区县名称的识别错误。最后我们再使用整理得到的数据绘制一幅区县地图。

首先我们使用上次课的代码处理“整理结果3.xlsx”：

# 处理 “整理结果3.xlsx”
library(tidyverse)
readxl::read_xlsx("整理结果3.xlsx", col_names = LETTERS[1:12]) -> df 
df %>% 
  fill(A) %>% 
  dplyr::filter(!(is.na(D) & is.na(E) & is.na(F) 
                  & is.na(G) & is.na(H) & is.na(I)
                  & is.na(J) & is.na(K) & is.na(L))) %>% 
  mutate(B = str_remove_all(B, " ")) %>% 
  mutate(B = str_remove_all(B, ","),
         B = str_remove_all(B, "-"),
         B = str_remove_all(B, "\\."),
         B = str_remove_all(B, "~"),
         B = str_remove_all(B, "・"),
         B = str_remove_all(B, "，"),
         B = str_remove_all(B, "、")) %>% 
  dplyr::filter(!is.na(B) & B != "一、基本情况行政区域面积" & B != "一基本情况行政区域面积") %>% 
  mutate(B = if_else(B %in% c("提供住宿的民政服务机构床位数",
                              "提供住宿的民谢艮务机构床位数",
                              "提供住宿的民政0艮务机构床位数"),
                     "提供住宿的民政服务机构床位数", B),
         B = if_else(B %in% c("提供住宿的民呦艮务机构",
                              "提供住宿的民政^务机构",
                              "提供住宿的民斑艮务机构",
                              "提供住宿的民班艮务机构"),
                     "提供住宿的民政服务机构", B)) %>% 
  mutate(z = if_else(B == "指标", row.names(.), "")) %>% 
  mutate(z = as.numeric(z)) %>% 
  fill(z) %>% 
  select(-C) %>% 
  gather(D:L, key = "variable", value = "value")  %>% 
  spread(B, value) %>% 
  select(z, A, 指标, everything()) %>% 
  dplyr::filter(!is.na(指标)) %>% 
  select(-variable) %>% 
  type_convert() %>% 
  rename(省 = A, 县 = 指标) %>% 
  select(省, 县, 行政区域面积, 乡, 镇, 街道办事处, 户籍人口, 
         地区生产总值, 第一产业增加值, 第二产业增加值, 第三产业增加值, 
         地方一般公共预算收入, 地方一般公共预算支出, 住户存款余额, 
         年末金融机构各项贷款余额, 设施农业种植占地面积, 油料产量, 
         棉花产量, 规模以上工业企业, 固定电话用户, 普通中学在校学生, 
         小学在校学生, 医疗卫生机构床位, 提供住宿的民政服务机构, 
         提供住宿的民政服务机构床位数) -> df4 

df4 

#> # A tibble: 2,054 × 25
#>    省     县     行政区域面积    乡    镇 街道办事处 户籍人口 地区生产总值
#>    <chr>  <chr>         <dbl> <dbl> <dbl>      <dbl>    <dbl>        <dbl>
#>  1 安徽省 肥东县         2182     6    12         NA    108.       8114070
#>  2 安徽省 肥西县         1695     4     8         NA     85.6     10186781
#>  3 安徽省 庐江县         2344    NA    17         NA    120        5471895
#>  4 安徽省 巢湖市         2046    NA    12          5     85.4      5231019
#>  5 安徽省 湾社区          650    NA     5         NA     35.4      3831248
#>  6 安徽省 长丰县         1841     2    12         NA     81.1      7619440
#>  7 安徽省 繁昌区          585    NA     6         NA     27.2      3641111
#>  8 安徽省 南陵县         1264    NA     8         NA     54.3      3206882
#>  9 安徽省 无为市         2022    NA    20         NA    119.       5770037
#> 10 安徽省 怀远县         2192     1    17         NA    134.       3562589
#> # ℹ 2,044 more rows
#> # ℹ 17 more variables: 第一产业增加值 <dbl>, 第二产业增加值 <dbl>,
#> #   第三产业增加值 <dbl>, 地方一般公共预算收入 <dbl>,
#> #   地方一般公共预算支出 <dbl>, 住户存款余额 <dbl>,
#> #   年末金融机构各项贷款余额 <dbl>, 设施农业种植占地面积 <dbl>, 油料产量 <dbl>,
#> #   棉花产量 <dbl>, 规模以上工业企业 <dbl>, 固定电话用户 <dbl>,
#> #   普通中学在校学生 <dbl>, 小学在校学生 <dbl>, 医疗卫生机构床位 <dbl>, …

由于之前我分享的县域统计年鉴数据都是使用的 2020 年行政区划代码，所以这次我们也同样。

2020 年行政区划代码可以从地理矢量数据得到（为了方便绘制地图）：

library(sf)
read_sf("2020行政区划/县.shp") -> county 
# 删除地理信息
county %>% 
  st_drop_geometry() %>% 
  select(-contains("类型")) -> countydb 

# 查看 省-县 的组合有无重复的
countydb %>% 
  group_by(省, 县) %>% 
  mutate(n = n()) %>% 
  arrange(desc(n)) %>% 
  ungroup() -> countydb 

# 这些重复的可能会影响下一步的匹配，所以先删除了
countydb %>% 
  filter(n <= 1) -> countycode1 
countycode1 
#> # A tibble: 2,862 × 7
#>    省     省代码 市     市代码 县     县代码     n
#>    <chr>   <dbl> <chr>   <dbl> <chr>  <chr>  <int>
#>  1 安徽省 340000 安庆市 340800 大观区 340803     1
#>  2 安徽省 340000 安庆市 340800 怀宁县 340822     1
#>  3 安徽省 340000 安庆市 340800 潜山市 340882     1
#>  4 安徽省 340000 安庆市 340800 宿松县 340826     1
#>  5 安徽省 340000 安庆市 340800 太湖县 340825     1
#>  6 安徽省 340000 安庆市 340800 桐城市 340881     1
#>  7 安徽省 340000 安庆市 340800 望江县 340827     1
#>  8 安徽省 340000 安庆市 340800 宜秀区 340811     1
#>  9 安徽省 340000 安庆市 340800 迎江区 340802     1
#> 10 安徽省 340000 安庆市 340800 岳西县 340828     1
#> # ℹ 2,852 more rows

首先匹配下看看能成功多少：

# tidylog 包的 join 族函数可以显示匹配效果：
df4 %>% 
  tidylog::left_join(countycode1) 
#> # A tibble: 2,054 × 30
#>    省     县     行政区域面积    乡    镇 街道办事处 户籍人口 地区生产总值
#>    <chr>  <chr>         <dbl> <dbl> <dbl>      <dbl>    <dbl>        <dbl>
#>  1 安徽省 肥东县         2182     6    12         NA    108.       8114070
#>  2 安徽省 肥西县         1695     4     8         NA     85.6     10186781
#>  3 安徽省 庐江县         2344    NA    17         NA    120        5471895
#>  4 安徽省 巢湖市         2046    NA    12          5     85.4      5231019
#>  5 安徽省 湾社区          650    NA     5         NA     35.4      3831248
#>  6 安徽省 长丰县         1841     2    12         NA     81.1      7619440
#>  7 安徽省 繁昌区          585    NA     6         NA     27.2      3641111
#>  8 安徽省 南陵县         1264    NA     8         NA     54.3      3206882
#>  9 安徽省 无为市         2022    NA    20         NA    119.       5770037
#> 10 安徽省 怀远县         2192     1    17         NA    134.       3562589
#> # ℹ 2,044 more rows
#> # ℹ 22 more variables: 第一产业增加值 <dbl>, 第二产业增加值 <dbl>,
#> #   第三产业增加值 <dbl>, 地方一般公共预算收入 <dbl>,
#> #   地方一般公共预算支出 <dbl>, 住户存款余额 <dbl>,
#> #   年末金融机构各项贷款余额 <dbl>, 设施农业种植占地面积 <dbl>, 油料产量 <dbl>,
#> #   棉花产量 <dbl>, 规模以上工业企业 <dbl>, 固定电话用户 <dbl>,
#> #   普通中学在校学生 <dbl>, 小学在校学生 <dbl>, 医疗卫生机构床位 <dbl>, …

查看匹配失败的：

# 查看匹配失败的 
df4 %>% 
  tidylog::anti_join(countycode1) 
#> # A tibble: 168 × 25
#>    省     县        行政区域面积    乡    镇 街道办事处 户籍人口 地区生产总值
#>    <chr>  <chr>            <dbl> <dbl> <dbl>      <dbl>    <dbl>        <dbl>
#>  1 安徽省 湾社区             650    NA     5         NA     35.4      3831248
#>  2 安徽省 滩溪县            1982    NA    11         NA    114.       5407788
#>  3 安徽省 板阳县            1473     1    15         NA     78.7      1882758
#>  4 安徽省 夥县               857     3     5         NA      9.2       509979
#>  5 安徽省 扬山县            1197    NA    13         NA    101.       2544835
#>  6 重庆市 由日 忠县         2187     6    19          4     96.2      4885521
#>  7 重庆市 新都区             497    NA     2          7     85.6     10001115
#>  8 福建省 沙县区            1799     4     6          2     27.1      3544375
#>  9 福建省 龙海区            1320     2    11          1     91.2     12657788
#> 10 福建省 长泰区             900     1     4         NA     21.1      3752268
#> # ℹ 158 more rows
#> # ℹ 17 more variables: 第一产业增加值 <dbl>, 第二产业增加值 <dbl>,
#> #   第三产业增加值 <dbl>, 地方一般公共预算收入 <dbl>,
#> #   地方一般公共预算支出 <dbl>, 住户存款余额 <dbl>,
#> #   年末金融机构各项贷款余额 <dbl>, 设施农业种植占地面积 <dbl>, 油料产量 <dbl>,
#> #   棉花产量 <dbl>, 规模以上工业企业 <dbl>, 固定电话用户 <dbl>,
#> #   普通中学在校学生 <dbl>, 小学在校学生 <dbl>, 医疗卫生机构床位 <dbl>, …

可以看到很多是由于空格和杂乱字符导致的匹配失败，所以我们先去除：

df4 %>% 
  tidylog::mutate(县 = str_remove_all(县, "[\\s;:]")) -> df4 
然后再匹配：

df4 %>% 
  tidylog::anti_join(countycode1) 
#> # A tibble: 83 × 25
#>    省     县       行政区域面积    乡    镇 街道办事处 户籍人口 地区生产总值
#>    <chr>  <chr>           <dbl> <dbl> <dbl>      <dbl>    <dbl>        <dbl>
#>  1 安徽省 湾社区            650    NA     5         NA     35.4      3831248
#>  2 安徽省 滩溪县           1982    NA    11         NA    114.       5407788
#>  3 安徽省 板阳县           1473     1    15         NA     78.7      1882758
#>  4 安徽省 夥县              857     3     5         NA      9.2       509979
#>  5 安徽省 扬山县           1197    NA    13         NA    101.       2544835
#>  6 重庆市 由日忠县         2187     6    19          4     96.2      4885521
#>  7 重庆市 新都区            497    NA     2          7     85.6     10001115
#>  8 福建省 沙县区           1799     4     6          2     27.1      3544375
#>  9 福建省 龙海区           1320     2    11          1     91.2     12657788
#> 10 福建省 长泰区            900     1     4         NA     21.1      3752268
#> # ℹ 73 more rows
#> # ℹ 17 more variables: 第一产业增加值 <dbl>, 第二产业增加值 <dbl>,
#> #   第三产业增加值 <dbl>, 地方一般公共预算收入 <dbl>,
#> #   地方一般公共预算支出 <dbl>, 住户存款余额 <dbl>,
#> #   年末金融机构各项贷款余额 <dbl>, 设施农业种植占地面积 <dbl>, 油料产量 <dbl>,
#> #   棉花产量 <dbl>, 规模以上工业企业 <dbl>, 固定电话用户 <dbl>,
#> #   普通中学在校学生 <dbl>, 小学在校学生 <dbl>, 医疗卫生机构床位 <dbl>, …

可以看到这个时候匹配不成功的就不是很多了，下面我们需要结合百度和 countycode1.dta 来逐个检查修正：

DT::datatable(countycode1)

这里建议先保存成一个 xlsx 文件，然后在 Excel 里面进行更正：

df4 %>% 
  anti_join(countycode1) %>% 
  writexl::write_xlsx("待修正.xlsx") 

# 待修正2.xlsx 是我手动调整之后得到的结果
readxl::read_xlsx("待修正2.xlsx") -> dftemp 

dftemp %>% 
  tidylog::anti_join(countycode1)
#> # A tibble: 0 × 25
#> # ℹ 25 variables: 省 <chr>, 县 <chr>, 行政区域面积 <dbl>, 乡 <dbl>, 镇 <dbl>,
#> #   街道办事处 <dbl>, 户籍人口 <dbl>, 地区生产总值 <dbl>, 第一产业增加值 <dbl>,
#> #   第二产业增加值 <dbl>, 第三产业增加值 <dbl>, 地方一般公共预算收入 <dbl>,
#> #   地方一般公共预算支出 <dbl>, 住户存款余额 <dbl>,
#> #   年末金融机构各项贷款余额 <dbl>, 设施农业种植占地面积 <dbl>, 油料产量 <dbl>,
#> #   棉花产量 <dbl>, 规模以上工业企业 <dbl>, 固定电话用户 <dbl>,
#> #   普通中学在校学生 <dbl>, 小学在校学生 <dbl>, 医疗卫生机构床位 <dbl>, …

这个时候就没有不匹配的了：

df4 %>% 
  anti_join(
    df4 %>% 
      select(省, 县) %>% 
      anti_join(countycode1)
  ) %>% 
  bind_rows(dftemp) -> df5 

df5 %>% 
  tidylog::left_join(countycode1) %>% 
  select(省, 省代码, 市, 市代码, 县, 县代码, everything()) -> df5 

然后我们再对变量进行重命名（和之前年份的保持一致）：

df5 %>% 
  mutate(乡 = if_else(is.na(乡), 0, 乡),
         镇 = if_else(is.na(镇), 0, 镇)) %>% 
  mutate(乡镇个数_个 = 乡 + 镇) %>% 
  rename(
    乡_个 = 乡, 
    镇_个 = 镇, 
    行政区域土地面积_平方公里 = 行政区域面积, 
    街道办事处个数_个 = 街道办事处, 
    户籍人口_人 = 户籍人口, 
    国内生产总值_万元 = 地区生产总值,
    第一产业增加值_万元 = 第一产业增加值,
    第二产业增加值_万元 = 第二产业增加值, 
    第三产业增加值_万元 = 第三产业增加值,
    一般公共预算收入_万元 = 地方一般公共预算收入,
    一般公共预算支出_万元 = 地方一般公共预算支出,
    住户储蓄存款余额_万元 = 住户存款余额, 
    年末各项贷款余额_万元 = 年末金融机构各项贷款余额, 
    设施农业占地面积_公顷 = 设施农业种植占地面积, 
    油料产量_吨 = 油料产量, 
    棉花产量_吨 = 棉花产量, 
    规模以上工业企业个数_个 = 规模以上工业企业, 
    固定电话用户_户 = 固定电话用户, 
    普通中学在校学生_人 = 普通中学在校学生, 
    小学在校学生数_人 = 小学在校学生, 
    医疗卫生机构床位_床 = 医疗卫生机构床位,
    提供住宿的社会工作机构_个 = 提供住宿的民政服务机构, 
    提供住宿的社会工作机构床位_床 = 提供住宿的民政服务机构床位数
  ) -> df6 

# 保存成 xlsx
df6 %>% 
  writexl::write_xlsx("2021年县市社会经济指标.xlsx") 

最后如果你想把该数据和之前年份的合并起来，只需要使用 bind_rows() 合并即可。

最后我们再使用该数据绘制一幅区县地图。这里使用的数据是我之前编辑过的一份 shp 数据。可以用于绘制带九段线小地图的中国地图。

library(ggspatial)
read_sf("chinacounty2020mini/chinacounty2020mini.shp") -> countymap 
read_sf("chinacounty2020mini/chinacounty2020mini_line.shp") -> countyline 

countymap %>% 
  filter(!is.na(objid)) -> countymap 

以地区生产总值为例：

df6 %>% 
  select(县代码, 国内生产总值_万元) %>% 
  mutate(v = 国内生产总值_万元 / 10000) -> df7 

缺失值使用所在市、所在省的均值填补，实在无法填补的设定为 -1：

countymap %>% 
  left_join(df7) %>% 
  group_by(省, 省代码, 市, 市代码) %>% 
  mutate(mean1 = mean(v, na.rm = T),
         v = if_else(is.na(v), mean1, v)) %>% 
  ungroup() %>% 
  group_by(省, 省代码) %>% 
  mutate(mean2 = mean(v, na.rm = T),
         v = if_else(is.na(v), mean2, v)) %>% 
  ungroup() %>% 
  mutate(v = if_else(is.na(v), -1, v)) -> df8 

下面我们将绘制两种地图，一种是使用连续变量绘制，另一种是使用分类变量绘制，为此，我们对地区生产总值变量进行分组：

# v 的范围
range(df8$v) 
#> [1]   -1.00 4748.06
quantile(df8$v, seq(from = 0, to = 1, length = 8)) %>% 
  `[`(2:7) %>% 
  round(digits = 2) -> cuts

df8 %>% 
  mutate(group = cut(v, breaks = c(-1, 0, cuts, 5000), 
                     include.lowest = T,
                     labels = c("无数据", "0 ~ 65.35",
                                "65.35 ~ 109.36", "109.36 ~ 154.74",
                                "154.74 ~ 229.96", "229.96 ~ 328.96",
                                "328.96 ~ 521.38", "> 521.38")))  -> df9 

countyline 变量还需要再处理下：

countyline %>% 
  tail(n = 9) %>% 
  select(class) %>% 
  filter(class %in% c("九段线", "海岸线", "小地图框格")) -> countyline 

连续变量的绘制：

# 绘制连续变量
ggplot(df9) + 
  geom_sf(aes(fill = v), linewidth = 0.001, color = "black") + 
  geom_sf(data = countyline, aes(color = class, linewidth = class), 
          show.legend = F) + 
  scale_color_manual(values = c(
    "九段线" = "black", 
    "海岸线" = "#0055AA", 
    "小地图框格"= "black"
  )) + 
  scale_linewidth_manual(values = c(
    "九段线" = 0.3, 
    "海岸线" = 0.3, 
    "小地图框格"= 0.2
  )) + 
  scale_fill_viridis_c(option = "A", trans = "log10") + 
  annotation_scale(
    width_hint = 0.2,
    text_family = cnfont 
  ) + 
  annotation_north_arrow(
    location = "tr", 
    width = unit(2, "cm"),
    height = unit(2, "cm"),
    which_north = "false",
    pad_x = unit(0.5, "cm"),
    pad_y = unit(0.5, "cm"),
    style = north_arrow_nautical(
      text_family = cnfont
    )
  ) + 
  guides(fill = guide_colorbar(title = "地区生产总值（亿元）")) + 
  theme_ipsum(base_family = cnfont, grid = F) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = c(0.12, 0.2),
        plot.background = element_rect(fill = "white", color = "white")) + 
  labs(title = "2021 年中国各县地区生产总值（亿元）",
       subtitle = "数据整理 & 绘制：微信公众号 RStata",
       caption = "数据来源：2022 年中国县域统计年鉴") -> p 

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

分类变量的绘制：

# 绘制分类变量
ggplot(df9) + 
  geom_sf(aes(fill = group), linewidth = 0.001, color = "black") + 
  geom_sf(data = countyline, aes(color = class, linewidth = class), 
          show.legend = F) + 
  scale_color_manual(values = c(
    "九段线" = "black", 
    "海岸线" = "#0055AA", 
    "小地图框格"= "black"
  )) + 
  scale_linewidth_manual(values = c(
    "九段线" = 0.3, 
    "海岸线" = 0.3, 
    "小地图框格"= 0.2
  )) + 
  scale_fill_viridis_d(option = "A") + 
  annotation_scale(
    width_hint = 0.2,
    text_family = cnfont 
  ) + 
  annotation_north_arrow(
    location = "tr", 
    width = unit(2, "cm"),
    height = unit(2, "cm"),
    which_north = "false",
    pad_x = unit(0.5, "cm"),
    pad_y = unit(0.5, "cm"),
    style = north_arrow_nautical(
      text_family = cnfont
    )
  ) + 
  guides(fill = guide_legend(title = "地区生产总值（亿元）",
                             ncol = 2)) + 
  theme_ipsum(base_family = cnfont, grid = F) + 
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(), 
        legend.position = c(0.15, 0.2), 
        plot.background = element_rect(fill = "white", color = "white")) + 
  labs(title = "2021 年中国各县地区生产总值（亿元）",
       subtitle = "数据整理 & 绘制：微信公众号 RStata",
       caption = "数据来源：2022 年中国县域统计年鉴") -> p 

ggsave("pic2.png", width = 10, height = 8, device = png)