直接使用 HDF 文件（不是 tif 文件）循环堆叠光栅答案

【问题标题】：Stack raster in loop working with HDF files directly (not tif files)直接使用 HDF 文件（不是 tif 文件）循环堆叠光栅
【发布时间】：2022-01-26 12:05:56
【问题描述】：

我正在处理 1.460 个 HDF 文件（4 年每日数据）。我有兴趣从所有文件中获取 MEAN AOD。使用以下代码，我只从最后一个文件中获取信息，而不是我正在使用的所有文件的组合，我不确定为什么这不起作用。（我对在此过程中获取 TIF 文件不感兴趣，但我不确定这是否有必要获得我想要的）

read_hdf = function(file, n) {
  sds = get_subdatasets(file)
  f = tempfile(fileext = ".tif")
  gdal_translate(sds[1], f)
  brick(f)
}

files <- dir(pattern = ".hdf")

for(f in files) {
  sds = get_subdatasets(f)
  Optical_Depth_047 = read_hdf(f, grep("grid1km:Optical_Depth_047", sds))
  names(Optical_Depth_047) = paste0("Optical_Depth_047.", letters[1:nlayers(Optical_Depth_047)])
  r = stack(Optical_Depth_047)
} 

meanAOD <- stackApply(r, indices =  rep(1,nlayers(r)), fun = "mean", na.rm = T)

【问题讨论】：

标签： r gis raster gdal rgdal

【解决方案1】：

欢迎来到 Stack Overflow M. Fernanda Valle Seijo！为了将来参考，请尝试发布一个可重现的示例。您可以查看他们的指南here。这些有助于回答您问题的人更好地诊断问题。

我认为您的代码在正确的轨道上，但看起来您只获得最后一项的原因是您在循环中包含了 stack() 函数。这将在循环的每次迭代中覆盖r。此外，您可能需要先将Optical_Depth_047 设置为列表，并将循环的每次迭代保存为该列表的元素。尝试像这样运行您的代码：

read_hdf = function(file, n) {
  sds = get_subdatasets(file)
  f = tempfile(fileext = ".tif")
  gdal_translate(sds[1], f)
  brick(f)
}

files <- dir(pattern = ".hdf")

Optical_Depth_047<-list()
for(f in files) {
  sds = get_subdatasets(f)
  Optical_Depth_047[[f]] = read_hdf(f, grep("grid1km:Optical_Depth_047", sds))
  names(Optical_Depth_047)[f] = paste0("Optical_Depth_047.", letters[f])
} 
r = stack(Optical_Depth_047)
meanAOD <- stackApply(r, indices =  rep(1,nlayers(r)), fun = "mean", na.rm = T)

【讨论】：

【解决方案2】：

使用terra，您可以采取一些捷径，因为您可以直接读取 HDF 文件。你可以这样做

library(terra)
r <- rast(f, "grid1km:Optical_Depth_047")

也许你所追求的可以这么简单：

x <- rast(files, "grid1km:Optical_Depth_047")
m <- mean(x)

【讨论】：