#' Linear interpolation between pairs of rasters
#'
#' Linearly interpolate cell values between pairs of raster grids, at specified
#' locations (e.g. "times"). Extrapolation to times beyond the last raster in
#' the series is also possible.
#'
#' @param s A \code{RasterStack} or \code{RasterBrick} containing the rasters
#' between which interpolation will be performed.
#' @param z_in A numeric vector indicating the positions (e.g., times)
#' associated with layers in \code{s} (in the same order as the layers of
#' \code{s} - see \code{names(s)}).
#' @param z_out a numeric vector that indicates the times to be interpolated to
#' (if \code{z_out} extends beyond the latest time slice in \code{z_in}, it
#' will be extrapolated using the rate of change from the period between the
#' last and second-to-last time in \code{z_in}).
#' @param outdir The directory to which files will be written (recursively
#' created if not already in existence) (character string).
#' @param prefix The output files will have pattern \code{prefix_x.tif}, where
#' \code{x} is the position (timestep) (potentially multiple digits), and
#' \code{prefix} is a string that you specify here (character string).
#' @param progress Logical. Show progress bar?
#' @param write_change Logical. Should the change grids (which define the change
#' in cells' values per timestep between each pair of time slices) be written
#' to \code{outdir}.
#' @param return_stack Logical. Should the interpolated grids (at timesteps
#' \code{z_out}) be returned as a \code{RasterStack}?
#' @param ... Additional arguments passed to writeRaster.
#' @return Unless \code{return_stack} is \code{TRUE}, \code{interpolate}
#' creates files in \code{outdir}, but returns \code{NULL}, invisibly, to R.
#' If \code{return_stack} is \code{TRUE}, a \code{RasterStack} of
#' interpolated \code{Raster} objects is returned to R.
#' @details Extrapolation below the earliest grid in \code{s} is currently not
#' implemented. Interpolation is linear between pairs of time slices existing
#' in \code{z_in}, and extrapolation beyond the last slice in \code{s}, if
#' requested, follows the trend observed between the last and second-to-last
#' slice in \code{z_in}.
#' @importFrom raster writeRaster nlayers stack overlay
#' @export
#' @examples
#' library(raster)
#' r1 <- raster(matrix(rnorm(100), ncol=10))
#' r2 <- raster(matrix(rnorm(100, mean=11), ncol=10))
#' s <- stack(r1, r2)
#' s_interp <- interpolate_linear(s, c(2001, 2011), 2001:2020, tempdir(),
#' 'example', return_stack=TRUE)
#' plot(s_interp, zlim=range(values(s_interp)))
interpolate_linear <-
function(s, z_in, z_out, outdir, prefix, progress=TRUE, write_change=TRUE,
return_stack=FALSE, ...) {
if(raster::nlayers(s) != length(z_in))
stop('Length of z_in must equal the number of layers in s.')
if(raster::nlayers(s) < 2)
stop('stack s must have at least 2 layers.')
if(!all(findInterval(z_out, range(z_in), rightmost.closed=TRUE) == 1)) {
if(any(z_out < min(z_in))) {
stop('This function does not extrapolate backwards (i.e. below the ',
'earliest element in z_in). All elements of z_out must be greater ',
'than min(z_in).')
} else {
warning('Some values of z_out require extrapolation beyond the range of ',
'z_in.\nThe rate of change for extrapolation is assumed to be ',
'equal to that for the period between the last and ',
'second-to-last elements of z_in (after temporal ordering).')
}
}
outdir <- normalizePath(sub('/$|\\\\$', '', outdir), winslash='/',
mustWork=FALSE)
if(!file.exists(outdir)) dir.create(outdir, recursive=TRUE)
z_out <- unique(z_out)
if(is.unsorted(z_in)) {
s <- s[[order(z_in)]]
z_in <- sort(z_in)
}
len <- diff(z_in)
base <- findInterval(z_out, z_in)
lower <- unique(base[base < nlayers(s)])
s.change <- raster::stack(
sapply(if(length(lower) > 0) lower else raster::nlayers(s) - 1,
function(x) {
message(
sprintf('Calculating change grid for %s to %s.',
z_in[x], z_in[x+1])
)
raster::overlay(
s[[x]], s[[x+1]], fun=function(x1, x2) (x2-x1)/len[x],
filename=ifelse(
write_change,
file.path(outdir, sprintf('changegrid_%s_%s',
z_in[x], z_in[x+1])),
''), recycle=FALSE, format='GTiff', ...)
})
)
multi <- z_out - z_in[base]
chg.ind <- ifelse(base > raster::nlayers(s.change),
raster::nlayers(s.change), base)
message('Calculating grids for years specified in z_out...')
if(progress) pb <- txtProgressBar(0, length(z_out), style=3)
invisible(sapply(seq_along(z_out), function(x) {
out.rast <- if(z_out[x] %in% z_in) {
s[[base[x]]]
} else {
raster::overlay(s[[base[x]]], s.change[[chg.ind[x]]],
fun=function(x1, x2) x1 + (x2*multi[x]))
}
raster::writeRaster(
out.rast,
filename=file.path(outdir, sprintf('%s_%s', prefix, z_out[x])),
format='GTiff', ...)
if(progress) setTxtProgressBar(pb, x)
}))
if(isTRUE(return_stack)) {
f <- file.path(outdir, paste0(prefix, '_', z_out, '.tif'))
return(raster::stack(
f[order(as.numeric(sub('.*_(\\d+)\\.tif$', '\\1', f)))]
))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.