将小时级栅格堆叠数据聚合为6小时均值层级

不用纠结索引整除问题，直接通过时间标签分组来处理是更可靠的方案，不管最后一天数据是否完整，都能精准匹配时段并计算均值。以下是基于R语言raster包的完整实现步骤：

步骤1：加载依赖包并确认栅格时间信息 #

首先确保你有raster和lubridate包，前者处理栅格数据，后者处理时间格式：

library(raster)
library(lubridate)

假设你的小时级SST栅格堆叠名为sst_stack，先提取或生成时间戳：

# 优先从栅格的Z维度读取时间（如果数据自带时间属性）
time_stamps <- getZ(sst_stack)

# 如果没有Z维度，从图层名称提取时间（假设图层名是类似"2021013100"的格式）
if (is.null(time_stamps)) {
  time_stamps <- ymd_h(names(sst_stack))
}

步骤2：定义时段分组与命名规则 #

根据你的需求，把每个小时数据映射到对应的6小时时段，并生成符合要求的图层名称：

# 提取每个时间戳的小时数
hour_of_day <- hour(time_stamps)

# 匹配时段对应的结束小时
period_end_hour <- case_when(
  hour_of_day >= 0 & hour_of_day < 6 ~ 6,
  hour_of_day >= 6 & hour_of_day < 12 ~ 12,
  hour_of_day >= 12 & hour_of_day < 18 ~ 18,
  hour_of_day >= 18 & hour_of_day <= 23 ~ 0
)

# 生成时段对应的日期（18:00-24:00时段仍用当天日期）
period_date <- date(time_stamps)

# 生成符合要求的图层名称，比如"X2021.01.31.06"
layer_names <- paste0(
  "X", 
  format(period_date, "%Y.%m.%d"), 
  ".", 
  sprintf("%02d", period_end_hour)  # 确保小时数是两位格式
)

# 生成分组索引：相同名称的图层属于同一时段
group_indices <- match(layer_names, unique(layer_names))

步骤3：计算时段均值并生成新栅格堆叠 #

用stackApply函数批量计算每个时段的均值，这是raster包专门优化的分组计算工具，效率比循环高：

# 计算各时段均值，na.rm=TRUE处理缺失值（根据你的数据情况调整）
sst_6hr_stack <- stackApply(
  x = sst_stack, 
  indices = group_indices, 
  fun = mean, 
  na.rm = TRUE
)

# 给新栅格堆叠设置正确的图层名称
names(sst_6hr_stack) <- unique(layer_names)

步骤4：处理不完整时段（可选） #

如果你的数据最后一天时段不足6小时，可选择过滤掉这些不完整的时段：

# 统计每个时段包含的小时数
group_counts <- table(group_indices)

# 筛选出包含至少6个小时的完整时段
valid_groups <- as.integer(names(group_counts)[group_counts >= 6])

# 生成仅包含完整时段的栅格堆叠
sst_6hr_stack_valid <- subset(sst_6hr_stack, valid_groups)

关键说明 #

• 用时间标签分组替代索引计算，彻底避免了"无法被6整除"的问题，自动适配所有时段
• case_when的规则完全匹配你的时段划分：00-06结束于06，18-24结束于00（当天）
• 如果不需要过滤不完整时段，直接保留sst_6hr_stack即可，它会用所有可用的小时数据计算均值

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