things: Mostly spatial things

#' Cell-wise raster standard deviation
#'
#' Calculate a raster layer whose cells' values are the standard deviations
#' across the bands of a single raster, or across multiple rasters.
#'
#' @param infile A \code{RasterStack} or \code{RasterBrick}, a file path to a
#'   single, multiband raster file, or a vector of paths to multiple raster
#'   files. If a vector of paths is provided, then if any of the rasters have
#'   multiple bands, the first band will be used. A maximum of 26 bands or files
#'   can be provided to \code{infile}.
#' @param outfile A path to the desired output raster file, which will be
#'   created if it doesn't already exist (though the containing directory must
#'   exist). If not provided and \code{return_raster} is \code{TRUE}, a 
#'   temporary output file will be created.
#' @param return_raster Logical. Should the output raster be read back in to R
#'   as a \code{Raster} object?
#' @param quiet Logical. Should GDAL output messages be suppressed?
#' @return Unless \code{return_raster} is \code{TRUE}, \code{gdal_sd} creates
#'   the file specified at \code{outfile} but returns \code{NULL}, invisibly, to
#'   R. If \code{return_raster} is \code{TRUE}, the standard deviation
#'   \code{Raster} object is returned to R.
#' @details The quantity calculated is the corrected sample standard deviation,
#'   calculated as \eqn{sqrt(sum((x - mean(x))^2)/(length(x)-1))}.
#' @importFrom raster raster writeRaster
#' @importFrom rgdal GDALinfo
#' @export
#' @examples
#' library(raster)
#' s <- stack(replicate(5, raster(matrix(rnorm(100), ncol=10))))
#' 
#' # Input is a RasterStack
#' gdal_sd(s)
#' 
#' # Input is multiple raster files
#' writeRaster(s, ff <- paste0(tempfile(), 1:5, '.tif'), bylayer=TRUE)
#' gdal_sd(ff)
#' 
#' # Input is a raster file with multiple bands 
#' writeRaster(s, f <- tempfile(fileext='.tif'))
#' gdal_sd(f)
#' 
#' # Input is a RasterBrick with multiple bands
#' b <- brick(f)
#' gdal_sd(b)
gdal_sd <- function(infile, outfile, return_raster=TRUE, quiet=TRUE) {
  gdal_calc <- Sys.which('gdal_calc.py')
  if(gdal_calc=='') stop('gdal_calc.py not found on system.')
  if(missing(outfile) && return_raster) outfile <- tempfile(fileext='.tif')
  if(file.exists(outfile)) stop('outfile already exists.')
  if(is(infile, 'RasterStack') | is(infile, 'RasterBrick'))
    writeRaster(infile, 
                infile <- paste0(tempfile(), seq_len(nlayers(infile)), '.tif'),
                bylayer=TRUE)
  if(any(!file.exists(infile))) stop('Check that all files listed at infile exist.')
  nbands <- sapply(infile, function(x) nrow(attr(GDALinfo(x), 'df')))
  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 standard deviation and writing to ', basename(outfile))
  cmd <- 'python %s %s --outfile=%s --calc="std([%s], 0, ddof=1)"'
  out <- system(
    sprintf(cmd, gdal_calc, inputs, outfile, 
            paste0(LETTERS[seq_len(n)], collapse=',')),
    show.output.on.console=!quiet, intern=TRUE
  )
  if(any(grepl('Error', out))) stop(out, call.=FALSE)
  if(return_raster) {
    r <- raster::raster(outfile)
    names(r) <- 'sd'
    r
  } else {
    invisible(NULL)
  }
}