从栅格列表中计算平均值并将其保存为不同的名称答案

【问题标题】：calculate mean from list of rasters and and save it in a different name从栅格列表中计算平均值并将其保存为不同的名称
【发布时间】：2019-11-15 15:38:04
【问题描述】：

我有多年的栅格列表（.tif 格式）。这是来自 landsat 的 16 天 NDVI，我想制作每月 NDVI（两个连续栅格的平均值）并将其作为每月平均值保存在相同或不同的目录中

我已经列出了栅格栅格并将其堆叠起来，后来我使用 stackApply 来计算平均值，但它会产生空栅格。我一年有 23 张图像，我想平均它并制作 12 个月。这就是我的光栅文件的样子

 "landsatNDVISC05SLC2000001.tif" "landsatNDVISC05SLC2000017.tif"
 "landsatNDVISC05SLC2000033.tif" "landsatNDVISC05SLC2000049.tif"
 "landsatNDVISC05SLC2000065.tif" "landsatNDVISC05SLC2000081.tif"
 "landsatNDVISC05SLC2000097.tif" "landsatNDVISC05SLC2000113.tif"
 "landsatNDVISC05SLC2000129.tif" "landsatNDVISC05SLC2000145.tif"
 "landsatNDVISC05SLC2000161.tif" "landsatNDVISC05SLC2000177.tif"
 "landsatNDVISC05SLC2000193.tif" "landsatNDVISC05SLC2000209.tif"
 "landsatNDVISC05SLC2000225.tif" "landsatNDVISC05SLC2000241.tif"
 "landsatNDVISC05SLC2000257.tif" "landsatNDVISC05SLC2000273.tif"
 "landsatNDVISC05SLC2000289.tif" "landsatNDVISC05SLC2000305.tif"
 "landsatNDVISC05SLC2000321.tif" "landsatNDVISC05SLC2000337.tif"
 "landsatNDVISC05SLC2000353.tif

此代码有效，但会产生超过 12 个空栅格，我还想将栅格砖保存为单个子集每月栅格

library(raster)
lrast<-list.files("G:/LANDSAT-NDVI/testAverage")
layers<-paste("landsatNDVISC05SLC2000", seq(from=001, to=353,by=16))
stak<-stack(lrast)
raster<-stackApply(stak, layers, fun = mean)

我想将 landsatNDVISC05SLC2000001.tif 和 landsatNDVISC05SLC2000017.tif 的月平均值设为 landsatNDVISC05SLC2000M1.tif。同样，33,49，由于我只有 23 个栅格，我想将 landsatNDVISC05SLC2000353.tif 保留为 landsatNDVISC05SLC2000M12.tif

块引用

【问题讨论】：

标签： r raster

【解决方案1】：

不确定 stackapply 是如何工作的，但这样的东西应该可以满足需要。

library(raster)
files <- list.files(path = "...", full.names = T, pattern = ".tif")

stk <- stack()

for (i in files){
  print(i)
  as <- raster(files[i])
  stk <- addLayer(stk, as)
}

jday <-c("landsatNDVISC05SLC2000017.tif","landsatNDVISC05SLC2000033.tif",
"landsatNDVISC05SLC2000049.tif","landsatNDVISC05SLC2000065.tif","landsatNDVISC05SLC2000081.tif",
"landsatNDVISC05SLC2000097.tif","landsatNDVISC05SLC2000113.tif","landsatNDVISC05SLC2000129.tif",
"landsatNDVISC05SLC2000145.tif","landsatNDVISC05SLC2000161.tif","landsatNDVISC05SLC2000177.tif",
"landsatNDVISC05SLC2000193.tif","landsatNDVISC05SLC2000209.tif","landsatNDVISC05SLC2000225.tif",
"landsatNDVISC05SLC2000241.tif","landsatNDVISC05SLC2000257.tif","landsatNDVISC05SLC2000273.tif",
"landsatNDVISC05SLC2000289.tif","landsatNDVISC05SLC2000305.tif","landsatNDVISC05SLC2000321.tif",
"landsatNDVISC05SLC2000337.tif","landsatNDVISC05SLC2000353.tif")

jday <- as.numeric(substr(jday, 24, 25)) #substract the julien days (which I think these number represent before .tif; or you can substract the names from the 'files' vector)

dates <- as.Date(jday, origin=as.Date("2000-01-01")) # create a Date vector

stk <- setZ(stk, dates) # assign the date vector to the raster stack

raster <- zApply(stk, by = format(dates,"%Y-%m"), fun = mean, na.rm = T) # create the monthly stack

【讨论】：

感谢 Andrei Nita，当我运行 setZ 函数时，它给出“setZ(stak, dates) 中的错误：length(z) == nlayers(x) is not TRUE”，是因为当我打印日期，它从 2000-01-02 开始，即使原点是 2000-01-01？堆叠层和日期的长度相似（在我的例子中是 23 层）。
length(z) == nlayers(x) 是真的吗？日期以 2000-01-02 开头的事实应该不是问题，只要它的长度与堆栈中的层数相同。
由于某种原因，这些单波段 tiff 图像被识别为两个波段 tiff，所以当我查询 nlayers 时，它显示了 46 层，但实际上它只有 23 层。仍然无法弄清楚问题所在。感谢安德烈的努力
我更新了答案。看起来您阅读了两次相同的光栅。但是，不看数据很难给出答案:)