Shift Operator¶

This tutorial aims to show how to handle the Shift Operator as used in SPDE algorithms.

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import matplotlib.pyplot as plt

We demonstrate the use of the Shift Operator while performing a Kriging operation using SPDE. The output data base is a regular 2-D grid (with square unit mesh of 1 unit) while the Input data base is reduced to a single point.

Using a toy example¶

In [2]:
ncell = 7
center = int(ncell / 2)
db_out_small = gl.DbGrid.create([ncell, ncell])
db_out_small.display()

Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 3
Total number of samples      = 49

Grid characteristics:
---------------------
Origin :       0.000      0.000
Mesh   :       1.000      1.000
Number :           7          7

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
In [3]:
db_in_small = gl.Db.createFromOnePoint([center, center])
db_in_small.addColumns([1], "z", gl.ELoc.Z)
db_in_small.display()

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      = 1

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x-1 - Locator = x1
Column = 2 - Name = x-2 - Locator = x2
Column = 3 - Name = z - Locator = z1

Constructing a Model composed of a single Matern Covariance, with range = 1 and param=0.5

In [4]:
range = 1
param = 0.5
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
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
---------------
Matern (Third Parameter = 0.5)
- Sill         =       1.000
- Range        =       1.000
- Theo. Range  =       0.408
Total Sill     =       1.000
Known Mean(s)      0.000

We need to create a Mesh (using the Turbo algorithm based on the output grid). Then we should create an array of meshes whose dimension must match the number of covariances in the Model (that we will keep equal to 1 in this tutorial).

In [5]:
mesh_small = gl.MeshETurbo.createFromGrid(db_out_small)
mesh_small.display()
meshes_small = gl.VectorMeshes([mesh_small])

Turbo Meshing
=============

Grid characteristics:
---------------------
Origin :       0.000      0.000
Mesh   :       1.000      1.000
Number :           7          7
Euclidean Geometry
Space Dimension           = 2
Number of Apices per Mesh = 3
Number of Meshes          = 72
Number of Apices          = 49

Bounding Box Extension
----------------------
Dim #1 - Min:0 - Max:6
Dim #2 - Min:0 - Max:6

Now let us concentrate on the Shift Operator. It may be performed:

  • using the Standard algorithm
  • or using the Stencil algorithm (all similar tiles)

This information is provided in the SPDEParam field

In [6]:
params = gl.SPDEParam()

We perform a first Kriging operation not using the Stencil lagorithm

In [7]:
params.setUseStencil(False)
err = gl.krigingSPDE(
    db_in_small,
    db_out_small,
    model,
    True,
    False,
    False,
    meshes_small,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-No-Stencil"),
)
result = gl.MatrixDense(ncell, ncell)
result.setValues(db_out_small.getColumn("K-No-Stencil.z.estim"), False)
result.display()

- Number of rows    = 7
- Number of columns = 7
                 [,  0]     [,  1]     [,  2]     [,  3]     [,  4]     [,  5]     [,  6]
      [  0,]      0.001      0.002      0.005      0.010      0.005      0.002      0.001
      [  1,]      0.002      0.004      0.014      0.035      0.014      0.004      0.002
      [  2,]      0.005      0.014      0.064      0.208      0.064      0.014      0.005
      [  3,]      0.010      0.035      0.208      0.977      0.208      0.035      0.010
      [  4,]      0.005      0.014      0.064      0.208      0.064      0.014      0.005
      [  5,]      0.002      0.004      0.014      0.035      0.014      0.004      0.002
      [  6,]      0.001      0.002      0.005      0.010      0.005      0.002      0.001

We perform a similar Kriging operation, but using the Stencil option this time

In [8]:
params.setUseStencil(True)
err = gl.krigingSPDE(
    db_in_small,
    db_out_small,
    model,
    True,
    False,
    False,
    meshes_small,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-With-Stencil"),
)
result = gl.MatrixDense(ncell, ncell)
result.setValues(db_out_small.getColumn("K-With-Stencil.z.estim"), False)
result.display()

- Number of rows    = 7
- Number of columns = 7
                 [,  0]     [,  1]     [,  2]     [,  3]     [,  4]     [,  5]     [,  6]
      [  0,]      0.000      0.000      0.000      0.000      0.000      0.000      0.000
      [  1,]      0.000      0.003      0.013      0.033      0.013      0.003      0.000
      [  2,]      0.000      0.013      0.064      0.208      0.064      0.013      0.000
      [  3,]      0.000      0.033      0.208      0.977      0.208      0.033      0.000
      [  4,]      0.000      0.013      0.064      0.208      0.064      0.013      0.000
      [  5,]      0.000      0.003      0.013      0.033      0.013      0.003      0.000
      [  6,]      0.000      0.000      0.000      0.000      0.000      0.000      0.000

The estimation results are not very similar:

  • this is essentially linked to the fact that the Stencil version required that a grid node can only be treated if the stencil (centered on the target node) is completely included within the grid. This implies that the border of the grid cannot be processed and is set to 0
  • the non-peripheral nodes are more similar but their count is small in this examples.

Using a reasonable grid size¶

The calls for using a much larger grid as it is done in the subsequent paragraph.

In [9]:
ncell = 51
center = int(ncell / 2)
db_out = gl.DbGrid.create([ncell, ncell])

db_in = gl.Db.createFromOnePoint([center, center])
db_in.addColumns([1], "z", gl.ELoc.Z)

mesh = gl.MeshETurbo.createFromGrid(db_out)
meshes = gl.VectorMeshes([mesh])
In [10]:
params.setUseStencil(False)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-No-Stencil"),
)
In [11]:
params.setUseStencil(True)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-With-Stencil"),
)
In [12]:
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
ax[0].raster(db_out, name="K-No-Stencil.z.estim", flagLegend=True)
ax[1].raster(db_out, name="K-With-Stencil.z.estim", flagLegend=True)
plt.show()
No description has been provided for this image

Clearly, with this choice of parameters, the results of the options are very close.

Checking the impact of the range¶

We use the Stencil option (quicker calculations) to visualize the impact of the range of the model. We recall that the grid dimension is 50 x 50.

In [13]:
params.setUseStencil(True)
param = 0.5
In [14]:
range = 1
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Range-1"),
)
In [15]:
range = 5
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Range-5"),
)
In [16]:
range = 10
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Range-10"),
)
In [17]:
range = 25
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Range-25"),
)
In [18]:
fig, ax = plt.subplots(2, 2, figsize=(11, 11))
ax[0, 0].raster(db_out, name="K-Range-1.z.estim", flagLegend=True)
ax[0, 1].raster(db_out, name="K-Range-5.z.estim", flagLegend=True)
ax[1, 0].raster(db_out, name="K-Range-10.z.estim", flagLegend=True)
ax[1, 1].raster(db_out, name="K-Range-25.z.estim", flagLegend=True)
plt.show()
No description has been provided for this image

Checking the impact of the regularity parameter¶

Still using the Stencil option (quicker calculations), we check the impact of the regularity term (param). The range is left constant equal to 2.

In [19]:
params.setUseStencil(True)
range = 10
In [20]:
param = 0.5
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Param-0.5"),
)
In [21]:
param = 2
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Param-2"),
)
In [22]:
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
ax[0].raster(db_out, name="K-Param-0.5.z.estim", flagLegend=True)
ax[1].raster(db_out, name="K-Param-2.z.estim", flagLegend=True)
plt.show()
No description has been provided for this image

Clearly the impact of the regularity parameter is not obvious on an estimation case: it would be much more visible in a simulation case.

Use of a Selection on the output grid¶

Let us consider that the output grid is partially masked off (using a selection for example) and compare the behavior of Kriging when use the Stencil option or not.

We first return to the small grid to visualize the impact of the selection on the values estimated at the nodes of the grid. We now create a selection in the upper left corner (smallet indices for X and Y).

In [23]:
ncell = 7
db_out_small = gl.DbGrid.create([ncell, ncell])

iuid = db_out_small.addColumnsByConstant(1, 1.0, "sel", gl.ELoc.SEL)
db_out_small.setValue("sel", db_out_small.indiceToRank([0, 0]), 0)
db_out_small.setValue("sel", db_out_small.indiceToRank([0, 1]), 0)
db_out_small.setValue("sel", db_out_small.indiceToRank([1, 0]), 0)
db_out_small.setValue("sel", db_out_small.indiceToRank([2, 0]), 0)
db_out_small.setValue("sel", db_out_small.indiceToRank([1, 1]), 0)
db_out_small.setValue("sel", db_out_small.indiceToRank([0, 2]), 0)

result = gl.MatrixDense(ncell, ncell)
result.setValues(db_out_small.getColumn("sel"), False)
result.display()

- Number of rows    = 7
- Number of columns = 7
                 [,  0]     [,  1]     [,  2]     [,  3]     [,  4]     [,  5]     [,  6]
      [  0,]      0.000      0.000      0.000      1.000      1.000      1.000      1.000
      [  1,]      0.000      0.000      1.000      1.000      1.000      1.000      1.000
      [  2,]      0.000      1.000      1.000      1.000      1.000      1.000      1.000
      [  3,]      1.000      1.000      1.000      1.000      1.000      1.000      1.000
      [  4,]      1.000      1.000      1.000      1.000      1.000      1.000      1.000
      [  5,]      1.000      1.000      1.000      1.000      1.000      1.000      1.000
      [  6,]      1.000      1.000      1.000      1.000      1.000      1.000      1.000
In [24]:
range = 1
param = 0.5
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
In [25]:
params.setUseStencil(False)
err = gl.krigingSPDE(
    db_in_small,
    db_out_small,
    model,
    True,
    False,
    False,
    meshes_small,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-No-Stencil"),
)
result = gl.MatrixDense(ncell, ncell)
result.setValues(db_out_small.getColumn("K-No-Stencil.z.estim"), False)
result.display()

- Number of rows    = 7
- Number of columns = 7
                 [,  0]     [,  1]     [,  2]     [,  3]     [,  4]     [,  5]     [,  6]
      [  0,]      0.000      0.000      0.000      0.010      0.005      0.002      0.001
      [  1,]      0.000      0.000      0.014      0.035      0.014      0.004      0.002
      [  2,]      0.000      0.014      0.064      0.208      0.064      0.014      0.005
      [  3,]      0.010      0.035      0.208      0.977      0.208      0.035      0.010
      [  4,]      0.005      0.014      0.064      0.208      0.064      0.014      0.005
      [  5,]      0.002      0.004      0.014      0.035      0.014      0.004      0.002
      [  6,]      0.001      0.002      0.005      0.010      0.005      0.002      0.001
In [26]:
params.setUseStencil(True)
err = gl.krigingSPDE(
    db_in_small,
    db_out_small,
    model,
    True,
    False,
    False,
    meshes_small,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-With-Stencil"),
)
result = gl.MatrixDense(ncell, ncell)
result.setValues(db_out_small.getColumn("K-With-Stencil.z.estim"), False)
result.display()

- Number of rows    = 7
- Number of columns = 7
                 [,  0]     [,  1]     [,  2]     [,  3]     [,  4]     [,  5]     [,  6]
      [  0,]      0.000      0.000      0.000      0.000      0.000      0.000      0.000
      [  1,]      0.000      0.000      0.013      0.033      0.013      0.003      0.000
      [  2,]      0.000      0.013      0.064      0.208      0.064      0.013      0.000
      [  3,]      0.000      0.033      0.208      0.977      0.208      0.033      0.000
      [  4,]      0.000      0.013      0.064      0.208      0.064      0.013      0.000
      [  5,]      0.000      0.003      0.013      0.033      0.013      0.003      0.000
      [  6,]      0.000      0.000      0.000      0.000      0.000      0.000      0.000

Once more, we can check the quality of the Stencil option:

  • it produces a result very close to the Standard option (although in a much faster way)
  • the guard zone is maintained: e.g. the masked area and the edge of the field

Using the large grid¶

We return to more decent grid size in order to visualize the impact of the selection. We add a selection based on a single cell.

In [27]:
iuid = db_out.addColumnsByConstant(1, 1.0, "sel", gl.ELoc.SEL)
db_out.setValue("sel", db_out.indiceToRank([20, 15]), 0)

We return to a long range model (range = 10).

In [28]:
range = 10
param = 0.5
model = gl.Model.createFromParam(gl.ECov.MATERN, range, 1.0, param)
In [29]:
params.setUseStencil(False)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Sel-No-Stencil"),
)
In [30]:
params.setUseStencil(True)
err = gl.krigingSPDE(
    db_in,
    db_out,
    model,
    True,
    False,
    False,
    meshes,
    None,
    None,
    None,
    None,
    None,
    params,
    False,
    namconv=gl.NamingConvention("K-Sel-With-Stencil"),
)
In [31]:
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
ax[0].raster(db_out, name="K-Sel-No-Stencil.z.estim", flagLegend=True)
ax[1].raster(db_out, name="K-Sel-With-Stencil.z.estim", flagLegend=True)
plt.show()
No description has been provided for this image

Finally, in order to demonstrate that the Stencil Option does not impair the quality of the result, we display a scatter plot between the two results.

In [32]:
fig, ax = plt.subplots()
ax.correlation(
    db_out, "K-Sel-No-Stencil.z.estim", "K-Sel-With-Stencil.z.estim", bins=100, cmin=1
)
ax.decoration(title="Scatter Plot between the two options")
plt.show()
No description has been provided for this image