I have a simple equation to calculate the "CL" from the sum of four different grids as:
CL = A + B + C + D
I am attempting to use the mc2d package to express uncertainty in the CL results based on uncertainties in each of the input grids.
The uncertainty in each input grid is expressed based on uniform distributions with uncertainty ranges of:
+/- 6 for raster A
+/- 10 for raster B
+/- 12 for raster C
+/- 5 for raster D
I would like to run a Monte Carlo simulation whereby the input grids are summed multiple times (e.g. n = 1000) based on values for each grid cell that are randomly selected from within these uncertainty ranges
The code below fails with:
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
cannot use this function
I am not sure how to address this error. Any suggestions would be much appreciated.
library(raster)
library(rgdal)
library(mc2d) #load package
# Input grid files
A <- raster(matrix(c(20, 35, 40, 60), nrow=2))
B <- raster(matrix(c(6, 10, 13, 14), nrow=2))
C <- raster(matrix(c(6, 8, 12, 14), nrow=2))
D <- raster(matrix(c(35, 40, 50, 60), nrow=2))
# combine the RasterLayer objects into a RasterStack
s <- stack(A, B, C, D)
# Uncertainty distance for each raster used to establish the
# min/max for the uniform distributions in the function just below
uA <- 6
uB <- 10
uC <- 12
uD <- 5
# Function for generating CL
ndunc(1000)
fun <- function(x) {
# Convert each raster to mcnode object
dA <- mcdata(x[1], type="0")
dB <- mcdata(x[2], type="0")
dC <- mcdata(x[3], type="0")
dD <- mcdata(x[4], type="0")
# Define uniform distirbutions for each raster
stA <- mcstoc(runif, type="U", min=dA-uA, max=dA+uA, rtrunc=TRUE, linf=0)
stB <- mcstoc(runif, type="U", min=dA-uB, max=dA+uB, rtrunc=TRUE, linf=0)
stC <- mcstoc(runif, type="U", min=dA-uC, max=dA+uC, rtrunc=TRUE, linf=0)
stD <- mcstoc(runif, type="U", min=dA-uD, max=dA+uD, rtrunc=TRUE, linf=0)
# Apply Monte Carlo to this equation
CL <- stA + stB + stC + stD
# Extract the min, mean, and max CL from the 1000 iterations
c(min(CL), mean(CL), max(CL))
}
# Need to create rasters for min(CL), mean(CL), and max(CL)
x <- calc(s, fun)
names(x) = c('min', 'mean', 'max')
plot(x)
The error in the function was caused by a typo (missing equals sign) in rtrunc=TRUE