使用 R 语言计算秦岭—淮河线两侧的省市区县（含结果）

最近有小伙伴问到如何判断每个区县位于秦岭淮河线的北侧、南侧还是被穿过。这个可以使用 R 语言计算，也就是一种空间拓扑计算，今天就让我们一起看一下如何操作。

附件中我准备了如下地理矢量数据：

1. 秦岭淮河线（新）.geojson
2. 2019行政区划：各省市区县行政区划矢量数据；
3. 秦岭淮河线北方区域.geojson： 也是地理矢量数据
4. 秦岭淮河线南方区域.geojson： 也是地理矢量数据

其中后面两个文件是使用这个应用制作的：https://tidyfriday.cn/geojson/

这里使用的秦岭淮河线仅供参考，大家要是感觉不准，也可以替换成适合自己研究的秦岭淮河线数据。

首先加载所需的 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("秦岭淮河线南方区域.geojson") %>%
  st_cast("POLYGON") %>%
  st_transform(mycrs) -> south
read_sf("秦岭淮河线北方区域.geojson") %>%
  st_cast("POLYGON") %>%
  st_transform(mycrs) -> north
read_sf("秦岭淮河线（新）.geojson") %>%
  rename(FID = Id) %>%
  st_transform(mycrs) -> qh
read_sf("2019行政区划/县.shp") %>%
  st_transform(mycrs) %>%
  rename(县代码 = PAC, 县 = NAME) %>%
  st_make_valid() -> county

这里的 st_make_valid() 用以纠正可能存在的拓扑错误。其中 north 和 south 就是南北方区域的范围了，使用 st_intersection() 就可以提取省市区县中位于分别位于南北方以及被秦岭淮河线穿过的了：

# 秦岭淮河线上的区县
county %>%
  st_make_valid() %>%
  st_intersection(qh) %>%
  st_drop_geometry() %>%
  select(-contains("类型"), -FID) -> county1

# 秦岭淮河线北方的区县
county %>%
  st_intersection(north) %>%
  st_drop_geometry() %>%
  select(-Id, -contains("类型")) %>%
  anti_join(county1)-> county2

# 秦岭淮河线南方的区县
county %>%
  st_intersection(south) %>%
  st_drop_geometry() %>%
  select(-Id, -contains("类型")) %>%
  anti_join(county1)-> county3

bind_rows(
  county1 %>% mutate(类别 = "中部区县"),
  county2 %>% mutate(类别 = "北方区县"),
  county3 %>% mutate(类别 = "南方区县")
) -> resdf

然后再计算质心和距离：

county %>%
  left_join(resdf) %>%
  st_centroid() %>%
  mutate(距离 = st_distance(., qh)[,1]) %>%
  mutate(距离 = as.numeric(距离) / 1000) %>%
  st_drop_geometry() %>%
  dplyr::select(-contains("类型")) %>%
  rename(距离秦岭淮河线_km = 距离) -> resdf

resdf %>%
  mutate(类别 = case_when(
    str_sub(省, 1, 2) == "新疆" ~ "北方区县",
    str_sub(省, 1, 2) == "青海" ~ "北方区县",
    str_sub(省, 1, 2) == "甘肃" ~ "北方区县",
    str_sub(省, 1, 2) == "西藏" ~ "南方区县",
    str_sub(省, 1, 2) == "四川" ~ "南方区县",
    str_sub(省, 1, 2) == "湖北" ~ "南方区县",
    T ~ 类别
  )) -> 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 北京市 北方区县              703.
#>  2 110102 西城区   110000 北京市 110000 北京市 北方区县              705.
#>  3 110105 朝阳区   110000 北京市 110000 北京市 北方区县              703.
#>  4 110106 丰台区   110000 北京市 110000 北京市 北方区县              702.
#>  5 110107 石景山区 110000 北京市 110000 北京市 北方区县              715.
#>  6 110108 海淀区   110000 北京市 110000 北京市 北方区县              721.
#>  7 110109 门头沟区 110000 北京市 110000 北京市 北方区县              737.
#>  8 110111 房山区   110000 北京市 110000 北京市 北方区县              708.
#>  9 110112 通州区   110000 北京市 110000 北京市 北方区县              680.
#> 10 110113 顺义区   110000 北京市 110000 北京市 北方区县              715.
#> # … with 2,890 more rows

城市和省份的操作类似，这里就不再展示代码了，经过计算可以得到下面两个文件：

• 2019年各省份南北方分类.xlsx
• 2019年各城市南北方分类.xlsx

下面我们再绘图展示下计算的结果：

library(ggspatial)
library(showtext)
showtext_auto(enable = TRUE)
font_add("cnfont", regular = "Tiejili Regular.ttf")
cnfont <- "cnfont"
source("theme.R")

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("海岸线/海岸线.shp") %>%
  st_transform(mycrs) -> hax
read_sf("九段线/九段线.shp") %>%
  st_transform(mycrs) -> jdx
  read_sf("秦岭淮河线（新）.geojson") %>%
  st_transform(mycrs) -> qh
  read_sf("2019行政区划/县.shp") %>%
  st_transform(mycrs) %>%
  rename(县代码 = PAC, 县 = NAME) %>%
  st_make_valid() -> county

readxl::read_xlsx("区县结果/2019年各区县南北方分类及其质心与秦岭淮河线的距离.xlsx") -> dfcounty

county %>%
  left_join(dfcounty) -> dfcounty

dfcounty %>%
  st_centroid() -> dfcounty2

ggplot() +
  geom_sf(data = dfcounty, aes(fill = 类别), size = 0.01, color = NA) +
  geom_sf(data = hax, color = "#0055AA", size = 0.3) +
  geom_sf(data = jdx, color = "black", size = 0.6) +
  geom_sf(data = qh, color = "#eac862", size = 2) +
  geom_sf(data = dfcounty2, aes(size = 距离秦岭淮河线_km),
          color = "black", shape = 21) +
  scale_fill_manual(values = c("#0055aa", "#c40003", "#00c19b")) +
  annotation_scale(location = "bl",
                   width_hint = 0.3,
                   text_family = cnfont) +
  labs(title = "秦岭淮河线北方、南方和穿过的区县",
       caption = "绘制：微信公众号 RStata",
       fill = "类别") +
  theme(legend.position = c(0.2, 0.2)) +
  annotation_north_arrow(
    location = "tr",
    which_north = "false",
    pad_y = unit(0.1, "cm"),
    style = north_arrow_fancy_orienteering(
      text_family = cnfont
    )
  ) -> p

ggsave("秦岭淮河线北方、南方和穿过的区县.png", device = png, width = 10, height = 10)
knitr::plot_crop("秦岭淮河线北方、南方和穿过的区县.png")

省市的绘图结果如下：

)

另外之前给大家分享过 1949～2021 年历年的省市区县矢量数据，这里正好我写了个循环计算了历年的结果。

×

另外可以结合这个课程掌握本文的代码：

×

相关代码和结果都在附件中～

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言计算秦岭—淮河线两侧的省市区县（含结果）