У меня есть простое уравнение для вычисления "CL" из суммы четырех различных сеток:
CL = A + B + C + D
Я пытаюсь использовать пакет mc2d
, чтобы выразить неопределенность результатов CL на основе неопределенностей в каждой из входных сеток.
Неопределенность в каждой входной сетке выражается на основе равномерных распределений с диапазонами неопределенности:
+/- 6 для растра A
+/- 10 для растра B
+/- 12 для растра C
+/- 5 для растра D
Я хотел бы запустить моделирование Монте-Карло, при котором входные сетки суммируются несколько раз (например, n = 1000) на основе значений для каждой ячейки сетки, которые случайным образом выбираются из этих диапазонов неопределенности.
Код ниже не работает:
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
cannot use this function
Я не знаю, как исправить эту ошибку. Любые предложения будут высоко ценится.
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)