The Grid containing the information is downloaded from the distribution

fileNF = loadData("FKA", "Image.ascii")
grid = DbGrid_createFromNF(fileNF)

ndim = 2
defineDefaultSpace(ESpaceType_RN(), ndim)
## NULL

The loaded file (called grid ) contains 3 variables:

dbfmt = DbStringFormat()
dbfmt$setFlags(flag_resume=FALSE,flag_vars=FALSE,flag_stats=TRUE, names="P")
## NULL
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
## NULL

Note that some pixels are not informed for variable P.

dbfmt$setFlags(flag_resume=FALSE,flag_vars=FALSE,flag_stats=TRUE, names="Ni")
## NULL
grid$display(dbfmt)
## 
## Data Base Grid Characteristics
## ==============================
## 
## Data Base Statistics
## --------------------
## 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
## NULL
dbfmt$setFlags(flag_resume=FALSE,flag_vars=FALSE,flag_stats=TRUE, names="Cr")
## NULL
grid$display(dbfmt)
## 
## Data Base Grid Characteristics
## ==============================
## 
## 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
## NULL
plot.init() + plot.correlation(grid,namex="Cr",namey="P", bins=100)

plot.init() + plot.correlation(grid,namex="Ni",namey="P", bins=100)

plot.init() + plot.correlation(grid,namex="Ni",namey="Cr", bins=100)

Using inverse square distance

grid$setLocator("P",ELoc_Z())
## NULL
err = DbHelper_dbgrid_filling(grid)

We concentrate on the variable of interest P (completed) in the next operations

p = plot.init()
p = p + plot(grid)
p = p + plot.decoration(title="P after completion")
plot.end(p)

varioparam = VarioParam_createMultipleFromGrid(grid,nlag=100)
varioP = Vario(varioparam)
err = varioP$compute(grid)
modelP = Model()
types = ECov_fromKeys(c("NUGGET", "SPHERICAL"))
err = modelP$fit(varioP, types=types, optvar=Option_VarioFit(TRUE, FALSE))
modelP$setDriftIRF(0,0)
## NULL
modelP
## 
## 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.014
## Spherical
## - Sill         =      75.975
## - Range        =       5.705
## Total Sill     =     431.989
## 
## Drift Part
## ----------
## Universality_Condition
plot.init() + plot.varmod(varioP,modelP)

We must define the Neighborhood

neigh = NeighImage(c(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.

modelP$setCovFiltered(0, TRUE)
## NULL
means = dbStatisticsMono(grid,"Fill.P",EStatOption_fromKeys("MEAN"))$getValues()
modelP$setMeans(means)
## NULL
modelP
## 
## 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.014
##   (This component is Filtered)
## Spherical
## - Sill         =      75.975
## - Range        =       5.705
## Total Sill     =     431.989
## 
## Drift Part
## ----------
## Universality_Condition
err = krimage(grid,modelP,neigh,flagFFT=TRUE, verbose=TRUE, namconv=NamingConvention("Mono"))
## 
## Convolution Pattern
## -------------------
##                          Minimum     Maximum
## Weight of Z1 for Z*1      -0.004       0.099
p = plot.init()
p = p + plot(grid, "Mono*.P")
p = p + plot.decoration(title="P denoised (monovariate)")
plot.end(p)

Correlation for P variable between Initial image and its Filtered version (monovariate FKA)

p = plot.init()
p = p + plot.correlation(grid,namex="Fill.P",namey="Mono.Fill.P", bins=100)
p = p + plot.decoration(xlab="P Filled",ylab="P Filtered")
plot.end(p)

varnames = c("Fill.P", "Cr", "Ni")
grid$setLocators(varnames, ELoc_Z())
## NULL
varioM = Vario(varioparam)
err = varioM$compute(grid)
modelM = Model()
err = modelM$fit(varioM, types= types, optvar=Option_VarioFit(TRUE, FALSE))
modelM$setDriftIRF(0,0)
## NULL
modelM
## 
## 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.141    235.177   -300.219
##       [  1,]    235.177 167267.806  12741.118
##       [  2,]   -300.219  12741.118 119349.758
## Spherical
## - Sill matrix:
##                  [,  0]     [,  1]     [,  2]
##       [  0,]     61.993   5144.087  -4483.493
##       [  1,]   5144.087 651989.758-594845.183
##       [  2,]  -4483.493-594845.183 633142.577
## - Range        =      10.595
## Total Sill
##                  [,  0]     [,  1]     [,  2]
##       [  0,]    436.134   5379.264  -4783.712
##       [  1,]   5379.264 819257.564-582104.065
##       [  2,]  -4783.712-582104.065 752492.335
## 
## 
## Drift Part
## ----------
## Universality_Condition
multi.varmod(varioM,modelM)

modelM$setCovFiltered(0, TRUE)
## NULL
means = dbStatisticsMono(grid,flagIso=TRUE, names=varnames,
                         EStatOption_fromKeys("MEAN"))$getValues()
modelM$setMeans(means)
## NULL
err  = krimage(grid,modelM,neigh,flagFFT=TRUE, verbose=TRUE, namconv=NamingConvention("Multi"))
## 
## Convolution Pattern
## -------------------
##                          Minimum     Maximum
## Weight of Z1 for Z*1      -0.001       0.032
## Weight of Z2 for Z*1       0.000       0.002
## Weight of Z3 for Z*1      -0.002       0.000
## Weight of Z1 for Z*2      -0.011       0.522
## Weight of Z2 for Z*2      -0.003       0.281
## Weight of Z3 for Z*2      -0.274       0.011
## Weight of Z1 for Z*3      -0.151       0.043
## Weight of Z2 for Z*3      -0.202       0.008
## Weight of Z3 for Z*3      -0.003       0.354

Note that, using the same neigh as in monovariate, the dimension of the Kriging System is now \(3 * 441 = 1323\)

p = plot.init()
p = p + plot(grid,"Multi*.P", zlim=c(0,150))
## Warning in geom_raster(data = df, mapping = aes(x = x, y = y, fill = data), :
## Ignoring unknown parameters: `zlim`
p = p + plot.decoration(title="P denoised (multivariate)")
plot.end(p)

Correlation for P variable between Initial image and its Filtered version (multivariate FKA)

p = plot.init()
p = p + plot.correlation(grid,namex="Fill.P",namey="Multi.Fill.P", bins=100)
p = p + plot.decoration(xlab="P Filled",ylab="P Filtered (Multi)")
plot.end(p)

Correlation for P filtered variable between the Monovariate and the Multivariate case

p = plot.init()
p = p + plot.correlation(grid, namex="Mono.Fill.P", namey="Multi.Fill.P", bins=100)
p = p + plot.decoration(xlab="P Filtered (Mono)",ylab="P Filtered (Multi)")
plot.end(p)