Interpolating rasterstack time series with spline in R

Question

I am trying to fill the gaps in a time series of ndvi images using spline.

I have created a raster stack with the ndvi images I have and some layers with only NA for the timesteps that I don't have. The raster stack is available as a Geotiff here: raster stack download

I have written the following function to interpolate the missing values:

f.int.spline <- function(x) {    # x is a raster stack or brick
              v=as.vector(x)     # transform x in a vector for easier manipulation
              z=which(v %in% NA) # find which pixel values are to be interpolated
              # interpolation with spline
              interp <- spline(x=c(1:NROW(v)), y = v,  
                 xout = z,       # missing values to be interpolated
                 method = "natural") 

x[z] <-interp$y # including the missing values in the original vector

}

The function works if I use it with one pixel (e.g x[ 50, 200 ]), but if I run it with calc(x, f.int.spline) it returns a generic error:

> error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

If I run it using f.int.spline(x) it returns an error related to memory usage:

> Error in as.vector(t((matrix(rep(i, rec2), nrow = rec2, byrow = TRUE)) +  : 
  error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error in t((matrix(rep(i, rec2), nrow = rec2, byrow = TRUE)) + add) : 
  error in evaluating the argument 'x' in selecting a method for function 't': Error: cannot allocate vector of size 4.9 Gb

1) Do you see any flaws or have any workarounds on how to make it work?

2) I can't understand exactly how the calc() function works for raster stacks: Does it take the values of each pixel in all the layers?

3) as suggested by Jeffrey Evans I am also looking for other interpolation functions that are better suited for the job. Any idea ?

Robert Hijmans Robert Hijmans · Accepted Answer · 2016-06-15T21:54:07

First create a function that works on a vector, including on some corner cases (that may or may not be relevant to you)

 f <- function(x) {  
    z <- which(is.na(x))
    nz <- length(z)
    nx <- length(x)
    if (nz > 0 & nz < nx) { 
        x[z] <- spline(x=1:nx, y=x, xout=z, method="natural")$y
    }
    x
}

Test the function

f(c(1,2,3,NA,5,NA,7))
##[1] 1 2 3 4 5 6 7
f(c(NA,NA,5,NA,NA))
##[1] 5 5 5 5 5
f(rep(NA, 8))
##[1] NA NA NA NA NA NA NA NA
f(rep(1, 8))
##[1] 1 1 1 1 1 1 1 1

Then use calc on a RasterStack or RasterBrick

Example data

r <- raster(ncols=5, nrows=5)
r1 <- setValues(r, runif(ncell(r)))
r2 <- setValues(r, runif(ncell(r)))
r3 <- setValues(r, runif(ncell(r)))
r4 <- setValues(r, runif(ncell(r)))
r5 <- setValues(r, NA)
r6 <- setValues(r, runif(ncell(r)))
r1[6:10] <- NA
r2[5:15] <- NA
r3[8:25] <- NA
s <- stack(r1,r2,r3,r4,r5,r6)
s[1:5] <- NA

Use the function

x <- calc(s, f)

Alternatively, you can use approxNA

x <- approxNA(s)

Interpolating rasterstack time series with spline in R

1 Answers