R语言选取rastlist前31个MODIS遥感影像堆叠的代码验证及实现

原代码问题说明 #

1. 效率问题：你直接将全部90个栅格文件读入堆叠，是运行耗时的核心原因，完全没必要提前加载不需要的后续59个文件。
2. 逻辑错误：
   • stackApply的第二个参数需要与栅格堆叠的图层数长度一致，你给90层的堆叠传入只有31个元素的indices[1:31]，参数长度不匹配，无法得到正确结果。
   • 就算无报错运行，该代码也没有实现仅筛选前31天数据的需求，而是基于全部90层数据做计算，结果不符合预期。

正确实现方式 #

方案1：确认rastlist已按年积日升序排列的场景 #

你提到rastlist已按年积日001到090排序，直接截取前31个路径读取即可，代码更简洁：

library(raster)
# 仅读取前31个文件堆叠
jan_rasts <- stack(rastlist[1:31])
# 计算所有图层的均值，na.rm=TRUE避免空值导致结果异常
Jan_NDVI <- mean(jan_rasts, na.rm = TRUE)

方案2：不确定rastlist排序的保险方案 #

如果担心文件排序错乱，可先从文件名提取年积日筛选后再读取：

library(raster)
library(stringr)

# 从路径中提取年积日并转数值
doy <- as.numeric(str_extract(rastlist, "(?<=MOD09Q1_2013_)\\d{3}"))
# 筛选年积日≤31的文件路径
jan_rast_paths <- rastlist[doy <= 31]
# 堆叠后计算均值
jan_rasts <- stack(jan_rast_paths)
Jan_NDVI <- mean(jan_rasts, na.rm = TRUE)

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