使用data.table语法计算数据表中两行间的地理距离

没问题，这事儿用data.table的分组和shift()函数就能高效搞定，完全不用写循环。下面是具体的实现方案：

1. 准备工作 #

首先确保你加载了data.table包，如果需要更便捷的地理距离计算，还可以用geosphere包（当然也可以自己实现公式）。先构造你的示例数据：

library(data.table)
library(geosphere)

# 构造示例data.table
dt <- data.table(
  Lat = c(52.21808, 52.25882, 52.24347, 52.21935, 52.25808, 52.24541),
  Lon = c(20.96675, 20.89850, 21.08460, 20.97186, 20.89790, 21.08522),
  Time = as.POSIXct(c("2018-04-20 21:27:26", "2018-04-20 21:27:23", "2018-04-20 21:27:27",
                      "2018-04-20 21:28:31", "2018-04-20 21:28:32", "2018-04-20 21:28:36")),
  Bus = c(3,8,1,3,8,1)
)

2. 核心操作：分组排序+计算相邻距离 #

首先要对每个Bus分组内的记录按Time排序，这样相邻的记录才是时间顺序的；然后用shift()函数获取每组内上一行的经纬度，最后计算距离：

# 先按Bus分组，每组内按Time排序（确保时间顺序正确）
dt <- dt[order(Bus, Time)]

# 计算同一Bus下相邻记录的地理距离（单位：米）
dt[, distance := distHaversine(
  cbind(Lon, Lat),  # 当前行的经纬度（注意geosphere函数要求先Lon后Lat）
  cbind(shift(Lon), shift(Lat))  # 上一行的经纬度
), by = Bus]

结果说明 #

运行后你会得到这样的结果：

Lat      Lon                Time Bus  distance
1: 52.24347 21.08460 2018-04-20 21:27:27   1        NA
2: 52.24541 21.08522 2018-04-20 21:28:36   1 223.62354
3: 52.21808 20.96675 2018-04-20 21:27:26   3        NA
4: 52.21935 20.97186 2018-04-20 21:28:31   3  63.76674
5: 52.25882 20.89850 2018-04-20 21:27:23   8        NA
6: 52.25808 20.89790 2018-04-20 21:28:32   8  87.16864

• 每组的第一条记录因为没有上一行，distance会是NA，这是合理的
• distHaversine()返回的距离单位是米，如果需要公里，直接除以1000即可：dt[, distance := distance / 1000]

3. 不用geosphere包？自己实现Haversine公式 #

如果你不想额外加载包，也可以自己写Haversine公式来计算：

# 自定义Haversine距离函数（单位：米）
haversine <- function(lon1, lat1, lon2, lat2) {
  # 转换为弧度
  lon1 <- lon1 * pi / 180
  lat1 <- lat1 * pi / 180
  lon2 <- lon2 * pi / 180
  lat2 <- lat2 * pi / 180
  
  dlon <- lon2 - lon1
  dlat <- lat2 - lat1
  
  a <- sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  
  6371000 * c  # 地球平均半径，单位米
}

# 用自定义函数计算距离
dt[, distance := haversine(Lon, Lat, shift(Lon), shift(Lat)), by = Bus]

这样也能得到和之前完全一致的结果。

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