使用 R 语言计算各个区县距离其所在城市的质心的距离

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

如何使用 R 语言计算每个区县的质心距离其所属城市的质心的距离？显然这也是一个地理计算的问题，我们可以通过下面的程序来实现。

如何使用 R 语言计算每个区县的质心距离其所属城市的质心的距离？显然这也是一个地理计算的问题，我们可以通过下面的程序来实现。

首先我们加载所需的 R 包，读取城市和区县的矢量数据：

library(tidyverse)
library(sf)

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

计算城市和区县的质心可以使用 st_centroid() 函数：

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

以广州为例，广州的各个区的质心与广州市的质心距离可以使用下面的代码计算：

city_centroid %>%
  dplyr::filter(市代码 == 440100) -> city_temp

county_centroid %>%
  dplyr::filter(市代码 == 440100) -> county_temp

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

#> Simple feature collection with 11 features and 8 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 113.2087 ymin: 22.77021 xmax: 113.7588 ymax: 23.67957
#> geographic CRS: CGCS_2000
#> # A tibble: 11 x 9
#>       PAC NAME   省代码 省     市代码 市     类型              geometry  dist
#>  *  <dbl> <chr>   <dbl> <chr>   <dbl> <chr>  <chr>          <POINT [°]> <dbl>
#>  1 440103 荔湾区 440000 广东省 440100 广州市 市辖区 (113.2224 23.09039)  43.0
#>  2 440104 越秀区 440000 广东省 440100 广州市 市辖区  (113.276 23.13553)  35.5
#>  3 440105 海珠区 440000 广东省 440100 广州市 市辖区 (113.3216 23.08441)  36.5
#>  4 440106 天河区 440000 广东省 440100 广州市 市辖区 (113.3706 23.16218)  26.7
#>  5 440111 白云区 440000 广东省 440100 广州市 市辖区 (113.3184 23.29228)  23.1
#>  6 440112 黄埔区 440000 广东省 440100 广州市 市辖区 (113.5026 23.21921)  14.7
#>  7 440113 番禺区 440000 广东省 440100 广州市 市辖区 (113.4012 22.97682)  43.4
#>  8 440114 花都区 440000 广东省 440100 广州市 市辖区 (113.2087 23.44502)  35.1
#>  9 440115 南沙区 440000 广东省 440100 广州市 市辖区  (113.504 22.77021)  64.1
#> 10 440117 从化区 440000 广东省 440100 广州市 市辖区 (113.6752 23.67957)  39.4
#> 11 440118 增城区 440000 广东省 440100 广州市 市辖区 (113.7588 23.34319)  22.8

使用类似的方法我们就可以计算所有城市的了：

lapply(unique(city$市代码), function(x){
  city_centroid %>%
    dplyr::filter(市代码 == x) -> city_temp
  county_centroid %>%
    dplyr::filter(市代码 == x) -> county_temp
  county_temp %>%
    mutate(dist = st_distance(county_temp,
                              city_temp)[,1],
           dist = units::set_units(dist, km),
           dist = as.numeric(dist))
}) %>%
  bind_rows() -> ccdist

ccdist
#> Simple feature collection with 2877 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,877 x 9
#>       PAC NAME     省代码 省     市代码 市     类型              geometry  dist
#>  *  <dbl> <chr>     <dbl> <chr>   <dbl> <chr>  <chr>          <POINT [°]> <dbl>
#>  1 110101 东城区   110000 北京市 110000 北京市 市辖区 (116.4106 39.91194)  30.4
#>  2 110102 西城区   110000 北京市 110000 北京市 市辖区 (116.3596 39.91068)  30.8
#>  3 110105 朝阳区   110000 北京市 110000 北京市 市辖区 (116.5094 39.95239)  27.2
#>  4 110106 丰台区   110000 北京市 110000 北京市 市辖区  (116.244 39.83434)  41.6
#>  5 110107 石景山区 110000 北京市 110000 北京市 市辖区 (116.1702 39.93219)  34.9
#>  6 110108 海淀区   110000 北京市 110000 北京市 市辖区  (116.2271 40.0256)  23.8
#>  7 110109 门头沟区 110000 北京市 110000 北京市 市辖区 (115.7856 39.99295)  57.6
#>  8 110111 房山区   110000 北京市 110000 北京市 市辖区 (115.8476 39.71777)  70.9
#>  9 110112 通州区   110000 北京市 110000 北京市 市辖区 (116.7267 39.80217)  50.3
#> 10 110113 顺义区   110000 北京市 110000 北京市 市辖区 (116.7179 40.14944)  26.3
#> # … with 2,867 more rows

把 ccdist 输出为 xlsx 文件：

bind_cols(
  st_drop_geometry(ccdist),
  st_coordinates(ccdist) %>%
    as_tibble() %>%
    set_names(c("经度", "纬度"))
) %>%
  writexl::write_xlsx("各个区县距离其所在城市的质心距离(单位：km).xlsx")

最后我们再绘制一幅图展示下：

# 读取中国省级地图数据
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 = ccdist,
          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
                         ))
各区县质心与其所属城市质心的距离

直播讲解

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

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

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

相关文章

评论