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
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import random
import math
gdoc.setNoScroll()
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 spirale.
model = gl.Model.createFromParam(gl.ECov.BESSEL_K, 1., 1., 1., [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 = gl.FunctionalSpirale(0., -1.4, 1., 1., 50., 50.)
nostat = gl.NoStatFunctional(spirale)
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 Known Mean(s) 0.000
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.,4.])
ax = gp.modelOnGrid(model, coarse, scale=2000)
Creating a output grid
grid = gl.DbGrid.create([101,101],[1.,1.])
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.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
fig = plt.figure(figsize=(16,12))
vmin = -4
vmax = +4
ax1 = fig.add_subplot(2,2,1)
ax1.raster(grid,name="SimuSPDE.1", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
ax2 = fig.add_subplot(2,2,2)
ax2.raster(grid,name="SimuSPDE.2", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
ax3 = fig.add_subplot(2,2,3)
ax3.raster(grid,name="SimuSPDE.3", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
ax4 = fig.add_subplot(2,2,4)
ax4.raster(grid,name="SimuSPDE.4", useSel=False, flagLegend = 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.1"])
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
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.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
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.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 = CondSimu.data.1 - Locator = z1 Column = 9 - Name = CondSimu.data.2 - Locator = z2 Column = 10 - Name = CondSimu.data.3 - Locator = z3 Column = 11 - Name = CondSimu.data.4 - Locator = z4
Representing the conditional simulations
fig = plt.figure(figsize=(16,12))
vmin = -4
vmax = +4
ax1 = fig.add_subplot(2,2,1)
ax1.raster(grid,name="CondSimu.*.1", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
ax2 = fig.add_subplot(2,2,2)
ax2.raster(grid,name="CondSimu.*.2", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
ax3 = fig.add_subplot(2,2,3)
ax3.raster(grid,name="CondSimu.*.3", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
ax4 = fig.add_subplot(2,2,4)
ax4.raster(grid,name="CondSimu.*.4", useSel=False, flagLegend = False, vmin=vmin, vmax=vmax)
fig.subplots_adjust(right=0.7)
