Stack raster in loop working with HDF files directly (not tif files)-CodePudding

I'm working with 1.460 HDF files (4 years daily data). I'm interested in obtaining MEAN AOD from all the files. With the following code I only obtain information from the last file, not the combination of all the files I'm working with and I'm not sure why this is not working. (I'm not interested in getting TIF files in the process, but I'm not sure if this is necessary to get what I want)

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)

CodePudding user response：

Welcome to Stack Overflow M. Fernanda Valle Seijo! For future reference, please try to post a reproducible example. You can see guidelines for them here. These help people answering your questions better diagnose the issue.

I think your code is on the right track, but it looks like the reason you are getting just the last item is that you have included the stack() function in you loop. This will overwrite r in each iteration of the loop. Also, you may need to setup Optical_Depth_047 as a list first and save each iteration of the loop as an element of that list. Try running your code like this:

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)

CodePudding user response：

with terra you can take some shortcuts, as you can read the HDF files directly. You can do

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

Perhaps what you are after can be as simple as this:

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