R 语言：如何计算各国各地区质心经纬度、质心间距离矩阵及邻接矩阵

1149 字 约 5 分钟 课程 距离与边界地形与地理R 语言编程技巧

最近有小伙伴需要各国、各地区间的距离矩阵和邻接矩阵。今天我们一起使用 R 语言来计算一下～

计算各地区质心经纬度

world_high_resolution_mill.geo.json 文件是我根据国内的境界要求调整的全球地图矢量数据。不过是米勒投影的：

library(sf)
library(tidyverse)
read_sf("world_high_resolution_mill.geo.json") -> wdmp

wdmp
#> Simple feature collection with 215 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -20015110 ymin: -7379023 xmax: 20015110 ymax: 12641180
#> Projected CRS: World_Miller_Cylindrical
#> # A tibble: 215 × 4
#>    name       name_en       code                                        geometry
#>    <chr>      <chr>         <chr>                             <MULTIPOLYGON [m]>
#>  1 阿尔及利亚 ALGERIA       DZA   (((-245324.1 4068138, -223332.1 4066332, -196…
#>  2 列支敦士登 Liechtenstein LIE   (((1067880 5656162, 1059093 5662487, 1055562 …
#>  3 埃及       EGYPT         EGY   (((3961141 2619833, 3915419 2577702, 3886419 …
#>  4 孟加拉国   BANGLADESH    BGD   (((10297641 2482880, 10299627 2472235, 102967…
#>  5 尼日尔     Niger         NER   (((400866.4 1306467, 394566.5 1311739, 385678…
#>  6 卡塔尔     QATAR         QAT   (((5699872 2791895, 5682040 2775874, 5663103 …
#>  7 纳米比亚   NAMIBIA       NAM   (((2223927 -2808902, 2223888 -2883519, 222377…
#>  8 保加利亚   Bulgaria      BGR   (((2486312 5005203, 2497552 5009048, 2505639 …
#>  9 玻利维亚   Bolivia       BOL   (((-6966045 -2514366, -6982458 -2504665, -698…
#> 10 加纳       Ghana         GHA   (((-299106.5 1058207, -309033.4 1069084, -306…
#> # ℹ 205 more rows
plot(wdmp[1])

使用 st_centroid() 计算质心经纬度再转换成 WGS84 坐标系即可：

# 质心
wdmp %>%
  st_make_valid() %>%
  st_centroid() -> wdmpcen

wdmpcen %>%
  st_transform(4326) -> wdmpcen

wdmpcen
#> Simple feature collection with 215 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -175.2887 ymin: -54.50234 xmax: 172.3443 ymax: 78.88678
#> Geodetic CRS:  WGS 84
#> # A tibble: 215 × 4
#>    name       name_en       code               geometry
#>  * <chr>      <chr>         <chr>           <POINT [°]>
#>  1 阿尔及利亚 ALGERIA       DZA     (2.644491 28.32924)
#>  2 列支敦士登 Liechtenstein LIE     (9.556793 47.16157)
#>  3 埃及       EGYPT         EGY     (29.77654 26.62068)
#>  4 孟加拉国   BANGLADESH    BGD     (90.26146 23.87647)
#>  5 尼日尔     Niger         NER     (9.418163 17.46307)
#>  6 卡塔尔     QATAR         QAT     (51.19399 25.29545)
#>  7 纳米比亚   NAMIBIA       NAM     (17.2231 -22.21198)
#>  8 保加利亚   Bulgaria      BGR     (25.23909 42.77278)
#>  9 玻利维亚   Bolivia       BOL   (-64.65832 -16.75428)
#> 10 加纳       Ghana         GHA    (-1.208645 7.968293)
#> # ℹ 205 more rows

再保存成 dta 文件：

bind_cols(
  wdmpcen %>%
    st_drop_geometry(),
  wdmpcen %>%
    st_coordinates() %>%
    as_tibble() %>%
    set_names("lon", "lat")
) %>%
  haven::write_dta("各国各地区质心经纬度.dta")

计算质心间距离矩阵

使用 st_distance() 即可计算质心间距离矩阵：

wdmpcen %>%
  st_distance(wdmpcen) -> distmat

把单位转换成 km 再保存成 dta 文件：

as_tibble(distmat/1000) %>%
  mutate_all(as.numeric) %>%
  set_names(wdmpcen$code) %>%
  mutate(code = wdmpcen$code) %>%
  select(code, everything()) -> distdf

distdf
#> # A tibble: 215 × 216
#>    code    DZA    LIE    EGY    BGD   NER    QAT   NAM    BGR    BOL   GHA
#>    <chr> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl> <dbl>
#>  1 DZA      0   2178.  2678.  8558. 1392.  4798. 5834.  2585.  8800. 2300.
#>  2 LIE   2178.     0   2889.  7405. 3302.  4381. 7752.  1324. 10226. 4476.
#>  3 EGY   2678.  2889.     0   6033. 2328.  2144. 5596.  1843. 11260. 3870.
#>  4 BGD   8558.  7405.  6033.     0  8330.  3939. 9408.  6234. 17287. 9798.
#>  5 NER   1392.  3302.  2328.  8330.    0   4394. 4493.  3187.  8957. 1562.
#>  6 QAT   4798.  4381.  2144.  3939. 4394.     0  6429.  3059. 13348. 5864.
#>  7 NAM   5834.  7752.  5596.  9408. 4493.  6429.    0   7273.  8502. 3911.
#>  8 BGR   2585.  1324.  1843.  6234. 3187.  3059. 7273.     0  11255. 4652.
#>  9 BOL   8800. 10226. 11260. 17287. 8957. 13348. 8502. 11255.     0  7497.
#> 10 GHA   2300.  4476.  3870.  9798. 1562.  5864. 3911.  4652.  7497.    0
#> # ℹ 205 more rows
#> # ℹ 205 more variables: PAK <dbl>, JOR <dbl>, LBR <dbl>, LBY <dbl>, MYS <dbl>,
#> #   PRI <dbl>, MYT <dbl>, PRK <dbl>, TZA <dbl>, PRT <dbl>, TMP <dbl>,
#> #   SPM <dbl>, PRY <dbl>, SAU <dbl>, LBN <dbl>, SVN <dbl>, BFA <dbl>,
#> #   SVK <dbl>, MRT <dbl>, HRV <dbl>, CHL <dbl>, CHN <dbl>, KNA <dbl>,
#> #   JAM <dbl>, DJI <dbl>, GIN <dbl>, FIN <dbl>, URY <dbl>, NPL <dbl>,
#> #   MAR <dbl>, YEM <dbl>, ZAF <dbl>, NIC <dbl>, PHL <dbl>, TGO <dbl>, …
distdf %>%
  haven::write_dta("各国间距离矩阵.dta")

计算邻接矩阵

邻接矩阵的计算通常应该使用 st_touches()，不过由于地图数据通常不是非常精确的，所以使用 st_intersects() 也是一样的：

wdmpcen %>%
  st_intersects(wdmpcen, sparse = F) -> touchmat

# 对角线替换成 FALSE
diag(touchmat) <- F
as_tibble(touchmat) %>%
  mutate_all(as.numeric) %>%
  set_names(wdmpcen$code) %>%
  mutate(code = wdmpcen$code) %>%
  select(code, everything()) -> touchdf

touchdf
#> # A tibble: 215 × 216
#>    code    DZA   LIE   EGY   BGD   NER   QAT   NAM   BGR   BOL   GHA   PAK   JOR
#>    <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#>  1 DZA       0     0     0     0     0     0     0     0     0     0     0     0
#>  2 LIE       0     0     0     0     0     0     0     0     0     0     0     0
#>  3 EGY       0     0     0     0     0     0     0     0     0     0     0     0
#>  4 BGD       0     0     0     0     0     0     0     0     0     0     0     0
#>  5 NER       0     0     0     0     0     0     0     0     0     0     0     0
#>  6 QAT       0     0     0     0     0     0     0     0     0     0     0     0
#>  7 NAM       0     0     0     0     0     0     0     0     0     0     0     0
#>  8 BGR       0     0     0     0     0     0     0     0     0     0     0     0
#>  9 BOL       0     0     0     0     0     0     0     0     0     0     0     0
#> 10 GHA       0     0     0     0     0     0     0     0     0     0     0     0
#> # ℹ 205 more rows
#> # ℹ 203 more variables: LBR <dbl>, LBY <dbl>, MYS <dbl>, PRI <dbl>, MYT <dbl>,
#> #   PRK <dbl>, TZA <dbl>, PRT <dbl>, TMP <dbl>, SPM <dbl>, PRY <dbl>,
#> #   SAU <dbl>, LBN <dbl>, SVN <dbl>, BFA <dbl>, SVK <dbl>, MRT <dbl>,
#> #   HRV <dbl>, CHL <dbl>, CHN <dbl>, KNA <dbl>, JAM <dbl>, DJI <dbl>,
#> #   GIN <dbl>, FIN <dbl>, URY <dbl>, NPL <dbl>, MAR <dbl>, YEM <dbl>,
#> #   ZAF <dbl>, NIC <dbl>, PHL <dbl>, TGO <dbl>, VIR <dbl>, SYR <dbl>, …
touchdf %>%
  haven::write_dta("各国间邻矩阵.dta")

点击这里跳转到 RStata 短书平台获取附件：R 语言：如何计算各国各地区质心经纬度、质心间距离矩阵及邻接矩阵

相关文章

评论