创建一个循环以从栅格堆栈创建 NDVI 图像并根据文件名命名它们

Question

我一直在尝试编写一个循环来遍历 Sentinel 2 卫星图像（波段 4 和 5）的两个文件夹并获取每个日期的 NDVI。

为每个波段创建一个堆栈，进行一些裁剪和重新采样以最终进行 NDVI 计算。我在循环中集成 NDVI 计算和创建文件名时遇到了困难。

我只想让我的循环为 x 个日期生成 x 个文件，然后为每个 NDVI 图像提供日期作为从文件名中提取的名称“YYYY/MM/DD.tif”。但我想不出办法，经过多次尝试和错误失败。

#list files
files4 <- list.files(path4, pattern = "jp2$", full.names = TRUE)
files5 <- list.files(path5, pattern = "jp2$", full.names = TRUE)

ms5 <- stack()
ms4 <- stack()

for (f in files4){
  # loading a raster
  r4 <- raster(f)
  proj4string(r4)
  proj4string(emprise)
  emprise <- spTransform(emprise, proj4string(r4)) 
  r4b <- crop(r4, emprise)
  ms4<- stack(ms4,r4b)
#copy the date from the file to give a name to the final NDVI image (I have to get ride of everything but the date
  x <- gsub("[A-z //.//(//)]", "", r4)
  y <- substr(x, 4, 11)
}


for (f in files5){
  # load the raster
  r5 <- raster(f)
  proj4string(r5)
  proj4string(emprise)
  emprise <- spTransform(emprise, proj4string(r5)) 
  r5b <- crop(r5, emprise)
  ms5<- stack(ms5,r5b)
}

#Resampling : setting the Band 5 to the same resolution as Band 4 
b5_resamp <- resample(ms5, ms4)

Answer 1

您是否考虑过循环日期而不是文件？没有示例数据，我无法给出更具体的建议，但这里是总体思路：

# List files
files4 <- list.files("./band4", pattern = ".tif", full.names = TRUE)
#> "band4/T31UDR_20170126T105321_B04.tif" "band4/T31UDR_20180126T105321_B04.tif"

files5 <- list.files("./band5", pattern = ".tif", full.names = TRUE)
#> "./band5/T31UDR_20170126T105321_B05.tif" "./band5/T31UDR_20180126T105321_B05.tif"

# Get dates
dates <- unique(gsub(pattern = ".*_(\d{8}).*", replacement = "\1", x = c(files4, files5)))
#> "20170126" "20180126"

# Define empty stacks
ms5 <- stack()
ms4 <- stack()

for(date in dates){
  
  ## Band 4
  
  f4 <- list.files("./band4", pattern = date, full.names = TRUE)
  
  # loading a raster
  r4 <- raster(f4)
  proj4string(r4)
  proj4string(emprise)
  emprise <- spTransform(emprise, proj4string(r4)) 
  r4b <- crop(r4, emprise)
  ms4 <- stack(ms4,r4b)
  
  ## Band 5
  
  f5 <- list.files("./band5", pattern = date, full.names = TRUE)
  
  # load the raster
  r5 <- raster(f5)
  proj4string(r5)
  proj4string(emprise)
  emprise <- spTransform(emprise, proj4string(r5)) 
  r5b <- crop(r5, emprise)
  ms5<- stack(ms5,r5b)
  
  ## Resampling : setting the Band 5 to the same resolution as Band 4 
  
  b5_resamp <- resample(ms5, ms4)

  ## Write to file
  
  writeRaster(b5_resamp, filename = paste0(date, ".tif"))
}

创建一个循环以从栅格堆栈创建 NDVI 图像并根据文件名命名它们

Creating a loop to create NDVI images from raster stacks and naming them from the file name

gis

r

raster