Gaussian simulations

Preamble

In this preamble, we load the gstlearn library and clean the workspace.

rm(list=ls())
library(gstlearn)
library(ggplot2)
library(ggpubr)
library(ggnewscale)

## Global plot option
plot.setDefaultGeographic(dims=c(8,8))

We load two data bases:

• a data base dat containing point observations of two variables across Scotland: the elevation (Elevation) and the temperature (January_temp)
• a data base target containing a grid of points covering Scotland with a selection variable (inshore) selecting the points that are on land, and a variable (Elevation) giving the elevation at every point on land

## Data points
fileNF = loadData("Scotland", "Scotland_Temperatures.NF")
dat = Db_createFromNF(fileNF)

## Target grid
fileNF = loadData("Scotland", "Scotland_Elevations.NF")
target = DbGrid_createFromNF(fileNF)

We also compute an experimental variogram on the observations and fit a model on it.

## Define and compute experimental variogram
varioparam = VarioParam_createOmniDirection(npas=40, dpas=10)
vario_raw2dir = Vario_create(varioparam)
err = vario_raw2dir$compute(dat)

## Fit model
fitmod = Model()
err = fitmod$fit(vario_raw2dir, 
                 types=ECov_fromKeys(c("NUGGET", "SPHERICAL", "CUBIC")))
fitmod$display()

## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 1
## Number of drift function(s)  = 0
## Number of drift equation(s)  = 0
## 
## Covariance Part
## ---------------
## Spherical
## - Sill         =      1.155
## - Range        =    135.133
## Total Sill     =      1.155
## Known Mean(s)     0.000

## NULL

neighU = NeighUnique_create()

ndim = 2
defineDefaultSpace(ESpaceType_RN(), ndim)

## NULL

---

Unconditional simulation

To generate unconditional simulations, we use the simtub function. This function generates samples from a Gaussian random field with a covariance model defined in a Model object, using the turning bands algorithm. We specify

• the data base containing the target points on which we want to simulate the model (argument dbout)
• the Model object defining the model we want to simulate (argument model)
• the number of samples we want to generate (argument nbsimu)
• the number of turning bands (argument nbtuba)

Optionally, we can specify a seed number for the simulation (to ensure reproducibility). The simtub function adds the simulated samples directly to the target data base specified in dbout (with a naming convention that can be set through the argument namconv). Note that the samples generated by this function have the same mean as the one specified in the model object. If this mean has not specified been specified (through the setMeans method), then zero-mean simulations are generated.

Let us generate a sample from the model fitmod we fitted on the observations. First, we simulate the model with a single turning band.

err = simtub(dbout=target, model=fitmod, 
             nbsimu=1,
             nbtuba=1, seed=12454,
             namconv=NamingConvention("Simu1"))
target$display()

## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 5
## Total number of samples      = 11097
## Number of active samples     = 3092
## 
## Grid characteristics:
## ---------------------
## Origin :     65.000   535.000
## Mesh   :      4.938     4.963
## Number :         81       137
## 
## Variables
## ---------
## Column = 0 - Name = Longitude - Locator = x1
## Column = 1 - Name = Latitude - Locator = x2
## Column = 2 - Name = Elevation - Locator = f1
## Column = 3 - Name = inshore - Locator = sel
## Column = 4 - Name = Simu1 - Locator = z1

## NULL

p = ggDefaultGeographic()
p = p + plot.grid(target, nameRaster = "Simu1",
            flagLegendRaster=TRUE,palette="Spectral",
            legendNameRaster="Value")
p = p + plot.decoration(title="Simulation with 1 band")
ggPrint(p)

Let us now simulate the model using 10 turning bands.

err = simtub(dbout=target, model=fitmod, 
             nbsimu=1,
             nbtuba=10, seed=12454,
             namconv=NamingConvention("Simu10"))

p = ggDefaultGeographic()
p = p + plot.grid(target, nameRaster = "Simu10",
            flagLegendRaster=TRUE,palette="Spectral",legendNameRaster="Value")
p = p + plot.decoration(title="Simulation with 1 band")
ggPrint(p)