gdal_project.R

## Function to use GDAL to project coordinate reference system

# See http://www.gdal.org/gdalwarp.html for additional details 
# `resampling` can be 'near' (nearest neighbour), 'bilinear', 'cubic', or
# 'lanczos' (Lanczos windowed sinc resampling). 
# `extent` should be a bbox object or a vector of c(xmin, ymin, xmax, ymax)
# `of` is the output format (use GDAL short name as given by the name field of
# gdalDrivers(), or at http://www.gdal.org/formats_list.html)
# `extension` is the output extension corresponding to the primary file
# `ot` is the output type (see http://www.gdal.org/gdal_translate.html)
# 'NB: This requires a GDAL installation, and requires that the path to the GDAL
# binaries is included on the system PATH.

gdal_project <- function(infile, outdir, xres, yres=xres, s_srs, t_srs, 
                         resampling, extent, of='GTiff', extension='tif', 
                         ot='Float32') {
  outfile <- infile
  extension(outfile) <- sub('^\\.+', '.', paste0('.', extension))
  outfile <- file.path(normalizePath(outdir), basename(outfile))
  system(sprintf('gdalwarp -ot %s -s_srs "%s" -t_srs "%s" -r %s -multi %s -tr %s -dstnodata -9999 -of %s "%s" "%s"',
                 ot, s_srs, t_srs, resampling, 
                 if(!missing(extent)) paste(c('-te', extent), collapse=' ') else '',
                 paste(xres, yres), of, infile, outfile))
  system(sprintf('gdalinfo -stats %s', outfile), show.output.on.console=FALSE)
}