1

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
Danple
  • 111
  • 2
  • 6

1 Answers1

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

}
lbusett
  • 5,801
  • 2
  • 24
  • 47
  • 1
    Thanks! Turns out the function wasn't working because of the for loop. Turns out `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. – Danple Feb 20 '17 at 21:34
  • glad that you solved it. However, keep in mind the warning in the help: if you are trying to subtract the average from each layer and `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). – lbusett Feb 20 '17 at 22:30