---
title: "Cokriging"
author: "D. Renard, N. Desassis, X. Freulon"
date: "3 novembre 2021"
output: beamer_presentation
---

```{r setup, include=FALSE}
library(knitr)
library(RGeostats)
rm(list=ls())
constant.define("asp",1)
knitr::opts_chunk$set(align="center")
```

## Preamble

\footnotesize

```{r Preamble, message=FALSE, warning=FALSE, include=FALSE}
rg.load("Exdemo_Scotland_Temperatures","dat")
rg.load("Exdemo_Scotland_Elevations","target")
dat=db.rename(dat, "January_temp", "Temperature")
dat=db.locate(dat, names = "Temperature", loctype = "z")
unique.neigh=neigh.create(type=0, ndim = 2)
```

## Statistics

\footnotesize

```{r EDA1, echo=FALSE, message=FALSE, warning=FALSE}
tab <- db.stat.multi(dat, c("Elevation", "Temperature"))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on observations")

tab <- db.stat.multi(target, c("Elevation"))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on the grid")
```


## Observations

```{r maps-observations, echo=FALSE}
plot(db.locate(dat, "El*", "z"), col = "gray", pch=19)
plot(dat, col = "red", add = TRUE)
legend("topright", 
       legend=c("Temperature","Elevation"), pch=c(1,19), col=c("red", "gray"))
```

## Map of the elevations

```{r maps-elevation, echo=FALSE}
plot(target, pos.legend=1)
```

## Temperature vs. Elevation

\footnotesize

```{r temp_vs_ele, echo=FALSE}
dum = correlation(dat, name2 = "Temperature", name1 = "Elevation", 
            flag.aspoint = TRUE, flag.regr = TRUE,
            title = "Correlation between Temperature and Elevation")
```

## Temperature vs. Elevation (statistiques)

\footnotesize

```{r temp_vs_ele2, echo=TRUE}
tab <- db.stat.multi(
  db.sel(dat, sel=!is.na(Temperature)), 
  c("Temperature","Elevation"))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on the grid")
```

## Variograms and Non-stationarity

\footnotesize

Directional variograms (N-S and E-W): no drift or with first order global trend

\center

```{r Experimental_Directional_Variograms, echo=FALSE}
vario.raw2dir = vario.calc(dat,lag=10,nlag=30,dir=c(0,90))
vario.res2dir = vario.calc(dat,lag=10,nlag=30,uc=1,dir=c(0,90))
plot(vario.raw2dir, main="Temperature (°C)", xlab = "Distance (km)", ylab = "Variogram")
plot(vario.res2dir,add=TRUE,lty=2)
legend("topleft", 
       legend=c("Easting","Northing","Easting-res","Northing-res" ), 
       lty=c(1,1, 2, 2), col=c("black", "red", "black", "red"))

```

## Global Trend

\footnotesize

Find the coefficients of the global regression (first order polynomial)

```{r Regression}
regression(dat, uc=1, flag.mode=2, save.coeff=TRUE)
```

## Modelling Variogram of Raw Data with Anisotropy

\footnotesize

```{r Model_for_Raw_Variogram, echo=FALSE}
# fitmod=model.auto(vario.raw2dir,struct=c(1,2,4))
fitmod=model.auto(vario.raw2dir,struct=c(1,2,4), draw=FALSE)
plot(fitmod, vario = vario.raw2dir, main="Temperature (°C)", 
     xlab = "Distance (km)", ylab = "Variogram", lw=2)
plot(vario.raw2dir,add=TRUE,lty=1)
legend("topleft", 
       legend=c("Easting","Northing"), 
       lty=c(1,1), col=c("black", "red"))
```


## Ordinary Kriging - Cross-Validation

\footnotesize

Perform Cross-validation (using Ordinary Kriging) and calculate the Mean Squared Error.

```{r Cross-validation_Ordinary_Kriging}
res.cv = xvalid(dat,fitmod,unique.neigh,radix="OK",
                modify.target = FALSE)
```

```{r Cross-validation_Ordinary_Kriging_stat, echo=FALSE}
tab <- as.matrix(db.stat.multi(res.cv, names = (c("OK.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption= "Cross-validation with Ordinary Kriging")
```

## Ordinary Kriging - Estimation

\footnotesize

```{r Estimation_Ordinary_Kriging}
target=kriging(dat,target,model=fitmod,unique.neigh,radix="OK")
```

```{r Estimation_Ordinary_Kriging_plot, echo=FALSE}
plot(target,name="OK*estim",title="Temperature - Ordinary Kriging",pos.legend=1)
plot(dat,add=TRUE)
```


## Ordinary Kriging - Statistics

\scriptsize

```{r Estimation_Ordinary_Kriging_stat, echo = FALSE}
tab <- as.matrix(db.stat.multi(target, names = (c("OK.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on the Ordinary Kriging")
```


## Fitting Variogram of Residuals

\footnotesize

```{r Fitting_Residuals}
fitmodUK = model.auto(vario.res2dir,struct=c(2),auth.aniso=FALSE)
```

Note that the residuals seem isotropic, hence use isotropic option for Fitting.

## Universal Kriging - Cross-Validation

\footnotesize

Perform Cross-validation (using Universal Kriging) and calculate the Mean Squared Error.

```{r Cross-validation_Universal_ Kriging}
res.cv = xvalid(res.cv,fitmodUK,unique.neigh,radix="UK",uc=1,
                modify.target=F)
```

```{r Cross-validation_Universal_ Kriging_stat, echo=FALSE}
tab <- as.matrix(db.stat.multi(res.cv, names = (c("UK.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption= "Cross-validation with Ordinary Kriging")
```


## Universal Kriging - Estimation

\scriptsize

```{r Universal_Kriging}
target=kriging(dat,target,model=fitmodUK,unique.neigh,radix="UK",uc=1)
```
 
```{r Universal_Kriging_stat, echo = FALSE}
plot(target,name="UK*estim",pos.legend=1,
     title="Temperature - Universal Kriging")
plot(dat,add=T)
```

## Universal Kriging - Statistics

\scriptsize

```{r Estimation_Universal_Kriging_stat, echo = FALSE}
tab <- as.matrix(db.stat.multi(target, names = (c("UK.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on the Universal Kriging")
```


## Comparing Ordinary and Universal Krigings

\scriptsize

```{r Comparing_Ordinary_and_Universal_Krigings, results="hide"}
dum = correlation(target,name1="OK*estim", name2="UK*estim",flag.diag=T,asp=0,
            title="",xlab="Ordinary Kriging",ylab="Universal Kriging")
```


## Comparing Temperature and Elevation

\scriptsize

```{r Correlation_of_Temperature_and_Elevation}
dum = correlation(dat,name1="Elevation",name2="Temperature",flag.regr = T,
            title="",ylab="Temperature",flag.aspoint=TRUE,size0=0.05)
```

## Comparing Temperature and Elevation

\scriptsize

Average of Elevations of all data and average of Elevations of meteorological stations

```{r Stat_on_temperature, echo=TRUE}
tab <- rbind(
  db.stat.multi(db.sel(dat), names = "Elevation"),
  db.stat.multi(db.sel(dat, !is.na(Temperature)), names = "Elevation"),
  db.stat.multi(db.sel(dat, !is.na(Temperature)), names = "Temperature"))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption= "Statistics on temperatures and elevations")
```

```{r selection_Temp_Ele, echo=TRUE}
dat = db.sel(dat)
dat = db.locate(dat, names = c("Temperature", "Elevation"), loctype = "z")
```

## Bivariate Modelling

\scriptsize

```{r Bivariate_Modelling, fig.align='center'}
varioexp2var = vario.calc(dat,lag=10,nlag=30,dir=c(0,90))
fitmod2var=model.auto(varioexp2var,struct=c(1,2,4))
```

## Cokriging with elevation - Cross-Validation

\scriptsize

Most of the processes are more time-consuming in Unique Neighborhood.
We create a small neighborhood for demonstration.

```{r CK_neigh}
moving.neigh = neigh.create(ndim=2, type=2, radius = 1000, nmaxi = 10)
```

Perform Cross-validation (Bivariate Model) and calculate the Mean Squared Error.

```{r CK_Cross-validation, echo=TRUE}
res.cv = db.locate(res.cv,c("Temperature","Elevation"),"z")
res.cv = xvalid(res.cv,fitmod2var,moving.neigh,radix="COK",
                modify.target=F)
```


```{r CK_Cross-validation_stat, echo=FALSE}
tab <- as.matrix(db.stat.multi(res.cv, names = (c("COK.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption="Cross-validation with Cokriging")
```

## Cokriging with elevation - Estimate

\scriptsize

```{r Cokriging}
dat = db.locate(dat,c("Temperature","Elevation"),"z")
target=kriging(dat,target,model=fitmod2var, unique.neigh,radix="COK")
```

```{r Cokriging_plot, echo=FALSE}
plot(target,name="COK.T*estim",pos.legend=1,
     title="Temperature - Cokriging Kriging (with elevation)")
plot(dat,add=T)
```


## Cokriging with elevation - Statistics

\scriptsize

```{r Estimation_cokriging_stat, echo = FALSE}
tab <- as.matrix(db.stat.multi(target, names = (c("COK.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on the Cokriging")
```

## Comparing Kriging and CoKriging

\scriptsize

```{r Comparing_Kriging_and_CoKriging, results="hide", echo=FALSE}
dum = correlation(target, name1="OK.T*estim",name2="COK.T*estim",flag.diag=TRUE,
            xlab="Ordinary Kriging",ylab="Ordinary CoKriging",title="")
tab <- as.matrix(db.stat.multi(target, 
      names = (c("COK.T*estim", "OK.T*estim"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption="Statistics on Kriging and CoKriging estimates")
```

Note that CoKriging produces estimates which are mostly larger than Kriging estimates.

## Kriging the residuals

\scriptsize

$$
Z_2(s)=b + a Z_1(s) + R(s)
$$
```{r residual_init, echo=TRUE}
dat <- db.sel(dat, !is.na(Temperature))
temperature <- db.extract(dat,"Temperature", flag.compress = TRUE)
elevation   <- db.extract(dat,"Elevation", flag.compress = TRUE)
a <- cov(temperature, elevation) / cov(elevation,elevation)
b <- mean(temperature) - a * mean(elevation)
dat <- db.add(dat, Res = Temperature - b - a * Elevation)
```


```{r residual_init_stat, echo=FALSE}
tab <- as.matrix(db.stat.multi(dat, names = c("Res")))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption="Statistics on the residual")
```


## Kriging the residuals - Correlation

\scriptsize

```{r residual_correlation, echo=TRUE}
dum = correlation(dat, name1 = "Elevation", name2 = "Res", 
            flag.aspoint = TRUE, flag.regr = TRUE)
```

## Kriging the residuals - Variogram of the residual

\scriptsize

```{r residual_variogram, echo=TRUE}
dat = db.locate(dat, names = c("Res"), loctype = "z", flag.locnew = TRUE)
varioexpR = vario.calc(dat,lag=10,nlag=30,dir=c(0,90))
fitmodR=model.auto(varioexpR,struct=c(1,2,13), auth.aniso = TRUE)
```

## Kriging the residuals - Variogram of the residual

\scriptsize

```{r residual_kriging, echo=TRUE}
dat = db.locate(dat,c("Res"),"z", flag.locnew = TRUE)
target=kriging(dat,target,model=fitmodR, unique.neigh,radix="OK")
```

```{r residual_kriging_plot, echo=FALSE}
plot(target,name="OK.R*estim",pos.legend=1,
     title="Kriging Residuals")
plot(db.locate(dat,"Res","z"),add=T)
```

## Kriging the residuals - Computing the estimate

\scriptsize

$$
Z_2^{\star} = b + a Z_1(s) + R(s)^{OK}
$$

```{r residual_estimate, echo=TRUE}
target <- db.add(target, KR.Temperature.estim = b + a * Elevation + OK.Res.estim)
```

```{r residual_estimate_plot, echo=FALSE}
plot(target,name="KR.T*estim",pos.legend=1,
     title="Temperature - Kriging from Residuals")
plot(db.locate(dat, "Temperature", "z", flag.locnew = TRUE),add=T)
```

## Kriging the residuals - Correlation

\scriptsize

```{r residual_estimate_corr, echo=FALSE}
dum = correlation(target, name1="OK.T*estim",name2="KR.T*estim",flag.diag=TRUE,
            xlab="Ordinary Kriging",ylab="Kriging with residual",title="")
```

## Kriging the residuals - Correlation

\scriptsize

```{r residual_estimate_stat, echo=FALSE}
tab <- as.matrix(db.stat.multi(target, 
      names = (c("OK.R*estim", "KR.T*estim", "OK.T*estim"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption="Statistics on OK and Kriging of the resisual")
```

## Using Elevation Map as External Drift

\scriptsize

Preparing the data bases

```{r Preparation_for_KDE}
dat = db.locate(dat,"Temperature","z", flag.locnew = TRUE)
dat = db.locate(dat,"Elevation","f")
target = db.locate(target,"Elevation","f")
```

## Linear Regression of Temperature knowing Elevation on Data

\scriptsize

We perform the Linear Regression of Temperature knowing Elevation (specified as External Drift)

```{r Regression_Temperature_Elevation}
regression(dat, save.coeff=TRUE, flag.mode=1)
```

## Experimental Variogram of Residuals (External Drift)

\scriptsize

```{r Fitting_Variogram_of_Residuals_with_ED}
varioKED = vario.calc(dat,lag=10,nlag=30,dir=c(0,90),uc=c("1","f1"))
```

```{r Fitting_Variogram_of_Residuals_with_ED_fig, echo=FALSE}
plot(vario.raw2dir, main="Temperature (°C)", 
     xlab = "Distance (km)", ylab = "Variogram")
plot(varioKED,add=TRUE,lty=2)
legend("topleft", 
       legend=c("Easting","Northing","Easting-KED","Northing-KED" ), 
       lty=c(1,1, 2, 2), col=c("black", "red", "black", "red"))
```

## Model of Residuals (External Drift)

\scriptsize

```{r Model_of_Residuals_with_External_Drift}
modelKED = model.auto(varioKED,struct=c(1,4,3))
legend("topleft", legend=c("Easting","Northing"), lty=c(1,1), 
       col=c("black", "red"))
```

## Kriging with External Drift - Cross-Validation

\scriptsize

Perform Cross-validation (External Drift) and calculate the Mean Squared Error.

```{r Cross-Validation_External_Drift}
res.cv=db.locate(res.cv,"Elevation","f")
res.cv=xvalid(db = res.cv, model = modelKED,neigh = unique.neigh, 
              uc=c("1","f1"),radix="KED")
```

```{r Cross-Validation_External_Drift_stat, echo=FALSE}
tab <- as.matrix(db.stat.multi(res.cv, names = (c("KED.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption="Cross-validation with Kriging with External Drift")
```

## Kriging with External Drift - Estimate

\scriptsize

```{r Kriging_with_External_Drift}
target = kriging(dat,target,modelKED,unique.neigh,uc=c("1","f1"), radix="KED")
```

```{r Kriging_with_External_Drift_plot, echo=FALSE}
plot(target,name="KED*estim",pos.legend=1,
     title="Temperature - Kriging with External Drift")
plot(dat,add=TRUE)
```

## Kriging with External Drift - Statistics

\scriptsize

```{r Estimation_ked_stat, echo = FALSE}
tab <- as.matrix(db.stat.multi(target, names = (c("KED.T*"))))
tab <- as.matrix(round(tab, 3))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
knitr::kable(tab, caption = "Statistics on the Kriging with External Drift")
```

## Comparing OK with KED

\scriptsize

```{r Comparing_Ordinary_Kriging_with_Kriging_with_External_Drift, results="hide", echo=FALSE}
dum = correlation(target,name1="OK.T*estim",name2="KED.T*estim",
            flag.diag=T,flag.same=T,flag.iso=T,
            title="",xlab="Ordinary Kriging",ylab="Kriging with External Drift")
```

Note that negative Estimates are present when using External Drift.

## Summary of Cross-validation scores

\scriptsize

Comparing the cross-validation Mean Squared Errors

```{r Comparing_Cross-validations}
tab <- as.matrix(round(db.stat.multi(res.cv, names = (c("*.T*esterr"))), 
                       digits = 4))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
rownames(tab) <- c("Ordinary Kriging", "Universal Kriging", 
                   "Co-kriging","Kriging with External Drift")
knitr::kable(tab,
  caption = "Statistics of the cross-validation error for the temperature"
  )
```

## Statistics on the estimates

\scriptsize

Computing the statistics on the various estimates.

```{r Comparing_estimates_estim, echo = TRUE}
tab <- as.matrix(round(db.stat.multi(target, names = (c("*.T*estim"))), 
                       digits = 4))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
rownames(tab) <- c("Ordinary Kriging", "Universal Kriging", 
                   "Co-kriging","Kriging residual", "Kriging with External Drift")
knitr::kable(tab,
  caption = "Statistics of the estimates for the temperature"
  )
```

## Statistics on the estimates

\scriptsize

Computing the statistics on the various estimates.

```{r Comparing_estimates_stdev, echo = TRUE}
tab <- as.matrix(round(db.stat.multi(target, names = (c("*.T*stdev"))), 
                       digits = 4))
colnames(tab) <- c("Nb", "Min", "Max", "Mean", "STD")
rownames(tab) <- c("Ordinary Kriging", "Universal Kriging", 
                   "Co-kriging","Kriging with External Drift")
knitr::kable(tab,
  caption = "Statistics of the kriging Std. for the temperature"
  )
```