[R-sig-Geo] parallel raster processing with calc and mc2d monte carlo simulation

[spkearney]

Hello all, and thanks in advance for any and all help you can give on this:

I have set up a function to extract the 2.5%, 50% and 97.5% percentiles from
a monte carlo simulation on three rasters that is to be called up using
calc() in the raster package and it works great on a test-sized stack/brick,
thanks to suggestions at this post here:
http://grokbase.com/t/r/r-sig-geo/123cb3daaq/apply-monte-carlo-simulation-for-each-cell-in-a-matrix-originally-raster

My problem, is that I want to run this function on a much larger Raster
Brick that, as written, takes hours to process.  I need to do this multiple
times, so I am trying to speed up the processing using clusterR (or another
option such as rasterEngine with multi-core processing).  However, I can't
get it to work!   Here is the code that works on the test raster brick:

brick <- brick(BC_BA, BC_BA_SE, SlopePer)  ## Stack three rasters into one
Raster Brick
testbrick <- crop(brick, extent(299700, 300100, 1553550, 1553650)) ## Crop
brick to manageable size

ndunc(101)
fun.CROP_AGWBC <- function(x) {
  dBC_BA <- mcdata(x[[1]], type="0")
  dBC_BA_SE <- mcdata(x[[2]], type = "0")
  SlopePer <-x[[3]]
  stBA <- mcstoc(rnorm, type = "U", rtrunc = TRUE, 
                 mean = dBC_BA, sd = dBC_BA_SE, linf = 0, lhs = FALSE)
  BC_AGWBC <- lm.final$coefficients[1] + 
    lm.final$coefficients[2]*stBA + 
    lm.final$coefficients[3]*SlopePer
  AGWBC <- (lambda_DV * BC_AGWBC + 1)^(1/lambda_DV)-1
  quantile(AGWBC[], c(0.025, 0.5, 0.975), na.rm=TRUE)
}
 
CROP_AGWBC <- calc(teststack, fun.CROP_AGWBC) ##Run the calc function
CROP_AGWBC ##Check the result
plot(CROP_AGWBC) ##Plot the three-raster brick result

##Extract the individual raster layers
CROP_AGWBC_PRED <- CROP_AGWBC[[2]]
CROP_AGWBC_LWR <- CROP_AGWBC[[1]]
CROP_AGWBC_UPR <- CROP_AGWBC[[3]]

As I mentioned, the code works great on a small sample.  I tried to speed it
up using clusterR as follows, first testing it on the 'testbrick' Raster
Brick with hopes to use it on the whole Raster Brick:

beginCluster(8)
clusterR(x = testbrick, fun = fun.CROP_AGWBC)

and I get the following error:
[1] "data should be numeric or logical"
attr(,"class")
[1] "snow-try-error" "try-error"     
Error in clusterR(x = testbrick, fun = fun.CROP_AGWBC) : cluster error

It is interesting because, if I try to run it again, I get this error
instead:
Error in as.vector((x[, 1] - 1) * ncol(object) + x[, 2]) : 
  error in evaluating the argument 'x' in selecting a method for function
'as.vector': Error in x[, 2] : subscript out of bounds

I have tried this many different ways, including along the lines of: 
f <- function (x) calc(x, fun.CROP_AGWBC)
y <- clusterR(testbrick, f)

which gives me the same error (more or less) of:
Error in checkForRemoteErrors(lapply(cl, recvResult)) : 
  2 nodes produced errors; first error: data should be numeric or logical

And I have tried using the rasterEngine() function (first without parallel
processing) by changing up the code in two ways, the first being: 
ndunc(101)
fun.CROP_AGWBC <- function(x) {
  dBC_BA <- mcdata(x[[1]], type="0")
  dBC_BA_SE <- mcdata(x[[2]], type = "0")
  SlopePer <-x[[3]]
  stBA <- mcstoc(rnorm, type = "U", rtrunc = TRUE, 
                 mean = dBC_BA, sd = dBC_BA_SE, linf = 0, lhs = FALSE)
  BC_AGWBC <- lm.final$coefficients[1] + 
    lm.final$coefficients[2]*stBA + 
    lm.final$coefficients[3]*SlopePer
  AGWBC <- (lambda_DV * BC_AGWBC + 1)^(1/lambda_DV)-1
  output <- quantile(AGWBC[], c(0.025, 0.5, 0.975), na.rm=TRUE)
  output_array <- array(output,dim=c(dim(x)[1],dim(x)[2],3))
  return(output_array)
}
re <- rasterEngine(x = testbrick, fun = fun.CROP_AGWBC)

which runs but gives me a 3-layer Raster Brick all with NA's or Inf.  The
second thing I tried used the same fun.CROP_AGWBC function as above, but
with the following rasterEngine code to call up the calc formula:
f <- function(x) {reout <- calc(x, fun.CROP_AGWBC)
                  reout_array <- array(getValues(reout),
dim=c(dim(x)[1],dim(x)[2],3))
                  return(reout_array)
}
re <- rasterEngine(x = testbrick, fun = f, chunk_format = "raster")

which gives me the following error, even though I thought I converted the
output to an array:
chunk processing units require array vector outputs.  Please check your
function.
Error in focal_hpc_test(x, fun, window_center, window_dims, args,
layer_names,  :

So, my questions are as follows:
*1) Does anyone know why the clusterR does not work for the calc() function
in my first attempt?  I imagine it has something to do with the conversion
of rasters to mcnodes in the function, but can't figure it out!  Any
suggestions?

2) Any thoughts on why I can't get this to work with the rasterEngine()
function?  I am converting the outputs to arrays with the same dimensions as
the input file, but still no luck.*

Again, any help is much appreciated.  Any suggestions for improving this
question are welcome and I'll do my best to update it - this is my first
post!

Kind Regards,
sean





--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/parallel-raster-processing-with-calc-and-mc2d-monte-carlo-simulation-tp7587901.html
Sent from the R-sig-geo mailing list archive at Nabble.com.