R语言遍历列表调用geosphere::bearing计算方位角报错排查

geosphere::bearing()的入参要求为数值型结构：支持两列的经纬度矩阵（第一列为经度、第二列为纬度）或长度为2的经纬度数值向量，不支持直接传入sf类空间点对象。
代码中通过st_as_sf转换得到的空间点是列表结构的S3类对象，传入bearing()后，函数内部的.pointsToMatrix()逻辑无法将列表强制转换为double类型的数值矩阵，因此触发该报错。
方位角计算的核心逻辑是计算当前观测点指向下一个相邻观测点的方向夹角，因此每个ID分组内需要先按观测时间/采样顺序排序，再配对前后点坐标计算，每个分组的第一个点无前置参照点，计算结果会返回NA。

方案1：不拆分列表，直接分组计算（推荐，适配大数据量场景） #

无需提前用split拆分为列表，直接按ID分组批量计算即可，代码最简洁、运行效率最高：

library(geosphere)
library(dplyr)

# 替换为你自己的原始数据框名、排序字段名：如果没有观测时间字段，可删除arrange行
result <- df %>%
  group_by(ID) %>%
  arrange(obs_time, .by_group = TRUE) %>%
  mutate(
    # 提取下一个点的经纬度作为方位角终点
    next_x = lead(x),
    next_y = lead(y),
    # 传入经纬度矩阵计算，注意顺序：经度在前，纬度在后
    bearing_angle = bearing(cbind(x, y), cbind(next_x, next_y))
  ) %>%
  select(-next_x, -next_y) %>%
  ungroup()

方案2：沿用split拆分逻辑的lapply写法 #

如果你已经熟悉列表拆分的逻辑，可以直接在每个子数据块中提取数值型经纬度计算，无需转换sf对象：

library(geosphere)

# 按ID拆分为列表
df_list <- split(df, df$ID)

# 遍历每个分组计算
cal_list <- lapply(df_list, function(sub_df){
  # 按观测顺序排序，无时间字段可删除该行
  sub_df <- sub_df[order(sub_df$obs_time), ]
  # 构造前后点的经纬度矩阵
  p1 <- cbind(sub_df$x, sub_df$y)
  p2 <- rbind(p1[-1, ], c(NA, NA)) # 最后一行补NA，和lead逻辑对齐
  sub_df$bearing_angle <- bearing(p1, p2)
  return(sub_df)
})

# 合并所有分组结果为完整数据框
result <- do.call(rbind, cal_list)

方案3：purrr包遍历写法 #

如果习惯用purrr族函数，可以用map_dfr直接返回合并后的数据框，省去手动合并步骤：

library(geosphere)
library(purrr)
library(dplyr)

result <- df %>%
  split(.$ID) %>%
  map_dfr(~{
    # 按观测顺序排序
    .x <- .x %>% arrange(obs_time)
    p1 <- cbind(.x$x, .x$y)
    p2 <- lead(p1)
    .x$bearing_angle <- bearing(p1, p2)
    return(.x)
  })

注意事项：geosphere::bearing()返回的角度以正北为0度，顺时针旋转取值范围为0-360度；传入坐标时必须保证经度在前、纬度在后，如果你的x/y字段顺序为纬度、经度，需要调换cbind内的字段顺序，否则计算结果会完全错误。

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