---
title: "Factorial Kriging Analysis"
author: "D. Renard, N. Desassis"
date: "22 August 2017"
output:
  rmarkdown::html_vignette:
    fig_width: 5
    fig_height: 4
vignette: >
  %\VignetteIndexEntry{Factorial Kriging Analysis}             
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r Loading_library,include=FALSE,echo=FALSE}
library(RGeostats)
constant.define("asp",1)
rm(list=ls())
```

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.
```{r Loading Data}
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. Then, we upscale the grid (by dividing the number of cells) in order to reduce execution time.

```{r Statistics}
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.

```{r Correlation_Cr_P}
correlation(grid,"Cr","P",title="Correlation between P and Cr")
```

```{r Correlation_Ni_P}
correlation(grid,"Ni","P",title="Correlation between P and Ni")
```

```{r Correlation_Ni_Cr}
correlation(grid,"Ni","Cr",title="Correlation between Cr and Ni")
```

# 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.

```{r Completing P}
grid = db.locate(grid,"P","z")
grid = db.grid.fill(grid,radix="P.Fill")
```

We concentrate on the variable of interest **P**. 
```{r Locate 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).
```{r Monovariate_Variogram_calculation}
varioP = vario.grid(grid,nlag=40)
modelP = model.auto(varioP,struct=c(1,2,12),
                    auth.aniso = FALSE,flag.noreduce=TRUE)
print(modelP)
```

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).

```{r Visualize_P}
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.

```{r Neighborhood Definition}
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 

```{r Performing the monovariate FKA}
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).

```{r Visualize_P_after_FKA_monovariate}
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.
```{r Locate Multivariate}
grid = db.locate(grid,c("P.Fill*","Cr","Ni"),"z")
```

Calculate the multivariate experimental variogram benefiting from the grid organization and fit the model.
```{r Multivariate_Variogram_calculation}
varioM = vario.grid(grid,nlag=40)
modelM = model.auto(varioM,struct=c(1,2,12),
                    auth.aniso = FALSE,flag.noreduce=TRUE)
modelM
```

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
```{r Perform multivariate FKA}
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).

```{r Visualize_P_after_FKA_multivariate}
plot(multi,pos.legend=1,title="P denoised (multivariate)",
     zlim=c(0,150))
```