使用 R 语言计算胡焕庸线两侧的省市区县（含结果）

前不久给大家分享了如何使用 R 语言计算各省市区县位于秦岭淮河线的北侧、南侧还是被穿过，并且分享了计算好的数据，最近又有小伙伴问是否也可以计算下胡焕庸线的结果。

胡焕庸线也即黑河-腾冲线，是我国地理学家胡焕庸在 1935 年提出的划分我国人口密度的对比线，例如之前使用 R 语言绘制的这幅图人口密度图：

该图来源于课程：https://rstata.duanshu.com/#/brief/course/8f921af271a04b6db79c52ce48a04d49

下面我们来以区县为例，使用 R 语言提取胡焕庸线东西侧的区县。

首先加载相应的 R 包以及读取数据：

library(sf)
library(tidyverse)
# 在不同的坐标系里面，胡焕庸线可能穿过不同的区县，这里以下面的坐标系计算
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("2019行政区划/县.shp") %>%
  st_transform(mycrs) %>%
  rename(县代码 = PAC, 县 = NAME) %>%
  st_make_valid() -> county
read_sf("胡焕庸线东侧的区域.geojson") %>%
  st_cast("POLYGON") %>%
  st_transform(mycrs) %>%
  st_make_valid() -> east
read_sf("胡焕庸线西侧的区域.geojson") %>%
  st_cast("POLYGON") %>%
  st_transform(mycrs) %>%
  st_make_valid() -> west
read_sf("胡焕庸线/胡焕庸线.shp") %>%
  st_transform(mycrs) %>%
  st_make_valid() -> hhy

其中 “胡焕庸线东侧的区域.geojson” 和 “胡焕庸线西侧的区域.geojson” 数据是根据胡焕庸线绘制的：https://tidyfriday.cn/geojson/

然后就可以分别计算和胡焕庸线相交的以及两侧的区县了：

# 胡焕庸线上的区县
county %>%
  st_make_valid() %>%
  st_intersection(hhy) %>%
  st_drop_geometry() %>%
  select(-contains("类型"), -FID) -> county1

# 胡焕庸线西侧的区县
county %>%
  st_intersection(west) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  anti_join(county1)-> county2

# 胡焕庸线东侧的区县
county %>%
  st_intersection(east) %>%
  st_drop_geometry() %>%
  select(-contains("类型")) %>%
  anti_join(county1)-> county3

合并三部分结果：

bind_rows(
  county1 %>% mutate(类别 = "中间区县"),
  county2 %>% mutate(类别 = "西侧区县"),
  county3 %>% mutate(类别 = "东侧区县")
) -> resdf

最后我们再计算下各个区县的质心以及这些质心距离秦岭淮河线的距离：

county %>%
  left_join(resdf) %>%
  st_centroid() %>%
  mutate(距离 = st_distance(., hhy)[,1]) %>%
  mutate(距离 = as.numeric(距离) / 1000) %>%
  st_drop_geometry() %>%
  dplyr::select(-contains("类型")) %>%
  rename(距离胡焕庸线_km = 距离) -> resdf

resdf %>%
  writexl::write_xlsx("2019年各区县东西侧分类及其质心与胡焕庸线的距离.xlsx")

这个结果大概是这样的：

resdf
#> # A tibble: 2,900 × 8
#>    县代码 县       省代码 省     市代码 市     类别     距离胡焕庸线_km
#>     <dbl> <chr>     <dbl> <chr>   <dbl> <chr>  <chr>              <dbl>
#>  1 110101 东城区   110000 北京市 110000 北京市 东侧区县            237.
#>  2 110102 西城区   110000 北京市 110000 北京市 东侧区县            234.
#>  3 110105 朝阳区   110000 北京市 110000 北京市 东侧区县            240.
#>  4 110106 丰台区   110000 北京市 110000 北京市 东侧区县            233.
#>  5 110107 石景山区 110000 北京市 110000 北京市 东侧区县            221.
#>  6 110108 海淀区   110000 北京市 110000 北京市 东侧区县            217.
#>  7 110109 门头沟区 110000 北京市 110000 北京市 东侧区县            193.
#>  8 110111 房山区   110000 北京市 110000 北京市 东侧区县            219.
#>  9 110112 通州区   110000 北京市 110000 北京市 东侧区县            265.
#> 10 110113 顺义区   110000 北京市 110000 北京市 东侧区县            236.
#> # … with 2,890 more rows