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]:
ndiscs = [5,5]
err = gl.kriging(data, grid, model, neigh, calcul=gl.EKrigOpt.BLOCK, ndiscs=ndiscs,
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, ndiscs=ndiscs)
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, ndiscs=ndiscs,
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, ndiscs=ndiscs,
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