如何计算多个空间点间的距离并生成距离矩阵与数据框？

嗨，我来帮你搞定这个空间点距离计算的问题！你已经用sf包加载了北卡罗来纳州的空间数据，接下来我会分两种常见方法，帮你生成所有质心两两之间的距离矩阵和数据框，咱们一步步来：

前置准备（统一加载包与处理数据） #

先把基础工作做好，避免重复操作：

library(sf)  
library(tidyverse)  # 包含dplyr、expand.grid等实用工具

# 加载并预处理NC数据：转成WGS84坐标系，拆分多边形
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE) %>%
  st_transform(st_crs(4326)) %>%
  st_cast('POLYGON')

# 计算每个多边形的质心
centroids <- st_centroid(nc)

# 给每个区域生成唯一ID（避免重名带来的混淆）
nc$unique_id <- paste(nc$NAME, nc$FIPS, sep="_")

Option 1: 基于sf的常规空间距离 #

sf包的st_distance()会根据你的坐标系自动选择合适的距离计算逻辑，对于WGS84（EPSG:4326）经纬度数据，它会自动计算球面距离，非常省心。

生成距离数据框

# 计算两两质心的距离，转换为千米单位
dist_matrix <- st_distance(centroids, centroids) %>%
  units::set_units(km)  

# 生成所有"from-to"配对的数据框，并添加对应距离值
dist_df <- expand.grid(from = nc$NAME, 
                       to = nc$NAME) %>%
  mutate(distance = as.vector(dist_matrix))

小解释：expand.grid()会生成所有区域之间的配对组合，as.vector(dist_matrix)把矩阵转成向量，刚好能和配对一一对应上。

生成带标签的距离矩阵

# 把距离结果转成标准矩阵，设置行名和列名为唯一ID
dist_matrix <- as.matrix(dist_matrix)
rownames(dist_matrix) <- nc$unique_id
colnames(dist_matrix) <- nc$unique_id

这样矩阵的行和列就对应具体的区域，你可以直接查看任意两个区域的距离值。

Option 2: Haversine球面距离（经典球面距离计算） #

如果你想专门用经典的Haversine公式计算球面距离，可以用geosphere包的distHaversine()函数，它是专门针对经纬度点设计的，结果更贴合球面地球的实际距离：

library(geosphere)

# 把质心转成经纬度矩阵（geosphere包需要这种格式输入）
centroid_coords <- st_coordinates(centroids)

# 计算Haversine距离，默认单位是米，转成千米
haversine_matrix <- distHaversine(centroid_coords) / 1000
rownames(haversine_matrix) <- nc$unique_id
colnames(haversine_matrix) <- nc$unique_id

# 转成数据框格式
haversine_df <- expand.grid(from = nc$unique_id, 
                            to = nc$unique_id) %>%
  mutate(distance = as.vector(haversine_matrix))

小提示 #

• 如果你不需要重复计算，建议把距离矩阵或数据框用saveRDS()存起来，下次直接加载就行：saveRDS(dist_df, "nc_centroid_distances.rds")
• 当数据量很大时，st_distance()可能会慢一些，这时候可以考虑用并行计算或者简化数据。

备注：内容来源于stack exchange，提问作者stats_noob