为一个 sf 对象创建距离矩阵。

Question

我有一个名为sf的对象，它是一个划分为行政区的地图。我想计算每个行政区的重心（使用st_point_on_surface），然后计算每个重心之间的相对距离，就像一个距离矩阵，在此基础上可以执行计算（例如保留在某个半径内的距离）并获得一个带有每个行政区标识符和符合条件的列表的数据库。

抱歉，缺乏可重现代码。

最简单的方法是什么？

提前感谢。

- MCS

---

你好！您可能想要查看以下链接（也许是整本书）：1）质心的几何运算；2）空间距离。 - agila

我不确定 st_point_on_surface 是否会给你质心。你是在使用它来避免质心实际上位于凹多边形外部的情况吗？ - camille

Answer 1

如果你想留在sf，可以尝试以下方法：

library(sf)

# I got my Irak map from UN's OCHA, see: https://data.humdata.org/dataset/iraq-admin-level-1-boundaries
shp_irak01 <- st_read(dsn = "./irq-administrative-divisions-shapefiles/irq_admbnda_adm1_cso_20190603.shp")

#The centroids:
shp_centroid <- st_point_on_surface(x = shp_irak01)

#The euclidian distance matrix:
mtx_distance <- st_distance(shp_centroid, shp_centroid)
mtx_distance

Units: [m]
          [,1]     [,2]     [,3]     [,4]      [,5]     [,6]      [,7]      [,8]     [,9]    [,10]    [,11]    [,...until18]
 [1,]      0.0 665051.7 484956.1 295870.9 383276.14 449273.4 309551.40 277419.82 338721.4 458264.2 438887.2
 [2,] 665051.7      0.0 206781.0 408086.2 287295.37 624344.0 378896.62 451360.03 494695.5 860740.2 751176.7
 [3,] 484956.1 206781.0      0.0 205871.6 160931.72 570092.5 251440.27 332638.38 410021.2 766396.3 672788.5
 [4,] 295870.9 408086.2 205871.6      0.0 189849.79 526919.9 203679.02 260171.36 360077.9 661294.1 592026.7
 [5,] 383276.1 287295.4 160931.7 189849.8      0.00 410667.4  96030.18 175862.79 249100.0 607665.1 512029.2
 [6,] 449273.4 624344.0 570092.5 526919.9 410667.36      0.0 335212.34 266829.78 168229.3 261881.0 141057.5
 [7,] 309551.4 378896.6 251440.3 203679.0  96030.18 335212.3      0.00  81209.73 167546.7 514984.8 424189.9
 [8,] 277419.8 451360.0 332638.4 260171.4 175862.79 266829.8  81209.73      0.00 100377.7 433844.3 345162.4
 [9,] 338721.4 494695.5 410021.2 360077.9 249100.03 168229.3 167546.72 100377.65      0.0 366935.5 264039.4
[10,] 458264.2 860740.2 766396.3 661294.1 607665.14 261881.0 514984.82 433844.32 366935.5      0.0 120940.2
[11,] 438887.2 751176.7 672788.5 592026.7 512029.23 141057.5 424189.92 345162.43 264039.4 120940.2      0.0
[12,] 236774.0 434018.5 275796.9 160698.5 147862.67 367278.5  78884.14 102143.31 202413.5 509250.3 433156.4
[13,] 339045.0 649972.5 556320.0 469992.6 396524.64 122292.9 305687.02 225374.60 155304.3 212297.2 122084.4
[14,] 552966.6 202333.0 244743.0 370890.8 186402.59 424966.6 243458.18 294501.93 310734.5 669733.2 556160.8
[15,] 325331.4 801693.1 681103.2 551881.7 529148.71 276836.1 433248.41 353886.47 314515.2 135975.2 171322.4
[16,] 264036.5 613728.6 499989.6 396476.1 344132.89 185842.8 249137.50 168288.85 133007.6 269034.3 198147.7
[17,] 504238.7 168077.1 131432.6 277522.5 121527.11 479013.4 210879.43 283915.79 335117.1 701926.8 597103.5
[18,] 404859.5 326577.9 258506.4 284146.2 108878.33 315038.5 100818.86 138534.47 168707.2 535609.5 430059.8

- Nicolás Velásquez

---

@ Nicolas：非常好的回答！有没有办法将这个距离从“欧几里得距离”改为“哈弗辛距离”？ - stats_noob

我不确定如何使用 {sf} 完成此操作，但是通过使用 {geosphere} 包，您可以使用 geosphere::distm(x = st_coordinates(shp_centroid), y = st_coordinates(shp_centroid), fun = geosphere::distHaversine) 来实现。 - Nicolás Velasquez PhD

@ Nicolas：谢谢您的回复！这是否意味着我只需替换以下行？ - stats_noob

mtx_distance <- geosphere::distm(x = st_coordinates(shp_centroid), y = st_coordinates(shp_centroid), fun = geosphere::distHaversine) - stats_noob

是的！ 您还可以加载geosphere库，然后在不使用双引号调用的情况下运行distm()命令。 library(geosphere) mtx_distance_harvesine <- distm(x = st_coordinates(shp_centroid), y = st_coordinates(shp_centroid), fun = geosphere::distHaversine) - Nicolás Velasquez PhD