如何在R中循环使用pointDistance()按个体计算点位距离

问题描述 #

现有包含四列的数据框：个体标识（individual）、观测日期（date）、经度（longitude）、纬度（latitude），需要按个体分组，计算每个个体的第一个位置坐标与其余位置坐标之间的距离。

原尝试代码返回NULL，无法得到预期结果。

原尝试代码 #

library(raster)
names(spdf)
[1] "individual" "date" "utm.longitude" "utm.latitude" "longitude"

coordinates(spdf) <- c("utm.longitude", "utm.latitude")
proj4string(spdf) <- CRS("+proj=utm +zone=33 ellps=WGS84")
spdf$date <- as.POSIXct(strptime(spdf$date, format="%Y-%m-%d %H:%M"))

test <- for(i in unique(spdf$individual)) {
pointDistance(spdf, lonlat=FALSE)
}
test
# NULL

附数据结构 #

dput(spdf)
structure(list(individual = c(1L, 2L, 3L, 1L, 3L, 2L, 1L, 3L, 
1L, 2L, 1L), date = structure(c(1629141240, 1631819700, 1631992680, 
1650136980, 1651952340, 1654976400, 1655063160, 1659556080, 1661111520, 
1681914060, 1682519280), class = c("POSIXct", "POSIXt"), tzone = ""), 
    longitude = c(455531L, 455529L, 455528L, 455528L, 455548L, 
    455548L, 455541L, 455530L, 455526L, 455555L, 455524L), latitude = c(5557649L, 
    5557649L, 5557647L, 5557646L, 5557640L, 5557636L, 5557638L, 
    5557645L, 5557648L, 5557720L, 5557651L)), row.names = c(NA, 
11L), class = "data.frame")

str(spdf)
# 'data.frame':    18 obs. of  4 variables:
#  $ individual: int  1 2 3 1 3 2 1 3 1 2 ...
#  $ date      : POSIXct, format: "2021-08-16 21:14:00" "2021-09-16 21:15:00" "2021-09-18 21:18:00" "2022-04-16 21:23:00" ...
#  $ longitude : int  455531 455529 455528 455528 455548 455548 455541 455530 455526 455555 ...
#  $ latitude  : int  5557649 5557649 5557647 5557646 5557640 5557636 5557638 5557645 5557648 5557720 ..

问题分析 #

原代码存在三个核心问题：

• 循环中未针对个体子集计算，直接使用整个spdf，导致计算所有点间距离而非个体内首个位置与其他位置的距离。
• for循环默认不返回结果，未手动存储每次循环的计算值。
• 数据结构中不存在utm.longitude字段，原代码的坐标转换逻辑错误。

正确实现方案 #

方法1：基础R循环实现 #

library(raster)

# 加载数据
spdf <- structure(list(individual = c(1L, 2L, 3L, 1L, 3L, 2L, 1L, 3L, 
1L, 2L, 1L), date = structure(c(1629141240, 1631819700, 1631992680, 
1650136980, 1651952340, 1654976400, 1655063160, 1659556080, 1661111520, 
1681914060, 1682519280), class = c("POSIXct", "POSIXt"), tzone = ""), 
    longitude = c(455531L, 455529L, 455528L, 455528L, 455548L, 
    455548L, 455541L, 455530L, 455526L, 455555L, 455524L), latitude = c(5557649L, 
    5557649L, 5557647L, 5557646L, 5557640L, 5557636L, 5557638L, 
    5557645L, 5557648L, 5557720L, 5557651L)), row.names = c(NA, 
11L), class = "data.frame")

# 转换为空间数据框并设置UTM投影
coordinates(spdf) <- c("longitude", "latitude")
proj4string(spdf) <- CRS("+proj=utm +zone=33 ellps=WGS84")

# 创建空列表存储结果
result_list <- list()

# 遍历每个个体
for (i in unique(spdf$individual)) {
    # 提取当前个体的所有位置
    ind_subset <- spdf[spdf$individual == i, ]
    # 获取个体的第一个位置
    first_point <- ind_subset[1, ]
    # 计算首个位置与其余位置的距离（UTM投影设置lonlat=FALSE）
    distances <- pointDistance(first_point, ind_subset, lonlat = FALSE)
    # 整理结果并存入列表
    result_list[[as.character(i)]] <- data.frame(
        individual = i,
        date = ind_subset$date,
        distance_from_first = distances
    )
}

# 合并所有结果为一个数据框
final_result <- do.call(rbind, result_list)
print(final_result)

方法2：dplyr分组实现（更简洁） #

library(raster)
library(dplyr)

# 加载数据并转为普通数据框
df <- structure(list(individual = c(1L, 2L, 3L, 1L, 3L, 2L, 1L, 3L, 
1L, 2L, 1L), date = structure(c(1629141240, 1631819700, 1631992680, 
1650136980, 1651952340, 1654976400, 1655063160, 1659556080, 1661111520, 
1681914060, 1682519280), class = c("POSIXct", "POSIXt"), tzone = ""), 
    longitude = c(455531L, 455529L, 455528L, 455528L, 455548L, 
    455548L, 455541L, 455530L, 455526L, 455555L, 455524L), latitude = c(5557649L, 
    5557649L, 5557647L, 5557646L, 5557640L, 5557636L, 5557638L, 
    5557645L, 5557648L, 5557720L, 5557651L)), row.names = c(NA, 
11L), class = "data.frame")

# 分组计算距离
final_result <- df %>%
    group_by(individual) %>%
    mutate(
        # 获取当前个体的首个经纬度
        first_lon = first(longitude),
        first_lat = first(latitude),
        # 计算与首个位置的距离
        distance_from_first = pointDistance(
            cbind(first_lon, first_lat),
            cbind(longitude, latitude),
            lonlat = FALSE
        )
    ) %>%
    select(-first_lon, -first_lat) %>%
    ungroup()

print(final_result)

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