Kriging with an external covariance function

Didier Renard, Léa Pannecoucke and Xavier Freulon

November 2019

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.

knitr::opts_chunk$set(echo = TRUE, fig.align="center")
library(RGeostats)
## Loading required package: Rcpp
## Loading RGeostats - Version:12.0.2
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.

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:

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

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.

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:

Finally the function my_cov_func needs to access three variables:

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.

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.

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.

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