[R-sig-Geo] temporal interpolation with stackApply

[Michael Loranty]

Hi All,

I have a raster stack, that represents a time series of temperature  
observations, and would like to interpolate temporally to fill gaps in  
the data. That is, to apply a gap-filling function across all layers  
in a stack, pixel by pixel. My hope was to accomplish this using  
stackApply, however this doesn't seem to work. The reason seems to be  
associated with how the data are handled within stackApply which is  
based upon the type of function that is used.

Has anyone else encountered a similar problem or come up with a  
solution? It seems to me that it may not be too difficult to make a a  
function to do this based upon few minor adjustments to the current  
raster function. A simple example is below.

Best,
Mike


 > r <- raster(ncol=10, nrow=10)
 > r[]=1:ncell(r)
 > s <- brick(r,flip(r,direction = "x"),flip(r,direction =  
"y"),r,flip(r,direction = "y"),flip(r,direction = "x"))
 >
 > set.NA <- function(x, na.rm = T){x[x==18] <- NA; return(x)}
 >
 > s <- stackApply(s,fun = set.NA, indices = 1:nlayers(s))
 >
 > fill.NA <- function(x,na.rm = F){
+ 			fill <- approxfun(x, y = NULL, method = 'linear', rule = 2);
+ 			new <- fill(1:length(x));
+ 			rec <- which(is.na(x));
+ 			x[rec] <- new[rec];
+ 			return(x)
+ 			}
 >
 > s < - stackApply(s,fun=fill.NA,indices = 1:nlayers(s))
Error in approxfun(x, y = NULL, method = "linear", rule = 2) :
   need at least two non-NA values to interpolate
 >
 > sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-apple-darwin9.8.0

other attached packages:
[1] rgdal_0.6-27  raster_1.6-22 sp_0.9-72



Michael Loranty Ph.D.
Postdoctoral Fellow
The Woods Hole Research Center
149 Woods Hole Rd
Falmouth, MA 02540
508.444.1560
mloranty at whrc.org