[R-sig-Geo] temporal interpolation with stackApply
Michael Loranty
mloranty at whrc.org
Thu Dec 2 03:04:35 CET 2010
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
>>
>
More information about the R-sig-Geo
mailing list