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.getNSample()
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]:
gp.plot(grid, "Point_Kriging*estim")
gp.plot(data, c="white")
gp.decoration(title="Point Kriging")
In [9]:
gp.plot(grid,"Point_Kriging*stdev")
gp.plot(data, c="white")
gp.decoration(title="Error for Point Kriging")
To better visualize the differences amongst 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)
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
=====================
Dimension of the Covariance Matrix = 10
Rank 1 2 3 4 5
1 1.000 0.244 0.261 0.288 0.284
2 0.244 1.000 0.008 0.016 0.197
3 0.261 0.008 1.000 0.653 0.055
4 0.288 0.016 0.653 1.000 0.061
5 0.284 0.197 0.055 0.061 1.000
6 0.453 0.212 0.245 0.282 0.114
7 0.431 0.069 0.488 0.505 0.173
8 0.055 0.168 0.000 0.000 0.000
9 0.159 0.025 0.207 0.237 0.000
10 0.290 0.540 0.030 0.045 0.147
Rank 6 7 8 9 10
1 0.453 0.431 0.055 0.159 0.290
2 0.212 0.069 0.168 0.025 0.540
3 0.245 0.488 0.000 0.207 0.030
4 0.282 0.505 0.000 0.237 0.045
5 0.114 0.173 0.000 0.000 0.147
6 1.000 0.331 0.153 0.332 0.314
7 0.331 1.000 0.002 0.172 0.103
8 0.153 0.002 1.000 0.100 0.261
9 0.332 0.172 0.100 1.000 0.086
10 0.314 0.103 0.261 0.086 1.000
RHS of Kriging matrix
=====================
Number of active samples = 10
Total number of equations = 10
Number of right-hand sides = 1
Punctual Estimation
Rank 1
1 0.428
2 0.215
3 0.225
4 0.262
5 0.102
6 0.669
7 0.305
8 0.172
9 0.342
10 0.322
(Co-) Kriging weights
=====================
Rank Data Z1*
1 0.946 0.122
2 -0.605 0.016
3 0.801 0.004
4 0.447 0.027
5 0.932 -0.017
6 -0.248 0.513
7 0.166 0.038
8 1.070 0.048
9 -0.040 0.126
10 -1.474 0.091
Sum of weights 0.968
Drift or Mean Information
=========================
Mean for Variable Z1 = 0.000000
(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]:
ndiscs = [5,5]
krigopt = gl.KrigOpt()
krigopt.setOptionCalcul(gl.EKrigOpt.BLOCK, ndiscs)
err = gl.kriging(data, grid, model, neigh, krigopt=krigopt,
namconv=gl.NamingConvention("Block_Kriging"))
The output Db now contains the two additional variables Block_Kriging.Simu.estim and Block_Kriging.Simu.stdev.
Displaying the maps of the Point Estimation and the corresponding Standard Deviation of Estimation errors
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
In [13]:
gp.plot(grid, "Block_Kriging*estim")
gp.plot(data, c="white")
gp.decoration(title="Block Kriging")
In [14]:
gp.plot(grid, "Block_Kriging*stdev")
gp.plot(data, c="white")
gp.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]:
gp.correlation(grid, namex="Point_Kriging*estim", namey="Block_Kriging*estim", bins=100)
gp.decoration(title="Estimation", xlabel="Point", ylabel="Block")
In [16]:
gp.correlation(grid, namex="Point_Kriging*stdev", namey="Block_Kriging*stdev", bins=100)
gp.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]:
krigopt = gl.KrigOpt()
krigopt.setOptionCalcul(gl.EKrigOpt.BLOCK, ndiscs)
err = gl.krigtest(data, grid, model, neigh, node, krigopt)
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
=====================
Dimension of the Covariance Matrix = 10
Rank 1 2 3 4 5
1 1.000 0.244 0.261 0.288 0.284
2 0.244 1.000 0.008 0.016 0.197
3 0.261 0.008 1.000 0.653 0.055
4 0.288 0.016 0.653 1.000 0.061
5 0.284 0.197 0.055 0.061 1.000
6 0.453 0.212 0.245 0.282 0.114
7 0.431 0.069 0.488 0.505 0.173
8 0.055 0.168 0.000 0.000 0.000
9 0.159 0.025 0.207 0.237 0.000
10 0.290 0.540 0.030 0.045 0.147
Rank 6 7 8 9 10
1 0.453 0.431 0.055 0.159 0.290
2 0.212 0.069 0.168 0.025 0.540
3 0.245 0.488 0.000 0.207 0.030
4 0.282 0.505 0.000 0.237 0.045
5 0.114 0.173 0.000 0.000 0.147
6 1.000 0.331 0.153 0.332 0.314
7 0.331 1.000 0.002 0.172 0.103
8 0.153 0.002 1.000 0.100 0.261
9 0.332 0.172 0.100 1.000 0.086
10 0.314 0.103 0.261 0.086 1.000
RHS of Kriging matrix
=====================
Number of active samples = 10
Total number of equations = 10
Number of right-hand sides = 1
Block Estimation : Discretization = 5 x 5
Rank 1
1 0.420
2 0.214
3 0.224
4 0.259
5 0.103
6 0.598
7 0.302
8 0.172
9 0.338
10 0.318
(Co-) Kriging weights
=====================
Rank Data Z1*
1 0.946 0.143
2 -0.605 0.021
3 0.801 0.007
4 0.447 0.034
5 0.932 -0.015
6 -0.248 0.414
7 0.166 0.048
8 1.070 0.056
9 -0.040 0.144
10 -1.474 0.104
Sum of weights 0.956
Drift or Mean Information
=========================
Mean for Variable Z1 = 0.000000
(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]:
krigopt = gl.KrigOpt()
krigopt.setOptionCalcul(gl.EKrigOpt.BLOCK, ndiscs, True);
err = gl.krigcell(data, grid, model, neigh, krigopt = krigopt,
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]:
gp.plot(grid, "Irregular_Kriging*estim")
gp.plot(data, c="white")
gp.decoration(title="Block Kriging")
In [22]:
gp.plot(grid, "Irregular_Kriging*stdev")
gp.plot(data, c="white")
gp.decoration(title="Error for Block Kriging")
Comparison between standard Block Kriging and the Block Kriging with Irregular block size
In [23]:
gp.correlation(grid, namex="Block_Kriging*estim", namey="Irregular_Kriging*estim", bins=100)
gp.decoration(title="Estimation", xlabel="Regular Block", ylabel="Irregular Block")
In [24]:
gp.correlation(grid, namex="Block_Kriging*stdev", namey="Irregular_Kriging*stdev", bins=100)
gp.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]:
krigopt = gl.KrigOpt()
krigopt.setOptionCalcul(gl.EKrigOpt.BLOCK, ndiscs, True)
err = gl.krigtest(data, grid, model, neigh, node, krigopt)
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
=====================
Dimension of the Covariance Matrix = 10
Rank 1 2 3 4 5
1 1.000 0.244 0.261 0.288 0.284
2 0.244 1.000 0.008 0.016 0.197
3 0.261 0.008 1.000 0.653 0.055
4 0.288 0.016 0.653 1.000 0.061
5 0.284 0.197 0.055 0.061 1.000
6 0.453 0.212 0.245 0.282 0.114
7 0.431 0.069 0.488 0.505 0.173
8 0.055 0.168 0.000 0.000 0.000
9 0.159 0.025 0.207 0.237 0.000
10 0.290 0.540 0.030 0.045 0.147
Rank 6 7 8 9 10
1 0.453 0.431 0.055 0.159 0.290
2 0.212 0.069 0.168 0.025 0.540
3 0.245 0.488 0.000 0.207 0.030
4 0.282 0.505 0.000 0.237 0.045
5 0.114 0.173 0.000 0.000 0.147
6 1.000 0.331 0.153 0.332 0.314
7 0.331 1.000 0.002 0.172 0.103
8 0.153 0.002 1.000 0.100 0.261
9 0.332 0.172 0.100 1.000 0.086
10 0.314 0.103 0.261 0.086 1.000
RHS of Kriging matrix
=====================
Number of active samples = 10
Total number of equations = 10
Number of right-hand sides = 1
Block Estimation : Discretization = 5 x 5
Rank 1
1 0.317
2 0.200
3 0.203
4 0.228
5 0.117
6 0.371
7 0.255
8 0.166
9 0.276
10 0.267
(Co-) Kriging weights
=====================
Rank Data Z1*
1 0.946 0.112
2 -0.605 0.052
3 0.801 0.031
4 0.447 0.052
5 0.932 0.025
6 -0.248 0.165
7 0.166 0.065
8 1.070 0.083
9 -0.040 0.155
10 -1.474 0.107
Sum of weights 0.846
Drift or Mean Information
=========================
Mean for Variable Z1 = 0.000000
(Co-) Kriging results
=====================
Target Sample = 156
Variable Z1
- Estimate = 0.040
- Std. Dev. = 0.223
- Variance = 0.050
- Cov(h=0) = 0.279