Generate sample (Gaussian distribution) and store simulations in raster stack in R

Question

I have a raster stack, one of the raster layers is the mean prediction of a property for that pixel, the second layer is the variance of the prediction.

I would like to write a function for generating a sample (n=50) from the Gaussian distribution, for each pixel, 50 simulations. And store the output in a raster brick or similar.

When I try directly with the function rnorm I get error message. When I try with overlay (raster), using rnorm as input function, also.

Does anyone have a clue of how doing this?

I did not calculate the spatial variograms, the soil property was predicted simply with regression.

For example:

r <- raster(nrow=10, ncol=10)
s1 <- setValues(r,runif(n = 100, min = 0, max=50))
s2 <- setValues(r,runif(n = 100, min = 1, max=5))
N <- setValues(r,50)

test <- overlay(N,s1,s2,fun=rnorm)
test
class       : RasterBrick 
dimensions  : 10, 10, 100, 3  (nrow, ncol, ncell, nlayers)
resolution  : 36, 18  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       :   layer.1,   layer.2,   layer.3 
min values  : -2.551580, -1.781912, -2.879981 
max values  :  2.533998,  2.245894,  2.828349

test <- rnorm(n=N, mean=s1,sd=s2)
Error in rnorm(n = N, mean = s1, sd = s2) : arguments incorrects

Answer 1

Based on your example, I am not entirely clear as to what you are after. It looks like you would like to generate a random number in a specified normal distribution using rasters to define what the mean and standard deviations are at each cell. This is fairly common in exploring uncertainty of spatial estimates. Perhaps something like this would work.

library(raster)
  r <- raster(nrow=10, ncol=10)
  s1 <- setValues(r,runif(n = 100, min = 0, max=50))
  s2 <- setValues(r,runif(n = 100, min = 1, max=5))

rand.norm <- calc(stack(s1, s2), fun = function(x) 
       { return(rnorm(1, mean = x[[1]], sd = x[[2]])) } )

If you would like multiple realizations, which I imagine would be required here, you could easily put this in a for loop and add the results to a stack or list object.

rand.norm <- stack(s1)
  for( i in 1:50) {
    rand.norm <- addLayer(rand.norm, calc(stack(s1, s2), fun = function(x) 
                 { return(rnorm(1, mean = x[[1]], sd = x[[2]])) } )) 
  }
rand.norm <- dropLayer(rand.norm, 1)  
( rand.var <- overlay(rand.norm, fun=var) )
  plot(rand.var)

You do not necessary have to produce the entire Gaussian distribution, for each cell, at once. An iterative approach, where each random number is generated at i point of the iteration, yields the same results.

You can easily proof this:

x <- rnorm(50)
y <- vector()  
for(i in 1:50) { y <- append(y, rnorm(1) ) }

par(mfrow=c(2,1))
  plot(density(x))
  plot(density(y))