使用ltraj分析轨迹得到异常低输出的原因排查
2026-6-25
鸟类追踪数据计算异常问题排查
问题说明
处理鸟类追踪数据时，使用adehabitatLT包计算轨迹参数（速度、距离），代码可正常运行，但结果远低于预期（研究对象为飞鸟，预期速度1-100km/h，最大距离2-3000km）。例如第1到第2条记录的dist预期约2099000米，实际输出仅1.76左右，数值异常。
数据结构
> head(springsub) X DateTime id lat lon 1 1 2022-05-03 19:39:00 N34246@096AFE12 45.7309083 13.3009833 2 2 2022-05-03 21:39:00 N34246@096AFE12 46.3943317 14.9383667 3 3 2022-05-03 23:39:00 N34246@096AFE12 47.3849317 16.5986833 4 4 2022-05-04 01:39:00 N34246@096AFE12 48.44073 18.0876833 5 5 2022-05-04 03:39:00 N34246@096AFE12 48.73939 18.5627 6 6 2022-05-04 05:39:00 N34246@096AFE12 48.739305 18.5627 > str(springsub) 'data.frame': 21053 obs. of 12 variables: $ X : int 1 2 3 4 5 6 7 8 9 10 ... $ DateTime : POSIXct, format: "2022-05-03 19:39:00" "2022-05-03 21:39:00" "2022-05-03 23:39:00" "2022-05-04 01:39:00" ... $ id : chr "N34246@096AFE12" "N34246@096AFE12" "N34246@096AFE12" "N34246@096AFE12" ... $ lat : chr "45.7309083" "46.3943317" "47.3849317" "48.44073" ... $ lon : chr "13.3009833" "14.9383667" "16.5986833" "18.0876833" ...
运行代码
library(adehabitatLT) library(sp) library(sf) spring_sf <- st_as_sf(springsub, coords = c("lon", "lat"), crs = 4326) individual_ids <- unique(springsub$id) # get unique individual IDs ltraj_list <- list() # stores ltraj obj for each id # loop plus some debugging that was needed for it to work properly for (id in individual_ids) { cat("Processing individual:", id, "\n") individual_data <- spring_sf[spring_sf$id == id, ] cat("Number of unique dates:", length(unique(individual_data$DateTime)), "\n") if (anyDuplicated(individual_data$DateTime) > 0) { individual_data$DateTime <- individual_data$DateTime + seq(0, length.out = nrow(individual_data), by = 1) } individual_coords <- st_coordinates(individual_data) tryCatch({ individual_ltraj <- as.ltraj(individual_coords, date = individual_data$DateTime, id = individual_data$id, typeII = TRUE) ltraj_list[[id]] <- individual_ltraj}, error = function(e) { cat("Error for individual:", id, "\n") print(e) }) } ld_results <- list() for (i in seq_along(ltraj_list)) { individual_ltraj <- ltraj_list[[i]] ld_result <- ld(individual_ltraj) ld_results[[i]] <- ld_result } ltraj_metrics <- do.call(rbind, ld_results) # adding relevant metrics # velocity ltraj_metrics$velocity <- ltraj_metrics$dist / ltraj_metrics$dt # euclidean distance and total distance euclidean_distance <- aggregate(dist ~ id, data = ltraj_metrics, FUN = sum) total_distance <- aggregate(dist ~ id, data = ltraj_metrics, FUN = sum) # straightness sum_R2n_per_ID <- aggregate(R2n ~ id, data = ltraj_metrics, FUN = sum) sum_R2n_per_ID <- aggregate(R2n ~ id, data = ltraj_metrics, FUN = sum) R2n_data <- merge(sum_R2n_per_ID, total_distance, by = "id") R2n_data$straightness <- R2n_data$R2n / R2n_data$dist # turning angle on x (abs.angle) ltraj_metrics$turning_angle_x <- c(NA, diff(ltraj_metrics$abs.angle)) # turning angle on y (rel.angle) ltraj_metrics$turning_angle_y <- c(NA, diff(ltraj_metrics$rel.angle))
异常输出示例
> head(ltraj_metrics) x y date dx dy dist dt R2n abs.angle rel.angle id 1 13.30098 45.73091 2022-05-03 19:39:00 1.6373834 0.6634234 1.7666790898 7200 0.000000 0.3849578 NA N34246@096AFE12 2 14.93837 46.39433 2022-05-03 21:39:00 1.6603166 0.9906000 1.9333751763 7200 3.121155 0.5379402 0.15298242 N34246@096AFE12 3 16.59868 47.38493 2022-05-03 23:39:00 1.4890000 1.0557983 1.8253303948 7200 13.610619 0.6167842 0.07884398 N34246@096AFE12 4 18.08768 48.44073 2022-05-04 01:39:00 0.4750167 0.2986600 0.5611048573 7200 30.255631 0.5612812 -0.05550300 N34246@096AFE12 5 18.56270 48.73939 2022-05-04 03:39:00 0.0000000 -0.0000850 0.0000850000 7200 36.736625 -1.5707963 -2.13207754 N34246@096AFE12 6 18.56270 48.73931 2022-05-04 05:39:00 0.0000833 0.0001100 0.0001379815 7200 36.736113 0.9226569 2.49345322 N34246@096AFE12 burst pkey velocity turning_angle_x turning_angle_y 1 N34246@096AFE12 N34246@096AFE12.2022-05-03 19:39:00 2.453721e-04 NA NA 2 N34246@096AFE12 N34246@096AFE12.2022-05-03 21:39:00 2.685243e-04 0.15298242 NA 3 N34246@096AFE12 N34246@096AFE12.2022-05-03 23:39:00 2.535181e-04 0.07884398 -0.07413844 4 N34246@096AFE12 N34246@096AFE12.2022-05-04 01:39:00 7.793123e-05 -0.05550300 -0.13434698 5 N34246@096AFE12 N34246@096AFE12.2022-05-04 03:39:00 1.180556e-08 -2.13207754 -2.07657455 6 N34246@096AFE12 N34246@096AFE12.2022-05-04 05:39:00 1.916410e-08 2.49345322 4.62553077
问题原因及修正方案
核心错误点
- 经纬度数据类型错误：从
str(springsub)输出可见，lat和lon列是字符型（chr），而非数值型。直接使用字符型数据计算会导致距离计算逻辑错误，得到的结果并非实际地理距离。 - 未使用平面坐标系计算：代码中使用WGS84地理坐标系（EPSG:4326）的坐标直接创建
ltraj对象，adehabitatLT默认用欧几里得距离计算，此时单位是度，而非米。1度纬度约等于111km，经度则随纬度变化，直接得到的dist是度的差值，并非实际米数。
修正步骤
1. 转换经纬度为数值型
# 将字符型经纬度转为数值型 springsub$lat <- as.numeric(springsub$lat) springsub$lon <- as.numeric(springsub$lon)
2. 投影到平面坐标系
将地理坐标转换为以米为单位的平面坐标系（如UTM投影，需根据数据所在区域选择对应带号，示例中欧洲东部区域使用EPSG:32633）：
# 创建sf对象并转换投影 spring_sf <- st_as_sf(springsub, coords = c("lon", "lat"), crs = 4326) %>% st_transform(crs = 32633) # 替换为数据所在区域的UTM带号EPSG代码
3. 调整速度计算单位
dt的单位是秒，dist转换为米后，计算得到的速度是m/s，如需转换为预期的km/h，需乘以3.6：
ltraj_metrics$velocity_kmh <- (ltraj_metrics$dist / ltraj_metrics$dt) * 3.6
修正后代码示例
library(adehabitatLT) library(sf) # 修正数据类型 springsub$lat <- as.numeric(springsub$lat) springsub$lon <- as.numeric(springsub$lon) # 转换为平面坐标系（UTM示例） spring_sf <- st_as_sf(springsub, coords = c("lon", "lat"), crs = 4326) %>% st_transform(crs = 32633) individual_ids <- unique(springsub$id) ltraj_list <- list() for (id in individual_ids) { cat("Processing individual:", id, "\n") individual_data <- spring_sf[spring_sf$id == id, ] cat("Number of unique dates:", length(unique(individual_data$DateTime)), "\n") # 处理重复时间戳（保留原逻辑） if (anyDuplicated(individual_data$DateTime) > 0) { individual_data$DateTime <- individual_data$DateTime + seq(0, length.out = nrow(individual_data), by = 1) } individual_coords <- st_coordinates(individual_data) tryCatch({ individual_ltraj <- as.ltraj(individual_coords, date = individual_data$DateTime, id = individual_data$id, typeII = TRUE) ltraj_list[[id]] <- individual_ltraj }, error = function(e) { cat("Error for individual:", id, "\n") print(e) }) } # 计算轨迹参数 ld_results <- lapply(ltraj_list, ld) ltraj_metrics <- do.call(rbind, ld_results) # 计算km/h单位的速度 ltraj_metrics$velocity_kmh <- (ltraj_metrics$dist / ltraj_metrics$dt) * 3.6 # 其他指标计算（保留原逻辑） total_distance <- aggregate(dist ~ id, data = ltraj_metrics, FUN = sum) sum_R2n_per_ID <- aggregate(R2n ~ id, data = ltraj_metrics, FUN = sum) R2n_data <- merge(sum_R2n_per_ID, total_distance, by = "id") R2n_data$straightness <- R2n_data$R2n / R2n_data$dist ltraj_metrics$turning_angle_x <- c(NA, diff(ltraj_metrics$abs.angle)) ltraj_metrics$turning_angle_y <- c(NA, diff(ltraj_metrics$rel.angle))
内容的提问来源于stack exchange，提问作者Devon
相关产品推荐
相关产品推荐
