Second order stationary random functions

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

# type_cov <- "GAUSSIAN"; nu = 10
type_cov <- "EXPONENTIAL"; nu = 1/2
# type_cov <- "BESSEL_K"; nu <- 1

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

1D simulation of a random function

# Definition of the Space
err = defineDefaultSpace(ESpaceType_RN(), 1)

# specification of the model (theoretical range is provided in 'a')
a <- c(25)
mod  = Model_createFromParam(type = ECov_fromKey(type_cov), 
                             ranges = a, param = nu, flagRange=FALSE)

# target data base
nx <- c(10000)
dx <- c(1.)
grid = DbGrid_create(nx = nx, dx = dx)

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

sims = SimuSpectral(mod)
for (n in 1:ns) {
  sims$simulate(nb = nb, seed = seed)
  sims$compute(grid, namconv = NamingConvention(paste("S", n, sep = "."))) 
  seed = 0 # The seed for the next simulations is not given anymore
}
# statistics
knitr::kable(dbStatisticsMono(db = grid, names = "S.*", opers = opers)$toTL(),
             digits = 3, caption = title)
EXPONENTIAL model
Number Minimum Maximum Mean St. Dev.
S.1 10000 -3.695 4.277 0.050 1.083
S.2 10000 -3.750 3.538 0.046 1.057
S.3 10000 -3.402 3.222 0.160 0.977
S.4 10000 -3.185 3.479 -0.088 0.997
S.5 10000 -3.337 3.258 0.041 0.987
S.6 10000 -3.817 4.587 0.170 1.012
S.7 10000 -4.207 4.039 0.047 0.932
S.8 10000 -3.409 3.605 0.009 0.998
S.9 10000 -3.808 3.442 -0.050 1.054
S.10 10000 -3.814 3.461 0.023 0.997 
# histogram
p = ggDefault() + 
  plot.hist(db = grid, name = "S.1", bins = 25, col = "skyblue") +
  plot.decoration(title = title)
ggPrint(p)

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

# variograms on all simulations
vl = list()
for (s in 1:ns) {
  err = grid$setLocator(name = paste("S", s, sep = "."), 
                        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

err = defineDefaultSpace(ESpaceType_RN(), 2)

# specification of the model (theoretical range is provided in 'a')
a  = c(10, 20)
angles = c(30., 0)
mod  = Model_createFromParam(type = ECov_fromKey(type_cov), 
                             ranges = a, param = nu, angles=angles,
                             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))

# Initialization of the simulator
seed <- 3132
sims = SimuSpectral(mod)
sims$simulate(nb = nb, seed = seed)
## [1] 0
# computing the simulation on two intersecting grids
err = sims$compute(db = grid1, namconv = NamingConvention("S1"))
err = sims$compute(db = grid2, namconv = NamingConvention("S1"))

p = ggDefaultGeographic() +
  plot.grid(grid1, nameRaster = "S1", palette = "Spectral") +
  plot.grid(grid2, nameRaster = "S1", palette = "Spectral") + 
  plot.decoration(xlab = "Easting", ylab = "Northing", title = title)
## Scale for colour_new is already present.
## Adding another scale for colour_new, which will replace the existing scale.
ggPrint(p)