Fast calculation of cellwise mean across a raster stack, using gdal_calc.py

gdal_weightedmean.R

gdal_weightedmean <- function(infile, outfile, weights, return_raster=FALSE, overwrite=FALSE) {
  # Be aware that the outfile type will be the same as the infile type
  require(rgdal)
  if(return_raster) require(raster)
  # infile:        The multiband raster file (or a vector of paths to multiple 
  #                 raster files) for which to calculate cell mean.
  # weights:       The weights to apply to each layer. If missing, equal weights
  #                 assumed.
  # outfile:       Path to raster output file.
  # return_raster: (logical) Should the output raster be read back into R?
  # overwrite:     (logical) Should outfile be overwritten if it exists?
  gdal_calc <- Sys.which('gdal_calc.py')
  if(gdal_calc=='') stop('gdal_calc.py not found on system.')
  if(file.exists(outfile) & !overwrite) 
    stop("'outfile' already exists. Use overwrite=TRUE to overwrite.", 
         call.=FALSE)
  nbands <- sapply(infile, function(x) nrow(attr(GDALinfo(x), 'df')))
  if(missing(weights)) weights <- rep(1, nbands)
  if(length(infile) > 26 || nbands > 26) stop('Maximum number of inputs is 26.')
  if(length(nbands) > 1 & any(nbands > 1)) 
    warning('One or more rasters have multiple bands. First band used.')
  
  if(length(infile)==1) {
    inputs <- paste0('-', LETTERS[seq_len(nbands)], ' ', infile, ' --', 
                     LETTERS[seq_len(nbands)], '_band ', seq_len(nbands), collapse=' ')
    n <- nbands
  } else {
    inputs <- paste0('-', LETTERS[seq_along(nbands)], ' ', infile, ' --', 
                     LETTERS[seq_along(nbands)], '_band 1', collapse=' ')
    n <- length(infile)
  }
  
  message('Calculating mean and writing to ', basename(outfile))
  system2('python',
          args=c(gdal_calc, inputs,
                 sprintf("--outfile=%s", outfile),
                 sprintf('--calc="average([%s], axis=0, weights=[%s])"', 
                         paste0(LETTERS[seq_len(n)], collapse=','),
                         paste(weights, collapse=', ')),
                 '--co="COMPRESS=LZW"',
                 #'--type="Float32"', 
                 if(overwrite) '--overwrite'),
          stdout=FALSE
  )
  if(return_raster) raster(outfile) else invisible(NULL)
}