最近有培训班小伙伴问到了这样一个问题，如何计算每个城市距离所在省的省会城市的距离。这个问题比较简单，我们就一起来解决下。
首先我们读取 2019 年中国各城市的行政区划矢量数据：
library(tidyverse)
library(sf)
read_sf("2019行政区划/市.shp") -> city
city
#> Simple feature collection with 371 features and 5 fields
#> geometry type: MULTIPOLYGON
#> dimension: XY
#> bbox: xmin: 73.50114 ymin: 6.323421 xmax: 135.0885 ymax: 53.5609
#> geographic CRS: CGCS_2000
#> # A tibble: 371 x 6
#> 省代码 省 市代码 市 类型 geometry
#> <dbl> <chr> <dbl> <chr> <chr> <MULTIPOLYGON [°]>
#> 1 110000 北京市 110000 北京市 直辖市… (((116.6753 41.0401, 116.6762 41.04006, 11…
#> 2 120000 天津市 120000 天津市 直辖市… (((117.4438 40.25101, 117.4561 40.24615, 1…
#> 3 130000 河北省 130100 石家庄市… 地级市… (((113.8242 38.75805, 113.8312 38.74815, 1…
#> 4 130000 河北省 130200 唐山市 地级市… (((118.8539 39.10692, 118.8493 39.10679, 1…
#> 5 130000 河北省 130300 秦皇岛市… 地级市… (((119.1521 40.6128, 119.1517 40.60917, 11…
#> 6 130000 河北省 130400 邯郸市 地级市… (((113.8711 37.01219, 113.8724 37.01182, 1…
#> 7 130000 河北省 130500 邢台市 地级市… (((115.1259 37.79847, 115.1287 37.79843, 1…
#> 8 130000 河北省 130600 保定市 地级市… (((115.4378 39.95016, 115.4435 39.9472, 11…
#> 9 130000 河北省 130700 张家口市… 地级市… (((114.8005 42.14749, 114.8045 42.14733, 1…
#> 10 130000 河北省 130800 承德市 地级市… (((117.7998 42.6137, 117.8 42.61273, 117.7…
#> # … with 361 more rows
计算城市的质心可以使用 st_centroid() 函数：
city %>%
st_centroid() -> city_centroid
# 展示质心的位置
library(leaflet)
library(leafem)
# devtools::install_git('https://gitee.com/tidyfriday/rstatatools.git')
library(rstatatools)
library(mapview)
leaflet() %>%
geoqmap(attribution = "绘制：微信公众号 RStata") %>%
addScaleBar() %>%
addLogo(img = "https://mdniceczx.oss-cn-beijing.aliyuncs.com/image_20201220175301.png") -> map
mapview(city_centroid, map = map,
zcol = "省", layer.name = "省份",
col.regions = paletteer::paletteer_d("ggsci::default_igv", 35))
各个城市质心的位置
省会城市有下面这些：
city_centroid %>%
dplyr::filter(市 %in% c("哈尔滨市", "长春市",
"沈阳市", "呼和浩特市",
"北京市", "天津市",
"石家庄市", "济南市",
"太原市", "西安市",
"兰州市", "乌鲁木齐市",
"拉萨市", "成都市",
"郑州市", "南京市",
"合肥市", "重庆市",
"贵阳市", "长沙市",
"上海市", "南昌市",
"杭州市", "福州市",
"昆明市", "南宁市",
"广州市", "海口市",
"香港特别行政区", "澳门特别行政区", "银川市", "武汉市", "西宁市")) -> sh_centroid
计算各个城市距离其所在省的省会城市的距离
各个城市距离其所在省的省会城市的距离
我们先以山东省为例，计算山东省的各个城市距离其省会城市济南市的距离：
city_centroid %>%
dplyr::filter(省 == "山东省") -> city_temp
sh_centroid %>%
dplyr::filter(省 == "山东省") -> sh_temp
st_distance(city_temp, sh_temp) -> mat
city_temp %>%
mutate(dist = mat[,1] %>%
units::set_units(km),
dist = as.numeric(dist))
#> Simple feature collection with 16 features and 6 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 115.6924 ymin: 34.91596 xmax: 121.9831 ymax: 37.62223
#> geographic CRS: CGCS_2000
#> # A tibble: 16 x 7
#> 省代码 省 市代码 市 类型 geometry dist
#> * <dbl> <chr> <dbl> <chr> <chr> <POINT [°]> <dbl>
#> 1 370000 山东省 370100 济南市 副省级市 (117.2148 36.6371) 0
#> 2 370000 山东省 370200 青岛市 副省级市 (120.1426 36.4532) 263.
#> 3 370000 山东省 370300 淄博市 地级市 (118.0533 36.60773) 75.1
#> 4 370000 山东省 370400 枣庄市 地级市 (117.3933 34.91596) 192.
#> 5 370000 山东省 370500 东营市 地级市 (118.6016 37.62223) 165.
#> 6 370000 山东省 370600 烟台市 地级市 (120.8006 37.24081) 326.
#> 7 370000 山东省 370700 潍坊市 地级市 (119.0711 36.54565) 166.
#> 8 370000 山东省 370800 济宁市 地级市 (116.735 35.3685) 147.
#> 9 370000 山东省 370900 泰安市 地级市 (117.0258 36.00098) 72.6
#> 10 370000 山东省 371000 威海市 地级市 (121.9831 37.12033) 428.
#> 11 370000 山东省 371100 日照市 地级市 (119.1351 35.58296) 209.
#> 12 370000 山东省 371300 临沂市 地级市 (118.2805 35.30909) 176.
#> 13 370000 山东省 371400 德州市 地级市 (116.6464 37.24772) 84.6
#> 14 370000 山东省 371500 聊城市 地级市 (115.8806 36.45852) 121.
#> 15 370000 山东省 371600 滨州市 地级市 (117.8377 37.53163) 114.
#> 16 370000 山东省 371700 菏泽市 地级市 (115.6924 35.15043) 215.
然后我们就可以循环计算每个城市距离其所在省省会城市的距离了：
# 首先准备一个空数据框
city_centroid %>%
slice(0) %>%
mutate(dist = 0) -> total
for (i in unique(city$省)) {
if(!i %in% c("中朝共有", "台湾省")){
city_centroid %>%
dplyr::filter(省 == i) -> city_temp
sh_centroid %>%
dplyr::filter(省 == i) -> sh_temp
st_distance(city_temp, sh_temp) -> mat
city_temp %>%
mutate(dist = mat[,1] %>%
units::set_units(km),
dist = as.numeric(dist)) -> temp
bind_rows(total, temp) -> total
}
}
bind_cols(
st_drop_geometry(total),
st_coordinates(total) %>%
as_tibble() %>%
set_names(c("经度", "纬度"))
) %>%
writexl::write_xlsx("各个城市距离其所在省的省会城市的距离(单位：km).xlsx")
计算各个城市距离最近的省会城市的距离
各个城市距离最近的省会城市的距离
那么另外如果想计算每个城市距离其最近的省会城市（不论是不是本省的）的距离怎么办呢？
这个就很简单了：
st_distance(city_centroid, sh_centroid) -> distmat3
city_centroid %>%
mutate(dist = distmat3[,1] %>%
units::set_units(km),
dist = as.numeric(dist)) -> total2
bind_cols(
st_drop_geometry(total2),
st_coordinates(total2) %>%
as_tibble() %>%
set_names(c("经度", "纬度"))
) %>%
writexl::write_xlsx("各个城市距离最近省会城市的距离(单位：km).xlsx")
计算各个区县距离其所在省的省会城市的距离
各个区县距离其所在省的省会城市的距离
这个问题也就很简单了：
read_sf("2019行政区划/县.shp") -> county
county %>%
st_centroid() -> county_centroid
# 首先准备一个空数据框
county_centroid %>%
slice(0) %>%
mutate(dist = 0) -> total3
for (i in unique(county$省)) {
if(!i %in% c("中朝共有", "台湾省") & !is.na(i)){
county_centroid %>%
dplyr::filter(省 == i) -> county_temp
sh_centroid %>%
dplyr::filter(省 == i) -> sh_temp
st_distance(county_temp, sh_temp) -> mat
county_temp %>%
mutate(dist = mat[,1] %>%
units::set_units(km),
dist = as.numeric(dist)) -> temp
bind_rows(total3, temp) -> total3
}
}
bind_cols(
st_drop_geometry(total3),
st_coordinates(total3) %>%
as_tibble() %>%
set_names(c("经度", "纬度"))
) %>%
rename(县 = NAME, 县代码 = PAC) %>%
writexl::write_xlsx("各个区县距离其所在省的省会城市的距离(单位：km).xlsx")
计算各个区县距离最近的省会城市的距离
各个区县距离最近的省会城市的距离
同样计算各个区县距离最近的省会城市的距离：
st_distance(county_centroid, sh_centroid) -> distmat4
county_centroid %>%
mutate(dist = distmat4[,1] %>%
units::set_units(km),
dist = as.numeric(dist)) -> total4
bind_cols(
st_drop_geometry(total4),
st_coordinates(total4) %>%
as_tibble() %>%
set_names(c("经度", "纬度"))
) %>%
writexl::write_xlsx("各个区县距离最近省会城市的距离(单位：km).xlsx")
这样我们就获得这四份数据啦：
- 各个城市距离其所在省的省会城市的距离(单位：km).xlsx
- 各个区县距离其所在省的省会城市的距离(单位：km).xlsx
- 各个城市距离最近省会城市的距离(单位：km).xlsx
- 各个区县距离最近省会城市的距离(单位：km).xlsx
