Setting raster projections, extent, and computing annual summaries

HDF_to_TIF.R

library(rgdal)
library(rgeos)
library(raster)
library(gdalUtils)


# gdalinfo("GFED40_MQ_199506_BA.hdf")
# sds <- get_subdatasets("GFED40_MQ_199506_BA.hdf")
# sds

files <- list.files('HDF', pattern = ".hdf", full.names = TRUE)

filename <- substr(files, 15, 20)
filename <- paste0("BA", filename, ".tif")
n_files <- length(files)

for (i in seq_along(files)){
  sds <- get_subdatasets(files[i])
  gdal_translate(sds[1], dst_dataset = filename[i])
}

rstack <- stack(filename)
projection(rstack) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") 
bb <- c(-180, 180, -90, 90)
rstack <- setExtent(rstack, bb)

# compute a summary function for each year
years <- substring(names(rstack), 3, 6)
year_index <- as.numeric(as.factor(years))

annual_summaries <- stackApply(rstack, indices = year_index, fun = mean)