[R-sig-Geo] temporal interpolation with stackApply

[Robert J. Hijmans]

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
>