Factorial Kriging Analysis¶
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import matplotlib.pyplot as plt
import numpy as np
import os
gdoc.setNoScroll()
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]:
ax = gp.correlation(grid, namex="Cr", namey="P", bins=100)
In [6]:
ax = gp.correlation(grid, namex="Ni", namey="P", bins=100)
In [7]:
ax = gp.correlation(grid, namex="Ni", namey="Cr", bins=100)
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]:
gp.setDefaultGeographic(dims=[8,8])
ax = grid.plot("Fill.P")
Variogram Calculation along Grid main axes
In [10]:
varioparam = gl.VarioParam.createMultipleFromGrid(grid, npas=100)
varioP = gl.Vario(varioparam)
err = varioP.compute(grid)
modelP = gl.Model()
err = modelP.fit(varioP, types=[gl.ECov.NUGGET, gl.ECov.SPHERICAL, gl.ECov.POWER],
optvar=gl.Option_VarioFit(True, False))
modelP.setDriftIRF(0,0)
In [11]:
ax = gp.varmod(varioP, modelP)
We must define the Neighborhood
In [12]:
neigh = gl.NeighImage([10,10])
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]:
modelP.setCovaFiltered(0, True)
means = gl.dbStatisticsMono(grid,["Fill.P"],[gl.EStatOption.MEAN]).getValues()
modelP.setMeans(means)
modelP
Out[13]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 3 Number of drift function(s) = 1 Number of drift equation(s) = 1 Covariance Part --------------- Nugget Effect - Sill = 273.033 (This component is Filtered) Spherical - Sill = 56.100 - Range = 3.598 Power (Third Parameter = 0.0879341) - Slope = 1.187 Drift Part ---------- Universality_Condition
In [14]:
err = gl.krimage(grid,modelP,neigh,namconv=gl.NamingConvention("Mono"))
In [15]:
ax = grid.plot("Mono*.P")
ax.decoration(title="P denoised (monovariate)")