Second order stationary random functions

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

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

# 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

# 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

err = simuSpectral(NULL, grid, mod, ns, nb, seed=seed , namconv = NamingConvention("S"))

# 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))

# computing the simulation on two intersecting grids
seed <- 3132
err = simuSpectral(NULL, grid1, mod, 1, nb, seed=seed , namconv = NamingConvention("S1"))
err = simuSpectral(NULL, grid2, mod, 1, nb, seed=seed , 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)
ggPrint(p)

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

# computing the simulation on scattered points (extracted from grid1)
grid1["sel"] <- VectorHelper_simulateBernoulli(grid1$getSampleNumber(), 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["grid1.S1"] <- grid1$getColumn(name = "S1", useSel = TRUE)
err = db$setLocators(c("x1", "x2"), ELoc_X(), cleanSameLocator = TRUE)
err = grid1$clearSelection()


# computing the simulation on scattered points
err = simuSpectral(NULL, db, mod, 1, nb, seed=seed , namconv = NamingConvention("S1"))


p1 = ggDefaultGeographic() +
  plot.grid(grid1, nameRaster = "S1", palette = "Spectral") +
  plot.point(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 = ggDefault() +
  plot.correlation(db, namex = "S1", namey = "grid1.S1", 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
knitr::kable(rbind(
  dbStatisticsMono(db = grid1, names = "S1", opers = opers)$toTL(),
  dbStatisticsMono(db = grid2, names = c("S1"), opers = opers, 
                   flagIso = TRUE)$toTL(),
  dbStatisticsMono(db = db, names = c("S1"), opers = opers, 
                   flagIso = TRUE)$toTL()),
  digits = 3, caption = title)

EXPONENTIAL model
	Number	Minimum	Maximum	Mean	St. Dev.
S1	65536	-3.999	3.579	-0.114	0.964
S11	10000	-3.066	4.201	0.295	0.988
S12	357	-2.565	2.434	-0.196	0.912

# variogram on the simulations
nlag = 25
coldir = c("green", "blue", "red")
err = grid1$setLocator("S1", ELoc_Z(), cleanSameLocator = TRUE)
varioPar = VarioParam_createMultipleFromGrid(grid1, npas=nlag)

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

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

3D simulation of a random function

err = defineDefaultSpace(ESpaceType_RN(), 3)

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

# statistics
knitr::kable(dbStatisticsMono(db = grid, names = "S1", opers = opers)$toTL(),
             digits = 3, caption = title)

EXPONENTIAL model
	Number	Minimum	Maximum	Mean	St. Dev.
S1	2e+05	-3.525	3.987	0.018	1.03

# histogram
p = ggDefault() + 
  plot.hist(db = grid, name = "S1", bins = 25, col = "skyblue")
ggPrint(p)

# map
p = ggDefaultGeographic() + 
  plot.grid(grid, nameRaster = "S1", palette = "Spectral") +
  plot.decoration(xlab = "Easting", ylab = "Northing", title = title)
ggPrint(p)

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

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

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

Spectral simulation of random functions

gstlearn Core Team

2024-02-13

Second order stationary random functions

1D simulation of a random function

2D simulation of a random function

3D simulation of a random function