R如何计算多个经纬度坐标点两两之间的公里距离

方案说明 #

我们使用球面坐标系下的Haversine公式计算两点距离，支持两种实现方式：使用成熟第三方包快速实现，或无依赖手动实现，输出结果单位均为公里。
注意：经纬度距离计算要求输入坐标顺序为经度在前，纬度在后，你的数据框列顺序已经符合要求，无需调整。

方式1：使用geosphere包实现（推荐，不易出错） #

步骤1：安装加载依赖包 #

# 首次使用先安装包，安装过可跳过此行
install.packages("geosphere")
# 加载包
library(geosphere)

步骤2：构造坐标数据 #

list1 <- data.frame(longitude = c(-3.0000,
                                   -1.3333,
                                   -6.5667), 
                     latitude = c(59.0000, 
                                  60.3333,  
                                  62.2500),
                      name= c('Orkney', 
                            'Shetlands', 
                               'Faroe'))

步骤3：计算距离 #

输出距离矩阵（快速查表）

矩阵行和列均为岛屿名称，交叉位置的数值就是对应两个地点的距离：

# 提取经纬度列
coord <- list1[, c("longitude", "latitude")]
# 计算距离矩阵，米转公里
dist_matrix <- distHaversine(coord) / 1000
# 绑定岛屿名称
rownames(dist_matrix) <- list1$name
colnames(dist_matrix) <- list1$name
# 打印结果
print(dist_matrix)

输出成对结果表（直观展示所有两两组合）

# 生成所有不重复的岛屿两两组合
pair_comb <- t(combn(list1$name, 2))
# 逐对计算距离
pair_dist <- apply(pair_comb, 1, function(x){
  p1_coord <- list1[list1$name == x[1], c("longitude", "latitude")]
  p2_coord <- list1[list1$name == x[2], c("longitude", "latitude")]
  return(distHaversine(p1_coord, p2_coord)/1000)
})
# 生成结果表，保留2位小数
result_df <- data.frame(
  地点A = pair_comb[,1],
  地点B = pair_comb[,2],
  距离_公里 = round(pair_dist, 2)
)
# 打印结果
print(result_df)

输出示例：

地点A     地点B 距离_公里
1    Orkney Shetlands    169.45
2    Orkney     Faroe    393.01
3 Shetlands     Faroe    286.27

方式2：无第三方依赖手动实现 #

不想安装额外包可以直接使用自定义Haversine公式函数计算：

# 自定义Haversine距离函数，输出单位公里
haversine_km <- function(lon1, lat1, lon2, lat2){
  R <- 6371 # 地球平均半径，单位公里
  d_lat <- (lat2 - lat1) * pi / 180
  d_lon <- (lon2 - lon1) * pi / 180
  a <- sin(d_lat/2)^2 + cos(lat1 * pi / 180) * cos(lat2 * pi / 180) * sin(d_lon/2)^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  return(R * c)
}

# 生成两两组合结果
pair_comb <- t(combn(list1$name, 2))
pair_dist <- apply(pair_comb, 1, function(x){
  p1 <- list1[list1$name == x[1], ]
  p2 <- list1[list1$name == x[2], ]
  return(haversine_km(p1$longitude, p1$latitude, p2$longitude, p2$latitude))
})
result_df <- data.frame(
  地点A = pair_comb[,1],
  地点B = pair_comb[,2],
  距离_公里 = round(pair_dist, 2)
)
print(result_df)

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