This paper is meant to demonstrate the capability of Filtering by Kriging using Factorial Kriging Analysis (FKA), performed on data gathered on a regular grid.
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:
- P (phosphorus) which is the variable of interest
- Cr (chromium) is an auxiliary variable
- Ni (nickel) another auxiliary variable
The first task is to get some statistics on these variables where we can check that only P is not defined everywhere.
print(grid,flag.stats=TRUE,names=c("P","Cr","Ni"))
##
## Data Base Characteristics
## =========================
##
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension = 2
## Number of fields = 6
## Maximum Number of attributes = 6
## Total number of samples = 262144
## Grid characteristics:
## Origin : 0.000 0.000
## Mesh : 1.000 1.000
## Number : 512 512
## Angles : 0.000 0.000
##
## 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
## 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
##
## Variables
## ---------
## Field = 1 - Name = rank - Locator = NA
## Field = 2 - Name = x1 - Locator = x1
## Field = 3 - Name = x2 - Locator = x2
## Field = 4 - Name = Cr - Locator = NA
## Field = 5 - Name = Ni - Locator = NA
## Field = 6 - Name = P - Locator = NA
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.2771063
correlation(grid,"Ni","P",title="Correlation between P and Ni")
## [1] -0.2545205
correlation(grid,"Ni","Cr",title="Correlation between Cr and Ni")
## [1] -0.7099362
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=100)
modelP = model.auto(varioP,struct=c(1,2,12),
auth.aniso = FALSE,flag.noreduce=TRUE)
##
## 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 = 334.273
## - Exponential
## Range = 5.179
## Theo. Range = 1.729
## Sill = 91.519
## - Linear
## Range = 99.000
## Slope = 0.238
##
## 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))
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=100)
modelM = model.auto(varioM,struct=c(1,2,12),
auth.aniso = FALSE,flag.noreduce=TRUE)
##
## 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 =
## [,1] [,2] [,3]
## [1,] 390.687 -96.118 -58.567
## [2,] -96.118 99696.440 75185.550
## [3,] -58.567 75185.550 56701.285
## - Exponential
## Range = 12.375
## Theo. Range = 4.131
## Sill matrix =
## [,1] [,2] [,3]
## [1,] 41.609 5420.531 -4860.501
## [2,] 5420.531710096.587-654766.91
## [3,] -4860.501-654766.91685816.033
## - Linear
## Range = 99.000
## Sill matrix =
## [,1] [,2] [,3]
## [1,] 10.952 279.392 372.939
## [2,] 279.392 52544.620-28142.980
## [3,] 372.939-28142.980 52396.200
##
## 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))
