Introduction

In this tutorial, we show how to simulated two dimensional stationary Gaussian Random Fields.

Its objectives are two folds:

All the fields are simulated on the domain \([0, 1]^2\) discretized on a regular grid.

To evaluate the conditioning algorithms, a realization of the random field, a non conditional simulation, is observed at points, uniformly scattered on the domain. These observations are used to condition a simulation.

Initialization and parameters

# simulation parameters
seed = 123
set.seed (seed)
# number of simulations for Monte-Carlo estimators
nsim = 10

# number of observations
ndat = 50

# parameters of the grid on [0,1]^2
nx = 101*c(1,1)
dx = 1/(nx-1)
x0 = c(0,0)

# parameters for the simulation algorithms
ntb = 1000 # number of band for the Turning Bands method (simtub)
tau = 0.01 # std. of the model observation process for SPDE (not used as default value is used instead)
matern_param = 1.0

Definition of the grid and the observation

# observations
dat = Db_create()
err = dat$setColumn(tab = runif(n=ndat), name = "x")
err = dat$setColumn(tab = runif(n=ndat), name = "y")
err = dat$setLocators(names = c("x", "y"), locatorType = ELoc_X(), cleanSameLocator = TRUE)
err = dat$setColumn(tab = rnorm(ndat), name = "Z", locatorType = ELoc_Z())
# grid for the simulations
grid    = DbGrid_create(nx,dx,x0)

# limits
limits <- list(XP = 1.0*c(0, 1, 1, 0, 0),
               YP = 1.0*c(0, 0, 1, 1, 0))
p_lim   = PolyElem(limits$XP, limits$YP)
pol_lim = Polygons()
err     = pol_lim$addPolyElem(p_lim)

# displays
p = ggDefaultGeographic() + 
    plot.grid(grid) +
    plot.polygon(poly = pol_lim, color = "orange", fill = NA) +
    plot.point(dat, color = "red") +
    plot.decoration(xlab = "Easting", ylab = "Northing", 
                    title = "Limits of the simulations")
ggPrint(p)

Auxiliary functions

# ------------------------------------------------------------
# QC of simulations
# ------------------------------------------------------------
qc_simulations <- function(grid, names, var_model, title, stat_keys) {
nsim = length(names)

# statistics
tab = dbStatisticsMono(grid, names = names, opers = EStatOption_fromKeys(stat_keys))
tab$setTitle(paste0(title, ": statistics"))
tab$display()

# histogram of simulation  #1
i = 1
hist(grid[names[i]], xlab = names[i], main = title, nc = 25, col = "orange", bg = "gray")
# p = ggplot() + 
#     plot.hist(grid, name = names[i], useSel = TRUE, col = "gray", fill = "orange") +
#     plot.decoration(xlab = names[i], title = title)
# ggPrint(p)

# base map
p = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = names[i], flagLegendRaster = TRUE, legendNameRaster = "Y",
              palette = "Spectral") +
    # plot.point(dat, color = "red") +
    plot.decoration(xlab = "Easting", ylab = "Northing", 
                  title = title)
ggPrint(p)

# variogram computations
nlags = grid$getNXs()[1]/2
varioparam = VarioParam_createMultipleFromGrid(grid, npas = nlags)
vario_list = list()
for (i in 1:nsim) {
  vario = Vario(varioparam, grid)
  err   = grid$setLocator(name = names[i], locatorType = ELoc_Z(), cleanSameLocator = TRUE)
  err   = vario$compute()
  vario_list[[1 + length(vario_list)]] <- vario
}

# plot of mean variogram per direction
err = qc_variogram(var_list = vario_list, var_param = varioparam, var_model = var_model, titre = title)

NULL
}

# ------------------------------------------------------------
# plot of mean variogram per direction
# ------------------------------------------------------------
qc_variogram <- function(var_list, var_param, var_model, titre) {
nsim = length(var_list)
ndir = var_param$getDirectionNumber()
# variogram matrices
res_gg = list()
for (idir in 1:ndir) {
  nlags = var_param$getDirParam(idir-1)$getLagNumber()
  gg_sim = matrix(NaN, nrow = nlags-1, ncol = nsim)
  for (i in 1:nsim) {
    gg_sim[,i] = var_list[[i]]$getGgVec(idir = idir-1, ivar = 0, jvar = 0)
  }
  res_gg[[1 + length(res_gg)]] <- gg_sim
}

# plot for each direction
for (idir in 1:ndir) {
  
  nlags = var_param$getDirParam(idir-1)$getLagNumber()
  hh = dx[idir]*(1:(nlags-1))
  v_mean = apply(X = res_gg[[idir]], MARGIN = 1, FUN = mean)
  v_sd   = apply(X = res_gg[[idir]], MARGIN = 1, FUN = sd)
  
  # initial plot
  dir = var_param$getDirParam(idir - 1)$getCodirs()
  sdir = sprintf("Direction = [%1d, %1d]", dir[1], dir[2])
  plot(NULL, NULL, xlim = c(0, dx[idir]*nlags), 
       ylim = 1.2*c(0, var_model$getTotalSill(ivar = 0, jvar = 0)),
       xlab = "Distance", ylab = "Variogram", main = paste0(titre, "\n", sdir),
       xaxs="i", yaxs="i")
  abline(h = 0, lty = 1, col = "black")
  abline(v = 0, lty = 1, col = "black")
  
  # experimental variogram
  lines(hh, v_mean, col = "orange", lw = 2)
  lines(hh, v_mean + 2 * v_sd, col = "orange", lty = 2)
  lines(hh, v_mean - 2 * v_sd, col = "orange", lty = 2)
  for (s in 1:nsim) {
    lines(hh, res_gg[[idir]][,s] , col = "gray", lty = 3)
  }
  
  # variogram model
  abline(h  = var_model$getTotalSill(ivar = 0, jvar = 0), lty = 2, col = "skyblue")
  lines(hh, var_model$getTotalSill(ivar = 0, jvar = 0) - var_model$sample(hh = hh, ivar = 0, jvar = 0, 
        codir = var_param$getDirParam(idir - 1)$getCodirs()), 
        lty = 1, lw = 2, col = "skyblue")
  
  # legend TODO
  legend("bottomright", 
         legend = c("model", "mean variogram", "+/- 2 x Std.", "empirical variograms"),
         lty = c(1, 1, 2, 3), lw = c(2, 2, 1, 1), col = c("skyblue", "orange", "orange", "gray")
  )
} # loop over the directions
  invisible()
}

A single structure

  err = grid$deleteColumns(names = c("*.nc.simu.*","kriging.Z*", 
                                     "SPDE.Z.kriging", "SPDE.Z.condSimu.*", "SPDE.Z.*"))
# a simple Matern's covariance model
  model = Model()
  err = model$addCovFromParam(type = ECov_fromKey("BESSEL_K"), sill = 1, ranges = 1*c(0.1, 0.3), 
                              angles = c(30, 0), param = 1.0)

Non-conditional simulations

  # non conditional simulations
  err = simtub(dbin = NULL, dbout = grid, model = model, nbsimu = nsim, seed=seed, 
               neigh = NULL, nbtuba = ntb, namconv = NamingConvention("TB.nc.simu"))
  err = simulateSPDE(dbin = NULL, dbout = grid, model = model, nbsimu = nsim, seed=seed, 
               namconv = NamingConvention("SPDE.nc.simu"))
  list_title = c("Using Turning Bands", "Using SPDE approach")
  list_key   = c("TB", "SPDE")
  for (i in seq_along(list_title)) {
    err = qc_simulations(grid = grid, 
                         names = paste(list_key[i], "nc.simu", 1:nsim, sep = "."), 
                         var_model = model, title = list_title[i], stat_keys = stat_keys)
  }
##                   Number    Minimum    Maximum       Mean   St. Dev.
##  TB.nc.simu.1  10201.000     -3.536      3.309     -0.182      1.059
##  TB.nc.simu.2  10201.000     -2.902      3.549      0.044      1.005
##  TB.nc.simu.3  10201.000     -3.613      3.127      0.141      0.993
##  TB.nc.simu.4  10201.000     -2.744      2.571     -0.106      0.889
##  TB.nc.simu.5  10201.000     -3.337      3.464      0.119      1.027
##  TB.nc.simu.6  10201.000     -2.537      3.722      0.085      0.950
##  TB.nc.simu.7  10201.000     -2.765      3.772      0.228      1.026
##  TB.nc.simu.8  10201.000     -3.201      3.629     -0.029      0.993
##  TB.nc.simu.9  10201.000     -3.220      2.962     -0.128      0.831
## TB.nc.simu.10  10201.000     -4.238      3.951     -0.124      1.045

##                     Number    Minimum    Maximum       Mean   St. Dev.
##  SPDE.nc.simu.1  10201.000     -3.449      3.703     -0.264      1.071
##  SPDE.nc.simu.2  10201.000     -3.149      2.559     -0.139      0.880
##  SPDE.nc.simu.3  10201.000     -3.304      3.076     -0.018      1.003
##  SPDE.nc.simu.4  10201.000     -3.679      2.686     -0.185      0.976
##  SPDE.nc.simu.5  10201.000     -3.176      3.336     -0.025      1.041
##  SPDE.nc.simu.6  10201.000     -3.555      3.391     -0.244      0.950
##  SPDE.nc.simu.7  10201.000     -2.535      3.233      0.238      0.930
##  SPDE.nc.simu.8  10201.000     -3.462      2.805     -0.308      1.064
##  SPDE.nc.simu.9  10201.000     -3.501      3.082     -0.060      1.013
## SPDE.nc.simu.10  10201.000     -3.554      3.561      0.068      1.043

Kriging

# simple kriging
err = model$delDrift(rank = 0)
err = model$setMean(mean = 0.0)
neigh = NeighUnique()
# kriging in unique neighborhood
  err = kriging(dbin = dat, dbout = grid, model = model,
               neigh = neigh, 
               flag_est = TRUE, flag_std = TRUE, flag_varz = TRUE, 
               namconv = NamingConvention("kriging"))

  err = krigingSPDE(dbin = dat, dbout = grid, model = model,
               flag_est = TRUE, flag_std = FALSE, flag_varz = FALSE, 
               namconv = NamingConvention("SPDE"))

  err = simulateSPDE(dbin = dat, dbout = grid, model = model,
                     nbsimu = nsim, seed = seed, namconv = NamingConvention("SPDE"))

  # cross plot for estimates (SPDE vs. kriging)
  p = ggplot() + 
    plot.correlation(grid, namex = "kriging.Z.estim", namey = "SPDE.Z.estim", 
                     flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE, bins=100) +
    plot.decoration(xlab = "Simple Kriging", ylab = "SPDE estimate", title = "Mono structure")
  ggPrint(p)

  # cross plot of kriging std. (SPDE vs. kriging)
  nm_sims = grid$getNames(names = "SPDE.Z.*")
  tab = matrix(grid$getColumns(names = nm_sims) , nrow = grid$getSampleNumber(), ncol = length(nm_sims))
  estim = apply(X = tab, MARGIN = 1, FUN = mean)
  stdev = apply(X = tab, MARGIN = 1, FUN = sd)
  grid["SPDE.ZMC.estim"] = estim
  grid["SPDE.ZMC.stdev"] = stdev

  # conditional simulations
  p1 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.1", palette = "Spectral")
  p2 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.2", palette = "Spectral")
  p3 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.3", palette = "Spectral")
  p4 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.4", palette = "Spectral")
  ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2)

  # estimate by MC
  p = ggplot() + 
    plot.correlation(grid, namex = "SPDE.Z.estim", namey = "SPDE.ZMC.estim", 
                     flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE, bins=100) +
    plot.decoration(xlab = "SPDE Kriging", ylab = "SPDE conditional exp.", title = "Mono structure")
  ggPrint(p)

  # stdev
  p = ggplot() + 
    plot.correlation(grid, namex = "kriging.Z.stdev", namey = "SPDE.ZMC.stdev", 
                     flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE, bins=100) +
    plot.decoration(xlab = "Simple Kriging Std.", ylab = "SPDE conditional Std.", 
                    title = "Mono structure")
  ggPrint(p)

  # estimator
  p1 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "kriging.Z.estim", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "kriging") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  p2 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "SPDE.Z.estim", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "SPDE") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  ggarrange(p1, p2, nrow = 1, ncol = 2)

  # standard deviation
  p1 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "kriging.Z.stdev", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "kriging") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  p2 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "SPDE.ZMC.stdev", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "SPDE") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  ggarrange(p1, p2, nrow = 1, ncol = 2)

Linear model of covariances

  err = grid$deleteColumns(names = c("*.nc.simu.*","kriging.Z*", 
                                     "SPDE.Z.kriging", "SPDE.Z.condSimu.*", "SPDE.Z.*"))
# a covariance model defined as a sum of two Matern's covariance functions and a nugget effect 
  model = Model()
  err = model$addCovFromParam(type = ECov_fromKey("BESSEL_K"), sill = 0.4, ranges = 1*c(0.1, 0.3), 
                              angles = c(30, 0), param = 1.0)
  err = model$addCovFromParam(type = ECov_fromKey("BESSEL_K"), sill = 0.5, ranges = 2*c(0.1, 0.3), 
                              angles = c(30, 0), param = 2.0)
  err = model$addCovFromParam(ECov_NUGGET(), sill=0.1)

Non-conditional simulations

  # non conditional simulations
  err = simtub(dbin = NULL, dbout = grid, model = model, nbsimu = nsim, seed=seed, 
               neigh = NULL, nbtuba = ntb, namconv = NamingConvention("TB.nc.simu"))
  err = simulateSPDE(dbin = NULL, dbout = grid, model = model, nbsimu = nsim, seed=seed, 
               namconv = NamingConvention("SPDE.nc.simu"))
  list_title = c("Using Turning Bands", "Using SPDE approach")
  list_key   = c("TB", "SPDE")
  for (i in seq_along(list_title)) {
    err = qc_simulations(grid = grid, 
                         names = paste(list_key[i], "nc.simu", 1:nsim, sep = "."), 
                         var_model = model, title = list_title[i], stat_keys = stat_keys)
  }
##                   Number    Minimum    Maximum       Mean   St. Dev.
##  TB.nc.simu.1  10201.000     -3.736      3.410     -0.052      1.000
##  TB.nc.simu.2  10201.000     -2.932      4.197      0.615      1.124
##  TB.nc.simu.3  10201.000     -4.448      3.034     -0.074      1.038
##  TB.nc.simu.4  10201.000     -2.385      2.967      0.125      0.864
##  TB.nc.simu.5  10201.000     -2.922      3.147      0.205      0.951
##  TB.nc.simu.6  10201.000     -3.565      3.197     -0.141      0.912
##  TB.nc.simu.7  10201.000     -3.410      3.452      0.298      1.001
##  TB.nc.simu.8  10201.000     -3.822      3.094     -0.206      1.051
##  TB.nc.simu.9  10201.000     -3.700      3.307     -0.042      1.006
## TB.nc.simu.10  10201.000     -3.265      3.688      0.007      0.903

##                     Number    Minimum    Maximum       Mean   St. Dev.
##  SPDE.nc.simu.1  10201.000     -3.343      3.215     -0.165      0.945
##  SPDE.nc.simu.2  10201.000     -3.605      2.413     -0.473      0.917
##  SPDE.nc.simu.3  10201.000     -3.843      3.687      0.096      1.133
##  SPDE.nc.simu.4  10201.000     -2.950      2.697     -0.070      0.938
##  SPDE.nc.simu.5  10201.000     -3.525      3.507      0.222      0.900
##  SPDE.nc.simu.6  10201.000     -2.883      4.069      0.107      0.897
##  SPDE.nc.simu.7  10201.000     -3.826      3.046     -0.289      1.056
##  SPDE.nc.simu.8  10201.000     -3.169      3.941     -0.068      1.046
##  SPDE.nc.simu.9  10201.000     -3.439      3.117      0.005      1.088
## SPDE.nc.simu.10  10201.000     -2.945      3.454     -0.068      0.836

Kriging

# simple kriging
err = model$delDrift(rank = 0)
err = model$setMean(mean = 0.0)
neigh = NeighUnique()
# kriging in unique neighbourhood
  err = kriging(dbin = dat, dbout = grid, model = model,
               neigh = neigh, 
               flag_est = TRUE, flag_std = TRUE, flag_varz = TRUE, 
               namconv = NamingConvention("kriging"))

  err = krigingSPDE(dbin = dat, dbout = grid, model = model,
               flag_est = TRUE, flag_std = FALSE, flag_varz = FALSE, 
               namconv = NamingConvention("SPDE"))

  err = simulateSPDE(dbin = dat, dbout = grid, model = model,
                     nbsimu = nsim, seed = seed, namconv = NamingConvention("SPDE"))

  # cross plot for estimates (SPDE vs. kriging)
  p = ggplot() + 
    plot.correlation(grid, namex = "kriging.Z.estim", namey = "SPDE.Z.estim", 
                     flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE, bins=100) +
    plot.decoration(xlab = "Simple Kriging", ylab = "SPDE estimate", title = "Multi structure")
  ggPrint(p)

  # cross plot of kriging std. (SPDE vs. kriging)
  nm_sims = grid$getNames(names = "SPDE.Z.*")
  tab = matrix(grid$getColumns(names = nm_sims) , nrow = grid$getSampleNumber(), ncol = length(nm_sims))
  estim = apply(X = tab, MARGIN = 1, FUN = mean)
  stdev = apply(X = tab, MARGIN = 1, FUN = sd)
  grid["SPDE.ZMC.estim"] = estim
  grid["SPDE.ZMC.stdev"] = stdev

  # conditional simulations
  p1 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.1", palette = "Spectral")
  p2 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.2", palette = "Spectral")
  p3 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.3", palette = "Spectral")
  p4 = ggDefaultGeographic() + plot.grid(grid, nameRaster = "SPDE.Z.4", palette = "Spectral")
  ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2)

  # estimate by MC
  p = ggplot() + 
    plot.correlation(grid, namex = "SPDE.Z.estim", namey = "SPDE.ZMC.estim", 
                     flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE, bins=100) +
    plot.decoration(xlab = "SPDE Kriging", ylab = "SPDE conditional exp.", title = "Mono structure")
  ggPrint(p)

  # stdev
  p = ggplot() + 
    plot.correlation(grid, namex = "kriging.Z.stdev", namey = "SPDE.ZMC.stdev", 
                     flagDiag = TRUE, diagColor = "red", flagSameAxes = TRUE, bins=100) +
    plot.decoration(xlab = "Simple Kriging Std.", ylab = "SPDE conditional Std.", 
                    title = "Multi structure")
  ggPrint(p)

  # estimator
  p1 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "kriging.Z.estim", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "kriging") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  p2 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "SPDE.ZMC.estim", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "SPDE") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  ggarrange(p1, p2, nrow = 1, ncol = 2)

  # standard deviation
  p1 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "kriging.Z.stdev", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "kriging") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  p2 = ggDefaultGeographic() + 
    plot.grid(grid, nameRaster = "SPDE.ZMC.stdev", 
              palette = "Spectral", flagLegendRaster = TRUE, legendNameRaster = "SPDE") +
    plot.decoration(xlab = "Easting", ylab = "Northing")
  ggarrange(p1, p2, nrow = 1, ncol = 2)

References