[R-sig-Geo] temporal interpolation with stackApply

[Robert J. Hijmans]

Michael,

I think that is because of a difference in versions. The error should
go away with a more recent version of 'raster'. You can get one from
R-Forge:  install.packages("raster",
repos="http://R-Forge.R-project.org")

The current version on CRAN, 1.6-22, has some glitches, and I just
uploaded version 1.7-6 to CRAN. This version should be available for
install from there within 24 hrs.

Robert


On Wed, Dec 1, 2010 at 6:04 PM, Michael Loranty <mloranty at whrc.org> wrote:
> 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
>>>
>>
>
>
>