如何按小时计算风向矢量平均值？解决算术平均偏差问题

核心原理 #

风向属于环形数据，算术平均会出现类似350°和10°平均得到180°的错误。正确做法是将风向转换为正弦/余弦分量，计算分量的平均值后再转回角度；若要标准化为1m/s，直接对分量取算术平均即可（无需风速加权）。

方法1：基于Base R的自定义聚合函数 #

先定义矢量平均计算函数，再替换aggregate中的默认均值函数：

# 角度与弧度转换工具函数
deg2rad <- function(x) x * pi / 180
rad2deg <- function(x) x * 180 / pi

# 计算标准化为1m/s的矢量平均风向
vector_mean_wd <- function(wd) {
  wd_rad <- deg2rad(wd)
  # 计算正弦/余弦分量
  sin_comp <- sin(wd_rad)
  cos_comp <- cos(wd_rad)
  # 求分量均值（忽略NA）
  mean_sin <- mean(sin_comp, na.rm = TRUE)
  mean_cos <- mean(cos_comp, na.rm = TRUE)
  
  # 处理全NA的情况
  if (is.na(mean_sin) || is.na(mean_cos)) return(NA)
  
  # 转回角度并修正到0-360°范围
  avg_rad <- atan2(mean_sin, mean_cos)
  avg_deg <- rad2deg(avg_rad)
  ifelse(avg_deg < 0, avg_deg + 360, avg_deg)
}

# 自定义聚合函数：风向用矢量平均，其余变量用算术平均
aggregate_weather_stats <- function(group_data) {
  data.frame(
    WD = vector_mean_wd(group_data$WD),
    WS = mean(group_data$WS, na.rm = TRUE),
    temp = mean(group_data$temp, na.rm = TRUE),
    RH = mean(group_data$RH, na.rm = TRUE),
    press = mean(group_data$press, na.rm = TRUE),
    stringsAsFactors = FALSE
  )
}

# 执行聚合计算均值
<name>_mean <- aggregate(
  . ~ hour + day + month + year,
  data = <name>[c("hour", "day", "month", "year", "WD", "WS", "temp", "RH", "press")],
  FUN = aggregate_weather_stats
)

# 展开嵌套的结果数据框
<name>_mean <- do.call(cbind, c(<name>_mean[, 1:4], <name>_mean$x))

方法2：用dplyr实现（更简洁直观） #

dplyr的分组聚合语法更清晰，适合同时处理多种统计量：

# 安装并加载dplyr
install.packages("dplyr")
library(dplyr)

# 角度转换工具函数
deg2rad <- function(x) x * pi / 180
rad2deg <- function(x) x * 180 / pi

# 标准化为1m/s的矢量平均风向函数
vector_mean_wd <- function(wd) {
  wd_rad <- deg2rad(wd)
  mean_sin <- mean(sin(wd_rad), na.rm = TRUE)
  mean_cos <- mean(cos(wd_rad), na.rm = TRUE)
  
  if (is.na(mean_sin) || is.na(mean_cos)) return(NA)
  
  avg_rad <- atan2(mean_sin, mean_cos)
  avg_deg <- rad2deg(avg_rad)
  ifelse(avg_deg < 0, avg_deg + 360, avg_deg)
}

# 环形标准差函数（风向专用，替代普通sd）
circ_sd_wd <- function(wd) {
  wd_rad <- deg2rad(wd)
  mean_sin <- mean(sin(wd_rad), na.rm = TRUE)
  mean_cos <- mean(cos(wd_rad), na.rm = TRUE)
  rho <- sqrt(mean_sin^2 + mean_cos^2)
  # 环形标准差公式，转回度数
  rad2deg(sqrt(-2 * log(rho)))
}

# 计算均值
<name>_mean <- <name> %>%
  group_by(hour, day, month, year) %>%
  summarize(
    WD = vector_mean_wd(WD),
    WS = mean(WS, na.rm = TRUE),
    temp = mean(temp, na.rm = TRUE),
    RH = mean(RH, na.rm = TRUE),
    press = mean(press, na.rm = TRUE),
    .groups = "drop"
  )

# 计算标准差（风向用环形标准差）
<name>_sd <- <name> %>%
  group_by(hour, day, month, year) %>%
  summarize(
    WD_sd = circ_sd_wd(WD),
    WS_sd = sd(WS, na.rm = TRUE),
    temp_sd = sd(temp, na.rm = TRUE),
    RH_sd = sd(RH, na.rm = TRUE),
    press_sd = sd(press, na.rm = TRUE),
    .groups = "drop"
  )

关于circular包的补充说明 #

如果之前用circular包失败，可能是未正确转换为环形对象，正确用法如下：

install.packages("circular")
library(circular)

# 将风向转为环形对象（地理坐标系，0°为正北）
<name>$WD_circ <- circular(<name>$WD, units = "degrees", template = "geographics")

# 单独聚合风向的矢量平均
wd_mean <- aggregate(WD_circ ~ hour + day + month + year, data = <name>, FUN = mean.circular)
# 转回普通数值
wd_mean$WD <- as.numeric(wd_mean$WD_circ)

# 再与其他变量的聚合结果合并即可

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