gdal_system2.R

# 2018-07-03
# An example of doing some slightly fancy stuff with GDAL on the command line 
# from R using OSGeo4W GDAL, which is on system PATH already

# Task: rescale all input covariate rasters to the last BM run to [0, 1], see 
# if it makes any difference to model outputs (spoiler: no)

library(raster)
library(tidyverse)
options(stringsAsFactors = FALSE)

# get rasters from in here
last_run <- file.path(getwd(), 'dsmart_results', 'BM_20171007')

# list rasters to rescale
covs <- list.files(file.path(last_run, 'inputs'), 
                   pattern = '^rs_.*\\.tif',
                   full.names = TRUE)

# output rasters will go into
if(!dir.exists(file.path(getwd(), 'covariates_out', 'rescaled'))) {
  dir.create(file.path(getwd(), 'covariates_out', 'rescaled'))
}
dest <- file.path(getwd(), 'covariates_out', 'rescaled')

# using GDAL
# first get a table of raster min/maxes and check they're ok
raster_stats <- map_df(covs, function(cov) {
  info <- system2('gdalinfo', 
                  args = c('-mm', covs), stdout = TRUE)
  mmstring <- substr(info[37], 22, nchar(info[37]))
  data.frame('name' = tools::file_path_sans_ext(basename(cov)),
             'min'  = as.numeric(stringr::str_split(mmstring, ',')[[1]][1]),
             'max'  = as.numeric(stringr::str_split(mmstring, ',')[[1]][2]),
             stringsAsFactors = FALSE)
})

# then rescale
pmap(list(covs, raster_stats$min, raster_stats$max),
     function(covs, min, max) {
  outname <- gsub('^rs_', 'SC01_', basename(covs))
  system2('gdal_translate', args = c(
          '-scale', min, max, 0, 1,
          '-ot Float32', # important if input is integer - categ data layers
          covs,
          file.path(dest, outname)
          ))
}
)

# job done \o/