This lecture demonstrate the capability of Filtering by Kriging using Factorial Kriging Analysis (FKA), performed on data gathered on a regular grid.

1 Loading Data

The data are stored in a 2-D grid (of 512 by 512 pixels) which is downloaded from the distribution kit and renamed for easier manipulation.

rg.load("Exdemo_Filter","grid")

The file contains 3 variables:

The first task is to get some statistics on these variables where we can check that only P is not defined everywhere. Then, we upscale the grid (by dividing the number of cells) in order to reduce execution time.

subdiv = 2
grid = db.grid.refine(grid, nmult = c(subdiv, subdiv), flag.refine = F, flag.copy = T)
plot(grid,name="Cr")

The different scatter plots show a weak correlation between P and the other variables, but a much stronger relation between the two auxiliary variables.

correlation(grid,"Cr","P",title="Correlation between P and Cr")

## [1] 0.267981
correlation(grid,"Ni","P",title="Correlation between P and Ni")

## [1] -0.24477
correlation(grid,"Ni","Cr",title="Correlation between Cr and Ni")

## [1] -0.7105469

2 Monovariate Factorial Kriging Analysis

The variable P is not completely defined. In order to use efficiently the grid organization of the information, we may use FKA on an image. This requires that P be completed beforehand.

grid = db.locate(grid,"P","z")
grid = db.grid.fill(grid,radix="P.Fill")

We concentrate on the variable of interest P.

grid = db.locate(grid,"P.Fill*","z")

We calculate its experimental variogram benefiting from the grid organization (i.e. along the two main grid axes) and fit it using an automatic procedure (where the anisotropy is discarded).

varioP = vario.grid(grid,nlag=40)
modelP = model.auto(varioP,struct=c(1,2,12),
                    auth.aniso = FALSE,flag.noreduce=TRUE)

print(modelP)
## 
## 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         =    412.308
## Exponential
## - Sill         =     47.104
## - Range        =     15.971
## - Theo. Range  =      5.331
## Linear
## - Slope        =      0.121
## 
## Drift Part
## ----------
## Universality Condition

The model shows a large nugget component (around 80% of the total sill). The principle of the Factorial Kriging Analysis performed in the monovariate case is to postulate that P data contains the part of interest which corresponds to the exponential structure whereas the nugget effect comes from the instrumental noise. FKA is used to extract the structured part.

It is now time to visualize the P variable (after completion).

plot(grid,pos.legend=1,title="P after completion")

This step requires a neighborhood definition. The neighborhood is defined a squared area centered around the target pixel. Note that an extension of \(N=10\) leads to neighborhood composed of \((2 * N + 1) ^2\) samples. In order to reduce the dimension of the kriging system, it would be possible to use a skipping ratio.

neigh = neigh.create(ndim=2,type=3,radimg=c(10,10))

We use the krimage function which performs FKA starting from values defined on the grid and returning the filtered result on the same grid. The extracted component is the one which corresponds to the second terms

means = db.stat(grid,flag.mono=TRUE,names="P.Fill*")
mono = krimage(grid,modelP,neigh,cov.extract = c(2,3),uc="",
               mean = means)

We can visualize the denoised P variable (monovariate FKA).

plot(mono,pos.legend=1,title="P denoised (monovariate)",
     zlim=c(0,150))

3 Multivariate approach

As we have already mentioned the correlation between the main variable (P) and two other auxiliary variables (Cr and Ni), we now wish to apply FKA on a multivariate basis.

We need to define the various variables in the grid file.

grid = db.locate(grid,c("P.Fill*","Cr","Ni"),"z")

Calculate the multivariate experimental variogram benefiting from the grid organization and fit the model.

varioM = vario.grid(grid,nlag=40)
modelM = model.auto(varioM,struct=c(1,2,12),
                    auth.aniso = FALSE,flag.noreduce=TRUE)

modelM
## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 3
## Number of basic structure(s) = 3
## Number of drift function(s)  = 1
## Number of drift equation(s)  = 3
## 
## Covariance Part
## ---------------
## Nugget Effect
## - Sill matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]   427.055  2063.838 -2014.163
##      [  1,]  2063.838384705.973-188976.917
##      [  2,] -2014.163-188976.917329236.126
## Exponential
## - Sill matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]    33.769  3736.507 -3218.138
##      [  1,]  3736.507507048.529-461662.728
##      [  2,] -3218.138-461662.728492881.395
## - Range        =     34.125
## - Theo. Range  =     11.391
## Linear
## - Slope matrix:
##                [,  0]    [,  1]    [,  2]
##      [  0,]     0.167     2.480     2.481
##      [  1,]     2.480    62.132    62.126
##      [  2,]     2.481    62.126    62.120
## 
## 
## Drift Part
## ----------
## Universality Condition

Note the nugget effect disappears on the cross-variograms, which prooves that using the auxiliary variables for filtering the target one will be beneficial.

We can now perform the multivariate FKA operation

means = db.stat(grid,flag.mono=TRUE,names=c("P.Fill*","Cr","Ni"))
multi = krimage(grid,modelM,neigh,cov.extract = c(2,3),uc="",
                mean = means)

We can visualize the denoised P variable (multivariate FKA).

plot(multi,pos.legend=1,title="P denoised (multivariate)",
     zlim=c(0,150))