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)
  err   = grid$setLocator(name = names[i], locatorType = ELoc_Z(), cleanSameLocator = TRUE)
  err   = vario$compute(grid)
  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
  c00 = var_model$getTotalSill()
  abline(h  = c00, lty = 2, col = "skyblue")
  lines(hh, c00 - var_model$sample(h = hh, 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("MATERN"), 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"))
  lsr = law_set_random_seed(seed)
  err = simulateSPDE(dbin = NULL, dbout = grid, model = model, nbsimu = nsim,
               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.702      3.364     -0.121      1.085
##  SPDE.nc.simu.2  10201.000     -3.227      3.025     -0.108      0.949
##  SPDE.nc.simu.3  10201.000     -2.944      3.242      0.105      0.967
##  SPDE.nc.simu.4  10201.000     -3.139      3.242      0.015      0.951
##  SPDE.nc.simu.5  10201.000     -2.875      3.646      0.357      1.050
##  SPDE.nc.simu.6  10201.000     -3.115      2.953     -0.131      0.924
##  SPDE.nc.simu.7  10201.000     -2.877      3.815      0.454      1.001
##  SPDE.nc.simu.8  10201.000     -3.656      3.082     -0.094      0.892
##  SPDE.nc.simu.9  10201.000     -2.719      3.429     -0.078      0.983
## SPDE.nc.simu.10  10201.000     -3.706      2.940      0.193      0.933

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, nbMC = 0,
               namconv = NamingConvention("SPDE"))
  lsr = law_set_random_seed(seed)
  err = simulateSPDE(dbin = dat, dbout = grid, model = model,
                     nbsimu = nsim, 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)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

  # 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)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

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("MATERN"), sill = 0.4, ranges = 1*c(0.1, 0.3), 
                              angles = c(30, 0), param = 1.0)
  err = model$addCovFromParam(type = ECov_fromKey("MATERN"), 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"))
  lsr = law_set_random_seed(seed)
  err = simulateSPDE(dbin = NULL, dbout = grid, model = model, nbsimu = nsim,
               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.608      3.712      0.136      1.013
##  SPDE.nc.simu.2  10201.000     -3.172      2.867     -0.248      0.807
##  SPDE.nc.simu.3  10201.000     -3.097      3.630      0.092      0.979
##  SPDE.nc.simu.4  10201.000     -3.397      3.069     -0.207      0.975
##  SPDE.nc.simu.5  10201.000     -2.891      3.252     -0.012      0.857
##  SPDE.nc.simu.6  10201.000     -2.777      3.767      0.476      0.845
##  SPDE.nc.simu.7  10201.000     -3.454      2.861      0.023      0.888
##  SPDE.nc.simu.8  10201.000     -3.265      3.219      0.250      0.954
##  SPDE.nc.simu.9  10201.000     -2.389      3.520      0.105      0.854
## SPDE.nc.simu.10  10201.000     -2.702      3.200     -0.003      0.904

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, nbMC = 0,
               namconv = NamingConvention("SPDE"))
  lsr = law_set_random_seed(seed)
  err = simulateSPDE(dbin = dat, dbout = grid, model = model,
                     nbsimu = nsim, 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)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

  # 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)
## Warning: `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.
## `guide_colourbar()` cannot be used for fill_ggnewscale_1.
## ℹ Use one of colour, color, or fill instead.

References