In this tutorial, we show how the use of SPDE for Varying Anisotropy when this Anisotropy must follow a Spiral shape (defined as an external function)

Defining the Model as a single Bessel structure. This function is defined as anisotropic: we clearly specify the extension of the ranges in the two main directions. The angle does not have to be defined here: it will be overwritten later as the non-stationary parameter. Note that it is essential to define the short range of the anisotropy ellipsoid first (for the definition of angle as defined in the Spiral function used as a function)… otherwise future results will represent the shape otabined as the orthogonal of the spiral.

model = Model_createFromParam(ECov_BESSEL_K(), 1., 1., 1., c(4.,45.))

A Spiral function is defined and attached to the Model: this is a manner to update the Model by transforming the anisotropy angle as the unique non-stationary parameter.

spirale = FunctionalSpirale(0., -1.4, 1., 1., 50., 50.)
nostat = NoStatFunctional(spirale)
err = model$addNoStat(nostat)
model$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
## ---------------
## K-Bessel (Third Parameter = 1)
## - Sill         =      1.000
## - Ranges       =      4.000    45.000
## - Theo. Ranges =      1.155    12.990
## Total Sill     =      1.000
## 
## Non-Stationary Parameters
## -------------------------
## Angle      : GRF=1 Str=1 V#1=1
## Functional
## NULL

Creating a output grid

grid = DbGrid_create(c(101,101), c(1.,1.))

Perform several non-conditional simulations on the grid, using the Model and the non-stationarity.

nbsimu = 4
iuid = simulateSPDE(NULL,grid,model,nbsimu)
grid
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 7
## Total number of samples      = 10201
## 
## Grid characteristics:
## ---------------------
## Origin :      0.000     0.000
## Mesh   :      1.000     1.000
## Number :        101       101
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = SimuSPDE.1 - Locator = z1
## Column = 4 - Name = SimuSPDE.2 - Locator = z2
## Column = 5 - Name = SimuSPDE.3 - Locator = z3
## Column = 6 - Name = SimuSPDE.4 - Locator = z4

We represent the non-conditional simulations

p1 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.1")
p2 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.2")
p3 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.3")
p4 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.4")
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)

Extracting a set of nodes randomly located in order to create a data file which will serve as conditioning. The data is extracted from the first non-conditional simulation.

data = Db_createSamplingDb(grid, number=100, names=c("x1", "x2", "SimuSPDE.1"))
err = data$setName("*SimuSPDE.1", "data")
data
## 
## Data Base Characteristics
## =========================
## 
## Data Base Summary
## -----------------
## File is organized as a set of isolated points
## Space dimension              = 2
## Number of Columns            = 4
## Total number of samples      = 100
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = data - Locator = z1
ggDefaultGeographic() + plot.point(data, nameColor="data")
## Warning: Using size for a discrete variable is not advised.

Use the previous data set (and the non-stationary Model) in order to perform an estimation

iuid = krigingSPDE(data,grid,model)
grid
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 8
## Total number of samples      = 10201
## 
## Grid characteristics:
## ---------------------
## Origin :      0.000     0.000
## Mesh   :      1.000     1.000
## Number :        101       101
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = SimuSPDE.1 - Locator = NA
## Column = 4 - Name = SimuSPDE.2 - Locator = NA
## Column = 5 - Name = SimuSPDE.3 - Locator = NA
## Column = 6 - Name = SimuSPDE.4 - Locator = NA
## Column = 7 - Name = KrigingSPDE.data.estim - Locator = z1

Representing the Estimation obtained on the Grid

ggDefaultGeographic() + plot.grid(grid, nameRaster="KrigingSPDE.data.estim")

Performing several conditional simulation

nbsimu = 4
iuid = simulateSPDE(data,grid,model,nbsimu)
grid
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of Columns            = 12
## Total number of samples      = 10201
## 
## Grid characteristics:
## ---------------------
## Origin :      0.000     0.000
## Mesh   :      1.000     1.000
## Number :        101       101
## 
## Variables
## ---------
## Column = 0 - Name = rank - Locator = NA
## Column = 1 - Name = x1 - Locator = x1
## Column = 2 - Name = x2 - Locator = x2
## Column = 3 - Name = SimuSPDE.1 - Locator = NA
## Column = 4 - Name = SimuSPDE.2 - Locator = NA
## Column = 5 - Name = SimuSPDE.3 - Locator = NA
## Column = 6 - Name = SimuSPDE.4 - Locator = NA
## Column = 7 - Name = KrigingSPDE.data.estim - Locator = NA
## Column = 8 - Name = SimuSPDE.data.1 - Locator = z1
## Column = 9 - Name = SimuSPDE.data.2 - Locator = z2
## Column = 10 - Name = SimuSPDE.data.3 - Locator = z3
## Column = 11 - Name = SimuSPDE.data.4 - Locator = z4

Representing the conditional simulations

p1 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.data.1")
p2 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.data.2")
p3 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.data.3")
p4 = ggDefaultGeographic() + plot.grid(grid, nameRaster="SimuSPDE.data.4")
ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2)