Search code examples
rrasterzoo

Efficient way to interpolate time series of raster bricks in R

I have some large raster bricks that represent time series (each layer is one time point). I wish to increase the time resolution by interpolating between layers (in the example below interpolation is using zoo::na.spline, although I am open to other interpolation methods if they help). I can't use na.spline directly on the bricks - it's not designed for that and just gives Error in as.array.default(X) : 'dimnames' applied to non-array.

At the moment I am using a series of loops to do this (1 loop to generate a brick with more time layers, a 2nd loop to insert the original brick into the correct layers of this, and a 3rd loop to use na.spline on each cell). However, this is horribly slow on large bricks and seems to be a rather inefficent approach.

A minimal reproducible example.

First lets make an intial raster brick representing the low time resolution data. Note that the top left cell is always NA, as any solution must be robust for cells that only contain NA:

library(raster); library(rasterVis); library(zoo)
r1 = raster(matrix(c(NA, rep(1,8)), 3, 3))
r2 = raster(matrix(c(NA, rep(2,8)), 3, 3))
r3 = raster(matrix(c(NA, rep(3,8)), 3, 3))    
b1 = brick(list(r1,r2,r3))
levelplot(b1)

enter image description here

Now, lets create an empty raster brick in which to insert the values:

b2 = brick(lapply(1:9, function(x) raster(matrix(NA, 3, 3))))

Next we insert the layers of b1 into the appropriate layers of b2. Note that I even resort to a loop here because subassignment of selected layers into raster bricks does not work.

old.layers = c(1,4,7)
for (i in 1:nlayers(b1)) {
  b2[[old.layers[i]]] = b1[[i]]
} 
levelplot(b2, layout=c(9,1))

enter image description here

Finally we loop over each cell to do the time series interpolation and insert the results back into b2:

for (cell in 1:ncell(b2)) {
  if (all(is.na(as.vector(b2[cell])))) {
    # na.spline wil fail if all are NA
    new.values = as.numeric(rep(NA, dim(b2)[3])) 
  } else {
    new.values = na.spline(as.vector(b2[cell]), na.rm = F) 
  }
  b2[cell][] = new.values
}
levelplot(b2, layout=c(9,1))

enter image description here

This works but seems to be extremely inefficient and inelegant. Is there a faster (and preferably also more elegant) way to do this?

Solution

  • You can use raster::calc with a function like this:

    yy <- rep(NA, 9)
fun <- function(y) { 
  if (all(is.na(y)))  yy else na.spline((yy[c(1,4,7)] <- y), xout=1:9)
}
b2 <- calc(b1,  fun)

    For reference, raster::approxNA comes close to what you want (but it does not add layers, and uses linear interpolation).