Evaluation of the simulation functions (Turning bands vs. SPDE)
gstlearn Team
2023-09-05
Introduction
In this tutorial, we show how to simulated two dimensional stationary Gaussian Random Fields.
Its objectives are two folds:
illustrate how to use the simulation functions with the Turning Bands methods and the SPDE approach
evaluate the reproduction of the first and second statistics, mean and covariance function, defining these fields
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.0Definition 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("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
