[R-sig-Geo] temporal interpolation with stackApply - seg fault
Robert J. Hijmans
r.hijmans at gmail.com
Fri Dec 3 06:45:55 CET 2010
Could be a memory problem, but it does not look like something that R
would say. Rather looks like something caused by the C routine of
approxfun. Perhaps it is related to the values in your rasters as I
was able to run the below (on a PC with ample of RAM).
r <- raster(nrow=508, ncol=3825, xmn=-175.4283, xmx=-100.0055,
ymn=59.98306, ymx=70)
s <- list()
for (i in 1:12) {
r[] = runif(ncell(r)) * 65535
n <- unique(round(runif(0.5 * ncell(r)) * ncell(r) ))
r[n] <- NA
s[i] <- r
}
s <- stack(s)
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(s, fun=fill.LST) #, filename = "day.LST.h12v02.2005.tif")
Perhaps you can debug this by doing the equivalent:
b <- brick(s, values=F)
s <- getValues(s)
s <- fill.LST(s)
# and if there is no error:
b <- setValues(b, s)
Best, Robert
On Thu, Dec 2, 2010 at 2:24 PM, Michael Loranty <mloranty at whrc.org> wrote:
> 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
>>>>>
>>>>
>>>
>>>
>>>
>>
>
> _______________________________________________
> 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