使用R的geosphere包计算经纬度点对距离时无法识别点对标签

样本数据 #

GridID  Longitude   Latitude
A13 -79.320029  44.033275
A15 -79.18051   44.056613
A17 -79.057251  44.026426
A8  -79.6297686 44.0395171
B15 -79.180559  44.017748
B9  -79.565778  44.007081
C20 -78.8883705 43.935996

期望输出 #

需要包含点对标识和以米为单位的距离，格式示例：

Point1  Point2  Distance
A13     A13     1452.3
A13     A15     562.1
A13     A17     1423
A13     A8      432
A13     B15     673.23
A13     B9      2345
A13     C20     123

现有问题 #

使用geosphere包的distm函数计算得到的距离矩阵，提取下三角后仅得到数值，无法对应具体点对：

data <- read.csv("samplepoints.csv", header=TRUE)
points <- data %>% select(GridID, Longitude, Latitude)

# 原代码存在语法错误：缺少一个闭合括号
geoDist <- distm(as.matrix(points[2:3]), fun=distGeo)
df <- data.frame(Distance = geoDist[lower.tri(geoDist)])

输出仅包含距离数值，无点对标签：

Distance
1     11478.59
2     21079.59
3     24837.51
4     11313.88
5     19917.70
6     36278.19

解决方案 #

方法1：生成全点对组合后计算距离 #

通过生成所有点的配对组合，再逐对计算距离，直接得到带标签的结果：

library(geosphere)
library(dplyr)

# 读取数据
data <- read.csv("samplepoints.csv", header=TRUE)
points <- data %>% select(GridID, Longitude, Latitude)

# 生成所有点对的笛卡尔积
point_pairs <- expand.grid(Point1 = points$GridID, Point2 = points$GridID, stringsAsFactors = FALSE)

# 合并两组点的坐标信息
point_pairs <- point_pairs %>%
  left_join(points, by = c("Point1" = "GridID")) %>%
  rename(Longitude1 = Longitude, Latitude1 = Latitude) %>%
  left_join(points, by = c("Point2" = "GridID")) %>%
  rename(Longitude2 = Longitude, Latitude2 = Latitude)

# 计算每对点的球面距离（单位：米）
point_pairs$Distance <- distGeo(
  point_pairs %>% select(Longitude1, Latitude1) %>% as.matrix(),
  point_pairs %>% select(Longitude2, Latitude2) %>% as.matrix()
)

# 保留目标列
result <- point_pairs %>% select(Point1, Point2, Distance)

# 查看结果
print(result)

方法2：将距离矩阵转换为带标签的长格式数据 #

如果已经生成了距离矩阵，可通过tidyr工具将宽矩阵转换为带点对标签的长格式：

library(geosphere)
library(dplyr)
library(tidyr)

data <- read.csv("samplepoints.csv", header=TRUE)
points <- data %>% select(GridID, Longitude, Latitude)

# 生成距离矩阵并设置行列名为GridID
geoDist <- distm(as.matrix(points[2:3]), fun=distGeo)
rownames(geoDist) <- points$GridID
colnames(geoDist) <- points$GridID

# 转换为长格式数据框
result <- as.data.frame(geoDist) %>%
  rownames_to_column(var = "Point1") %>%
  pivot_longer(cols = -Point1, names_to = "Point2", values_to = "Distance")

# 查看结果
print(result)

可选：过滤非重复点对 #

如果不需要自身配对和重复配对（如A13-A15与A15-A13只保留一个），可添加过滤逻辑：

result_unique <- result %>%
  filter(Point1 != Point2, Point1 < Point2)

内容的提问来源于stack exchange，提问作者Katherine Chau