直接处理HDF文件的循环中的堆栈光栅(不是tif文件)

Question

我正在处理1.460份HDF文件(4年的每日数据)。我有兴趣从所有的文件中获得平均的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)

Answer 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)

Answer 2

使用terra，您可以选择一些快捷方式，因为您可以直接读取HDF文件。你能做到的

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

也许你想要的东西可以这么简单：

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