The Grid containing the information is downloaded from the distribution

fileNF = "Image.ascii"
if(flagInternetAvailable){
  download.file(paste0("https://soft.minesparis.psl.eu/gstlearn/data/FKA/",fileNF), fileNF)
}
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
ggplot() + plot.correlation(grid,"Cr","P", bins=100)

ggplot() + plot.correlation(grid,"Ni","P", bins=100)

ggplot() + plot.correlation(grid,"Ni","Cr", bins=100)

Using inverse square distance

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

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

p = ggDefaultGeographic()
p = p + plot(grid)
p = p + plot.decoration(title="P after completion")
ggPrint(p)

varioparam = VarioParam_createMultipleFromGrid(grid,npas=100)
varioP = Vario(varioparam, grid)
err = varioP$compute()
modelP = Model()
err = modelP$fit(varioP, 
                 types=ECov_fromKeys(c("NUGGET", "SPHERICAL", "POWER")),
                 optvar=Option_VarioFit(TRUE, FALSE))
modelP
## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 3
## Number of drift function(s)  = 0
## Number of drift equation(s)  = 0
## 
## Covariance Part
## ---------------
## Nugget Effect
## - Sill         =    273.080
## Spherical
## - Sill         =     56.089
## - Range        =      3.600
## Power (Third Parameter = 0.0879341)
## - Slope        =      1.187
ggplot() + 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$setCovaFiltered(0, TRUE)
## NULL
means = dbStatisticsMono(grid,"Fill.P",EStatOption_fromKeys("MEAN"))
modelP$setMeans(means)
## NULL
modelP
## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 3
## Number of drift function(s)  = 0
## Number of drift equation(s)  = 0
## 
## Covariance Part
## ---------------
## Nugget Effect
## - Sill         =    273.080
##   (This component is Filtered)
## Spherical
## - Sill         =     56.089
## - Range        =      3.600
## Power (Third Parameter = 0.0879341)
## - Slope        =      1.187
err = krimage(grid,modelP,neigh,namconv=NamingConvention("Mono"))
p = ggDefaultGeographic()
p = p + plot(grid, "Mono*.P")
p = p + plot.decoration(title="P denoised (monovariate)")
ggPrint(p)

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

p = ggplot()
p = p + plot.correlation(grid,"Fill.P","Mono.Fill.P", bins=100)
p = p + plot.decoration(xlab="P Filled",ylab="P Filtered")
ggPrint(p)

grid$setLocators(c("Fill.P", "Cr", "Ni"), ELoc_Z())
## NULL
varioM = Vario(varioparam, grid)
err = varioM$compute()
modelM = Model()
err = modelM$fit(varioM, 
                 types= ECov_fromKeys(c("NUGGET", "SPHERICAL", "POWER")),
                 optvar=Option_VarioFit(TRUE, FALSE))
modelM
## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 3
## Number of basic structure(s) = 3
## Number of drift function(s)  = 0
## Number of drift equation(s)  = 0
## 
## Covariance Part
## ---------------
## Nugget Effect
## - Sill matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]   376.850   452.534  -476.811
##      [  1,]   452.534194188.109-11524.845
##      [  2,]  -476.811-11524.845145939.572
## Spherical
## - Sill matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]    57.513  5031.290 -4489.149
##      [  1,]  5031.290636076.559-583291.704
##      [  2,] -4489.149-583291.704616673.997
## - Range        =     12.375
## Power (Third Parameter = 1.99)
## - Slope matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]     0.263    -0.414     6.133
##      [  1,]    -0.414   145.976    44.478
##      [  2,]     6.133    44.478   163.446
multi.varmod(varioM,modelM)

modelM$setCovaFiltered(0, TRUE)
## NULL
means = dbStatisticsMono(grid,flagIso=TRUE, 
                         names=c("Fill.P","Cr","Ni"),
                         EStatOption_fromKeys("MEAN"))
modelM$setMeans(means)
## NULL
err  = krimage(grid,modelM,neigh,namconv=NamingConvention("Multi"))

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

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

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

p = ggplot()
p = p + plot.correlation(grid,"Fill.P","Multi.Fill.P", bins=100)
p = p + plot.decoration(xlab="P Filled",ylab="P Filtered (Multi)")
ggPrint(p)

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

p = ggplot()
p = p + plot.correlation(grid, "Mono.Fill.P", "Multi.Fill.P", bins=100)
p = p + plot.decoration(xlab="P Filtered (Mono)",ylab="P Filtered (Multi)")
ggPrint(p)