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,]   389.529  1513.801 -1479.786
##      [  1,]  1513.801325227.614-133555.049
##      [  2,] -1479.786-133555.049273226.914
## Spherical
## - Sill matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]    52.325  4340.380 -3648.747
##      [  1,]  4340.380558573.145-505121.489
##      [  2,] -3648.747-505121.489541429.895
## - Range        =     24.750
## Total Sill
##                [,  0]    [,  1]    [,  2]
##      [  0,]   441.854  5854.181 -5128.533
##      [  1,]  5854.181883800.759-638676.538
##      [  2,] -5128.533-638676.538814656.809
## 
## 
## 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.000      0.015
## Weight of Z2 for Z*1      0.000      0.001
## Weight of Z3 for Z*1      0.000      0.000
## Weight of Z1 for Z*2     -0.007      0.178
## Weight of Z2 for Z*2     -0.001      0.098
## Weight of Z3 for Z*2     -0.076      0.003
## Weight of Z1 for Z*3     -0.004      0.048
## Weight of Z2 for Z*3     -0.056      0.002
## Weight of Z3 for Z*3      0.000      0.119

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)