R语言rasterbrick按年提取均值及年最小3个月均值方法咨询

基于现有raster工具链即可实现全流程栅格层面计算，不需要转换为data.frame，运算效率远高于nc读取+提取值+数据框分组的工作流。

前置准备：数据读取与目标时段子集 #

library(R.utils)
library(raster)
library(tidyverse)

# 下载并解压示例海冰数据
download.file("https://www.metoffice.gov.uk/hadobs/hadisst/data/HadISST_ice.nc.gz","HadISST_ice.nc.gz")
gunzip("HadISST_ice.nc.gz", ext="gz", FUN=gzfile)

# 读取NetCDF为栅格时间序列对象，注意文件路径匹配本地存储位置
hadISST <- brick('HadISST_ice.nc') 

# 子集目标研究时段（示例取1900-1909共10年月度数据）
hadISST_sub <- hadISST %>% 
  subset(., which(getZ(.) >= as.Date("1900-01-01") & getZ(.) <= as.Date("1909-12-31")))

问题1：逐像元计算逐年均值 #

直接使用stackApply()按年份分组做逐像元计算即可，不需要展开数据：

# 从Z维度时间戳提取每个栅格层对应的年份，作为分组索引
time_vec <- getZ(hadISST_sub)
year_index <- format(time_vec, "%Y") %>% as.numeric()

# 逐像元计算逐年均值，输出结果为多波段栅格，每个波段对应1年的均值结果
hadISST_yearly_mean <- stackApply(
  x = hadISST_sub,
  indices = year_index,
  fun = mean,
  na.rm = TRUE
)
# 给结果波段命名为对应年份，方便后续识别
names(hadISST_yearly_mean) <- paste0("year_", unique(year_index))

问题2：每年取数值最小的3个月份计算年度均值 #

自定义栅格计算函数，同样传入stackApply()按年分组计算，全程保持栅格对象结构：

# 自定义计算函数：对单个像元的年度月度序列，取最小3个有效值求均值
calc_min3_mean <- function(x, na.rm = TRUE){
  if(na.rm) x <- x[!is.na(x)]
  # 当年有效观测不足3个时返回NA，避免无效计算
  if(length(x) < 3) return(NA_real_)
  # 排序取最小3个值计算均值
  sort(x, decreasing = FALSE)[1:3] %>% mean()
}

# 按年分组调用自定义函数计算
hadISST_yearly_min3mean <- stackApply(
  x = hadISST_sub,
  indices = year_index,
  fun = calc_min3_mean,
  na.rm = TRUE
)
names(hadISST_yearly_min3mean) <- paste0("year_min3_", unique(year_index))

方案优势 #

• 无冗余转换开销：全程操作栅格原生对象，不需要调用ncvar_get()读取全量数组，也不需要extract()把栅格转成数据框
• 内存效率高：raster默认支持分块处理，GB级大文件也不会直接占满运行内存
• 扩展性强：如果需要换其他统计逻辑（比如取最大3个月、计算分位数等），只需要修改传入stackApply()的自定义函数即可

提示：如果处理的是高分辨率、长时间序列的大体积NetCDF文件，可以将raster替换为terra包，对应使用rast()读取数据、tapp()做分组计算，接口逻辑和raster基本一致，运算速度可提升3-10倍，内存占用更低。

内容的提问来源于stack exchange，提问作者marine-ecologist