Factorial Kriging Analysis¶
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import numpy as np
import os
gdoc.setNoScroll()
gl.OptCst.define(gl.ECst.NTROW, -1)
gl.OptCst.define(gl.ECst.NTCOL, -1)
gl.OptCst.define(gl.ECst.NTCAR, 15)
figsize=(8,8)
The Grid containing the information is downloaded from the distribution.
The loaded file (called **grid **) contains 3 variables:
- P (phosphorus) which is the variable of interest
- Cr (chromium) is an auxiliary variable
- Ni (nickel) another auxiliary variable
In [2]:
fileNF = gdoc.loadData("FKA", "Image.ascii")
grid = gl.DbGrid.createFromNF(fileNF)
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
In [3]:
dbfmt = gl.DbStringFormat()
dbfmt.setFlags(flag_resume=False,flag_vars=False,flag_stats=True, names="P")
grid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Statistics -------------------- 6 - Name P - Locator NA Nb of data = 262144 Nb of active values = 242306 Minimum value = 0.000 Maximum value = 314.000 Mean value = 31.767 Standard Deviation = 21.759 Variance = 473.457
Note that some pixels are not informed for variable P.
Statistics on auxiliary variables
In [4]:
dbfmt.setFlags(flag_vars=False, flag_resume=True, flag_stats=True, names=["Cr", "Ni"])
grid.display(dbfmt)
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 6 Total number of samples = 262144 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 1.000 1.000 Number : 512 512 Data Base Statistics -------------------- 4 - Name Cr - Locator NA Nb of data = 262144 Nb of active values = 262144 Minimum value = 2591.000 Maximum value = 24982.000 Mean value = 16800.231 Standard Deviation = 936.213 Variance = 876495.558 5 - Name Ni - Locator NA Nb of data = 262144 Nb of active values = 262144 Minimum value = 1840.000 Maximum value = 12593.000 Mean value = 10111.444 Standard Deviation = 884.996 Variance = 783217.898
Correlation between variables
In [5]:
res = gp.correlation(grid, namex="Cr", namey="P", bins=100, cmin=1)
In [6]:
res = gp.correlation(grid, namex="Ni", namey="P", bins=100, cmin=1)
In [7]:
res = gp.correlation(grid, namex="Ni", namey="Cr", bins=100, cmin=1)
Using inverse square distance for completing the variable P
In [8]:
grid.setLocator("P", gl.ELoc.Z)
err = gl.DbHelper.dbgrid_filling(grid,0,13432,1)
We concentrate on the variable of interest P completed (Fill.P) in the next operations
In [9]:
fig, ax = gp.init(figsize=figsize, flagEqual=True)
ax.raster(grid,"Fill.P")
gp.close()
Variogram Calculation along Grid main axes
In [10]:
varnames = ["Fill.P"]
varioparam = gl.VarioParam.createMultipleFromGrid(grid, nlag=100)
varioP = gl.Vario(varioparam)
err = varioP.compute(grid)
modelP = gl.Model()
types = [gl.ECov.NUGGET, gl.ECov.SPHERICAL]
err = modelP.fit(varioP, types=types, optvar=gl.Option_VarioFit(True, False))
modelP.setDriftIRF(0,0)
modelP.setCovFiltered(0, True)
means = gl.dbStatisticsMono(grid,varnames,[gl.EStatOption.MEAN]).getValues()
modelP.setMeans(means)
modelP
Out[10]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 2 Number of drift function(s) = 1 Number of drift equation(s) = 1 Covariance Part --------------- Nugget Effect - Sill = 356.005 (This component is Filtered) Spherical - Sill = 75.980 - Range = 5.703 Total Sill = 431.985 Drift Part ---------- Universality_Condition
In [11]:
res = gp.varmod(varioP, modelP)
We must define the Neighborhood
In [12]:
neigh = gl.NeighImage([5,5])
The image neighborhood is based on (2∗10+1)2=441 pixels (centered on the target pixel).
During the estimation, only the contribution of second and third basic structures are kept (Nugget Effect is filtered out): Factorial Kriging Analysis.
In [13]:
err = gl.krimage(grid,modelP,neigh,flagFFT=True, verbose=True,namconv=gl.NamingConvention("Mono"))
Convolution Pattern
-------------------
Minimum Maximum
Weight of Z1 for Z*1 -0.002 0.100
In [14]:
fig, ax = gp.init(1,2, figsize=figsize, flagEqual=True)
ax[0].raster(grid,"Fill.P")
ax[1].raster(grid,"Mono*.P")
fig.decoration(title="Filtering P (monovariate)")
Correlation for P variable between Initial image (completed) and its Filtered version (monovariate FKA)
In [15]:
gp.correlation(grid, namex="Fill.P", namey="Mono.Fill.P", bins=100, cmin=1)
gp.decoration(xlabel="P Filled",ylabel="P Filtered (Mono)")
Multivariate approach¶
Having a look at the two covariables
In [16]:
fig, ax = gp.init(1,2, figsize=figsize, flagEqual=True)
ax[0].raster(grid,"Cr")
ax[1].raster(grid,"Ni")
fig.decoration(title="Covariables")
In [17]:
varnames = ["Fill.P", "Ni", "Cr"]
grid.setLocators(varnames, gl.ELoc.Z)
varioM = gl.Vario(varioparam)
err = varioM.compute(grid)
modelM = gl.Model()
err = modelM.fit(varioM, types=types, optvar=gl.Option_VarioFit(True, False))
modelM.setDriftIRF(0,0)
modelM.setCovFiltered(0, True)
means = gl.dbStatisticsMono(grid,varnames,[gl.EStatOption.MEAN]).getValues()
modelM.setMeans(means)
modelM
Out[17]:
Model characteristics
=====================
Space dimension = 2
Number of variable(s) = 3
Number of basic structure(s) = 2
Number of drift function(s) = 1
Number of drift equation(s) = 3
Covariance Part
---------------
Nugget Effect
- Sill matrix:
[, 0] [, 1] [, 2]
[ 0,] 391.152 -1595.520 1651.095
[ 1,] -1595.520 290378.614 -149556.582
[ 2,] 1651.095 -149556.582 342922.388
(This component is Filtered)
Spherical
- Sill matrix:
[, 0] [, 1] [, 2]
[ 0,] 50.294 -3504.361 4168.490
[ 1,] -3504.361 519976.264 -485106.506
[ 2,] 4168.490 -485106.506 536440.228
- Range = 24.750
Total Sill
[, 0] [, 1] [, 2]
[ 0,] 441.446 -5099.881 5819.585
[ 1,] -5099.881 810354.878 -634663.088
[ 2,] 5819.585 -634663.088 879362.616
Drift Part
----------
Universality_Condition
In [18]:
res = gp.varmod(varioM, modelM)
Multivariable Factorial Kriging Analysis
In [19]:
err = gl.krimage(grid,modelM,neigh,flagFFT=True,verbose=True,namconv=gl.NamingConvention("Multi"))
Convolution Pattern
-------------------
Minimum Maximum
Weight of Z1 for Z*1 0.006 0.017
Weight of Z2 for Z*1 0.000 0.000
Weight of Z3 for Z*1 0.000 0.001
Weight of Z1 for Z*2 -0.016 0.043
Weight of Z2 for Z*2 -0.001 0.111
Weight of Z3 for Z*2 -0.052 0.002
Weight of Z1 for Z*3 -0.027 0.161
Weight of Z2 for Z*3 -0.070 0.003
Weight of Z3 for Z*3 0.000 0.092
Note that, using the same neigh as in monovariate, the dimension of the Kriging System is now 3∗441=1323
In [20]:
fig, ax = gp.init(1,2, figsize=figsize, flagEqual=True)
ax[0].raster(grid,"Mono*.P")
ax[1].raster(grid,"Multi*.P")
fig.decoration(title="P filtered")
Correlation for P variable between Initial image and its Filtered version (multivariate FKA)
In [21]:
gp.correlation(grid, namex="Fill.P", namey="Multi.Fill.P", bins=100, cmin=1)
gp.decoration(xlabel="P Filled",ylabel="P Filtered (Multi)")
Correlation for P filtered variable between the Monovariate and the Multivariate case
In [22]:
gp.correlation(grid, namex="Mono.Fill.P", namey="Multi.Fill.P", bins=100, cmin=1)
gp.decoration(xlabel="P Filtered (Mono)",ylabel="P Filtered (Multi)")
