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