trendSurface04.rmd

---
##-- arkriger: August 2023
title: "Trend Surface Analysis with R (part IV)"
output:
  html_notebook: default
  html_document:
    df_print: paged
---
<style> 
p.comment {
background-color: #DBDBDB;
padding: 10px;
border: 1px solid black;
margin-left: 25px;
border-radius: 5px;
font-style: italic;
}

</style>

We end our journey through *higher order* **_Trend Surface analysis with R_** with this Notebook. Here we cover a lot of ground.

1. Generalized Least Squares (GLS)
	a) Variogram Modeling and
	b) Ordinary Kriging
2. GLS-Regression Kriging
3. Universal Kriging

<div class="alert alert-danger">
  <strong>REQUIRED!</strong> 
  
You are required to insert your outputs and any comment into this document. The document you submit should therefore contain the existing text in addition to:

 - Plots and other outputs from executing the code chunks
 - Discussion of your plots and other outputs as well as conclusions reached.
 - This should also include any hypotheses and assumptions made as well as factors that may affect your conclusions.
</div>

At the close of the previous Notebook we talked about **_spatial autocorrelation_**. Patterns and particularly residuals (errors) have a spatial structure / are connected.

What we mean by this is; that even though an estimate _(prediction)_ is unbiased the **structure of the underlying data can be**. A general error term applied to our entire study area might not be ideal.

To overcome the challenge; we use Generalized Least Squares and introduce a **covariance structure between residuals**. Our `linear model`, from [trendSurface02](https://github.com/AdrianKriger/r-trend-surfaces/blob/main/trendSurface02.rmd) then becomes:

$$
    𝜷_{gls} = (X^TC^{-1}X)^{-1}X^TC^{-1} . y
$$

We see our new term includes $C$; the covariance matrix of the spatially correlated residuals.

The obvious question is. _'OK; so what are the correlated residuals?'_
And the answer is: _'We dont know'_.

What we do is fit the covariance structure (based on imprecise observations) at the same time we calculate the trend surface coefficients.

```{r install }
options(prompt="> ", continue="+ ", digits=3, width=70,  show.signif.stars=F, repr.plot.width=7, repr.plot.height=7)
rm(list=ls())

# Install necessary packages: You only need to run this part once
##- install.packages(c("sf", "gstat", "ggplot2", "gridExtra","units", "terra","mgcv","fields", "nlme"))

library(sf) # 'simple features' representations of spatial objects
library(gstat) # geostatistics
library(ggplot2) # geostatistics
library(gridExtra)
library(units) # units of measure
library(terra) # gridded data structures ("rasters")
library(mgcv)
library(fields)
library(nlme)
#- 3d surface
library(plotly)
library(reshape2)
```

We've already covered creating a 500m grid and 2nd-order prediction and interpolation with previous Notebooks and move through the introduction quickly.

```{r read-dataset }
#-- read
file = 'cptFlatsAquifer_watertable-Aug2023.txt'
```

```{r import-grid-2nd-ols-gam }
#-- import
cfaq <- read.csv2(file, header = 1, sep = ';', dec = ',')
#- set crs
cfaq.sf <- st_as_sf(cfaq, coords=c("long", "lat"), crs = 4326) #wgs84
#- transform to local crs
cfaq.sf <- st_transform(cfaq.sf, crs = 32734) #utm 34s
cfaq$X <- st_coordinates(cfaq.sf)[, "X"]
cfaq$Y <- st_coordinates(cfaq.sf)[, "Y"]

model.ts2 <- lm(waterLevel ~ X + Y + I(X^2) + I(Y^2) + I(X*Y), data=drop_units(cfaq))

#--  create grid
(n.col <- length(seq.e <- seq(min.x <- floor(min(cfaq$X)/1000)*1000,
                              max.x <- ceiling(max(cfaq$X)/1000)*1000, by=1000)))
(n.row <- length(seq.n <- seq(min.y <- floor(min(cfaq$Y)/1000)*1000,
                              max.y <- ceiling(max(cfaq$Y)/1000)*1000, by=1000)))

#we want a XXXXm grid
grid <- rast(nrows = n.row, ncols = n.col,
             xmin=min.x, xmax=max.x,
             ymin=min.y, ymax=max.y, crs = st_crs(cfaq.sf)$proj4string,
             resolution = 500, names="waterLevel")

values(grid) <- NA_real_

grid.df <- as.data.frame(grid, xy = TRUE, na.rm = FALSE)
names(grid.df)[1:2] <- c("X", "Y") # match the names of the point dataset
summary(grid.df)

#--
pred.ts2 <- predict.lm(model.ts2,
                       newdata = grid.df,
                       interval = "prediction", level = 0.95)

#-- add the three prediction fields (fit, lower, upper) to the data
grid.df[, 3:5] <- pred.ts2
names(grid.df)[3:5] <- c("ts2.fit", "ts2.lwr", "ts2.upr")
values(grid) <- pred.ts2[,"fit"]

#-- gam
model.gam <- gam(waterLevel ~ s(X, Y, k = 29), data=drop_units(cfaq))
#-- predict aquifer elevation, standard error of prediction using the fitted GAM
tmp <- predict.gam(object=model.gam,
                   newdata=grid.df,
                   se.fit=TRUE)

grid.df$pred.gam <- as.numeric(tmp$fit)
grid.df$pred.gam.se <- as.numeric(tmp$se.fit)

grid.gam <- grid
values(grid.gam) <- grid.df$pred.gam
```

```{r look }
#-- look
head(cfaq ,3)
```

To compute the trend and covariance at the same time we use the `gls` function in the the `nlme` package.

```{r gls-model-class }
#- gls with coefficient
model.ts2.gls <- gls(
  model = waterLevel ~ Y + X + I(Y^2) + I(X^2) + I(X * Y),
  data = drop_units(cfaq),
  method = "ML",
  correlation = corExp(form = ~X + Y,
                       nugget = FALSE,
                       value = 10000) # initial value of the range parameter
)
#model.ts2.gls
class(model.ts2.gls)
```

```{r gls-summary }
summary(model.ts2.gls)
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 1.**  What is the range of spatial correlation of the exponential model, as estimated by GLS?

<p class="comment">
[ Answer 1. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

**Look at how the GLS and OLS coefficients differ.**

We work with the **Mixed GAM Computation Vehicle** `mgcv` library default smoothing setting and explore **thin plate splines**, as an alternate smoothing function later.

```{r summarize-gls-uncertainty }
##-- absolute values
#- generic coef method extracts coefficients from model objects
coef(model.ts2.gls) - coef(model.ts2)
```

```{r summarize-gls-uncertainty-percent }
#- percentage
round(100*(coef(model.ts2.gls) - coef(model.ts2))/coef(model.ts2),1)
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 2.** Why are the GLS coefficients different than the OLS coefficients?

<p class="comment">
[ Answer 2. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

**Lets go deeper.**

Display the $90\%$ confidence intervals of the GLS model, calcualte the residuals of this surface and compare those to the OLS.

```{r gls-intervals }
#- generic intervals method has a specific method for a fitted GLS model
intervals(model.ts2.gls, level=0.90)
```

```{r gls-2nd-residuals }
summary(res.ts2.gls <- residuals(model.ts2.gls))
```


```{r gls-2nd-residuals-summary }
res.ts2 <- set_units(residuals(model.ts2), m)
summary(res.ts2)
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 3.** What are the main differences between these sets of residuals? Which surface, in this case, most closely fits the points? Why?

<p class="comment">
[ Answer 3. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

```{r predict-gls-summary }
#- predict.gls function (of the nlme package) specific method for a fitted GLS model
pred.ts2.gls <- predict(model.ts2.gls, newdata=grid.df)
summary(pred.ts2.gls)
```

```{r spat-rast }
#- SpatRast of predictions
grid.gls <- grid
values(grid.gls) <- pred.ts2.gls
summary(values(grid.gls))
```

```{r plot-gls }
#- plot
plot(grid.gls, col = rainbow(100), main = "GLS 2nd order predicted surface")
points(st_coordinates(cfaq.sf)[,2] ~ st_coordinates(cfaq.sf)[,1],
       pch=21,
       col=ifelse((res.ts2.gls < 0), "red", "green"),
       cex=2*abs(res.ts2.gls)/max(abs(res.ts2.gls)),
       bg = adjustcolor("black", alpha.f = 0.1),
       lwd = 0.5)
```

**Notice how closely this prediction resembles the OLS 2nd-order trend surface!**

Lets see how different?

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 4.** What are the main differences between these sets of residuals? Which surface, in this case, most closely fits the points? Why?

<p class="comment">
[ Answer 4. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>


```{r plot-diff-ols-gls }
#- difference between the OLS and GLS trendsurfaces
grid.gls.ols.diff <- (grid.gls - grid)
plot(grid.gls.ols.diff, col = topo.colors(64),
     main = "GLS - OLS 2nd order trend surfaces, m")
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 5.** Where are the largest differences between the OLS and GLS trend surfaces? Explain why.

<p class="comment">
[ Answer 5. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

Lets look at the residuals

```{r residuals-gls-summary }
#- residuals from the GLS trend surface and as a postplot
summary(res.ts2.gls)
```

```{r gls-fit-residuals }
#- plot  ??local spatial correlation??
plot(cfaq$Y ~ cfaq$X, cex=3*abs(res.ts2.gls)/max(abs(res.ts2.gls)),
     col=ifelse(res.ts2.gls > 0, "green", "red"),
     xlab="X", ylab="Y",
     main="Residuals from 2nd-order trend, GLS fit",
     sub="Positive: green; negative: red", asp=1)
```

<div class="alert alert-info">
  <strong>Look at how correlated the residuals are.</strong> Nearby data points have similar values
</div>

We now account for the spatial correlation of the residuals

```{r empirical-variogram }
#-  empirical variogram model residuals from the GLS trend surface model

#- extract the residuals into the point observations object
cfaq.sf$res.ts2.gls <- residuals(model.ts2.gls)
vr.gls <- variogram(res.ts2.gls ~ 1, loc=cfaq.sf, cutoff = 65000)
#- plot empirical variogram
plot(vr.gls, plot.numbers=T,
     main="Residuals from second-order GLS trend", cex=1,
     xlab = "separation [m]", ylab = "semivariance [m^2]")
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 6.** What are the approximate variogram parameters?

<p class="comment">
[ Answer 6. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

### 1. a) Variogram Modeling

In order to solve any Kriging formula we need to compute the semivariance at _any_ seperation distance. We want to know the values that will remove the _'overlarge influence'_ nearby datapoints have on a prediction.

In order to do so we fit a variogram function to the empirical variogram. The function represent the structure of spatial autocorrelation. There are many variogram functions: _'Exponential, Spherical, Gauss, Matern'_.  

We will not go into these here but remember to `effective-range / 3` when chooing an _'Exponential'_ model.

```{r variogram-model }
#-  specify variogram model and parameters
vr.gls.m <- vgm(psill=250, model="Mat", range=20000)#, nugget=25)
#vr.gls.m <- vgm(psill=40, model="Exp", range=20000/3, nugget=0)
(vr.gls.m.f <- fit.variogram(vr.gls, vr.gls.m))
```

```{r plot-variogram-model }
#- plot
plot(vr.gls, model=vr.gls.m.f, plot.numbers=T,
     xlab = "separation [m]", ylab = "semivariance [m^2]")
```

```{r print-model }
print(vr.gls.m.f)
```

```{r variogram-intervals }
intervals(model.ts2.gls)$corStruct[2]
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 7.** Does the range parameter of this fitted model agree with the estimate from the GLS fit?

<p class="comment">
[ Answer 7. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

### 1. b) The Ordinary Kriging system

Kriging is a form of linear prediction of a value at an unknown point as a weighted sum of values at known points.

This makes sense. We know the structure (spatial correlation) of the residuals, their values and their location. And we use this knowledge to refine our GLS prediction.  

But what really makes Kriging special _(and computationally expensive)_ is the system ensures that each prediction has the least possible prediction variance. The _**uncertainty**_ of the predictions are minimized.

```{r ok}
#-- ok
grid.sf <- st_as_sf(grid.df, coords = c("X", "Y"))
st_crs(grid.sf) <- st_crs(grid)
kr <- krige(res.ts2.gls ~ 1,
            loc = cfaq.sf,
            newdata = grid.sf,
            model=vr.gls.m.f)
```

```{r ok-summary }
#- summary
summary(kr)
```

```{r ok-class }
class(kr)
```

```{r ok-residuals }
#- ok residuals
plot(kr["var1.pred"], pch=15, nbreaks=24, main="Residuals from GLS trend, m")
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 8.** Which areas were most changed by interpolating the residuals? Why?

<p class="comment">
[ Answer 8. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

```{r ok-std-dev-summary }
#- prediction standard deviations and plot
kr$var1.sd <- sqrt(kr$var1.var)
summary(kr)
```

```{r plot-ok-std-dev-residuals }
#-
plot(kr["var1.sd"], pch=15, nbreaks=24, pal = heat.colors,
     main="Standard errors of residuals from GLS trend, m")
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 9.** Which areas have the most and least uncertainty? Why?

<p class="comment">
[ Answer 10. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>`

## 2. Regression Kriging

We now have the ingredients to predict a very good trend surface. The global function and the local deviations from it.

$$
    Z(s) = Z*(s) + ϵ(s) + ϵ'(s)
$$

We modeled a trend surface, $Z*(s) + ϵ(s)$, with an OLS polynomial and later derived a model for the trend's spatial structure, $ϵ(s) + ϵ'(s)$, by Ordinary Kriging. These were modeled together with GLS.

**We now add the OK prediction of the GLS residuals to the GLS surface. This is called Generalized Least Squares Trend-Regression Kriging (GLS-RK)**

```{r gls-rk-summary }
#-- gls-regression kriging
grid.df$kr <- kr$var1.pred
grid.df$pred.ts2.gls <- pred.ts2.gls
grid.df$rkgls <- grid.df$pred.ts2.gls + grid.df$kr
summary(grid.df)
```

```{r gls-rk-to-grid }
grid.rkgls <- grid
values(grid.rkgls) <- grid.df$rkgls
```

```{r plot gls-rk }
#- plot
plot(grid.rkgls, col = rainbow(100),
     main="GLS-RK prediction, aquifer elevation, m.a.s.l.")
```
Compare the GLS-RK with the GAM

```{r diff-gls-rk-minus-gam }
summary(grid.rkgls.gam <- (grid.rkgls - grid.gam))
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 10.** Where are the largest differences between these two trend surface predictions? Explain why, considering how the two surfaces are computed

<p class="comment">
[ Answer 10. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>

## 3. Universal Kriging

We've come to the last leg of our journey.

Universal Kriging (UK) is an extension of OK. UK does not follow a step-wise process of adding an OK prediction of GLS residuals to a GLS surface (like GLS-RK does).

UK bypasses the GLS and predicts a surface from a 2nd-order OLS. In other words the spatial structure of the residuals (the variogram) is not modeled from the GLS residuals but the 2nd-order model directly.

**We start with a clean slate**

```{r look-see }
#- what do we have thusfar?
names(grid.sf)
```

```{r str-grid }
str(grid.sf$geometry)
```

```{r cfaq-sf-names}
#-- 
names(cfaq.sf)
```

```{r str-cfaq-geom}
str(cfaq.sf$geometry)
```

**With UK we must include the coordinates (the predictor variable) in both the observation and prediction grid**

```{r str-grid-sf-coords }
#-- 
str(st_coordinates(cfaq.sf))
```

```{r grid-to-sf}
#- variables to grid
grid.sf$X <- st_coordinates(grid.sf)[ , "X"]
grid.sf$Y <- st_coordinates(grid.sf)[ , "Y"]
names(grid.sf)
```

```{r variables-to-cfaq }
#- variables to cfaq
cfaq.sf$X <- st_coordinates(cfaq.sf)[ , "X"]
cfaq.sf$Y <- st_coordinates(cfaq.sf)[ , "Y"]
names(cfaq.sf)
```
Variogram from the OLS
```{r ols-empirical-variogram-plot}
#- empirical variogram
vr <- variogram(waterLevel ~ X + Y + I(X^2) + I(Y^2) + I(X*Y), locations = cfaq.sf, 
                cutoff = 65000)
#- plot
plot(vr, plot.numbers = TRUE, main = "Residuals from 2nd-order OLS trend surface", xlab = "separation (m)", ylab = "semivariance (m^2)")
```

```{r ols-variogram-model}
#- variogram model 
#(vr.m.f <- fit.variogram(vr, vgm(35, "Exp", 22000/3, 0)))
vr.m.f <- fit.variogram(vr, vgm(psill=250, model="Mat", range=20000))#, nugget=25)
```

```{r plot-ols-variogram-model}
#- plot
plot(vr, plot.numbers=TRUE, xlab="separation (km)", ylab="semivariance (m^2)", model=vr.m.f, main="Fitted variogram model, residuals from 2nd-order OLS trend surface")
```


```{r uk }
#- uk
k.uk <- krige(waterLevel ~ X + Y + I(X^2) + I(Y^2) + I(X*Y), locations = cfaq.sf, 
              newdata = grid.sf, model=vr.m.f)
```

```{r plot-uk }
#- plot uk predictors
plot(k.uk["var1.pred"], pch=15, col = rainbow(100), nbreaks=24, main="UK predictions, m")
#plot(k.uk["var1.pred"], pch=15, col = rainbow(100), main="UK predictions, m")

```

```{r uk-std-dev-summary }
#- uk-prediction std-dev.
k.uk$var1.sd <- sqrt(k.uk$var1.var)
summary(k.uk)
```

```{r plot-uk-std-dev }
#- plot std-dev
plot(k.uk["var1.sd"], pch=15, nbreaks=24, pal = heat.colors, main="Standard errors of UK predictions, m")
```

Compare the UK standard deviations to the GLS-Regression Kriging

```{r uk-std }
#- summary uk
summary(k.uk$var1.sd)
```

```{r gls-ok-std }
#- summary uk
summary(kr$var1.sd)
```

**Note:** the UK should be larger because it includes the uncertainty of the trend surface.

Lets illustrate the UK and GLS-RK differences as a histogram

```{r diff-uk-minus-gls-rk-summary }
#- diff. uk and gls-rk summary
grid.uk <- grid
values(grid.uk) <- k.uk$var1.pred
summary(grid.diff.uk.rkgls <- (grid.uk - grid.rkgls))
```

```{r diff-uk-minus-gls-rk-hist }
#- diff. uk and gls-rk hist.
hist(grid.diff.uk.rkgls, main = "UK - GLS-RK prediction differences", freq = FALSE, xlab = "difference, UK - GLS-RK")
```

```{r plot-diff-uk-minus-gls-rk}
#- plot uk minus gls-rk diff surface
plot(grid.diff.uk.rkgls, sub="UK - GLS-RK predictions", main="difference, m", xlab="East", ylab="North", col = topo.colors(64))
```

<div class="alert alert-warning"> <strong>QUESTION!</strong> </div>

- **Question 11.** How large are differences between the UK and GLS-RK trend surface predictions? Where are the largest differences? Explain why there is a difference.

<p class="comment">
[ Answer 10. click in this cell and type your answer here. your answer must be between the outer [] brackets ]
</p>


<div class="alert alert-info">
  <strong>DISCUSSION</strong> </div>

- **Thoughts 1**. In this study area, which of the prediction methods would you recommend, and why?

<p class="comment">
[ Thoughts 1. click in this cell and type your comments here. your answer must be between the outer [] brackets ]
</p>

- **Thoughts 2**. In this study area, which of the prediction methods would you recommend, and why? For each method introduced, in what situations would you prefer it to the other methods?

<p class="comment">
[ Thoughts 2. click in this cell and type your comments here. your answer must be between the outer [] brackets ]
</p>

<div class="alert alert-success">
  <strong>This is an excellent opportunity to create a _fancy interactive_ 3D surface plot</strong> </div>


```{r plot-fancy }
#grid.df
grid.df_matrix <- acast(grid.df, X~Y, value.var="rkgls")
# More rows than column, and the x and y should have the same size
x <- as.numeric(rownames(grid.df_matrix))
y <- as.numeric(colnames(grid.df_matrix))

#- change z-axis labels (ticks)
axz <- list(
  nticks = 2,
  range = c(-20,120),
  title='Elevation (m)'
)

fig <- plot_ly()#showscale = FALSE)
fig <- fig %>% add_surface(z = grid.df_matrix, x = x, y = y, colors = rainbow(100), 
                           type = "surface") %>%
  colorbar(title = "Water Level") %>%
  layout(
    scene=list(
    xaxis=list(title='x (m)'),
    yaxis=list(title='y (m)'),
    zaxis=axz, #list(title='Elevation (m)'),
    aspectmode = 'manual', aspectratio = list(x=1, y=1, z=0.1)))  # set range of z-axis

htmlwidgets::saveWidget(as_widget(fig), "plotly.html")

fig
```