编写循环处理Sentinel 2影像生成NDVI并按文件名提取日期命名求助

问题解决思路 #

你的代码存在三个核心问题：

1. 日期提取逻辑错误：你对raster对象执行gsub替换，而非对原始文件名处理，导致无法正确提取日期
2. 文件配对逻辑缺失：未按日期匹配对应的4波段、5波段影像，直接堆叠全量影像容易出现日期错位
3. 缺少逐日期的NDVI计算与导出逻辑

提示：标准NDVI计算使用红波段（B4）和近红外波段（B8），如果你明确需要使用红边波段（B5）计算红边NDVI可直接使用下述代码，否则建议将5波段替换为8波段影像路径。
注意：操作系统不允许文件名中出现/，如果你要求的YYYY/MM/DD.tif是指分年月日三级目录存储最终影像，代码会自动创建对应目录；如果是文件名格式，请自行将分隔符替换为-或其他合法字符。

完整可运行代码 #

# 加载依赖包
library(raster)
library(rgdal)

# 配置路径（请替换为你的实际路径）
path4 <- "你的4波段文件夹路径"
path5 <- "你的5波段文件夹路径"
output_root <- "NDVI输出根目录路径"

# 读取文件列表
files4 <- list.files(path4, pattern = "jp2$", full.names = TRUE)
files5 <- list.files(path5, pattern = "jp2$", full.names = TRUE)

# 定义日期提取函数（适配Sentinel 2标准文件名格式）
get_s2_date <- function(filename){
  # 从文件名中提取8位日期字符串YYYYMMDD
  file_basename <- gsub(".*/", "", filename)
  date_str <- substr(file_basename, 12, 19)
  return(date_str)
}

# 提取所有文件的日期，按日期排序文件确保配对正确
dates4 <- sapply(files4, get_s2_date)
dates5 <- sapply(files5, get_s2_date)
files4 <- files4[order(dates4)]
files5 <- files5[order(dates5)]

# 校验两个波段的文件数量和日期是否匹配
if(length(files4) != length(files5) || !all(sort(dates4) == sort(dates5))){
  stop("4波段和5波段文件数量/日期不匹配，请检查文件夹内容")
}

# 逐日期处理
for(i in seq_along(files4)){
  # 读取当前日期的两个波段
  r4 <- raster(files4[i])
  r5 <- raster(files5[i])
  
  # 裁剪处理
  emprise <- spTransform(emprise, proj4string(r4))
  r4_crop <- crop(r4, emprise)
  r5_crop <- crop(r5, emprise)
  
  # 重采样5波段到4波段分辨率
  r5_resamp <- resample(r5_crop, r4_crop)
  
  # 计算NDVI
  ndvi <- (r5_resamp - r4_crop) / (r5_resamp + r4_crop)
  
  # 生成输出路径和文件名
  date_str <- get_s2_date(files4[i])
  y <- substr(date_str, 1,4)
  m <- substr(date_str, 5,6)
  d <- substr(date_str, 7,8)
  # 创建三级目录 YYYY/MM
  dir.create(file.path(output_root, y, m), recursive = TRUE, showWarnings = FALSE)
  output_path <- file.path(output_root, y, m, paste0(d, ".tif"))
  
  # 若不需要分目录存储，仅需要文件名格式为YYYY-MM-DD.tif，可替换为以下代码：
  # output_path <- file.path(output_root, paste0(y, "-", m, "-", d, ".tif"))
  
  # 导出NDVI影像
  writeRaster(ndvi, filename = output_path, format = "GTiff", overwrite = TRUE)
}

内容的提问来源于stack exchange，提问作者Perrin Remonté