I need to run calculations on large multi-band rasters and export a RasterBrick, and am trying to do so using the calc() function in the raster package, for the purpose of memory efficiency. The function runs fine on its own, but when I try to include it in calc(), I keep getting this error:
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
cannot use this function
How can I make this work?
Simplified code:
fn = system.file("external/test.grd", package="raster")
s = stack(fn, fn, fn, fn)
out = calc(s, fun = function(x){
for (i in 1:nlayers(x)){
x[[i]] = x[[i]] - cellStats(x[[i]], "min")
x[[i]] = x[[i]]* 5
}
list = unstack(x)
out = brick(list)
return(out)
}
)
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
cannot use this function
1 Answer 1
From calc help:
For large objects calc will compute values chunk by chunk. This means that for the result of fun to be correct it should not depend on having access to all values at once. For example, to scale the values of a Raster* object by subtracting its mean value (for each layer), you would not do, for Raster object x:
calc(x, function(x)scale(x, scale=FALSE))
Because the mean value of each chunk will likely be different. Rather do something like
m <- cellStats(x, 'mean')
x - m
therefore, your function, even if it worked, would probably give you incorrect results. I'm not entirely sure why the function doesn't work, however: maybe there ìs an internal check in calc to avoid the use of cellStats.
To do "your" computation, however, you could use simply:
out = s
for (i in 1:nlayers(s)) {
out [[i]] = (s [[i]] - cellStats(s[[i]], 'min', na.rm = T))*5
}
2 Comments
calc distinguishes bands and indices on its own. Your solution is what I had before, but I need calc because it's more memory efficient. As for cellStats, I thin having NA values was indeed the problem.calc works in chunks, you're going to have wrong outputs (i.e., the value that you subtract to the first group of lines will be different from the one you subtract from the last).