Second order stationary random functions

Selection of the covariance model. It can only be: - GAUSSIAN - EXPONENTIAL - MATERN

idx_type = 3
type_cov <- c("GAUSSIAN", "EXPONENTIAL", "MATERN")[idx_type]
nu <- c(10, 1/2, 3/4)[idx_type]
ranges = c(5, 10, 15)
angles = c(30., 0, 0)

# initialization of the simulator
nb <- 1000 # number of spectral components
ns <- 20   # number of simulations
seed <- 42 # Seed for the first simulation

#var_name = paste(prefix, "V1", "S1", sep = ".")
var_name = NC_getNameEncoded(prefix, NULL, 1, 1, 1, ns)

# Define the title (used for plots)
title = paste0(type_cov, " model")
if (type_cov == "MATERN") {title = paste0(title, " (nu = ", nu, ")")}

1D simulation of a random function

ndim = 1
# definition of the Space and creation of the grid
err = defineDefaultSpace(ESpaceType_RN(), ndim)
nx   = c(10000)
dx   = c(1.)
grid = DbGrid_create(nx = nx, dx = dx)

# specification of the model using Model
mod  = Model_createFromParam(type = ECov_fromKey(type_cov),
                             ranges = ranges[1:ndim], param = nu, flagRange=FALSE)
# using the *SimuSpectralTest* interface
sim = SimuSpectralRN(1, nb, 100, 123, TRUE)
err = sim$setModelGeneric(mod)
err = sim$simulateSpectralTest()
## Simulation of the spectrum
## - Space dimension   = R1
## - Number of variables  = 1
## - Number of spectral components = 1000
# computing the mathematical expression
gamma = sim$getSpectrum()$getGamma()$toTL()
omega = sim$getSpectrum()$getOmega()$toTL()
phi   = sim$getSpectrum()$getPhi()
coor  = grid$getColumnsAsMatrix(names = c("x1"), useSel = FALSE)$toTL()
val   = sqrt(2) * as.numeric(t(gamma) %*% cos(omega %*% t(coor) + phi))
# computing using the gstlearn functions
err = sim$computeSpectralTest(dbout = grid)
## Spectral Simulation on a set of Isolated Points
## - Number of samples = 10000
stopifnot (all(abs(val - grid["Simu"]) < 1e-6))

# using the simuSpectral interface
err  = simuSpectral(NULL, dbout = grid, model = mod, nbsimu = ns, ns = nb, seed=seed,
                    namconv = NamingConvention(prefix))
# statistics
knitr::kable(dbStatisticsMono(db = grid, names = paste(prefix, "*", sep = "."), opers = opers)$toTL(),
             digits = 3, caption = title)
MATERN model (nu = 0.75)
Number Minimum Maximum Mean St. Dev.
Simu 10000 -3.368 3.784 0.030 0.967
Simu.S1 10000 -3.342 3.376 0.002 0.971
Simu.S2 10000 -3.454 3.717 -0.009 0.947
Simu.S3 10000 -3.341 3.376 0.002 0.971
Simu.S4 10000 -3.150 3.162 0.010 0.961
Simu.S5 10000 -3.486 3.410 -0.012 0.973
Simu.S6 10000 -3.151 3.151 0.010 0.961
Simu.S7 10000 -3.486 3.410 -0.012 0.973
Simu.S8 10000 -3.331 3.696 0.007 0.979
Simu.S9 10000 -3.916 3.849 -0.003 0.938
Simu.S10 10000 -3.331 3.696 0.007 0.979
Simu.S11 10000 -3.916 3.849 -0.003 0.938
Simu.S12 10000 -3.331 3.696 0.007 0.979
Simu.S13 10000 -3.787 3.850 -0.002 0.935
Simu.S14 10000 -3.331 3.696 0.007 0.979
Simu.S15 10000 -3.287 3.623 0.014 0.934
Simu.S16 10000 -3.391 3.698 0.007 0.980
Simu.S17 10000 -3.287 3.623 0.014 0.934
Simu.S18 10000 -3.817 3.238 -0.028 0.948
Simu.S19 10000 -3.287 3.623 0.014 0.934
Simu.S20 10000 -3.850 3.279 -0.028 0.948 
# histogram
p = plot.init() +
  plot.hist(db = grid, name = var_name, bins = 25, col = "skyblue") +
  plot.decoration(title = title)
plot.end(p)

# experimental variogram for the first simulation and simulated model
nlag = 50
varioPar = VarioParam_createMultipleFromGrid(grid, nlag = nlag)
err = grid$setLocator(name = var_name, locatorType = ELoc_Z(),
                      cleanSameLocator = TRUE)
vario = Vario_computeFromDb(varioPar, grid)
p = plot.init() +
  plot.varmod(vario = vario, model = mod, flagLegend=TRUE) +
  plot.decoration(title = title)
plot.end(p)

# variograms on all simulations
vl = list()
for (s in 1:ns) {
  err = grid$setLocator(name = NC_getNameEncoded(prefix, NULL, 1, 1, s, ns),
                        locatorType = ELoc_Z(), cleanSameLocator = TRUE)
  vl[[1 + length(vl)]] <- Vario_computeFromDb(varioPar, grid)
}
err = compare_variograms(var_list = vl, title = title)
hh = as.matrix((0:nlag)*dx)
gg = mod$sample(hh, mode=CovCalcMode_create(asVario=TRUE))
lines(hh, gg, col = "skyblue" )

2D simulation of a random function

ndim = 2
err = defineDefaultSpace(ESpaceType_RN(), ndim)

# specification of the model (theoretical range is provided in 'a')
mod  = Model_createFromParam(type = ECov_fromKey(type_cov),
                             ranges = ranges[1:ndim], param = nu, angles=angles[1:ndim],
                             flagRange=FALSE)

# Two intersecting grids
grid1 = DbGrid_create(nx = c(256, 256))
grid2 = DbGrid_create(nx = c(100, 100), x0 = c(200, 200), dx = c(1,1))

# computing the simulation on two intersecting grids
err = simuSpectral(NULL, dbout=grid1, model=mod, ns=nb, seed=seed , namconv = NamingConvention(prefix))
err = simuSpectral(NULL, dbout=grid2, model=mod, ns=nb, seed=seed , namconv = NamingConvention(prefix))
var_name = NC_getNameEncoded(prefix, NULL)

p = plot.init(asp=1) +
  plot.raster(grid1, name = var_name, palette = "Spectral", limits = 3.5 * c(-1,1)) +
  plot.raster(grid2, name = var_name, palette = "Spectral", limits = 3.5 * c(-1,1)) +
  plot.decoration(xlab = "Easting", ylab = "Northing", title = title)
plot.end(p)

# Control that the values on the intersection of the two grids are identical
err = migrate(dbin = grid1, dbout = grid2, name = var_name,
              namconv = NamingConvention("grid1"))
p = plot.init() +
  plot.correlation(grid2, namex = var_name, namey = paste("grid1", var_name, sep = "."), asPoint = TRUE,
                   flagDiag = TRUE, flagSameAxes = TRUE) +
  plot.decoration(xlab = "simulation on grid1", ylab = "simulation on grid2",
                  title = title)
plot.end(p)

# computing the simulation on scattered points (extracted from grid1)
grid1["sel"] <- VectorHelper_simulateBernoulli(grid1$getNSample(), proba=0.005)
err = grid1$setLocator("sel", ELoc_SEL(), cleanSameLocator = TRUE)
db = Db()
db["x1"] <- grid1$getColumn(name = "x1", useSel = TRUE)
db["x2"] <- grid1$getColumn(name = "x2", useSel = TRUE)
db[paste("grid1", var_name, sep = ".")] <- grid1$getColumn(name = var_name, useSel = TRUE)
err = db$setLocators(c("x1", "x2"), ELoc_X(), cleanSameLocator = TRUE)
err = grid1$clearSelection()

# computing the simulation on scattered points
err = simuSpectral(NULL, dbout=db, model=mod,  ns=nb, seed=seed , namconv = NamingConvention(prefix))

p1 = plot.init(asp=1) +
  plot.raster(grid1, name = var_name, palette = "Spectral") +
  plot.symbol(db, color = "red") +
  plot.decoration(xlab = "Easting", ylab = "Norhting", title = title)

# control that the values on the intersection of the two grids are identical
p2 = plot.init() +
  plot.correlation(db, namex = var_name, namey = paste("grid1", var_name, sep = "."),
                   asPoint = TRUE,
                   flagDiag = TRUE, flagSameAxes = TRUE) +
  plot.decoration(xlab = "simulation on db", ylab = "simulation on grid1",
                  title = title)

ggarrange(p1, p2, nrow = 1, ncol = 2)

# statistics
tab = rbind(
  dbStatisticsMono(db = grid1, names = var_name, opers = opers)$toTL(),
  dbStatisticsMono(db = grid2, names = c(var_name), opers = opers,
                   flagIso = TRUE)$toTL(),
  dbStatisticsMono(db = db, names = c(var_name), opers = opers,
                   flagIso = TRUE)$toTL())
rownames(tab) = c("grid1", "grid2", "scattered points")
knitr::kable(tab,  digits = 3, caption = title)
MATERN model (nu = 0.75)
Number Minimum Maximum Mean St. Dev.
grid1 65536 -3.764 3.649 -0.048 1.028
grid2 10000 -3.764 3.495 -0.027 1.054
scattered points 329 -3.564 3.492 -0.108 1.107 
# variogram on the simulations
nlag = 25
coldir = c("green", "blue", "red")
err = grid1$setLocator(var_name, ELoc_Z(), cleanSameLocator = TRUE)
varioPar = VarioParam_createMultipleFromGrid(grid1, nlag=nlag)

# computing the experimental variogram
vario = Vario_computeFromDb(varioPar, grid1)

p = plot.init() +
  plot.varmod(vario = vario, model=mod, flagLegend=TRUE) +
  plot.decoration(title = title)
plot.end(p)

3D simulation of a random function

ndim = 3
err = defineDefaultSpace(ESpaceType_RN(), ndim)

# specification of the model (theoretical range is provided in 'a')
mod  = Model_createFromParam(type = ECov_fromKey(type_cov),
                             ranges = ranges[1:ndim], param = nu, angles=angles[1:ndim],
                             flagRange=FALSE)

# Two intersecting grids
nx <- c(100, 100, 20)
dx <- c(1., 1., 1.)
x0 <- c(10, 1000, 2)
grid = DbGrid_create(nx = nx, x0 = x0, dx = dx)

# computing the simulation on the grid
err = simuSpectral(NULL, dbout=grid, model=mod, ns=nb, seed=seed , namconv = NamingConvention(prefix))
# statistics
knitr::kable(dbStatisticsMono(db = grid, names = var_name, opers = opers)$toTL(),
             digits = 3, caption = title)
MATERN model (nu = 0.75)
Number Minimum Maximum Mean St. Dev.
Simu 2e+05 -3.889 4.616 0.064 1.053 
# histogram
p = plot.init() +
  plot.hist(db = grid, name = var_name, bins = 25, col = "skyblue")
plot.end(p)

# map
p = plot.init(asp=1) +
  plot.raster(grid, name = var_name, palette = "Spectral", limits = 3.5 * c(-1,1)) +
  plot.decoration(xlab = "Easting", ylab = "Northing", title = title)
plot.end(p)

# variogram on the simulations
nlag = 20
coldir = c("green", "red", "blue")
err = grid$setLocator(var_name, ELoc_Z(), cleanSameLocator = TRUE)
varioPar = VarioParam_createMultipleFromGrid(grid, nlag=nlag)

# computing the experimental variogram
vario = Vario_computeFromDb(varioPar, grid)

p = plot.init() +
  plot.varmod(vario = vario, model=mod, flagLegend=TRUE) +
  plot.decoration(title = title)
plot.end(p)