SPDE for Spiral Anisotropy¶

In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

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)

In [2]:
import gstlearn as gl
import gstlearn.plot as gp
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import random
import math

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.

In [3]:
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.

In [4]:
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

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.

In [5]:
coarse = gl.DbGrid.create([26,26],[4.,4.])
ax = gp.modelOnGrid(model, coarse, scale=2000)

Creating a output grid

In [6]:
grid = gl.DbGrid.create([101,101],[1.,1.])

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

In [7]:
nbsimu = 4
iuid = gl.simulateSPDE(None,grid,model,nbsimu)
grid
Out[7]:
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

In [8]:
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.

In [9]:
data = gl.Db.createSamplingDb(grid, number=100, names=["x1", "x2", "SimuSPDE.1"])
data.setName("SimuSPDE.1", "data")
data
Out[9]:
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
In [10]:
gp.plot(data, nameColor="data")

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

In [11]:
iuid = gl.krigingSPDE(data,grid,model)
grid
Out[11]:
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

In [12]:
gp.plot(grid, "KrigingSPDE.data.estim")

Performing several conditional simulation

In [13]:
nbsimu = 4
iuid = gl.simulateSPDE(data,grid,model,nbsimu, namconv=gl.NamingConvention("CondSimu"))
grid
Out[13]:
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

In [14]:
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)