[R-sig-Geo] temporal interpolation with stackApply

[Michael Loranty]

Thanks for the quick response, Robert.

I hadn't realized I was misusing the stackApply. For some reason I'm  
getting a similar error with calc.  I don't know how to look at the  
function definition for calc to diagnose further.

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)
+ 			}
 >
 > ss < - calc(s,fun=fill.NA)
Error in approxfun(x, y = NULL, method = "linear", rule = 2) :
   need at least two non-NA values to interpolate

 > s[13]
[1] 13 NA 83 13 83 NA
 > fill.NA(s[13])
[1] 13 48 83 13 83 83
 >

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








On Dec 1, 2010, at 7:34 PM, Robert J. Hijmans wrote:

> Michael,
>
> Thanks for sending an good self-contained example. If you use
> stackApply with layers=1:nlayers you are requesting that the function
> should be applied to each layer separately. That is OK for the first
> case (set.NA) although 'calc' would be more direct. But it does not
> for fill.NA, because that function needs the values of all layers in a
> single computations. In this case you must use calc. The below works:
>
>
> r <- raster(ncol=10, nrow=10)
> r[] <- 1:ncell(r)
>
> s <- stack(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 <- calc(s, fun=set.NA)
>
> fill.NA <- function(x,na.rm = F){
>    fill <- approxfun(x, y = NULL, method = 'linear', rule = 2)
>    rec <- which(is.na(x))
>    x[rec] <- fill(rec)
>    return(x)
> }
>
> ss <- calc(s,fun=fill.NA)
>
>
> Best, Robert
>
> On Wed, Dec 1, 2010 at 1:44 PM, Michael Loranty <mloranty at whrc.org>  
> wrote:
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>