1

I am trying to fill NA values of a raster stack with interpNA but the first raster stays empty. I wish to obtain similar results as in substituteNA in a sense that all rasters of the stack are filled.

How can I get around this?

#rasters 
set.seed(123)
library(raster)
library(timeSeries)

r1 <- raster(nrows = 1, ncols = 1, res = 0.5, xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5, vals = 0.3)
rr <- lapply(1:10, function(i) setValues(r1,seq(1,ncell(r1),1)))
for (i in 2:length(rr)){
  rr[[i]] <- rr[[1]]+i
}

#make and stack with NA
s <- stack(rr)
pb = txtProgressBar(min = 0, max = nlayers(s), initial = 0, style = 3) 
for (i in  seq(1,nlayers(s),2)){
  values(s[[i]])[values(s[[i]])%% 2 == 0] <-  NA
  setTxtProgressBar(pb,i)
}
plot(s)

##First raster will be empty:
fun <- function(x) {
  v=as.vector(x)
  z=interpNA(v, type="linear") #substituteNA(v, type="mean")
}
s1<-calc(s, fun)
plot(s1)
  • 1
    Where does interpNA come from? Its not in the raster package (v 2.8.19). Your code won't run as written. – Spacedman May 10 at 7:23
  • Uh, my mistake. I've edited the question. That is from timeSeries. Raster package is (v 2.8.19). – Majid May 10 at 7:29
2

calc works by passing each stack of pixels at each cell location to the function specified. So for the pixel at cell coordinate (1,2) it gets this:

> s[1,2,]
     layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 layer.7 layer.8 layer.9
[1,]      NA       4       5       6       7       8       9      10      11
     layer.10
[1,]       12

That gets fed to interpNA, which will not replace NA at the start or end of a vector:

> interpNA(c(NA,1,2,3,NA,5,6,NA), method="linear")
     [,1]
[1,]   NA
[2,]    1
[3,]    2
[4,]    3
[5,]    4
[6,]    5
[7,]    6
[8,]   NA

So any cell with an NA in layer 1 will end up with an NA in layer 1, and any cell with an NA in the last layer will also be NA - but you don't have any NAs in layer 10 so you don't notice this.

To fix this you need an interpolator that extrapolates to fill in NA values at the start and end. For example you could leverage the na.omit that the timeSeries package has for its timeSeries objects:

> na.omit(timeSeries(c(NA,2,NA,4,NA)),method="ie", interp="linear")

     SS.1
[1,]    2
[2,]    2
[3,]    3
[4,]    4
[5,]    4

Note that it doesn't return 1 for the first NA. Maybe you want to fit a linear regression using the non-missing data and use predictions to fill in the NAs? In any event, check your function passed to calc returns non-NA values when passed a simple vector with some NAs in it before applying it to your stack.

  • Well described answer. Much appreciated. A bit more trouble than I expected but I will smooth the whole stack later on so the fact that it doesn't return 1 for the first NA should be OK I suppose. – Majid May 12 at 3:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.