在R语言中求解运动向量平均轨迹的技术问询

问题背景 #

现有一批带罗盘方位角（0-360°，0为正北、顺时针递增）和速度的运动向量，需要计算其平均向量（平均方位角+平均速度）。当前已有零散实现，但需要更简便的通用方案，且需支持非原点起始向量，最终适配带CRS的地理空间数据。

通用解决方案 #

1. 封装通用计算函数 #

将角度转换、加权平均、速度计算逻辑整合为一个函数，直接输入包含起点坐标、罗盘角、速度的数据框，输出平均向量的关键参数：

library(dplyr)
library(circular)
library(DescTools)

calculate_mean_vector <- function(df, x_col = "x", y_col = "y", angle_col = "compassdegree", vel_col = "velocity") {
  # 转换罗盘角为数学弧度（0为正东、逆时针递增）
  df_processed <- df %>%
    mutate(
      radian_degree = ((-!!sym(angle_col)) + 450) %% 360,
      radian = DescTools::DegToRad(radian_degree),
      # 计算每个向量的终点坐标
      x_end = !!sym(x_col) + !!sym(vel_col) * cos(radian),
      y_end = !!sym(y_col) + !!sym(vel_col) * sin(radian)
    )
  
  # 计算平均终点（处理非原点/不同起点的情况）
  mean_point <- df_processed %>%
    summarize(
      mean_x = mean(!!sym(x_col)),
      mean_y = mean(!!sym(y_col)),
      mean_x_end = mean(x_end),
      mean_y_end = mean(y_end)
    )
  
  # 计算平均位移的角度和速度
  delta_x <- mean_point$mean_x_end - mean_point$mean_x
  delta_y <- mean_point$mean_y_end - mean_point$mean_y
  
  mean_velocity <- sqrt(delta_x^2 + delta_y^2)
  mean_radian <- atan2(delta_y, delta_x)
  # 转换回罗盘角
  mean_compass_degree <- ((-circular::deg(mean_radian)) + 450) %% 360
  
  # 返回结果
  tibble(
    mean_compass_degree = round(mean_compass_degree, 2),
    mean_velocity = round(mean_velocity, 2),
    mean_start_x = mean_point$mean_x,
    mean_start_y = mean_point$mean_y,
    mean_end_x = mean_point$mean_x_end,
    mean_end_y = mean_point$mean_y_end
  )
}

2. 测试原点起始向量 #

用示例数据测试函数：

# 示例数据（原点起始）
df <- data.frame(
  x=0, 
  y=0,
  compassdegree = c(270,275,277,280,285,330, 40),
  velocity = c(2,2,2,2,1,1,1)
)

# 计算平均向量
mean_vec <- calculate_mean_vector(df)
print(mean_vec)
# 输出结果：mean_compass_degree = 286.85, mean_velocity = 1.86

3. 扩展到非原点/多起点向量 #

如果向量起点不同，函数会自动计算所有起点的均值，再结合终点均值得到平均位移向量：

# 非原点示例数据
df_non_origin <- data.frame(
  x = c(1, 2, 3),
  y = c(1, 2, 3),
  compassdegree = c(90, 90, 90),
  velocity = c(2, 2, 2)
)

mean_vec_non_origin <- calculate_mean_vector(df_non_origin)
print(mean_vec_non_origin)
# 平均起点为(2,2)，平均终点为(4,2)，罗盘角90°，速度2km/h

空间数据（带CRS）适配说明 #

核心要求：必须使用等距投影坐标系 #

平面向量的角度和距离计算依赖于等距属性，因此处理地理空间数据时：

1. 若数据是地理坐标系（如WGS84，EPSG:4326），需先转换为等距投影（如UTM投影，根据数据所在区域选择对应EPSG代码）
2. 在投影坐标系下完成向量平均计算
3. 如需返回地理坐标，再将结果转换回原地理坐标系

示例（使用sf包） #

library(sf)

# 创建带CRS的空间向量数据（WGS84）
sf_df <- st_as_sf(
  data.frame(
    lon = c(116.397, 116.400),
    lat = c(39.904, 39.905),
    compassdegree = c(90, 90),
    velocity = c(5, 5) # 单位：km
  ),
  coords = c("lon", "lat"),
  crs = 4326
)

# 转换为UTM投影（以北京附近为例，EPSG:32650）
sf_df_utm <- st_transform(sf_df, crs = 32650)

# 提取坐标到数据框
df_utm <- sf_df_utm %>%
  mutate(x = st_coordinates(.)[,1], y = st_coordinates(.)[,2]) %>%
  st_drop_geometry()

# 计算平均向量
mean_vec_utm <- calculate_mean_vector(df_utm)

# 将结果转换回空间对象并转WGS84
mean_sf <- st_as_sf(
  mean_vec_utm,
  coords = c("mean_start_x", "mean_start_y"),
  crs = 32650
) %>%
  st_transform(crs = 4326)

print(mean_sf)

关键注意事项 #

• 罗盘角与数学角的转换逻辑需严格对应：函数中((-compassdegree) + 450) %% 360确保0°正北转换为数学上的90°（正东），符合R绘图函数的角度定义
• 加权平均逻辑：函数通过终点均值实现了速度加权（速度大的向量对终点的影响更大），与原方法的加权角度计算结果一致
• 空间数据中，避免直接在地理坐标系（经纬度）下计算向量，因为经纬度的距离和角度会随纬度变化而变形

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