## ----setup, include=FALSE, cache=FALSE------------------------------
# set global chunk options
## this code is only required to set up knitr for reproducible documentation
## it has nothing to do with the exercise itself, please ignore!
require("knitr")
opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'), 
              out.format='latex', use.highlight=TRUE)
opts_chunk$set(fig.path='kgraph/ex_TrendSurface_sf', fig.align='center', 
               fig.show='hold', prompt=FALSE,
               fig.width=6, fig.height=6, out.width='.7\\linewidth', size='scriptsize')
## options to be read by formatR
options(replace.assign=TRUE,width=70)


## ----echo=FALSE, results='hide'-------------------------------------
## make sure we have a clean environment when knitting
rm(list=ls())


## ----child-data, child='ex_TrendSurface_ex1-knitr_sf.Rnw'-----------

## ----ts1-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('ex_TrendSurface_sf.Rnw')


## ----load-packages--------------------------------------------------
require(sf)        # 'simple features' representations of spatial objects
require(gstat)     # geostatistics
require(ggplot2)   # Grammer of Graphics plots
require(gridExtra) # arrange multiple ggplot graphics on one figure
require(units)     # units of measure
require(terra)     # gridded data structures ("rasters")
require(mgcv)      # for Generalized Additive Models
require(fields)    # NCAR etc. approach to surfaces
require(nlme)      # Linear and Nonlinear Mixed Effects Models


## ----set-options----------------------------------------------------
options(show.signif.stars=FALSE)


## ----read-aq, eval=F------------------------------------------------
## aq <- read.table("AQUIFER.TXT", skip=1)


## ----Read-aq-2, echo=F----------------------------------------------
aq <- read.table("./ds/kansas/AQUIFER.TXT", skip=1)


## ----str-aq---------------------------------------------------------
str(aq)
names(aq) <- c("E", "N", "z")
str(aq)


## ----ft-to-m--------------------------------------------------------
aq$z <- set_units(aq$z, feet)
aq$zm <- set_units(aq$z, m)
summary(aq)


## ----convert-spdf---------------------------------------------------
aq.sf <- st_as_sf(aq,  coords=c("E","N"))
str(aq.sf)


## ----st-coords------------------------------------------------------
str(st_coordinates(aq.sf))
#
summary(st_coordinates(aq.sf))


## ----show-crs-na----------------------------------------------------
st_crs(aq.sf)


## ----set-crs--------------------------------------------------------
st_crs(aq.sf) <- 26914 # EPSG code
print(st_crs(aq.sf))


## ----save-aq--------------------------------------------------------
save(aq, aq.sf, file="aquifer.rda")


## -------------------------------------------------------------------
dim(aq.sf)
summary(aq.sf)


## ----compute_bbox---------------------------------------------------
st_bbox(aq.sf)


## ----compute_range--------------------------------------------------
range(aq$E)/1000
range(aq$N)/1000
print(diff(range(aq$E)) *diff(range(aq$N))/1000^2)


## -------------------------------------------------------------------
range(aq.sf$zm); diff(range(aq.sf$zm))


## ----aquifer-height-text, fig.height=12, fig.width=12---------------
plot(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], 
     pch=20, cex=0.2, col="blue", asp=1,
     xlab="UTM 14N E", ylab="UTM 14N N")
grid()
text(st_coordinates(aq.sf)[,1], st_coordinates(aq.sf)[,2], 
     round(aq$zm), adj=c(0.5,0.5))
# text(aq$E, aq$N, round(aq$zm), adj=c(0.5,0.5))
title("Elevation of aquifer, m")


## ----aquifer-height-size, tidy=FALSE--------------------------------
plot(aq$N ~ aq$E, 
     cex=1.8*aq$zm/max(aq$zm),
     col="blue", bg="red", pch=21, asp=1,
     xlab="UTM14N_E", ylab="UTM14N_N") 
grid()
title("Elevation of aquifer, m.a.s.l.")


## ----aquifer-height-color-size, tidy=FALSE--------------------------
plot(aq$N ~ aq$E, pch=21,
     xlab="UTM14N_E", ylab="UTM14N_N",
     bg=sp::bpy.colors(length(aq$zm))[rank(aq$zm)],
     cex=1.8*aq$zm/max(aq$zm), asp=1) 
grid()
title("Elevation of aquifer, m.a.s.l.")


## ----model-matrix---------------------------------------------------
head(model.matrix(~E + N, data=aq))


## ----fit-1st-order--------------------------------------------------
model.ts1 <- lm(zm ~ E + N, data=drop_units(aq))
summary(model.ts1)
model.ts1.sf <- lm(zm ~ st_coordinates(aq.sf), data=drop_units(aq.sf))
summary(model.ts1.sf)


## ----hist-res-ts1, tidy=FALSE, fig.width=8, fig.height=6------------
res.ts1 <- set_units(residuals(model.ts1), m); summary(res.ts1)
hist(res.ts1, breaks=16, 
     main="Residuals from 1st-order trend",
     xlab="residual elevation (m)")
rug(res.ts1)
range(res.ts1)
max(abs(res.ts1))/median(aq$zm)*100


## ----ts1-resid, tidy=FALSE, fig.width=10, fig.height=5, out.width='\\textwidth'----
par(mfrow=c(1,2))
plot(model.ts1, which=1:2)
par(mfrow=c(1,1))


## ----post-res-ts1, tidy=FALSE---------------------------------------
plot(aq$N ~ aq$E, cex=3*abs(res.ts1)/max(abs(res.ts1)),
     col=ifelse(res.ts1 > set_units(0, m), "green", "red"),
     xlab="E", ylab="N", 
     main="Residuals from 1st-order trend",
     sub="Positive: green; negative: red", asp=1)
grid()


## ----show-design-2--------------------------------------------------
head(model.matrix(~ E + N + I(E^2) + I(N^2) + I(E*N),
                data=aq))


## ----fit.ts2--------------------------------------------------------
model.ts2 <- lm(zm ~ E + N + I(N^2) + I(E^2) + I(E*N),
                data=drop_units(aq))
summary(model.ts2)


## ----anova-ts1-ts2--------------------------------------------------
anova(model.ts1, model.ts2)


## ----hist-res-ts2, fig.width=8, fig.height=6------------------------
res.ts2 <- set_units(residuals(model.ts2), m)
summary(res.ts2)
hist(res.ts2, breaks=16, 
     main="Residuals from 2nd-order trend",
     xlab="residual elevation (m)")
rug(res.ts2)
max(abs(res.ts2))/median(aq$zm)


## ----ts2-resid, tidy=FALSE, fig.width=10, fig.height=5, out.width='\\textwidth'----
par(mfrow=c(1,2))
plot(model.ts2, which=1:2)
par(mfrow=c(1,1))


## ----find-largest-resid-ts2-2---------------------------------------
summary(sres.ts2 <- rstandard(model.ts2))
(ix <- which(sres.ts2 < -2))
(cbind(actual=aq[ix, "zm"], fitted=fitted(model.ts2)[ix],
            residual=res.ts2[ix],
            std.res <- sres.ts2[ix]))


## ----post-res-ts2, tidy=FALSE, fig.width=8, fig.height=8, out.width='.9\\linewidth'----
plot(aq$N ~ aq$E, cex=3*abs(res.ts2)/max(abs(res.ts2)),
     col=ifelse(res.ts2 > set_units(0, m), "green", "red"),
     xlab="E", ylab="N", asp=1,
     main="Residuals from 2nd-order trend",
     sub="Positive: green; negative: red; black dots: severe over-predictions")
points(aq[ix, "N"] ~ aq[ix, "E"], pch=20)
text(aq[ix, "E"], aq[ix, "N"], round(sres.ts2[ix], 2), pos=4)
grid()


## ----find-bbox------------------------------------------------------
range(aq$E); range(aq$N)
range(st_coordinates(aq.sf)[,"X"]); range(st_coordinates(aq.sf)[,"Y"])


## ----make-grid-1----------------------------------------------------
(n.col <- length(seq.e <- seq(min.x <- floor(min(aq$E)/1000)*1000,
                              max.x <- ceiling(max(aq$E)/1000)*1000, by=1000)))
(n.row <- length(seq.n <- seq(min.y <- floor(min(aq$N)/1000)*1000, 
                              max.y <- ceiling(max(aq$N)/1000)*1000, by=1000)))


## ----make-rast------------------------------------------------------
grid1km <- rast(nrows = n.row, ncols = n.col, 
                xmin=min.x, xmax=max.x, 
                ymin=min.y, ymax=max.y, crs = st_crs(aq.sf)$proj4string, 
                resolution = 1000, names="z")
values(grid1km) <- NA_real_
class(grid1km)
dim(grid1km)
summary(grid1km)
st_crs(grid1km)$proj4string
st_bbox(grid1km)    # replaces `raster` package `bbox()`


## ----plot-rast------------------------------------------------------
plot(grid1km); grid()


## ----rast-to-df-----------------------------------------------------
grid1km.df <- as.data.frame(grid1km, xy = TRUE, na.rm = FALSE) 
names(grid1km.df)[1:2] <- c("E", "N")            # match the names of the point dataset
summary(grid1km.df)


## ----tidy=FALSE-----------------------------------------------------
pred.ts2 <- predict.lm(model.ts2, 
                       newdata = grid1km.df,
                       interval = "prediction", level = 0.95)
summary(pred.ts2)
class(pred.ts2)


## ----ts2-grid-------------------------------------------------------
summary(values(grid1km))
summary(pred.ts2[,"fit"])
values(grid1km) <- pred.ts2[,"fit"]
summary(values(grid1km))


## ----plot-ts2-------------------------------------------------------
plot(grid1km, main = "OLS 2nd order predicted surface")
points(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], pch=16, 
     col = ifelse(res.ts2 < set_units(0, m), "red", "green"),
     cex=2*abs(res.ts2)/max(abs(res.ts2))
     )


## ----ts2-to-df------------------------------------------------------
summary(pred.ts2)
summary(grid1km.df)
grid1km.df[, 3:5] <- pred.ts2
names(grid1km.df)[3:5] <- c("ts2.fit", "ts2.lwr", "ts2.upr")
summary(grid1km.df)


## ----summarize-ts2-uncert-------------------------------------------
summary(grid1km.df$ts2.diff.range <- grid1km.df$ts2.upr - grid1km.df$ts2.lwr)
summary(100*grid1km.df$ts2.diff.range/grid1km.df$ts2.fit)


## ----map-ols-diff-range, tidy=FALSE, width='0.8\textwidth'----------
grid1km.diff <- grid1km
values(grid1km.diff) <- grid1km.df$ts2.diff.range 
plot(grid1km.diff, col=cm.colors(64),
     main="Range of 95% prediction interval, 2nd-order trend, OLS fit")
grid()
points(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], pch=16, 
     col = "gray")



## ----child-data, child='ex_TrendSurface_ex2-knitr_sf.Rnw'-----------

## ----ts2-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('ex_TrendSurface_sf.Rnw')


## ----plot-marginal-loess, fig.width=10, fig.height=5, out.width='\\textwidth'----
g1 <- ggplot(drop_units(aq), aes(x=N, y=zm)) +
    geom_point() +
    geom_smooth(method="loess") +
    labs(y = "elevation [m]")
g2 <- ggplot(drop_units(aq), aes(x=E, y=zm)) +
    geom_point() +
    geom_smooth(method="loess")  +
    labs(y = "elevation [m]")
grid.arrange(g1, g2, ncol = 2)


## ----fit-gam--------------------------------------------------------
model.gam <- gam(zm ~ s(E, N), data=drop_units(aq))
summary(model.gam)


## ----compare-gam-ts, fig.width=10, fig.height=5, out.width='\\textwidth', tidy=FALSE----
resid.gam <- residuals(model.gam)
summary(resid.gam)
summary(residuals(model.ts2))
par(mfrow=c(1,2))
hist(resid.gam, xlim=c(-20, 20), 
     breaks=seq(-20, 20,by=4), main="Residuals from GAM")
rug(residuals(model.gam))
hist(residuals(model.ts2), xlim=c(-20, 20),
     breaks=seq(-20,20,by=4), main="Residuals from 2nd-order OLS trend")
rug(residuals(model.ts2))
par(mfrow=c(1,1))


## ----compare-gam-ts-gg, fig.width=10, fig.height=5, out.width='\\textwidth', tidy=FALSE----
g1 <- ggplot(data=as.data.frame(resid.gam), aes(resid.gam)) +
    geom_histogram(breaks=seq(-20,20,by=2), 
                   fill="lightblue", color="black", alpha=0.9) +
    geom_rug() +
    labs(title = "Residuals from GAM",
         x = expression(paste(Delta, m)))
g2 <- ggplot(data=as.data.frame(residuals(model.ts2)), aes(residuals(model.ts2))) +
    geom_histogram(breaks=seq(-20,20,by=2), 
                   fill="lightgreen", color="darkblue", alpha=0.9) +
    geom_rug() +
    labs(title = "Residuals from 2nd order polynomial trend surface", 
         x = expression(paste(Delta, m)))
grid.arrange(g1, g2, nrow=1)


## ----gam-resid-bubble, out.width='.6\\linewidth'--------------------
aq$resid.gam <- resid.gam
ggplot(data = drop_units(aq)) +
    aes(x=E, y=N, size = abs(resid.gam), 
        col = ifelse((resid.gam < 0), "red", "green")) +
    geom_point(alpha=0.7) +
    scale_size_continuous(name = expression(paste(plain("residual ["), 
                                                  reDelta, m, plain("]"))),
                          breaks=seq(0,12, by=2)) +
    scale_color_manual(name = "polarity",
                       labels = c("negative","positive"),
                       values = c("red","green","blue"))


## ----gam-resid-vgm, out.width='.6\\linewidth', fig.height=5, fig.width=7, eval=FALSE, echo=FALSE----
## aq.sf$resid.gam <- resid.gam
## vr <- variogram(resid.gam ~ 1, loc=aq.sf)
## plot(vr, pl=T)


## ----plot-gam-2D, out.width='\\linewidth'---------------------------
plot.gam(model.gam, rug = TRUE, se = TRUE, select=1,
         scheme=1, theta=30+130, phi=30)


## ----plot-gam-2D-se, out.width='\\linewidth'------------------------
vis.gam(model.gam, plot.type="persp", color="terrain", 
        theta=160, zlab="elevation", se=1.96)


## ----plot-gam-fits--------------------------------------------------
(rmse.gam <- sqrt(sum(residuals(model.gam)^2)/length(residuals(model.gam))))


## ----predict-gam----------------------------------------------------
tmp <- predict.gam(object=model.gam, 
                   newdata=grid1km.df, 
                   se.fit=TRUE)
summary(tmp$fit)
summary(tmp$se.fit)
str(tmp)


## ----spgrid-gam-----------------------------------------------------
grid1km.df$pred.gam <- as.numeric(tmp$fit)
grid1km.df$pred.gam.se <- as.numeric(tmp$se.fit)


## ----map-gam, tidy=FALSE, out.width='0.8\\textwidth'----------------
grid1km.gam <- grid1km
values(grid1km.gam) <- grid1km.df$pred.gam
plot(grid1km.gam, main="GAM prediction"); grid()
points(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], pch=16, 
     col=ifelse(aq$resid.gam < 0, "red", "green"),
     cex=2*abs(aq$resid.gam)/max(abs(aq$resid.gam)))


## ----map-gam-se, tidy=FALSE, out.width='0.8\\textwidth'-------------
grid1km.gam.se <- grid1km
values(grid1km.gam.se) <- grid1km.df$pred.gam.se 
plot(grid1km.gam.se, main="GAM prediction standard error",
     col=cm.colors(64))
grid()
points(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], pch=16, 
     col="grey")


## ----diff-gam-ts2, tidy=FALSE, out.width='0.7\\textwidth'-----------
summary(grid1km.df$diff.gam.ols <- 
            grid1km.df$pred.gam - grid1km.df$ts2.fit)
grid1km.diff.gam.ols <- grid1km
values(grid1km.diff.gam.ols) <- grid1km.df$diff.gam.ols
plot(grid1km.diff.gam.ols, 
     main="Difference (GAM - 2nd order trend surface) predictions, [m]",
     col=topo.colors(64))
grid()


## ----TPS-setup------------------------------------------------------
aq.tps <- aq[, c("E","N", "zm")]
aq.tps$coords <- matrix(c(aq.tps$E, aq.tps$N), byrow=F, ncol=2)
str(aq.tps$coords)


## ----TPS-fit--------------------------------------------------------
surf.1 <-Tps(aq.tps$coords, aq.tps$zm)
summary(surf.1)


## ----TPS-grid-------------------------------------------------------
grid.coords.m <- as.matrix(grid1km.df[, c("E", "N")], ncol=2)
str(grid.coords.m)


## ----TPS-predict----------------------------------------------------
surf.1.pred <- predict.Krig(surf.1, grid.coords.m)
summary(grid1km.df$pred.tps <- as.numeric(surf.1.pred))


## ----TPS-get-resid--------------------------------------------------
grid1km.tps <- grid1km 
values(grid1km.tps) <- surf.1.pred
tmp <- extract(grid1km.tps, st_coordinates(aq.sf))
names(tmp)
summary(aq.sf$resid.tps <- (drop_units(aq.sf$zm) - tmp$z))


## ----TPS-resid-hist, fig.width=10, fig.height=7, out.width='0.8\\textwidth'----
hist(aq.sf$resid.tps, main="Thin-plate spline residuals", breaks=16)
rug(aq.sf$resid.tps)


## ----TPS-pred-map, fig.width=10, fig.height=7, out.width='0.8\\textwidth'----
plot(grid1km.tps,
     main = "2D Thin-plate spline surface")
points(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], pch=16, 
     col=ifelse(aq.sf$resid.tps < 0, "red", "green"),
     cex=2*abs(aq.sf$resid.tps)/max(abs(aq.sf$resid.tps)))


## ----TPS-GAM-diff, out.width='0.8\\textwidth'-----------------------
grid1km.diff.tps.gam <- grid1km
grid1km.diff.tps.gam <- grid1km.tps - grid1km.gam
values(grid1km.diff.tps.gam) <- (values(grid1km.tps) - values(grid1km.gam))
plot(grid1km.diff.tps.gam,
     col = topo.colors(64),
     main = "Thin-plate spline less GAM fits, [m]",
     xlab = "East", ylab = "North")



## ----child-data, child='ex_TrendSurface_ex3-knitr_sf.Rnw'-----------

## ----ts3-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('ex_TrendSurface.Rnw')


## ----extract-resid--------------------------------------------------
aq$fit.ts2 <- fitted(model.ts2)
aq.sf$fit.ts2 <- fitted(model.ts2)
aq$res.ts2  <- residuals(model.ts2)
aq.sf$res.ts2  <- residuals(model.ts2)


## ----compute-v, fig.width=12, fig.height=6, out.width='\\textwidth'----
vr.c <- variogram(res.ts2 ~ 1, loc = aq.sf, cutoff = 40000, cloud = T)
vr <- variogram(res.ts2 ~ 1, loc = aq.sf, cutoff = 40000)
p1 <- plot(vr.c, col = "blue", pch = 20, cex = 0.5,
                xlab = "separation [m]", ylab = "semivariance [m^2]")
p2 <- plot(vr, plot.numbers = T, col = "blue", pch = 20, cex = 1.5,
                xlab = "separation [m]", ylab = "semivariance [m^2]")
print(p1, split = c(1,1,2,1), more = T)
print(p2, split = c(2,1,2,1), more = F)


## ----fit-gls2, tidy=FALSE-------------------------------------------
model.ts2.gls <- gls(
   model = drop_units(zm) ~ N + E + I(N^2) + I(E^2) + I(E * N), 
   data = aq,
   method = "ML",
   correlation = corExp(form = ~E + N, 
      nugget = FALSE,
      value = 10000)   # initial value of the range parameter
 )
class(model.ts2.gls)
summary(model.ts2.gls)


## ----compare-gls-ols-coef-------------------------------------------
coef(model.ts2.gls) - coef(model.ts2)
round(100*(coef(model.ts2.gls) - coef(model.ts2))
                   /coef(model.ts2),1)


## ----gls-show-intervals---------------------------------------------
intervals(model.ts2.gls, level=0.90)


## ----gls-ols-resid-compare------------------------------------------
summary(res.ts2.gls <- residuals(model.ts2.gls))
summary(res.ts2)


## ----tidy=FALSE-----------------------------------------------------
pred.ts2.gls <- predict(model.ts2.gls, newdata=grid1km.df)
summary(pred.ts2.gls)


## ----gls-grid-------------------------------------------------------
grid1km.gls <- grid1km
values(grid1km.gls) <- pred.ts2.gls
summary(values(grid1km.gls))


## ----plot-ts2-gls---------------------------------------------------
plot(grid1km.gls, main = "GLS 2nd order predicted surface")
     points(st_coordinates(aq.sf)[,2] ~ st_coordinates(aq.sf)[,1], pch=16, 
     col = ifelse((res.ts2.gls < 0), "red", "green"),
     cex=2*abs(res.ts2.gls)/max(abs(res.ts2.gls))
)


## ----gls-ols-diff---------------------------------------------------
grid1km.gls.ols.diff <- (grid1km.gls - grid1km)
plot(grid1km.gls.ols.diff, col = topo.colors(64),
     main = "GLS - OLS 2nd order trend surfaces, m")


## ----post-gls-resid, tidy=FALSE-------------------------------------
summary(res.ts2.gls)
plot(aq$N ~ aq$E, cex=3*abs(res.ts2.gls)/max(abs(res.ts2.gls)),
     col=ifelse(res.ts2.gls > 0, "green", "red"),
     xlab="E", ylab="N",
     main="Residuals from 2nd-order trend, GLS fit",
     sub="Positive: green; negative: red", asp=1)
grid()


## ----resid-v-gls, out.width='0.8\\textwidth', fig.height=8, fig.width=10----
aq.sf$res.ts2.gls <- residuals(model.ts2.gls)
vr.gls <- variogram(res.ts2.gls ~ 1, loc=aq.sf,
                    cutoff = 40000)
plot(vr.gls, plot.numbers=T,
     main="Residuals from second-order GLS trend",
     xlab = "separation [m]", ylab = "semivariance [m^2]")


## ----resid-v-gls-fitted, out.width='0.8\\textwidth', fig.height=8, fig.width=10----
vr.gls.m <- vgm(psill=40, model="Exp", range=22000/3, nugget=0)
(vr.gls.m.f <- fit.variogram(vr.gls, vr.gls.m))
plot(vr.gls, model=vr.gls.m.f, plot.numbers=T,
     xlab = "separation [m]", ylab = "semivariance [m^2]")


## ----compare-gls-vgm-gls-range--------------------------------------
print(vr.gls.m.f)
intervals(model.ts2.gls)$corStruct[2]


## ----ok-ts2-resids--------------------------------------------------
grid1km.sf <- st_as_sf(grid1km.df, coords = c("E", "N"))
st_crs(grid1km.sf) <- st_crs(grid1km)
kr <- krige(res.ts2.gls ~ 1,
            loc = aq.sf, 
            newdata = grid1km.sf,
            model=vr.gls.m.f)
summary(kr)
class(kr)


## ----plot-ok-resids, tidy=FALSE, out.width='0.7\\textwidth'---------
plot(kr["var1.pred"], pch=15, nbreaks=24,
     main="Residuals from GLS trend, m")


## ----plot-ok-resids-se, tidy=FALSE, out.width='0.7\\textwidth'------
kr$var1.sd <- sqrt(kr$var1.var)
summary(kr)
plot(kr["var1.sd"], pch=15, nbreaks=24, pal = heat.colors,
     main="Standard errors of residuals from GLS trend, m")


## ----make-rkgls-grid------------------------------------------------
grid1km.df$kr <- kr$var1.pred
grid1km.df$pred.ts2.gls <- pred.ts2.gls
grid1km.df$rkgls <- grid1km.df$pred.ts2.gls + grid1km.df$kr
summary(grid1km.df)
grid1km.rkgls <- grid1km
values(grid1km.rkgls) <- grid1km.df$rkgls


## ----plot-rk-pred,  tidy=FALSE, out.width='0.8\\textwidth'----------
plot(grid1km.rkgls,
     main="GLS-RK prediction, aquifer elevation, m.a.s.l.")


## ----diff-gam-gls, tidy=FALSE, out.width='0.7\\textwidth'-----------
summary(grid1km.rkgls.gam <- (grid1km.rkgls - grid1km.gam))
plot(grid1km.rkgls.gam,  main ="RK-GLS - GAM fits, difference, m",
             col = topo.colors(64))



## ----child-data, child='ex_TrendSurface_ex4-knitr_sf.Rnw'-----------

## ----ts4-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('ex_TrendSurface.Rnw')


## ----add-E-N-aq.sf--------------------------------------------------
str(st_coordinates(aq.sf))
names(aq.sf)
aq.sf$E <- st_coordinates(aq.sf)[ , "X"]
aq.sf$N <- st_coordinates(aq.sf)[ , "Y"]
names(aq.sf)


## ----uk-vgm, out.width='0.8\\textwidth', fig.height=8, fig.width=10,tidy=FALSE----
# summary(apply(st_coordinates(aq.sf), MARGIN = 1, FUN = prod))
vr <- variogram(zm ~ E + N + I(E^2) + I(N^2) + I(E*N), 
                locations = aq.sf,
                cutoff = 40000)
plot(vr, plot.numbers = TRUE,
     main = "Residuals from 2nd-order OLS trend surface",
     xlab = "separation (m)",
     ylab = "semivariance (m^2)")


## ----uk-vgm-model, out.width='0.8\\textwidth', fig.height=8, fig.width=10,tidy=FALSE----
(vr.m.f <- fit.variogram(vr, vgm(35, "Exp", 22000/3, 0)))
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")


## ----compare-vgms-ols-gls-------------------------------------------
print(vr.m.f)       # OLS trend residuals
print(vr.gls.m.f)   # GLS trend residuals


## ----internals-sf-objects-------------------------------------------
names(grid1km.sf)
str(grid1km.sf$geometry)
names(aq.sf)
str(aq.sf$geometry)


## ----add-E-N-grid1km.sf---------------------------------------------
str(st_coordinates(grid1km.sf))
grid1km.sf$E <- st_coordinates(grid1km.sf)[ , "X"]
grid1km.sf$N <- st_coordinates(grid1km.sf)[ , "Y"]
names(grid1km.sf)


## ----predict-uk-----------------------------------------------------
k.uk <- krige(zm ~ E + N + I(E^2) + I(N^2) + I(E*N),
              locations = aq.sf,
              newdata = grid1km.sf, 
              model=vr.m.f)
summary(k.uk)


## ----plot-uk-resids, tidy=FALSE, out.width='0.7\\textwidth'---------
plot(k.uk["var1.pred"], pch=15, nbreaks=24,
     main="UK predictions, m")


## ----plot-uk-resids-se, tidy=FALSE, out.width='0.7\\textwidth'------
k.uk$var1.sd <- sqrt(k.uk$var1.var)
summary(k.uk)
plot(k.uk["var1.sd"], pch=15, nbreaks=24, pal = heat.colors,
     main="Standard errors of UK predictions, m")


## ----compare-uk-kr-sd-----------------------------------------------
summary(k.uk$var1.sd)
summary(kr$var1.sd)


## -------------------------------------------------------------------
grid1km.uk <- grid1km
values(grid1km.uk) <- k.uk$var1.pred
summary(grid1km.diff.uk.rkgls <- (grid1km.uk - grid1km.rkgls))


## ----hist-uk-glsrk-diff, tidy=FALSE, fig.width=8, fig.height=6------
hist(grid1km.diff.uk.rkgls, main = "UK - GLS-RK prediction differences",
     freq = FALSE, xlab = "difference, UK - GLS-RK")


## ----plot-uk-glsrk-diff, tidy=FALSE, out.width='0.7\\textwidth'-----
plot(grid1km.diff.uk.rkgls,  sub="UK - GLS-RK predictions",
             main="difference, m", xlab="East", ylab="North",
             col = topo.colors(64))



## ----child-data, child='ex_TrendSurface_answers-knitr_sf.Rnw'-------

## ----tsa-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('ex_TrendSurface_sf.Rnw')