Point & Block Kriging¶

This chapter describes the difference between Point Kriging, Regular Block Kriging and Irregular Block Kriging. It is delonstrated on a simulated 2-D data set.

Import packages¶

In [1]:
import numpy as np
import pandas as pd
import sys
import os
import matplotlib.pyplot as plt
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc

gdoc.setNoScroll()

Creating the environmenet¶

We create a regular grid

In [2]:
grid = gl.DbGrid.create(nx=[20,15],dx=[10.,10.])
ncell = grid.getSampleNumber()
grid
Out[2]:
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      = 300

Grid characteristics:
---------------------
Origin :      0.000     0.000
Mesh   :     10.000    10.000
Number :         20        15

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2

We create a Data Base of randomly located samples within a field covering the grid extension.

In [3]:
coormin = grid.getCoorMinimum()
coormax = grid.getCoorMaximum()
nech = 100
data = gl.Db.createFromBox(nech, coormin, coormax)

We create a Model composed of a Nugget Effect and a Spherical components. This Model is used to simulate a variable at the sample locations.

In [4]:
model = gl.Model()
model.addCovFromParam(gl.ECov.NUGGET, sill = 0.3)
model.addCovFromParam(gl.ECov.SPHERICAL, range=40., sill=0.7)
err = gl.simtub(None, dbout=data, model=model)
data
Out[4]:
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 = x-1 - Locator = x1
Column = 2 - Name = x-2 - Locator = x2
Column = 3 - Name = Simu - Locator = z1

We define a Neighborhood. On purpose, this neighborhood is small (10 data maximum) to produce legible printouts.

In [5]:
nmini = 1
nmaxi = 10
radius = 30.
nsect = 8
nsmax = 3
neigh = gl.NeighMoving.create(nmaxi=nmaxi, radius=radius, nmini=nmini, nsect=nsect, nsmax=nsmax)
neigh
Out[5]:
Moving Neighborhood
===================
Minimum number of samples           = 1
Maximum number of samples           = 10
Number of angular sectors           = 8
Maximum number of points per sector = 3
Maximum horizontal distance         = 30

Point Kriging¶

In [6]:
err = gl.kriging(data, grid, model, neigh, 
                 namconv=gl.NamingConvention("Point_Kriging"))

Two newly created variables have been added to the output Data Base (grid):

  • one for the estimation (called Point_Kriging.Simu.estim)
  • one for the standard deviation of the estimation error (called Point_Kriging.simu.stdev)

These two variables are stored as a consequance of the argument flag_est and flag_std which are not explicitely specified, but which are both defaulted to True.

In [7]:
grid
Out[7]:
Data Base Grid Characteristics
==============================

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

Grid characteristics:
---------------------
Origin :      0.000     0.000
Mesh   :     10.000    10.000
Number :         20        15

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Point_Kriging.Simu.estim - Locator = z1
Column = 4 - Name = Point_Kriging.Simu.stdev - Locator = NA

Displaying the maps of the Point Estimation and the corresponding Standard Deviation of Estimation errors. Note that the estimation is performed at the nodes of the grid (center of each grid cell)... evan if the block raster is performed.

In [8]:
ax = grid.plot("Point_Kriging*estim")
ax = data.plot(color="white")
ax.decoration(title="Point Kriging")
In [9]:
ax = grid.plot("Point_Kriging*stdev")
ax = data.plot(color="white")
ax.decoration(title="Error for Point Kriging")

To better visualize the differences alongst the different kriging, we focus on a target grid cell (or node) and print out all the information used during the Kriging procedure. We focus on the node

In [10]:
node = 155
err = gl.krigtest(data, grid, model, neigh, node, calcul=gl.EKrigOpt.POINT)
 Target location
 ---------------
 Sample #156 (from 300)
 Coordinate #1 = 150.000000
 Coordinate #2 = 70.000000
 
 Data selected in neighborhood
 -----------------------------
        Rank      Sample          x1          x2      Sector 
           1           2     156.414      61.537           3 
           2           7     140.011      52.437           2 
           3          34     169.243      74.058           5 
           4          36     167.433      74.021           5 
           5          39     159.821      45.060           3 
           6          40     151.089      69.507           4 
           7          49     165.529      66.754           4 
           8          55     127.590      71.435           8 
           9          73     151.278      84.171           6 
          10          88     140.091      58.596           2 
 
 LHS of Kriging matrix (compressed)
 ==================================
 Number of active samples    = 10
 Total number of equations   = 10
 Reduced number of equations = 10
 
        Rank                       1           2           3           4           5 
                    Flag           1           2           3           4           5 
           1           1       1.000       0.244       0.261       0.288       0.284 
           2           2       0.244       1.000       0.008       0.016       0.197 
           3           3       0.261       0.008       1.000       0.653       0.055 
           4           4       0.288       0.016       0.653       1.000       0.061 
           5           5       0.284       0.197       0.055       0.061       1.000 
           6           6       0.453       0.212       0.245       0.282       0.114 
           7           7       0.431       0.069       0.488       0.505       0.173 
           8           8       0.055       0.168       0.000       0.000       0.000 
           9           9       0.159       0.025       0.207       0.237       0.000 
          10          10       0.290       0.540       0.030       0.045       0.147 
 
        Rank                       6           7           8           9          10 
                    Flag           6           7           8           9          10 
           1           1       0.453       0.431       0.055       0.159       0.290 
           2           2       0.212       0.069       0.168       0.025       0.540 
           3           3       0.245       0.488       0.000       0.207       0.030 
           4           4       0.282       0.505       0.000       0.237       0.045 
           5           5       0.114       0.173       0.000       0.000       0.147 
           6           6       1.000       0.331       0.153       0.332       0.314 
           7           7       0.331       1.000       0.002       0.172       0.103 
           8           8       0.153       0.002       1.000       0.100       0.261 
           9           9       0.332       0.172       0.100       1.000       0.086 
          10          10       0.314       0.103       0.261       0.086       1.000 
 
 RHS of Kriging matrix (compressed)
 ==================================
 Number of active samples    = 10
 Total number of equations   = 10
 Reduced number of equations = 10
 Number of right-hand sides  = 1
 Punctual Estimation
 
        Rank        Flag           1 
           1           1       0.428 
           2           2       0.215 
           3           3       0.225 
           4           4       0.262 
           5           5       0.102 
           6           6       0.669 
           7           7       0.305 
           8           8       0.172 
           9           9       0.342 
          10          10       0.322 
 
 (Co-) Kriging weights
 =====================
        Rank          x1          x2        Data         Z1* 
           1     156.414      61.537       0.946       0.122 
           2     140.011      52.437      -0.605       0.016 
           3     169.243      74.058       0.801       0.004 
           4     167.433      74.021       0.447       0.027 
           5     159.821      45.060       0.932      -0.017 
           6     151.089      69.507      -0.248       0.513 
           7     165.529      66.754       0.166       0.038 
           8     127.590      71.435       1.070       0.048 
           9     151.278      84.171      -0.040       0.126 
          10     140.091      58.596      -1.474       0.091 
 Sum of weights                                     0.968 
 
 (Co-) Kriging results
 =====================
 Target Sample = 156
 Variable Z1 
  - Estimate  =       -0.103 
  - Std. Dev. =        0.709 
  - Variance  =        0.503 
  - Cov(h=0)  =        1.000 
 

Standard Block Kriging¶

We perform an estimation of the average value of the variable over the cells (of the previous grid).

In [11]:
err = gl.kriging(data, grid, model, neigh, calcul=gl.EKrigOpt.BLOCK, ndisc=[5,5],
                namconv=gl.NamingConvention("Block_Kriging"))

The output Db now contains the two additional variables Block_Kriging.Simu.estim and Block_Kriging.Simu.stdev.

In [12]:
grid
Out[12]:
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      = 300

Grid characteristics:
---------------------
Origin :      0.000     0.000
Mesh   :     10.000    10.000
Number :         20        15

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Point_Kriging.Simu.estim - Locator = NA
Column = 4 - Name = Point_Kriging.Simu.stdev - Locator = NA
Column = 5 - Name = Block_Kriging.Simu.estim - Locator = z1
Column = 6 - Name = Block_Kriging.Simu.stdev - Locator = NA

Displaying the maps of the Point Estimation and the corresponding Standard Deviation of Estimation errors

In [13]:
ax = grid.plot("Block_Kriging*estim")
ax = data.plot(color="white")
ax.decoration(title="Block Kriging")
In [14]:
ax = grid.plot("Block_Kriging*stdev")
ax = data.plot(color="white")
ax.decoration(title="Error for Block Kriging")

We do not see too much difference between estimations, due to the small dimension of the grid cells, as it is demonstrated in the following scatter plot.

In [15]:
ax = gp.correlation(grid, namex="Point_Kriging*estim", namey="Block_Kriging*estim", bins=100)
ax.decoration(title="Estimation", xlabel="Point", ylabel="Block")
In [16]:
ax = gp.correlation(grid, namex="Point_Kriging*stdev", namey="Block_Kriging*stdev", bins=100)
ax.decoration(title="St. Dev.", xlabel="Point", ylabel="Block")

We check our reference target cell again. Remember that here the target is the average over the cell. In the final part of the printout, we clearly see the value of the cell extension i.e. 10 (important for future comparison).

In [17]:
err = gl.krigtest(data, grid, model, neigh, node, calcul=gl.EKrigOpt.BLOCK, ndisc=[5,5])
 Target location
 ---------------
 Sample #156 (from 300)
 Coordinate #1 = 150.000000
 Coordinate #2 = 70.000000
 
 Data selected in neighborhood
 -----------------------------
        Rank      Sample          x1          x2      Sector 
           1           2     156.414      61.537           3 
           2           7     140.011      52.437           2 
           3          34     169.243      74.058           5 
           4          36     167.433      74.021           5 
           5          39     159.821      45.060           3 
           6          40     151.089      69.507           4 
           7          49     165.529      66.754           4 
           8          55     127.590      71.435           8 
           9          73     151.278      84.171           6 
          10          88     140.091      58.596           2 
 
 LHS of Kriging matrix (compressed)
 ==================================
 Number of active samples    = 10
 Total number of equations   = 10
 Reduced number of equations = 10
 
        Rank                       1           2           3           4           5 
                    Flag           1           2           3           4           5 
           1           1       1.000       0.244       0.261       0.288       0.284 
           2           2       0.244       1.000       0.008       0.016       0.197 
           3           3       0.261       0.008       1.000       0.653       0.055 
           4           4       0.288       0.016       0.653       1.000       0.061 
           5           5       0.284       0.197       0.055       0.061       1.000 
           6           6       0.453       0.212       0.245       0.282       0.114 
           7           7       0.431       0.069       0.488       0.505       0.173 
           8           8       0.055       0.168       0.000       0.000       0.000 
           9           9       0.159       0.025       0.207       0.237       0.000 
          10          10       0.290       0.540       0.030       0.045       0.147 
 
        Rank                       6           7           8           9          10 
                    Flag           6           7           8           9          10 
           1           1       0.453       0.431       0.055       0.159       0.290 
           2           2       0.212       0.069       0.168       0.025       0.540 
           3           3       0.245       0.488       0.000       0.207       0.030 
           4           4       0.282       0.505       0.000       0.237       0.045 
           5           5       0.114       0.173       0.000       0.000       0.147 
           6           6       1.000       0.331       0.153       0.332       0.314 
           7           7       0.331       1.000       0.002       0.172       0.103 
           8           8       0.153       0.002       1.000       0.100       0.261 
           9           9       0.332       0.172       0.100       1.000       0.086 
          10          10       0.314       0.103       0.261       0.086       1.000 
 
 RHS of Kriging matrix (compressed)
 ==================================
 Number of active samples    = 10
 Total number of equations   = 10
 Reduced number of equations = 10
 Number of right-hand sides  = 1
 Block Estimation : Discretization =  5  x  5 
 
        Rank        Flag           1 
           1           1       0.420 
           2           2       0.214 
           3           3       0.224 
           4           4       0.259 
           5           5       0.103 
           6           6       0.598 
           7           7       0.302 
           8           8       0.172 
           9           9       0.338 
          10          10       0.318 
 
 (Co-) Kriging weights
 =====================
        Rank          x1          x2       Size1       Size2        Data         Z1* 
           1     156.414      61.537      10.000      10.000       0.946       0.143 
           2     140.011      52.437      10.000      10.000      -0.605       0.021 
           3     169.243      74.058      10.000      10.000       0.801       0.007 
           4     167.433      74.021      10.000      10.000       0.447       0.034 
           5     159.821      45.060      10.000      10.000       0.932      -0.015 
           6     151.089      69.507      10.000      10.000      -0.248       0.414 
           7     165.529      66.754      10.000      10.000       0.166       0.048 
           8     127.590      71.435      10.000      10.000       1.070       0.056 
           9     151.278      84.171      10.000      10.000      -0.040       0.144 
          10     140.091      58.596      10.000      10.000      -1.474       0.104 
 Sum of weights                                                           0.956 
 
 (Co-) Kriging results
 =====================
 Target Sample = 156
 Variable Z1 
  - Estimate  =       -0.065 
  - Std. Dev. =        0.372 
  - Variance  =        0.138 
  - Cov(h=0)  =        0.565 
 

Irregular Block Kriging¶

Now we add two vectors in the Target Grid data base, which will contain the cell extension (variable per cell). Nevertheless, here, the cell extension is set to a constant value (for simplicity sake): this value (35) is different from the grid mesh (10). The variable block extension is assigned the locator BLEX.

In [18]:
size = 35.
iuid = grid.addColumnsByConstant(1, size, "X-ext", gl.ELoc.BLEX, 0)
iuid = grid.addColumnsByConstant(1, size, "Y-ext", gl.ELoc.BLEX, 1)
In [19]:
err = gl.krigcell(data, grid, model, neigh, ndisc=[5,5],
                  namconv=gl.NamingConvention("Irregular_Kriging"))

The two newly created results are added to the output data base grid with the radix *Irregular_Kriging".

In [20]:
grid
Out[20]:
Data Base Grid Characteristics
==============================

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

Grid characteristics:
---------------------
Origin :      0.000     0.000
Mesh   :     10.000    10.000
Number :         20        15

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = Point_Kriging.Simu.estim - Locator = NA
Column = 4 - Name = Point_Kriging.Simu.stdev - Locator = NA
Column = 5 - Name = Block_Kriging.Simu.estim - Locator = NA
Column = 6 - Name = Block_Kriging.Simu.stdev - Locator = NA
Column = 7 - Name = X-ext - Locator = dblk1
Column = 8 - Name = Y-ext - Locator = dblk2
Column = 9 - Name = Irregular_Kriging.Simu.estim - Locator = z1
Column = 10 - Name = Irregular_Kriging.Simu.stdev - Locator = NA

Graphic display of the new resulting maps.

In [21]:
ax = grid.plot("Irregular_Kriging*estim")
ax = data.plot(color="white")
ax.decoration(title="Block Kriging")
In [22]:
ax = grid.plot("Irregular_Kriging*stdev")
ax = data.plot(color="white")
ax.decoration(title="Error for Block Kriging")

Comparison between standard Block Kriging and the Block Kriging with Irregular block size

In [23]:
ax = gp.correlation(grid, namex="Block_Kriging*estim", namey="Irregular_Kriging*estim", bins=100)
ax.decoration(title="Estimation", xlabel="Regular Block", ylabel="Irregular Block")
In [24]:
ax = gp.correlation(grid, namex="Block_Kriging*stdev", namey="Irregular_Kriging*stdev", bins=100)
ax.decoration(title="St. Dev.", xlabel="Regular Block", ylabel="Irregular Block")

We can now focus on our target block. Note that, unlike inn the case of regular block Kriging, the extension of the cell is now 35.

In [25]:
err = gl.krigtest(data, grid, model, neigh, node, calcul=gl.EKrigOpt.BLOCK, ndisc=[5,5],
                 flagPerCell = True)
 Target location
 ---------------
 Sample #156 (from 300)
 Coordinate #1 = 150.000000
 Coordinate #2 = 70.000000
 
 Data selected in neighborhood
 -----------------------------
        Rank      Sample          x1          x2      Sector 
           1           2     156.414      61.537           3 
           2           7     140.011      52.437           2 
           3          34     169.243      74.058           5 
           4          36     167.433      74.021           5 
           5          39     159.821      45.060           3 
           6          40     151.089      69.507           4 
           7          49     165.529      66.754           4 
           8          55     127.590      71.435           8 
           9          73     151.278      84.171           6 
          10          88     140.091      58.596           2 
 
 LHS of Kriging matrix (compressed)
 ==================================
 Number of active samples    = 10
 Total number of equations   = 10
 Reduced number of equations = 10
 
        Rank                       1           2           3           4           5 
                    Flag           1           2           3           4           5 
           1           1       1.000       0.244       0.261       0.288       0.284 
           2           2       0.244       1.000       0.008       0.016       0.197 
           3           3       0.261       0.008       1.000       0.653       0.055 
           4           4       0.288       0.016       0.653       1.000       0.061 
           5           5       0.284       0.197       0.055       0.061       1.000 
           6           6       0.453       0.212       0.245       0.282       0.114 
           7           7       0.431       0.069       0.488       0.505       0.173 
           8           8       0.055       0.168       0.000       0.000       0.000 
           9           9       0.159       0.025       0.207       0.237       0.000 
          10          10       0.290       0.540       0.030       0.045       0.147 
 
        Rank                       6           7           8           9          10 
                    Flag           6           7           8           9          10 
           1           1       0.453       0.431       0.055       0.159       0.290 
           2           2       0.212       0.069       0.168       0.025       0.540 
           3           3       0.245       0.488       0.000       0.207       0.030 
           4           4       0.282       0.505       0.000       0.237       0.045 
           5           5       0.114       0.173       0.000       0.000       0.147 
           6           6       1.000       0.331       0.153       0.332       0.314 
           7           7       0.331       1.000       0.002       0.172       0.103 
           8           8       0.153       0.002       1.000       0.100       0.261 
           9           9       0.332       0.172       0.100       1.000       0.086 
          10          10       0.314       0.103       0.261       0.086       1.000 
 
 RHS of Kriging matrix (compressed)
 ==================================
 Number of active samples    = 10
 Total number of equations   = 10
 Reduced number of equations = 10
 Number of right-hand sides  = 1
 Block Estimation : Discretization =  5  x  5 
 
        Rank        Flag           1 
           1           1       0.317 
           2           2       0.200 
           3           3       0.203 
           4           4       0.228 
           5           5       0.117 
           6           6       0.371 
           7           7       0.255 
           8           8       0.166 
           9           9       0.276 
          10          10       0.267 
 
 (Co-) Kriging weights
 =====================
        Rank          x1          x2       Size1       Size2        Data         Z1* 
           1     156.414      61.537      35.000      35.000       0.946       0.112 
           2     140.011      52.437      35.000      35.000      -0.605       0.052 
           3     169.243      74.058      35.000      35.000       0.801       0.031 
           4     167.433      74.021      35.000      35.000       0.447       0.052 
           5     159.821      45.060      35.000      35.000       0.932       0.025 
           6     151.089      69.507      35.000      35.000      -0.248       0.165 
           7     165.529      66.754      35.000      35.000       0.166       0.065 
           8     127.590      71.435      35.000      35.000       1.070       0.083 
           9     151.278      84.171      35.000      35.000      -0.040       0.155 
          10     140.091      58.596      35.000      35.000      -1.474       0.107 
 Sum of weights                                                           0.846 
 
 (Co-) Kriging results
 =====================
 Target Sample = 156
 Variable Z1 
  - Estimate  =        0.040 
  - Std. Dev. =        0.223 
  - Variance  =        0.050 
  - Cov(h=0)  =        0.279