Factorial Kriging Analysis¶
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
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,] 374.142 -300.529 235.236
[ 1,] -300.529 119376.532 12715.586
[ 2,] 235.236 12715.586 167294.878
(This component is Filtered)
Spherical
- Sill matrix:
[, 0] [, 1] [, 2]
[ 0,] 61.992 -4483.615 5144.024
[ 1,] -4483.615 633134.809 -594837.311
[ 2,] 5144.024 -594837.311 651982.517
- Range = 10.597
Total Sill
[, 0] [, 1] [, 2]
[ 0,] 436.135 -4784.145 5379.260
[ 1,] -4784.145 752511.341 -582121.725
[ 2,] 5379.260 -582121.725 819277.395
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.002 0.033
Weight of Z2 for Z*1 -0.002 0.000
Weight of Z3 for Z*1 0.000 0.002
Weight of Z1 for Z*2 -0.149 0.045
Weight of Z2 for Z*2 -0.002 0.356
Weight of Z3 for Z*2 -0.202 0.009
Weight of Z1 for Z*3 -0.026 0.527
Weight of Z2 for Z*3 -0.275 0.012
Weight of Z3 for Z*3 -0.002 0.283
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)")
