interpolateRasts: Interpolate values from a series of rasters

Description

This function returns a series of rasters interpolated from another series of rasters. For example, the input might represent rasters of a process measured at times t, t + 1, and t + 4. The rasters at t + 2 and t + 3 could be interpolated based on the values in the other rasters. Note that this function can take a lot of time and memory, even for relatively small rasters.

Usage

interpolateRasts(
  rasts,
  interpFrom,
  interpTo,
  type = "linear",
  onFail = NA,
  useRasts = FALSE,
  na.rm = TRUE,
  verbose = TRUE,
  ...
)

Value

A SpatRaster "stack" with one layer per element in interpTo.

Arguments

rasts

A "stack" of SpatRasterss.

interpFrom

Numeric vector, one value per raster in rasts. Values represent "distance" along the set of rasters rasters (e.g., time).

interpTo

Numeric vector, values of "distances" at which to interpolate the rasters.

type

Character. The type of model used to do the interpolation. Note that some of these (the first few) are guaranteed to go through every point being interpolated from. The second set, however, are effectively regressions so are not guaranteed to do through any of the points. Note that some methods cannot handle cases where at least some series of cells have < a given number of non-NA values (e.g., smooth splines will not work if there are < 4 cells with non-NA values).

linear: A model based on linear segments "fastened" at each value of interpFrom. The segments will intersect each value being interpolated from.
spline: A natural splines-based model. Splines will intersect each value being interpolated from.
gam: A generalized additive model. Note that the GAM is not guaranteed to intersect each value being interpolated from. Arguments to gam can be supplied via .... Especially note the family argument! You can use the onFail argument with this method since in some cases gam if there are too few data points.
glm: A generalized linear model. Note that the GLM is not guaranteed to intersect each value being interpolated from. Arguments to gam can be supplied via .... Especially note the family argument (the main reason for why you would use a GLM versus just linear interpolation)! You can use the onFail argument with this method since in some cases glm if there are too few data points.
ns: A natural splines model. Note that the NS is not guaranteed to intersect each value being interpolated from. Arguments to trainNS can be supplied via .... Especially note the family argument and the df argument! If df is not supplied, then the number of splines attempted will be equal to 1:(length(interpFrom) - 1). You can use the onFail argument with this method.
poly: A polynomial model. This method constructs an n-degree polynomial where n = length(interpFrom) - 1. The most parsimonious model is then selected from all possible subsets of models (including an intercept-only model) using AICc. This method is not guaranteed to intersect each value being interpolated from. Arguments to glm can be supplied via .... Especially note the family argument! If family is not supplied, then the response is assumed to have a Gaussian distribution. You can use the onFail argument with this method.
bs: A basis-spline model. This method constructs a series of models with n-degree basis-spline model where n ranges from 3 to length(interpFrom) - 1. The most parsimonious model is then selected from all possible subsets of models (including an intercept-only model) using AICc. This method is not guaranteed to intersect each value being interpolated from. Arguments to glm can be supplied via .... Especially note the family argument! If family is not supplied, then the response is assumed to have a Gaussian distribution. You can use the onFail argument with this method.
smooth.spline: A smooth-spline model (see smooth.spline). This method is not guaranteed to intersect each value being interpolated from. Arguments to smooth.spline can be supplied via .... Unlike some other methods, a family cannot be specified (Gaussian is assumed)! You can use the onFail argument with this method.

onFail

Either NA (default) or any one of 'linear', 'spline', or 'poly'. If a method specified by type fails (i.e., because there are fewer than the required number of values to interpolate from), this method is used in its place. If this is NA and the method fails, then an error occurs.

useRasts

Logical. If FALSE (default), then the calculations are done using arrays. This can be substantially faster than using rasters (when useRasts = TRUE), but also run into memory issues.

na.rm

Logical, if TRUE (default) then ignore cases where all values in the same cells across rasters from which interpolations are made are NA (i.e., do not throw an error). If FALSE, then throw an error when this occurs.

verbose

Logical. If TRUE (default), display progress.

...

Other arguments passed to approx or spline (do not include any of these arguments: x, y, or xout), or to glm, gam, or smooth.spline.

Details

This function can be very memory-intensive for large rasters. It may speed things up (and make them possible) to do interpolations piece by piece (e.g., instead of interpolating between times t0, t1, t2, t3, ..., interpolate between t0 and t1, then t1 and t2, etc. This may give results that differ from using the entire set, however, depending on what type of interpolation is used. Note that using linear and splines will often yield very similar results, except that in a small number of cases splines may produce very extreme interpolated values.

Examples

Run this code


# \donttest{
# NB: The example below can take a few minutes to run.

library(terra)

interpFrom <- c(1, 3, 4, 8, 10, 11, 15)
interpTo <- 1:15
rx <- rast(nrows=10, ncols=10)
r1 <- setValues(rx, rnorm(100, 1))
r3 <- setValues(rx, rnorm(100, 3))
r4 <- setValues(rx, rnorm(100, 5))
r8 <- setValues(rx, rnorm(100, 11))
r10 <- setValues(rx, rnorm(100, 3))
r11 <- setValues(rx, rnorm(100, 5))
r15 <- setValues(rx, rnorm(100, 13))
rasts <- c(r1, r3, r4, r8, r10, r11, r15)
names(rasts) <- paste0('rasts', interpFrom)

linear <- interpolateRasts(rasts, interpFrom, interpTo)
spline <- interpolateRasts(rasts, interpFrom, interpTo, type='spline')
gam <- interpolateRasts(rasts, interpFrom, interpTo, type='gam', onFail='linear')
ns <- interpolateRasts(rasts, interpFrom, interpTo, type='ns', onFail='linear', verbose=FALSE)
poly <- interpolateRasts(rasts, interpFrom, interpTo, type='poly', onFail='linear')
bs <- interpolateRasts(rasts, interpFrom, interpTo, type='bs', onFail='linear')
ss <- interpolateRasts(rasts, interpFrom, interpTo, type='smooth.spline', onFail='linear',
verbose=FALSE)

# examine trends for a particular point on the landscape
pts <- matrix(c(-9, 13), ncol = 2)
pts <- vect(pts)
linearExt <- unlist(terra::extract(linear, pts, ID=FALSE))
splineExt <- unlist(terra::extract(spline, pts, ID=FALSE))
gamExt <- unlist(terra::extract(gam, pts, ID=FALSE))
nsExt <- unlist(terra::extract(ns, pts, ID=FALSE))
polyExt <- unlist(terra::extract(poly, pts, ID=FALSE))
bsExt <- unlist(terra::extract(bs, pts, ID=FALSE))
ssExt <- unlist(terra::extract(ss, pts, ID=FALSE))

mins <- min(linearExt, splineExt, gamExt, nsExt, polyExt, bsExt, ssExt)
maxs <- max(linearExt, splineExt, gamExt, nsExt, polyExt, bsExt, ssExt)

plot(interpTo, linearExt, type='l', lwd=2, ylim=c(mins, maxs), ylab='Value')
lines(interpTo, splineExt, col='blue')
lines(interpTo, gamExt, col='green')
lines(interpTo, nsExt, col='orange')
lines(interpTo, polyExt, col='gray')
lines(interpTo, bsExt, col='magenta')
lines(interpTo, ssExt, col='cyan')

ext <- unlist(extract(rasts, pts, ID = FALSE))
points(interpFrom, ext)

legend('topleft', inset=0.01, lty=c(rep(1, 7), NA),
legend=c('linear', 'spline', 'GAM', 'NS', 'polynomial', 'B-spline',
'Smooth spline', 'Observed'), col=c('black', 'blue', 'green',
'orange', 'gray', 'magenta', 'cyan'), pch=c(rep(NA, 7), 1))

# }

Run the code above in your browser using DataLab