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)