[R-sig-Geo] temporal interpolation with stackApply - seg fault
Michael Loranty
mloranty at whrc.org
Thu Dec 2 23:24:33 CET 2010
Hello again...
I've successfully applied the previously discussed function over a
raster stack using calc. However, the process seems to choke on large
data sets. It ran successfully on a stack with nlayers: 13, nrow:
119, ncol 1048 but aborted with the stack below.
Could this be a memory issue? I'm not sure if this is on my end, or
something to do with the function.
Many thanks,
Mike
> day
class : RasterStack
nlayers : 13
nrow : 508
ncol : 3825
ncell : 1943100
projection : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
+towgs84=0,0,0
min value : 0 0 0 0 0 0 0 0 0 0 ...
max value : 65535 65535 65535 65535 65535 65535 65535 65535 65535
65535 ...
extent : -175.4283, -100.0055, 59.98306, 70 (xmin, xmax, ymin,
ymax)
resolution : 0.01971838, 0.01971838 (x, y)
> fill.LST <- function(x,na.rm = F){
+ x[x==0] <- NA
+ if(length(which(is.na(x)==F))>2){
+ fill <- approxfun(x, y=NULL, method='linear', rule=2)
+ rec <- which(is.na(x))
+ x[rec] <- fill(rec)
+ }
+ x <- (x*0.02)-273.15
+ return(x)
+ }
> day.fill <- calc(day, fun=fill.LST, filename = "LST/
day.LST.h12v02.2005.tif")
*** caught segfault ***
address 0x1020, cause 'memory not mapped'
Traceback:
1: .Call("RGDAL_PutRasterData", raster, rasterData,
as.integer(offset), PACKAGE = "rgdal")
2: putRasterData(raster at file@transient, t(v), band = i, c(0, 0))
3: .writeGDALall(x, filename = filename, format = filetype, ...)
4: .local(x, filename, ...)
5: writeRaster(x, filename, ...)
6: writeRaster(x, filename, ...)
7: .local(x, fun, ...)
8: calc(day, fun = fill.LST, filename = "LST/day.LST.h12v02.2005.tif")
9: calc(day, fun = fill.LST, filename = "LST/day.LST.h12v02.2005.tif")
> sessionInfo()
R version 2.10.0 (2009-10-26)
x86_64-redhat-linux-gnu
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spatial_7.3-1 raster_1.7-5 rgdal_0.6-28 sp_0.9-66
loaded via a namespace (and not attached):
[1] grid_2.10.0 lattice_0.17-26
>
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 10:37 PM, Robert J. Hijmans wrote:
> 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
>>>>
>>>
>>
>>
>>
>
More information about the R-sig-Geo
mailing list