在R中用MOD13Q1 EVI像素填补MYD13Q1 EVI的像素空缺

针对你需要用Terra（MOD13Q1）有效像素填补Aqua（MYD13Q1）NA像素的需求，推荐使用raster包原生的cover()函数——它专门用于栅格对象的NA值填补，比手动索引更高效，且适配批量处理场景。

核心步骤与代码 #

1. 加载依赖包

library(raster)
library(purrr) # 可选，用于简化批量处理逻辑

2. 规范加载数据（确保时间顺序匹配）

先修正数据加载流程，保证两个影像列表的文件严格按时间排序，避免跨月份匹配错误：

# 加载Terra MOD13Q1数据
setwd("C:/...sample/mod")
pathtif_evi <- "C:/...sample/mod/"
tiflist_evi <- list.files(path = pathtif_evi, pattern = "*.tif", full.names = TRUE)
tiflist_evi <- sort(tiflist_evi) # 按文件名排序，保证时间顺序一致
all_tiffs_evi <- lapply(tiflist_evi, raster)

# 加载Aqua MYD13Q1数据
setwd("C:/...sample/myd")
pathtif_myd <- "C:/...sample/myd/"
tiflist_myd <- list.files(path = pathtif_myd, pattern = "*.tif", full.names = TRUE)
tiflist_myd <- sort(tiflist_myd) # 与MOD列表同步排序
all_tiffs_myd <- lapply(tiflist_myd, raster)

3. 批量完成像素填补与保存

方法1：基础循环（适合新手理解）

# 创建结果存储目录（提前手动创建或用dir.create）
dir.create("C:/...sample/merged", showWarnings = FALSE)

# 循环处理每一对影像
for(i in 1:length(all_tiffs_evi)){
  # 用Terra的非NA值替换Aqua的NA像素
  merged_raster <- cover(all_tiffs_myd[[i]], all_tiffs_evi[[i]])
  # 保存融合结果
  out_path <- paste0("C:/...sample/merged/merged_", basename(tiflist_myd[i]))
  writeRaster(merged_raster, filename = out_path, format = "GTiff", overwrite = TRUE)
}

方法2：purrr批量处理（代码更简洁）

dir.create("C:/...sample/merged", showWarnings = FALSE)

# 遍历影像对，完成填补并保存
walk2(all_tiffs_myd, all_tiffs_evi, function(myd_raster, mod_raster){
  merged <- cover(myd_raster, mod_raster)
  out_name <- paste0("C:/...sample/merged/merged_", basename(filename(myd_raster)))
  writeRaster(merged, filename = out_name, format = "GTiff", overwrite = TRUE)
})

为什么你之前的代码无效？ #

你之前的代码将all_tiffs_evi和all_tiffs_myd当作向量处理，但它们是栅格对象的列表，无法直接用[]索引整个列表的像素。必须逐个提取列表中的栅格对象，再进行像素级操作。

注意事项 #

• 务必保证两个影像列表的文件顺序严格对应（按时间排序），否则会出现跨月份错误填补
• 处理22年月度数据（共264个影像）时，建议边处理边保存，避免内存溢出
• 确认所有影像的投影、分辨率、范围完全一致（你提到已在GEE完成统一处理，此步骤可跳过）

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