---
title: "Kriging with an external covariance function"
author: "Didier Renard, Léa Pannecoucke and Xavier Freulon"
date: "November 2019"
output:
  rmarkdown::html_vignette:
    fig_width: 5
    fig_height: 4
vignette: >
  %\VignetteIndexEntry{Kriging with an external covariance function}             
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

# Introduction

This vignette illustrates the use of an external covariance function with **RGeostats** in a kriging procedure. For some geostatistical methods, such as kriging with numerical covariances or variograms (see for example Pannecoucke et al., 2020. Combining geostatistics and simulations of flow and transport to
characterize contamination within the unsaturated zone. *Science of the Total Environment*), the use of an external covariance function is needed. 

In those applications, the classical covariance model is replaced by a covariance matrix externally built by the user.

```{r Loading_library,include=FALSE,echo=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align="center")
library(RGeostats)
constant.define("asp",1)
```

# Kriging with an external covariance function

## Grid and data
The target grid is a 50 $\times$ 50 square, with 5 observation points. In the following, the target database (`dbTarget`) and the observation database (`dbObs`) are created. Values and coordinates of observations are randomly chosen.

```{r}
nx = 50
x0 = 0
dx = 1
dbTarget = db.create(nx=rep(nx+1,2),x0=rep(x0,2),dx=rep(dx,2))
dbTarget = db.locate(dbTarget,2:3,"x")

nObs = 5
set.seed(99)
coordObs = cbind(sample(seq(nx),nObs,replace=F),sample(seq(nx),nObs,replace=F))
valObs = rnorm(nObs)
dbObs = db.create(cbind(coordObs,valObs))
dbObs = db.locate(dbObs,2:3,"x")
dbObs = db.locate(dbObs,4,"z")
plot(dbObs,pos.legend=1,title="Observations")
```

## Covariance matrices
Two covariance matrices and one vector are computed to perform kriging with an external covariance function:

* `cxx`: the observation-observation matrix,

* `cx0`: the observation-target matrix,

* `c00`: the $C(0)$ vector needed to compute kriging error variance.

In this example, the covariance matrices are computed from the exponential model, with range 10 and sill 1. 

```{r}
exponentialCov <- function(h,range,sill){
  sill*exp(-h/range)
}
a = 10
C = 1
```

In the following, we compute the distance matrices between pairs of observations and pairs of targets and observations, and then `cxx`, `cx0` and `c00`.

```{r}
distObs = sqrt(outer(coordObs[,1],coordObs[,1],"-")^2 + outer(coordObs[,2],coordObs[,2],"-")^2)
cxx = exponentialCov(distObs,range=a,sill=C)

distTar = sqrt(outer(coordObs[,1],dbTarget[,"x1"],"-")^2 + outer(coordObs[,2],dbTarget[,"x2"],"-")^2)
cx0 = exponentialCov(distTar,range=a,sill=C)

c00 = exponentialCov(rep(0,dbTarget$nech),range=a,sill=C)
```

## Kriging

In order to use `cxx` and `cx0` instead of a classical model, we must define the external covariance function `my_cov_func`. This function is called internally during the kriging procedure (here via the `def_cov` function). Its prototype contains the following arguments:

* `dist`: the distance between two points (not used),

* `db1`: the type of `db` (1 if `dbin` or 2 if `dbout`) of the first point,

* `iech1`: the index of the first point in `db1`,

* `db2`: the type of db of the second point,

* `iech2`: the index of the second point in `db2`,

* `incr`: not used,

* `x1`: the coordinate vector of the first point (not used),

* `x2`: the coordinate vector of the second point (not used).

Finally the function `my_cov_func` needs to access three variables:

* `cxx`

* `cx0`

* `c00`

```{r}
def_cov <- function(cxx,cx0,c00) {
  function(dist,db1=NA,iech1=NA,db2=NA,iech2=NA,incr=NA,x1=NA,x2=NA,...){
    result = 0
    if (db1 == 1 && db2 == 1)
        result = cxx[iech1,iech2] 
    else if (db1 == 1 && db2 == 2)
        result = cx0[iech1,iech2]
    else if (db1 == 2 && db2 == 1) 
        result = cx0[iech2,iech1]
    else if (db1 == 2 && db2 == 2)
        result = c00[iech1]
    else
        cat("db1=",db1," db2=",db2,"\n")
    result
  }
}

my_cov_func <- def_cov(cxx=cxx,cx0=cx0,c00=c00)
```

Finally, we perform ordinary kriging on a unique neighborhood.

```{r}
unique.neigh = neigh.create(type=0)
dbTarget = kriging(dbin=dbObs,dbout=dbTarget,model=my_cov_func,neigh=unique.neigh,uc="1")
plot(dbTarget,pos.legend=1,title="Kriging estimate")
plot(dbTarget,name="*.stdev",pos.legend=1,title="Standard deviation of kriging error")
```

## Comparison with classical procedure
We compare the results obtained previously with the results obtained with the classical procedure. We first define an exponential model with the same parameters and then perform ordinary kriging with this model. 

```{r} 
model = model.create(vartype="exponential",scale=a,sill=C)
dbTarget = kriging(dbin=dbObs,dbout=dbTarget,model=model,neigh=unique.neigh,uc="1",radix="OK")
```

The correlation plots show that the estimations and standard deviations are identical with the two methods.

```{r}
dummy = correlation(dbTarget,name1="Kriging.*.estim",name2="OK.*.estim",
                    flag.iso=T,flag.diag=T,
                    title="Estimations")
dummy = correlation(dbTarget,name1="Kriging.*.stdev",name2="OK.*.stdev",
                    flag.iso=T,flag.diag=T,
                    title="Standard deviations")
```