SPDE for Spiral Anisotropy¶
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)
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
gdoc.setNoScroll()
Defining the Model as a single Matern 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 spirale.
model = gl.Model.createFromParam(gl.ECov.MATERN, 1.0, 1.0, 1.0, [4.0, 45.0])
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 = gl.FunctionalSpirale(0.0, -1.4, 1.0, 1.0, 50.0, 50.0)
cova = model.getCovAniso(0)
cova.makeAngleNoStatFunctional(spirale)
A visualisation of the non-stationarity can be otanined in the following paragraph. The angle is represented at each node of a grid. For better legibility the grid is defined as a coarse grid.
coarse = gl.DbGrid.create([26, 26], [4.0, 4.0])
res = gp.covaOnGrid(cova, coarse, scale=2000)
Creating a output grid
grid = gl.DbGrid.create([101, 101], [1.0, 1.0])
Perform several non-conditional simulations on the grid, using the Model and the non-stationarity.
nbsimu = 4
iuid = gl.simulateSPDE(None, 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.S1 - Locator = z1 Column = 4 - Name = SimuSPDE.S2 - Locator = z2 Column = 5 - Name = SimuSPDE.S3 - Locator = z3 Column = 6 - Name = SimuSPDE.S4 - Locator = z4
We represent the non-conditional simulations
vmin = -4
vmax = +4
fig, ax = gp.init(2, 2, figsize=(16, 12))
ax[0, 0].raster(grid, name="SimuSPDE.S1", useSel=False, vmin=vmin, vmax=vmax)
ax[0, 1].raster(grid, name="SimuSPDE.S2", useSel=False, vmin=vmin, vmax=vmax)
ax[1, 0].raster(grid, name="SimuSPDE.S3", useSel=False, vmin=vmin, vmax=vmax)
ax[1, 1].raster(grid, name="SimuSPDE.S4", useSel=False, vmin=vmin, vmax=vmax)
fig.subplots_adjust(right=0.7)
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 = gl.Db.createSamplingDb(grid, number=100, names=["x1", "x2", "SimuSPDE.S1"])
data.setName("SimuSPDE.S1", "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
res = gp.plot(data, nameColor="data")
Use the previous data set (and the non-stationary Model) in order to perform an estimation
iuid = gl.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.S1 - Locator = NA Column = 4 - Name = SimuSPDE.S2 - Locator = NA Column = 5 - Name = SimuSPDE.S3 - Locator = NA Column = 6 - Name = SimuSPDE.S4 - Locator = NA Column = 7 - Name = KrigingSPDE.data.estim - Locator = z1
Representing the Estimation obtained on the Grid
res = gp.plot(grid, "KrigingSPDE.data.estim")
Performing several conditional simulation
nbsimu = 4
iuid = gl.simulateSPDE(
data, grid, model, nbsimu, namconv=gl.NamingConvention("CondSimu")
)
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.S1 - Locator = NA Column = 4 - Name = SimuSPDE.S2 - Locator = NA Column = 5 - Name = SimuSPDE.S3 - Locator = NA Column = 6 - Name = SimuSPDE.S4 - Locator = NA Column = 7 - Name = KrigingSPDE.data.estim - Locator = NA Column = 8 - Name = CondSimu.data.S1 - Locator = z1 Column = 9 - Name = CondSimu.data.S2 - Locator = z2 Column = 10 - Name = CondSimu.data.S3 - Locator = z3 Column = 11 - Name = CondSimu.data.S4 - Locator = z4
Representing the conditional simulations
vmin = -4
vmax = +4
fig, ax = gp.init(2, 2, figsize=(16, 12))
ax[0, 0].raster(grid, name="CondSimu.*.S1", useSel=False, vmin=vmin, vmax=vmax)
ax[0, 1].raster(grid, name="CondSimu.*.S2", useSel=False, vmin=vmin, vmax=vmax)
ax[1, 0].raster(grid, name="CondSimu.*.S3", useSel=False, vmin=vmin, vmax=vmax)
ax[1, 1].raster(grid, name="CondSimu.*.S4", useSel=False, vmin=vmin, vmax=vmax)
fig.subplots_adjust(right=0.7)
