使用 R 语言计算各个城市_区县的质心距离省界边界的距离

2138 字 约 8 分钟 课程 距离与边界地形与地理R 语言

这是培训班的一个小伙伴遇到的问题，使用一些地理计算就能解决了。

计算各个城市/区县的质心距离省界的距离主要用两个函数，st_centroid() 计算质心，st_distance() 计算距离。

首先加载相关 R 包然后读取省市区县的行政区划数据：

library(tidyverse)
library(sf)

# 读取数据
read_sf("2019行政区划/县.shp") -> county
read_sf("2019行政区划/市.shp") -> city
read_sf("2019行政区划/省.shp") -> prov

计算城市和区县的质心并且从 prov 里面提取省界（线条）：

county %>%
  st_centroid() -> county_centroid
city %>%
  st_centroid() -> city_centroid

# 提取省界
st_cast(prov, "MULTILINESTRING") -> provline

计算所有城市区县之前我们先以安徽省为例计算下：

# 以安徽省为例
provline %>%
  dplyr::filter(省 == "安徽省") -> prov_temp
county_centroid %>%
  dplyr::filter(省 == "安徽省") -> county_temp
city_centroid %>%
  dplyr::filter(省 == "安徽省") -> city_temp

county_temp %>%
  mutate(dist = st_distance(county_temp,
                            prov_temp)[,1],
         dist = units::set_units(dist, km),
         dist = as.numeric(dist))

#> Simple feature collection with 105 features and 8 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 115.2367 ymin: 29.66561 xmax: 119.3527 ymax: 34.44854
#> geographic CRS: CGCS_2000
#> # A tibble: 105 x 9
#>       PAC NAME   省代码 省     市代码 市     类型              geometry  dist
#>  *  <dbl> <chr>   <dbl> <chr>   <dbl> <chr>  <chr>          <POINT [°]> <dbl>
#>  1 340102 瑶海区 340000 安徽省 340100 合肥市 市辖区 (117.3475 31.93624)  95.1
#>  2 340103 庐阳区 340000 安徽省 340100 合肥市 市辖区 (117.1978 31.91693) 109.
#>  3 340104 蜀山区 340000 安徽省 340100 合肥市 市辖区 (117.0553 31.86933) 107.
#>  4 340111 包河区 340000 安徽省 340100 合肥市 市辖区  (117.326 31.74349)  99.6
#>  5 340121 长丰县 340000 安徽省 340100 合肥市 县     (117.1931 32.23517) 113.
#>  6 340122 肥东县 340000 安徽省 340100 合肥市 县     (117.5707 32.00112)  74.4
#>  7 340123 肥西县 340000 安徽省 340100 合肥市 县     (117.0026 31.68642) 105.
#>  8 340124 庐江县 340000 安徽省 340100 合肥市 县     (117.3216 31.26925) 122.
#>  9 340181 巢湖市 340000 安徽省 340100 合肥市 县级市 (117.7172 31.64657)  68.2
#> 10 340202 镜湖区 340000 安徽省 340200 芜湖市 市辖区 (118.4262 31.27212)  25.9
#> # … with 95 more rows

city_temp %>%
  mutate(dist = st_distance(city_temp,
                            prov_temp)[,1],
         dist = units::set_units(dist, km),
         dist = as.numeric(dist))

#> Simple feature collection with 16 features and 6 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 115.7038 ymin: 29.90664 xmax: 118.8522 ymax: 33.85868
#> geographic CRS: CGCS_2000
#> # A tibble: 16 x 7
#>    省代码 省     市代码 市       类型              geometry  dist
#>  *  <dbl> <chr>   <dbl> <chr>    <chr>          <POINT [°]> <dbl>
#>  1 340000 安徽省 340100 合肥市   地级市 (117.3558 31.76344)  96.4
#>  2 340000 安徽省 340200 芜湖市   地级市 (118.1346 31.16267)  55.7
#>  3 340000 安徽省 340300 蚌埠市   地级市 (117.3254 33.10857)  58.5
#>  4 340000 安徽省 340400 淮南市   地级市   (116.7684 32.472)  80.0
#>  5 340000 安徽省 340500 马鞍山市 地级市   (118.3651 31.637)  18.6
#>  6 340000 安徽省 340600 淮北市   地级市 (116.7435 33.72656)  20.1
#>  7 340000 安徽省 340700 铜陵市   地级市 (117.5561 30.88418) 111.
#>  8 340000 安徽省 340800 安庆市   地级市 (116.4867 30.57605)  55.0
#>  9 340000 安徽省 341000 黄山市   地级市 (118.0711 29.90664)  36.9
#> 10 340000 安徽省 341100 滁州市   地级市 (118.1026 32.54416)  31.4
#> 11 340000 安徽省 341200 阜阳市   地级市 (115.7038 32.91734)  41.2
#> 12 340000 安徽省 341300 宿州市   地级市 (117.2108 33.85868)  23.0
#> 13 340000 安徽省 341500 六安市   地级市    (116.2297 31.66)  34.2
#> 14 340000 安徽省 341600 亳州市   地级市 (116.1807 33.43546)  30.2
#> 15 340000 安徽省 341700 池州市   地级市 (117.3656 30.28377)  42.0
#> 16 340000 安徽省 341800 宣城市   地级市 (118.8522 30.68624)  37.1

然后我们就可以循环所有的省了，这里我没有使用循环，而是使用了 lapply()：

lapply(unique(provline$省), function(x){
  provline %>%
    dplyr::filter(省 == x) -> prov_temp
  county_centroid %>%
    dplyr::filter(省 == x) -> county_temp

  county_temp %>%
    mutate(dist = st_distance(county_temp,
                              prov_temp)[,1],
           dist = units::set_units(dist, km),
           dist = as.numeric(dist))
}) %>%
  bind_rows() -> df
df

#> Simple feature collection with 2899 features and 8 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 74.90689 ymin: 16.84913 xmax: 134.2809 ymax: 52.93326
#> geographic CRS: CGCS_2000
#> # A tibble: 2,899 x 9
#>       PAC NAME     省代码 省     市代码 市     类型              geometry  dist
#>  *  <dbl> <chr>     <dbl> <chr>   <dbl> <chr>  <chr>          <POINT [°]> <dbl>
#>  1 110101 东城区   110000 北京市 110000 北京市 市辖区 (116.4106 39.91194)  29.6
#>  2 110102 西城区   110000 北京市 110000 北京市 市辖区 (116.3596 39.91068)  33.9
#>  3 110105 朝阳区   110000 北京市 110000 北京市 市辖区 (116.5094 39.95239)  20.7
#>  4 110106 丰台区   110000 北京市 110000 北京市 市辖区  (116.244 39.83434)  27.7
#>  5 110107 石景山区 110000 北京市 110000 北京市 市辖区 (116.1702 39.93219)  36.5
#>  6 110108 海淀区   110000 北京市 110000 北京市 市辖区  (116.2271 40.0256)  34.5
#>  7 110109 门头沟区 110000 北京市 110000 北京市 市辖区 (115.7856 39.99295)  15.8
#>  8 110111 房山区   110000 北京市 110000 北京市 市辖区 (115.8476 39.71777)  14.0
#>  9 110112 通州区   110000 北京市 110000 北京市 市辖区 (116.7267 39.80217)  10.1
#> 10 110113 顺义区   110000 北京市 110000 北京市 市辖区 (116.7179 40.14944)  13.7
#> # … with 2,889 more rows

# 然后我们就可以把这个数据输出成 xlsx 文件了：
bind_cols(
  st_drop_geometry(df),
  st_coordinates(df) %>%
    as_tibble() %>%
    set_names(c("经度", "纬度"))
) %>%
  writexl::write_xlsx("各个区县距离其所在省份边界的距离(单位：km).xlsx")

类似的方法求各个城市的：

lapply(unique(provline$省), function(x){
  provline %>%
    dplyr::filter(省 == x) -> prov_temp
  city_centroid %>%
    dplyr::filter(省 == x) -> city_temp

  city_temp %>%
    mutate(dist = st_distance(city_temp,
                              prov_temp)[,1],
           dist = units::set_units(dist, km),
           dist = as.numeric(dist))
}) %>%
  bind_rows() -> df2

df2

#> Simple feature collection with 371 features and 6 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 75.94554 ymin: 16.84913 xmax: 132.5427 ymax: 52.32577
#> geographic CRS: CGCS_2000
#> # A tibble: 371 x 7
#>    省代码 省     市代码 市       类型              geometry  dist
#>  *  <dbl> <chr>   <dbl> <chr>    <chr>          <POINT [°]> <dbl>
#>  1 110000 北京市 110000 北京市   直辖市 (116.4123 40.18554)  35.1
#>  2 120000 天津市 120000 天津市   直辖市  (117.334 39.29343)  39.0
#>  3 130000 河北省 130100 石家庄市 地级市 (114.4399 38.13127)  49.6
#>  4 130000 河北省 130200 唐山市   地级市 (118.3359 39.72124)  38.7
#>  5 130000 河北省 130300 秦皇岛市 地级市  (119.185 40.08797)  37.1
#>  6 130000 河北省 130400 邯郸市   地级市   (114.5421 36.553)  37.7
#>  7 130000 河北省 130500 邢台市   地级市 (114.8159 37.21264)  74.3
#>  8 130000 河北省 130600 保定市   地级市 (115.1706 39.02127)  64.0
#>  9 130000 河北省 130700 张家口市 地级市 (115.0324 40.86512)  70.6
#> 10 130000 河北省 130800 承德市   地级市 (117.5486 41.34761)  61.3
#> # … with 361 more rows

bind_cols(
  st_drop_geometry(df2),
  st_coordinates(df2) %>%
    as_tibble() %>%
    set_names(c("经度", "纬度"))
) %>%
  writexl::write_xlsx("各个城市距离其所在省份边界的距离(单位：km).xlsx")

最后我们再使用 ggplot2 绘制两幅图展示下：

区县的：

# 读取中国省级地图数据
read_sf('china_prov_full_map.json') -> cn

# 绘图
library(ggspatial)
mybreaks <- seq(from = 0, to = 600, by = 100)
ggplot(cn) +
  geom_sf(fill = NA, size = 0.1, color = "black") +
  geom_sf(data = df,
          aes(size = dist,
              color = dist),
          alpha = 0.8) +
  scale_size_continuous(range = c(0.01, 3),
                        breaks = mybreaks) +
  scico::scale_color_scico(breaks = mybreaks, palette = "lisbon") +
  guides(color = guide_legend()) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "bottom") +
  labs(size = "距离省界(km)",
       color = "距离省界(km)",
       title = "各区县质心距离省界的距离",
       subtitle = "使用 2019 年的行政区划数据计算得到",
       caption = "数据来源：2019 年中国各区县行政区划") +
  annotation_scale(location = "bl", width_hint = 0.3,
                   text_family = cnfont) +
  annotation_north_arrow(location = "tr", which_north = "false",
                         pad_x = unit(0.9, "cm"),
                         pad_y = unit(0.4, "cm"),
                         style = north_arrow_fancy_orienteering(
                           text_family = cnfont
                         ))
各区县质心距离省界的距离

城市的：

ggplot(cn) +
  geom_sf(fill = NA, size = 0.1, color = "black") +
  geom_sf(data = df2,
          aes(size = dist,
              color = dist),
          alpha = 0.8) +
  scale_size_continuous(range = c(0.01, 4),
                        breaks = mybreaks) +
  scico::scale_color_scico(breaks = mybreaks, palette = "lisbon") +
  guides(color = guide_legend()) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "bottom") +
  labs(size = "距离省界(km)",
       color = "距离省界(km)",
       title = "各城市质心距离省界的距离",
       subtitle = "使用 2019 年的行政区划数据计算得到",
       caption = "数据来源：2019 年中国各城市行政区划") +
  annotation_scale(location = "bl", width_hint = 0.3,
                   text_family = cnfont) +
  annotation_north_arrow(location = "tr", which_north = "false",
                         pad_x = unit(0.9, "cm"),
                         pad_y = unit(0.4, "cm"),
                         style = north_arrow_fancy_orienteering(
                           text_family = cnfont
                         ))
各城市质心距离省界的距离

直播讲解

为了让大家更好的理解上面的代码，欢迎各位培训班的小伙伴参加明晚 8 点的直播课「使用 R 语言计算各个城市/区县距离省界的距离」

  1. 直播时间：2021 年 3 月 19 日晚上 8 点；
  2. 直播地址：腾讯会议（需要报名培训班参加）
  3. 如何报名 RStata 培训班，可以阅读这篇推文了解：听说这里又能学习 R 语言和 Stata，又能答疑！零基础也不用担心！

点击这里跳转到 RStata 短书平台获取附件：使用 R 语言计算各个城市/区县的质心距离省界边界的距离

相关文章

评论