Spectral simulation of random functions
gstlearn Core Team
2024-02-13
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
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)| 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)
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 = sims$compute(db = db, 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)| 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)
# Initialization of the simulator
sims = SimuSpectral(mod)
sims$simulate(nb = nb, seed = seed)## [1] 0# computing the simulation on the grid
err = sims$compute(dbout = grid, namconv = NamingConvention("S1"))# statistics
knitr::kable(dbStatisticsMono(db = grid, names = "S1", opers = opers)$toTL(),
digits = 3, caption = title)| 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)