## ----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!
opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'), out.format='latex', use.highlight=TRUE)
opts_chunk$set(fig.path='kgraph/gsintro_sf', fig.align='center', fig.show='hold', prompt=FALSE, comment=NA, 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'-------------------------------------
rm(list=ls())


## ----child-gs_load, child='gs_load_sf.Rnw'--------------------------

## ----load-set-parent, echo=FALSE, cache=FALSE-----------------------
set_parent('gs_intro_sf.Rnw')


## ----eval=FALSE-----------------------------------------------------
## install.packages(c("sf", "sp", "gstat", "ggplot2"), dependencies=TRUE)


## ----load-first-expr------------------------------------------------
2*pi/360


## ----load-load-libraries--------------------------------------------
library(sf)
library(gstat)


## -------------------------------------------------------------------
search()


## ----load-show-datasets, eval=FALSE---------------------------------
## data()


## ----load-show-datasets-2, eval=FALSE-------------------------------
## data(package="sp")


## -------------------------------------------------------------------
ls()
data("meuse", package = "sp")
ls()


## ----load-str-df----------------------------------------------------
str(meuse)


## ----load-show-df-as-matrix-----------------------------------------
dim(meuse)
dim(as.matrix(meuse))
str(as.matrix(meuse))
head(as.matrix(meuse))


## ----eval=FALSE-----------------------------------------------------
## help(meuse)


## ----load-summary-meuse---------------------------------------------
summary(meuse)


## ----child-gs_break, child='gs_break_sf.Rnw'------------------------

## ----break-set-parent, echo=FALSE, cache=FALSE----------------------
set_parent('gs_intro_sf.Rnw')


## ----eval=FALSE-----------------------------------------------------
## q()


## ----reload-pkgs----------------------------------------------------
library(sf)
library(gstat)


## ----child-gs_uni, child='gs_uni_sf.Rnw'----------------------------

## ----uni-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('gs_intro_sf.Rnw')


## -------------------------------------------------------------------
meuse$zinc
sort(meuse$zinc)


## ----uni-show-hist-Zn, out.width='0.7\\textwidth'-------------------
hist(meuse$zinc, breaks=16, main="Meuse River soil pollution study",
     xlab="Zn [ppm]")
rug(meuse$zinc)
summary(meuse$zinc)


## ----uni-show-hist-logZn, out.width='0.7\\textwidth'----------------
meuse$logZn <- log10(meuse$zinc)
hist(meuse$logZn, breaks=16, main="Meuse River soil pollution study",
     xlab="Zn [log10(ppm)]")
rug(meuse$logZn)
summary(meuse$logZn)


## -------------------------------------------------------------------
head(meuse$zinc)
head(meuse$logZn)


## -------------------------------------------------------------------
meuse$logCu <- log10(meuse$copper)
str(meuse)


## ----child-gs_bi, child='gs_bi_sf.Rnw'------------------------------

## ----bi-set-parent, echo=FALSE, cache=FALSE-------------------------
set_parent('gs_intro_sf.Rnw')


## ----bi-plot-logZnlogCu---------------------------------------------
plot(meuse$logZn ~ meuse$logCu,
     xlab = "Cu [log10(ppm)]",
     ylab = "Zn [log10(ppm)]",
     main = "Meuse River soil pollution study")


## ----bi-find-unusual-1----------------------------------------------
which(meuse$logZn < 2.6)
which(meuse$logCu > 1.6)


## ----bi-find-unusual-2----------------------------------------------
which((meuse$logZn < 2.6) & (meuse$logCu > 1.6))


## ----bi-find-unusual-3----------------------------------------------
ix <- which((meuse$logZn < 2.6) & (meuse$logCu > 1.6))


## ----bi-show-unusual------------------------------------------------
meuse[ix, ]


## ----bi-row.names.meuse---------------------------------------------
print(row.names(meuse))


## ----m-lzn-lcu------------------------------------------------------
m.lzn.lcu <- lm(logZn ~ logCu, data=meuse)


## ----bi-m-lzn-lcu-matrix--------------------------------------------
model.matrix(m.lzn.lcu)[1:6,]
meuse$logCu[1:6]


## -------------------------------------------------------------------
ls()
class(meuse)
class(m.lzn.lcu)


## -------------------------------------------------------------------
summary(m.lzn.lcu)


## ----bi-plot-hist-LznLcu-resid, out.width='0.7\\textwidth'----------
hist(residuals(m.lzn.lcu),
     main = "Residuals from linear model log10(Zn) ~ log10(Cu)",
     xlab = "Zn (log10[ppm]")
rug(residuals(m.lzn.lcu))


## ----bi-plot-m-lzn-lcu, out.width='\\textwidth', fig.width=12, fig.height=4----
par(mfrow=c(1,3))
plot(m.lzn.lcu, which=c(1,2,5))
par(mfrow=c(1,1))


## ----find-largest-resids--------------------------------------------
(which.ix <- which(abs(residuals(m.lzn.lcu)) > 0.3))
# the records in the data frame that correspond to large absolute residuals
meuse[which.ix,]
# these large absolute residuals
residuals(m.lzn.lcu)[which.ix]
# the largest one
which.max(abs(residuals(m.lzn.lcu))[which.ix])


## ----bi-show.base.colours-------------------------------------------
palette()


## ----bi-plot-lZn-lCu-ffreq, out.width='0.75\\textwidth', keep.source=TRUE----
plot(meuse$logZn ~ meuse$logCu, asp=1, col=meuse$ffreq, pch=20,
     xlab="log10(Cu ppm)", ylab="log10(Zn ppm)")
abline(m.lzn.lcu)
legend("topleft", legend=c("2 years","10 years", "50 years"),
       pch=20, col=1:3)


## ----bi-ffreq.table-------------------------------------------------
table(meuse$ffreq)


## ----bi-boxplot-Zn-ffreq, out.width='0.6\\textwidth', keep.source=T----
boxplot(meuse$logZn ~ meuse$ffreq, xlab="Flood frequency class",
        ylab="log10-Zn ppm",
        main="Metal concentration per flood frequency class",
        boxwex=0.4, col="lightblue")


## ----bi-m-lzn-ff----------------------------------------------------
m.lzn.ff <- lm(logZn ~ ffreq, data=meuse)


## ----bi-m-lzn-ff-matrix---------------------------------------------
model.matrix(m.lzn.ff)[1:6,]


## ----bi-summary-lzn-ff----------------------------------------------
summary(m.lzn.ff)


## ----bi-means-------------------------------------------------------
pred.ff <- predict(m.lzn.ff, newdata=data.frame(ffreq=as.factor(c(1:3))), se.fit=TRUE,
                   interval="prediction")
pred.ff$fit
pred.ff$se.fit


## ----bi-m-mixed-----------------------------------------------------
m.lzn.ff.lcu <- lm(logZn ~ ffreq + logCu, data=meuse)
summary(m.lzn.ff.lcu)


## ----bi-m-mixed-anova-----------------------------------------------
anova(m.lzn.ff.lcu, m.lzn.lcu)


## ----echo=FALSE, results='hide'-------------------------------------
# save anova results for answers
tmp.anova.result <- anova(m.lzn.ff.lcu, m.lzn.lcu)
-tmp.anova.result$Df[2]
-round(tmp.anova.result$"Sum of Sq"[2],3)
round(tmp.anova.result$"Pr(>F)"[2],2)


## ----bi-m-mixed-i---------------------------------------------------
m.lzn.ff.lcu.i <- lm(logZn ~ ffreq * logCu, data=meuse)
summary(m.lzn.ff.lcu.i)
anova(m.lzn.ff.lcu.i, m.lzn.lcu)


## ----echo=FALSE, results='hide'-------------------------------------
# save anova results for answers
tmp.anova.result.i <- anova(m.lzn.ff.lcu.i, m.lzn.lcu)
-tmp.anova.result.i$Df[2]
-round(tmp.anova.result.i$"Sum of Sq"[2],3)
round(tmp.anova.result.i$"Pr(>F)"[2],2)


## ----bi-plot-mixed-i, keep.source=T---------------------------------
with(meuse, plot(logZn ~ logCu, col=ffreq, pch = 20,
                 xlab = "log10(Cu)", ylab = "log10(Zn)"))
legend("topleft", legend=levels(meuse$ffreq), pch=20,
       col=1:3)
title(main = "Relation between log10(Zn) and log10(Cu)")
title(sub = 
          "Interaction: solid lines; per class; single: dashed line")
abline(lm(logZn ~ logCu, data = meuse),col = "purple",
       lty=3, lwd=2.5)
abline(lm(logZn ~ logCu, data = meuse,
          subset=(meuse$ffreq==1)),col = 1)
abline(lm(logZn ~ logCu, data = meuse,
          subset=(meuse$ffreq==2)),col = 2)
abline(lm(logZn ~ logCu, data = meuse,
          subset=(meuse$ffreq==3)),col = 3)


## ----child-gs_sf, child='gs_sf.Rnw'---------------------------------

## ----sp-set-parent, echo=FALSE, cache=FALSE-------------------------
set_parent('gs_intro_sf.Rnw')


## ----sp-load-load-libraries-----------------------------------------
library(sf)
library(gstat)


## ----show-classes---------------------------------------------------
class(1); class("A"); class(TRUE)
class(list(1, "A", TRUE))
class(as.matrix(1:9, nrow=3))
class(meuse)
class(meuse$ffreq)


## ----assign-coords--------------------------------------------------
class(meuse)
meuse.sf <- st_as_sf(meuse, coords = c("x","y"))
class(meuse.sf)


## -------------------------------------------------------------------
str(meuse.sf)


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


## ----assign-crs-rdh-------------------------------------------------
st_crs(meuse.sf) <- 28992
print(st_crs(meuse.sf))


## ----load-meuse-riv-------------------------------------------------
data(meuse.riv, package="sp")
class(meuse.riv)
meuse.riv.sf <- st_sfc(st_linestring(meuse.riv, dim="XY"), crs = st_crs(meuse.sf))
class(meuse.riv.sf)
summary(meuse.riv.sf)


## ----spatvar-plot-meuse-riv, out.width='0.7\\textwidth'-------------
require(ggplot2)
ggplot(data = meuse.sf) +
    geom_sf(mapping = aes(size=zinc, color=dist.m)) +
    labs(x = "Longitude", y = "Latitude", title = "Meuse River (NL)",
         color="Distance to river [m]", size = "Zn concentration, ppm")


## ----convert-meuse-grid---------------------------------------------
data(meuse.grid, package="sp")
class(meuse.grid)
names(meuse.grid)
meuse.grid.sf <- st_as_sf(meuse.grid, coords = c("x", "y"))
st_crs(meuse.grid.sf) <- st_crs(meuse.sf)


## ----show-grid-attr-------------------------------------------------
summary(meuse.grid.sf)


## ----plot-ffreq-grid------------------------------------------------
plot(meuse.grid.sf["ffreq"], pch = 15,
     main = "Meuse River, flooding frequency classes")


## ----child-gs_trees, child='gs_trees_sf.Rnw'------------------------

## ----trees-set-parent, echo=FALSE, cache=FALSE----------------------
set_parent('gs_intro_sf.Rnw')


## ----trees-load-if-necess, echo=FALSE, results='hide'---------------
if (!exists("meuse")) data(meuse, package="sp")
if (!exists("meuse$logZn")) meuse$logZn <- log10(meuse$zinc)
if (!exists("meuse.sf")) { 
    require(sf)
    meuse.sf <- st_as_sf(meuse, coords = c("x", "y"))
    }
if (!exists("meuse.grid.sf")) { 
    data(meuse.grid, package="sp")
    meuse.grid.sf <- st_as_sf(meuse.grid, coords = c("x", "y"))
    }


## ----trees-load-rpart-----------------------------------------------
library(rpart)
library(rpart.plot)


## ----trees-set-seed-regstree, echo=FALSE----------------------------
# to get consistent results for comments
set.seed(6345789)


## ----trees-comp-rpart-----------------------------------------------
m.lzn.rp <- rpart(logZn ~ ffreq + dist + elev + soil,
                  data = meuse.sf,
                  minsplit=2,
                  cp=0.003)


## ----trees-print-unpruned-------------------------------------------
print(m.lzn.rp)


## ----trees-plot-unpruned, out.width='\\textwidth'-------------------
rpart.plot(m.lzn.rp, digits=3, type=4, extra=1)


## ----trees-count-leaves---------------------------------------------
sum(m.lzn.rp$frame$var == '<leaf>')
sum(m.lzn.rp$frame$var !="<leaf>")


## ----trees-import-rpart---------------------------------------------
x <- m.lzn.rp$variable.importance 
data.frame(variableImportance = 100 * x / sum(x))


## ----trees-plotcp-rp, out.width='\\textwidth'-----------------------
printcp(m.lzn.rp)
plotcp(m.lzn.rp)


## ----trees-find-min-xerror------------------------------------------
head(cp.table <- m.lzn.rp[["cptable"]], 12)
(cp.ix <- which.min(cp.table[,"xerror"]))
(xerror.min <- cp.table[cp.ix,"xerror"])
print(cp.table[cp.ix,])
cp.min <- cp.table[cp.ix,"CP"]


## ----trees-cp-table-min---------------------------------------------
(cp.min.plus.sd <- cp.table[cp.ix,"xerror"] + cp.table[cp.ix,"xstd"])
cp.ix.sd <- min(which(cp.table[,"xerror"] < cp.min.plus.sd))
print(cp.table[cp.ix.sd,])
cp.min.sd <- cp.table[cp.ix.sd,"CP"]


## ----trees-prune-rpart----------------------------------------------
(m.lzn.rpp <- prune(m.lzn.rp, cp=cp.min))


## ----trees-plot-rpart-----------------------------------------------
rpart.plot(m.lzn.rpp, digits=3, type=4, extra=1)


## ----trees-predict-rpart, out.width='0.8\\textwidth', tidy = FALSE----
p.rpp <- predict(m.lzn.rpp, newdata=meuse.sf)
length(unique(p.rpp))
summary(r.rpp <- meuse$logZn - p.rpp)
sqrt(sum(r.rpp^2)/length(r.rpp))
plot(meuse$logZn ~ p.rpp, asp=1, pch=20, xlab="fitted", ylab="actual",
     xlim=c(2,3.3), ylim=c(2,3.3),
     main="log10(Zn), Meuse topsoils, Regression Tree")
grid()
abline(0,1)


## ----show-meuse-grid-attr-------------------------------------------
names(meuse.grid.sf)


## ----gam-idw-elev, tidy = FALSE-------------------------------------
tmp <- gstat::idw(elev ~ 1, locations=meuse.sf,
           nmax=16, idp=2, newdata=meuse.grid.sf)
meuse.grid.sf$elev <- tmp$var1.pred; rm(tmp)
summary(meuse.sf$elev)
summary(meuse.grid.sf$elev)
plot(meuse.grid.sf["elev"], pch=15,
     nbreaks = 64,
     pal = terrain.colors,
     main = "Elevation [m.a.s.l.]")


## ----trees-predict-map----------------------------------------------
rt.grid <- predict(m.lzn.rpp, newdata = meuse.grid.sf)
meuse.grid.sf$rt.pred <- rt.grid; rm(rt.grid)
plot(meuse.grid.sf["rt.pred"], pch=15,
     nbreaks = 64,
     main = "Regression tree prediction, Log10(Zn [ppm])")


## ----trees-set-seed-classtree, echo=FALSE---------------------------
# to get consistent results for comments
set.seed(6345789)


## ----trees-build-classtree------------------------------------------
m.ff.rp <- rpart(ffreq ~ dist + elev + soil,
                  data=meuse.sf,
                  minsplit=2,
#                  method="class",
                  cp=0.01)


## ----trees-show-classtree, fig.width = 10, fig.height = 10, out.width='\\textwidth'----
rpart.plot(m.ff.rp, type=4, extra=1, tweak = 0.8)


## ----trees-show-classtree-cp, fig.width=10, fig.height=6------------
printcp(m.ff.rp)
plotcp(m.ff.rp)
cp.table.class <- m.ff.rp[["cptable"]]
# total variance explained is sum of the CP for all split
sum(cp.table.class[,"CP"])


## ----trees-root-node-error------------------------------------------
(n <- length(m.ff.rp$y))
(class.majority <- which.max(m.ff.rp$parms$prior))
(class.majority.proportion <- 
     m.ff.rp$parms$prior[class.majority])
(1 - (class.majority.proportion))


## ----trees-hide-classtree-cp, echo=FALSE,results='hide'-------------
# for the answers, others can see from the table
cp.ix.class <- which.min(cp.table.class[,"xerror"])
# print(cp.table.class[cp.ix.class,])
(cp.min.class <- cp.table.class[cp.ix.class,"CP"])


## ----trees-prune-classtree------------------------------------------
cp.ix.class <- which.min(cp.table.class[,"xerror"])
(cp.min.class <- cp.table.class[cp.ix.class,"CP"])
(m.ff.rpp <- prune(m.ff.rp, cp=cp.min.class))


## ----trees-show-pruned-classtree, fig.width = 10, fig.height = 10, out.width='\\textwidth'----
par(mfrow=c(1,2))
rpart.plot(m.ff.rpp, type=4, extra=1, tweak = 0.8)
rpart.plot(m.ff.rpp, type=4, extra=4, tweak = 0.8)
par(mfrow=c(1,1))


## ----trees-model-rf, tidy = FALSE-----------------------------------

library(ranger)
m.lzn.rf <- ranger(logZn ~ ffreq + dist+ elev + soil,
                   data = as.data.frame(meuse.sf),
                   importance = "permutation")
print(m.lzn.rf)
print(sqrt(m.lzn.rf$prediction.error))
print(m.lzn.rf$r.squared)


## ----trees-mse-rf, out.width='0.55\\linewidth'----------------------
nt <- seq(1, 500, by = 2)
oob_mse <- vector("numeric", length(nt))
for(i in 1:length(nt)) {
  rf <- ranger(logZn ~ ffreq + dist+ elev + soil, meuse, 
  num.trees = nt[i],
  write.forest = FALSE, importance = "none")
  oob_mse[i] <- rf$prediction.error
}
plot(x = nt, y = oob_mse, type = "l", 
  main = "Error rate, RF model of log10(Zn)",
  ylab = "OOB mean square error")
grid()


## ----trees-varimp-rf, out.width='.55\\linewidth'--------------------
print(sort(sqrt(importance(m.lzn.rf)), decreasing = TRUE))
print(round(sort(100*sqrt(importance(m.lzn.rf))/sqrt(rf$prediction.error),
                 decreasing = TRUE)))


## ----rfexplain-depth-plot, fig.width=8, fig.height=5, out.width='\\linewidth'----
require(randomForestExplainer)
min_depth_frame <- min_depth_distribution(m.lzn.rf)
str(min_depth_frame)  # has results for all the trees
plot_min_depth_distribution(min_depth_frame)


## ----rfexplain-importance-frame-------------------------------------
importance.frame <- measure_importance(m.lzn.rf)
print(importance.frame)


## ----fig.width=8, fig.height=7, out.width='\\linewidth'-------------
plot_predict_interaction(m.lzn.rf, meuse, "dist", "elev")


## ----fig.width=6, fig.height=4, out.width='0.85\\linewidth'---------
require(pdp)
partial.dist <- pdp::partial(m.lzn.rf, pred.var = "dist", pred.data = meuse.sf)
plotPartial(partial.dist)
partial.elev <- pdp::partial(m.lzn.rf, pred.var = "elev", pred.data = meuse.sf)
plotPartial(partial.elev)


## ----fig.width=8, fig.height=7, out.width='\\linewidth'-------------
partial.dist.elev <- pdp::partial(m.lzn.rf, pred.var = c("dist", "elev", "ffreq"), 
                                  pred.data = meuse.sf)
plotPartial(partial.dist.elev)


## ----trees-predict-rf, out.width='0.75\\textwidth', tidy = FALSE----
p.rf <- predict(m.lzn.rf, data=as.data.frame(meuse.sf))
str(p.rf)
length(unique(p.rf$predictions))
summary(r.rpp <- meuse$logZn - p.rf$predictions)
sqrt(sum(r.rpp^2)/length(r.rpp))
plot(meuse$logZn ~ p.rf$predictions, asp=1, pch=20, xlab="fitted values",
     ylab="actual values",
     xlim=c(2,3.3), ylim=c(2,3.3),
     main="log10(Zn), Meuse topsoils, Random Forest")
grid()
abline(0,1)


## ----trees-oob-rf, tidy = FALSE-------------------------------------
summary(r.rpp.oob <- meuse$logZn - m.lzn.rf$predictions)
sqrt(sum(r.rpp.oob^2)/length(r.rpp.oob))
plot(meuse$logZn ~ m.lzn.rf$predictions, asp=1, pch=20,
     xlab="Out-of-bag cross-validation estimates",
     ylab="actual values", xlim=c(2,3.3), ylim=c(2,3.3),
     main="log10(Zn), Meuse topsoils, Random Forest")
grid()
abline(0,1)


## ----rf-predict-map-------------------------------------------------
rf.grid <- predict(m.lzn.rf, data=as.data.frame(meuse.grid.sf))
meuse.grid.sf$rf.pred <- rf.grid$predictions; rm(rf.grid)
plot(meuse.grid.sf["rf.pred"], pch = 15,
     nbreaks = 64, 
     main="Random forest prediction, log10(Zn [ppm])")


## ----trees-build-classRF--------------------------------------------
m.ff.rf <- ranger(ffreq ~ dist + elev + soil,
                  data = as.data.frame(meuse.sf),
                  importance = "impurity",
                  write.forest = FALSE)
print(m.ff.rf)
str(m.ff.rf)
table(m.ff.rf$predictions)
print(m.ff.rf$prediction.error)
print(sort(importance(m.ff.rf, type=1), decreasing = TRUE))


## ----trees-oob-classRF----------------------------------------------
p.rf <- m.ff.rf$predictions
table(meuse$ffreq) # observed
table(p.rf)        # predicted
(p.rf.cm <- table(p.rf, meuse$ffreq, dnn=c("Predicted","Observed")))


## ----trees-oob-classRF-matrix---------------------------------------
(d <-diag(p.rf.cm))
(row.sums <- apply(p.rf.cm, MARGIN=1, "sum"))
(mu.purity <- d/row.sums)


## ----child-gs_spatvar, child='gs_spatvar_sf.Rnw'--------------------

## ----spatvar-set-parent, echo=FALSE, cache=FALSE--------------------
set_parent('gs_intro_sf.Rnw')


## ----spatvar-plot-meuse-Zn-riv--------------------------------------
plot(meuse.sf["zinc"],
     reset = FALSE,
     nbreaks = 64, pch = 20,
     cex=4*meuse.sf$zinc/max(meuse.sf$zinc),
     main = "Zn concentration [ppm]")
plot(meuse.riv.sf, add = TRUE)


## -------------------------------------------------------------------
n <- length(meuse.sf$logZn)
n*(n-1)/2


## -------------------------------------------------------------------
dim(st_coordinates(meuse.sf))
st_coordinates(meuse.sf)[1,]
st_coordinates(meuse.sf)[2,]
(sep <- dist(st_coordinates(meuse.sf)[1:2,]))
(gamma <- 0.5 * (meuse.sf$logZn[1] - meuse.sf$logZn[2])^2)


## ----show-vgm-cloud, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
head(vc <- variogram(logZn ~ 1, meuse.sf, cutoff=120, cloud=TRUE))
class(vc)
print(plot(vc, pch = 20, cex=2, xlim=c(0,120),
     ylab="semivariance, (log10(Zn)^2)", xlab="separation, m"))


## ----which-vgm-cloud-large------------------------------------------
order(vc$gamma, decreasing=TRUE)
order(vc$gamma, decreasing=TRUE)[1]
which.max(vc$gamma)
(unusual.pair <- vc[which.max(vc$gamma),])
meuse[138, "logZn"]; meuse[76, "logZn"]
(gamma.76.138 <- 0.5 * (meuse[138, "logZn"] - meuse[76, "logZn"])^2)


## ----spatvar-show-vgm-logZn, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
(v <- variogram(logZn ~ 1, meuse.sf, cutoff=1300, width=90))
print(plot(v, plot.numbers = TRUE, pch = 20,
           xlab = "separation", ylab = "semivariance [log10(Zn)^2]"))


## ----spatvar-showvgms, out.width='\\textwidth'----------------------
print(show.vgms())


## ----spatvar-plot-guess-vgm-model, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
vm <- vgm(psill=0.12,model="Sph",range=850,nugget=0.01)
print(plot(v, pl=T, model=vm, pch=20,
           xlab = "separation", ylab = "semivariance [log10(Zn)^2]"))


## ----spatvar-plot-fitted-vgm-model, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
(vmf <- fit.variogram(v, vm))
print(plot(v, pl = T, model = vmf, pch = 20,
      xlab = "separation", ylab = "semivariance [log10(Zn)^2]"))


## ----child-gs_ok, child='gs_ok_sf.Rnw'------------------------------

## ----ok-set-parent, echo=FALSE, cache=FALSE-------------------------
set_parent('gs_intro_sf.Rnw')


## -------------------------------------------------------------------
k40 <- krige(logZn ~ 1, locations = meuse.sf,
             newdata = meuse.grid.sf,
             model = vmf)


## -------------------------------------------------------------------
str(k40)


## ----ok-plot-k40pred, tidy = FALSE----------------------------------
plot(k40["var1.pred"], pch=15,
     nbreaks = 64,
     main="OK prediction [log10(Zn ppm)]")


## ----ok-bpy---------------------------------------------------------
head(sp::bpy.colors(64))
sp::bpy.colors(64)[32]


## ----ok-plot-k40var, tidy = FALSE-----------------------------------
plot(k40["var1.var"], pch=15,
     nbreaks = 64,
     pal = cm.colors,
     main="OK prediction variance [log10(Zn ppm)^2]")


## ----ok-plot-k40pred-pts, tidy = FALSE------------------------------
plot(k40["var1.pred"], pch=15,
     nbreaks = 64,
     reset = FALSE,  # allow further elements to be plotted
     main="OK prediction [log10(Zn ppm)]")
plot(meuse.sf["logZn"], add = TRUE, pch = 21,
     col = "white",
     cex=4*meuse.sf$zinc/max(meuse.sf$zinc))


## ----ok-plot-k40var-pts, tidy = FALSE-------------------------------
plot(k40["var1.var"], pch=15,
     nbreaks = 64,
     reset = FALSE,  # allow further elements to be plotted
     pal = cm.colors,
     main="OK prediction variance [log10(Zn ppm)^2]")
plot(meuse.sf["logZn"], add = TRUE, pch = 20,
     col = "darkgray")



## ----child-gs_ik, child='gs_ik_sf.Rnw'------------------------------

## ----ik-set-parent, echo=FALSE, cache=FALSE-------------------------
set_parent('gs_intro_sf.Rnw')


## ----ik-compute-indicator-------------------------------------------
meuse.sf$zn.i <- (meuse.sf$zinc < 150)
summary(meuse.sf$zn.i)


## ----ik-show-i-vi, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
vi <- variogram(zn.i ~ 1, location=meuse.sf, cutoff=1500)
print(plot(vi, pl=T, pch = 20,
           xlab = "separation", ylab = "semivariance of indicator"))


## ----ik-show-i-vimf, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
(vimf <- fit.variogram(vi, vgm(0.12, "Sph", 1300, 0)))
print(plot(vi, pl = T, model = vimf, pch = 20,
           xlab = "separation", ylab = "semivariance of indicator"))


## ----ik-spplot-k40i-pred, out.width='0.65\\textwidth', tidy = FALSE----
k40.i <- krige(zn.i ~ 1, loc = meuse.sf,
               newdata = meuse.grid.sf, model = vimf)
plot(k40.i["var1.pred"], pch = 15,
     nbreaks = 64,
     pal = heat.colors,
     main="Probability Zn < 150")


## ----ik-spplot-k40i-var, out.width='0.65\\textwidth', tidy = FALSE----
k40.i$var1.sd <- sqrt(k40.i$var1.var)
plot(k40.i["var1.sd"], pch = 15,
     nbreaks = 64,
     pal = cm.colors,
     main="standard deviation, probability Zn < 150")


## ----ik-safe-80pct--------------------------------------------------
k40.i$safe150.80pct <- ifelse((k40.i$var1.pred <= 0.8), TRUE, FALSE)
summary(k40.i$safe150.80pct)
plot(k40.i["safe150.80pct"], pch = 15, pal = c("darkgreen", "red"),
             main="(p >= 0.8) below critical level (Zn < 150)")


## ----child-gs_rk, child='gs_rk_sf.Rnw'------------------------------

## ----rk-set-parent, echo=FALSE, cache=FALSE-------------------------
set_parent('gs_intro_sf.Rnw')


## ----rk-summary-meuse-grid-ffreq, tidy = FALSE----------------------
summary(meuse.grid.sf$ffreq)
plot(meuse.grid.sf["ffreq"],
     pal = topo.colors(3), pch = 15,
     main="Flooding frequency class")


## ----rk-ggplot-kffreq-pred, tidy = FALSE----------------------------
k.ffreq <- krige(logZn ~ ffreq, locations=meuse.sf,
                 newdata=meuse.grid.sf, model=NULL)
plot(k.ffreq["var1.pred"],
     pal = sp::bpy.colors, pch = 15,
     main="prediction by flood frequency, log10(Zn ppm)")


## ----rk-ggplot-kffreq-var, tidy = FALSE-----------------------------
plot(k.ffreq["var1.var"],
     pal = cm.colors, pch = 15,
     main="prediction variance, log-ppm Zn^2")


## ----rk-show-ffreq-vr, fig.width=8, fig.height=6, out.width='0.65\\textwidth', tidy = FALSE----
(vr <- variogram(logZn ~ ffreq, location=meuse.sf,
                 cutoff=1300, width=90))
print(plot(vr, plot.numbers=T, pch = 20, 
           xlab = "separation", ylab = "semivariance of residuals [log10(Zn)^2]",
           main="Residuals, flood frequency co-variable"))


## -------------------------------------------------------------------
(vrmf <- fit.variogram(vr, 
                       vgm(psill=0.08, model="Sph",
                           range=800, nugget=0.01)))
print(vmf)


## ----rk-show-ffreq-vr-model, fig.width=8, fig.height=6, out.width='0.65\\textwidth', tidy = FALSE----
print(plot(vr, plot.numbers=T, pch = 20, model=vrmf,
           xlab = "separation", ylab = "semivariance of residuals [log10(Zn)^2]",
           main="Residuals, flood frequency co-variable"))


## ----tidy = FALSE---------------------------------------------------
kr40 <- krige(logZn ~ ffreq, locations=meuse.sf,
              newdata=meuse.grid.sf, model=vrmf)


## ----rk-ggplot-kr40, out.width='0.65\\textwidth', tidy = FALSE------
plot(kr40["var1.pred"],
     nbreaks = 64, pch = 15,
     pal = sp::bpy.colors,
     main="KED-ffreq prediction, log-ppm Zn")


## ----tidy = FALSE---------------------------------------------------
(zmax <- max(k40$var1.pred,kr40$var1.pred))
(zmin <- min(k40$var1.pred,kr40$var1.pred))
(zmax <- round(zmax, 1) + 0.1)
(zmin <- round(zmin, 1) - 0.1)
(ramp <- seq(from=zmin, to=zmax, by=.1))


## ----plot-compare-ok-ked, out.width='\\textwidth', fig.width=10, fig.height=5, tidy = FALSE----
g.ok <- ggplot() + 
    geom_sf(data = k40, aes(col = var1.pred)) +
    labs(title = "OK", col = "log10(Zn)") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
g.ked <-ggplot() + 
    geom_sf(data = kr40, aes(col = var1.pred))  +
    labs(title = "KED", col = "log10(Zn)") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
gridExtra::grid.arrange(g.ok, g.ked, nrow=1)


## ----rk-ggplot-kr40-var, out.width='0.65\\textwidth', tidy = FALSE----
plot(kr40["var1.var"],
     nbreaks = 64, pch = 15,
     pal = cm.colors,
     main="KED-ffreq prediction variance, log10-ppm Zn")


## -------------------------------------------------------------------
summary(kr40$var1.var)
summary(k40$var1.var)


## -------------------------------------------------------------------
zmax <- round(max(k40$var1.var,kr40$var1.var), 3) + 0.001
zmin <- round(min(k40$var1.var,kr40$var1.var), 3) - 0.001
(ramp <- seq(from=zmin, to=zmax, by=.005))


## ----plot-compare-ok-ked-var, out.width='\\textwidth', fig.width=10, fig.height=5, tidy = FALSE----
g.ok <- ggplot() + 
    geom_sf(data = k40, aes(col = var1.var)) +
    labs(title = "OK", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = cm.colors(length(ramp)))
g.ked <-ggplot() + 
    geom_sf(data = kr40, aes(col = var1.var))  +
    labs(title = "KED flood frequency", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = cm.colors(length(ramp)))
gridExtra::grid.arrange(g.ok, g.ked, nrow=1)


## -------------------------------------------------------------------
names(meuse.grid.sf)
intersect(names(meuse.sf), names(meuse.grid.sf))


## ----rk-ggplot-meuse-grid-dist, tidy = FALSE------------------------
plot(meuse.grid.sf["dist"],
     pal = topo.colors, pch = 15, nbreaks = 64,
     main = "Normalized distance to river")


## ----rk-plot-Zn-dist, out.width='0.65\\textwidth'-------------------
plot(logZn ~ dist, data=meuse.sf, 
     col=meuse.sf$ffreq, pch = 20)
legend(x=0.5, y=3.2, 
       legend=c("1 = once in 2 years",
                "2  = once in 10 years","3  = once in 50 years"), 
       pch=20, col=1:3)


## ----rk-m-lzn-dist--------------------------------------------------
m.lzn.dist <- lm(logZn ~ dist, data=meuse.sf)
summary(m.lzn.dist)


## ----rk-ggplot-kdist-pred, out.width='0.55\\textwidth', tidy = FALSE----
k.dist <- krige(logZn ~ dist, locations=meuse.sf,
                 newdata=meuse.grid.sf, model=NULL)


## ----rk-ggplot-kdist-pred-cf, out.width='\\textwidth', fig.width=10, fig.height=5, tidy = FALSE----
g.idw <- ggplot() + 
    geom_sf(data = k.dist, aes(col = var1.pred), pch = 15) +
    labs(title = "Prediction", col = "log10(Zn)") +
    scale_color_gradientn(colors = sp::bpy.colors(64))
g.idw.v <-ggplot() + 
    geom_sf(data = k.dist, aes(col = var1.var), pch = 15)  +
    labs(title = "Prediction variance", col = "log10(Zn)^2") +
    scale_color_gradientn(colors = cm.colors(64))
gridExtra::grid.arrange(g.idw, g.idw.v, nrow=1)


## ----rk-m-lzn-ff-dist, tidy = FALSE---------------------------------
m.lzn.ff.dist <- lm(logZn ~ ffreq + dist, data=meuse.sf)
m.lzn.ff.dist.i <- lm(logZn ~ ffreq * dist, data=meuse.sf)
anova(m.lzn.ff.dist.i, m.lzn.ff.dist, m.lzn.dist, m.lzn.ff)


## ----rk-AIC---------------------------------------------------------
AIC(m.lzn.dist, m.lzn.ff, m.lzn.ff.dist, m.lzn.ff.dist.i)


## -------------------------------------------------------------------
summary(m.lzn.ff.dist.i)


## ----rk-plot-m-lzn-ff-dist-i, out.width='\\textwidth', fig.width=12, fig.height=4----
par(mfrow=c(1,3))
plot(m.lzn.ff.dist.i, which=c(1,2,5))
par(mfrow=c(1,1))


## ----rk-find-worst-row-name-----------------------------------------
(ix <- which(row.names(meuse.sf) == "76"))
meuse.sf[ix,]


## ----rk-find-bad-rows-----------------------------------------------
# which row numbers correspond to the observations witH large residuals?
(bad.pt <- which(row.names(meuse.sf) %in% c("76","51","157")))
# where are they?
st_coordinates(meuse.sf)[bad.pt, ]
# make a logical vector of all rows, whether they have large residuals
#  or not
is.row.bad <- (row(meuse.sf)[,1] %in% bad.pt)


## ----rk-show-bad-pts, out.width='0.8\\textwidth', tidy = FALSE------
colours.ffreq = c("red","orange","green")
plot(st_coordinates(meuse.sf), asp=1,
     col=colours.ffreq[meuse.sf$ffreq],
     # select print character, large residual or not?
     # 20 = filled circle; 1 = open circle
     pch=ifelse(is.row.bad, 20, 1),
     # symbol size proportional to Zn concentration
     cex=4*meuse.sf$zinc/max(meuse.sf$zinc),
     main="Suspicious regression residuals (solid circles)",
     sub="Symbol size proportional to Zn concentration")
grid()
legend(178000, 333000, pch=1, col=colours.ffreq,
       legend=c("Frequent", "Occasional", "Rare"))
text(st_coordinates(meuse.sf)[bad.pt,],
    c("51","76","157"), pos=4)


## -------------------------------------------------------------------
rm(bad.pt, is.row.bad)


## ----rk-plot-vr2, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
(vr2 <- variogram(logZn ~ ffreq*dist, location=meuse.sf, 
                  cutoff=1300, width=90))
print(plot(vr2, plot.numbers=T, pch = 20,
           main="Residuals, ffreq*dist"))


## ----rk-vrm2f-------------------------------------------------------
(vrm2f <- fit.variogram(vr2, vgm(psill=0.04, model="Sph", 
                                 range=700, nugget=0.01)))
print(vrmf)
print(vmf)


## ----rk-plot-vrm2f, fig.width=8, fig.height=6, out.width='0.65\\textwidth'----
print(plot(vr2, plot.numbers=T, model=vrm2f, main="Residuals, ffreq*dist"))


## ----rk-kr240-------------------------------------------------------
kr240 <- krige(logZn ~ ffreq*dist, locations=meuse.sf, 
               newdata=meuse.grid.sf, model=vrm2f)


## ----rk-plot-kr240, tidy = FALSE------------------------------------
plot(kr240["var1.pred"], pch = 15, nbreaks = 64,
             pal = sp::bpy.colors,
             main="KED-ffreq*dist prediction, log-ppm Zn")


## ----rk-plot-k40kr40kr240-pred, out.width='\\textwidth', fig.width=12, fig.height=4----
zmax <- round(max(k40$var1.pred,
                  kr40$var1.pred,
                  kr240$var1.pred), 1) + 0.1
zmin <- round(min(k40$var1.pred,
                  kr40$var1.pred,
                  kr240$var1.pred), 1) - 0.1
ramp <- seq(from=zmin, to=zmax, by=.1)
g.ok <- ggplot() + 
    geom_sf(data = k40, aes(col = var1.pred)) +
    labs(title = "OK", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
g.ked <-ggplot() + 
    geom_sf(data = kr40, aes(col = var1.pred))  +
    labs(title = "KED flood frequency", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
g.ked.2 <-ggplot() + 
    geom_sf(data = kr240, aes(col = var1.pred))  +
    labs(title = "KED flood frequency  * distance", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
gridExtra::grid.arrange(g.ok, g.ked, g.ked.2, nrow=1)


## -------------------------------------------------------------------
summary(kr240$var1.var)
summary(kr40$var1.var)
summary(k40$var1.var)


## ----tidy = FALSE---------------------------------------------------
zmax <- round(max(k40$var1.var,
                  kr40$var1.var,
                  kr240$var1.var), 3) + 0.001
zmin <- round(min(k40$var1.var,
                  kr40$var1.var,
                  kr240$var1.var), 3) - 0.001
ramp <- seq(from=zmin, to=zmax, by=.005)


## ----rk-plot-k40kr40kr240-var, out.width='\\textwidth', fig.width=12, fig.height=4----
g.ok <- ggplot() + 
    geom_sf(data = k40, aes(col = var1.var)) +
    labs(title = "OK", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = cm.colors(length(ramp)))
g.ked <-ggplot() + 
    geom_sf(data = kr40, aes(col = var1.var))  +
    labs(title = "KED flood frequency", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = cm.colors(length(ramp)))
g.ked.2 <-ggplot() + 
    geom_sf(data = kr240, aes(col = var1.var))  +
    labs(title = "KED flood frequency*distance", col = "log10(Zn)^2") +
    scale_color_gradientn(limits = range(ramp),
                         colors = cm.colors(length(ramp)))
gridExtra::grid.arrange(g.ok, g.ked, g.ked.2, nrow=1)


## ----child-gs_rkgls, child='gs_rkgls_sf.Rnw'------------------------

## ----rkgls-set-parent, echo=FALSE, cache=FALSE----------------------
set_parent('gs_intro_sf.Rnw')


## ----rkgls-compute-prop-nugget--------------------------------------
vrm2f[2,"range"]
(prop.nugget <- vrm2f[1,"psill"]/sum(vrm2f[,"psill"]))


## ----rkgls-fit-reml, tidy = FALSE-----------------------------------
library(nlme)
m.gls <- gls(model=logZn ~ ffreq * dist,
             data=meuse,
             correlation=corSpher(
                 form=~x + y,
                 nugget=TRUE,
                 value=c(vrm2f[2,"range"], prop.nugget),
                 ))


## ----rkgls-summary-reml---------------------------------------------
summary(m.gls)


## ----rkgls-compare-ols-gls-coeff, tidy = FALSE----------------------
coef(m.gls); coef(m.lzn.ff.dist.i)
# percent change
round(100*(coefficients(m.gls)
           - coefficients(m.lzn.ff.dist.i))/
      coefficients(m.lzn.ff.dist.i),2)


## ----rkgls-ci-reml--------------------------------------------------
intervals(m.gls, level=0.95)$coef


## ----rkgls-summary-intervals----------------------------------------
round(intervals(m.gls, level=0.95)$corStruct,4)
print(paste("range:", round(vrm2f[2,"range"], 4), "m" ))
print(paste("proportional nugget", round(prop.nugget, 4)))


## ----rkgls-plot-gls-fits, tidy = FALSE------------------------------
plot(meuse$logZn ~ predict(m.gls),
     col=meuse$ffreq, pch=20, asp=1,
     xlab="Fitted by GLS",
     ylab="Actual",
     main="log10(Zn), ppm")
legend("topleft", levels(meuse$ffreq), pch=20, col=1:4)
grid()
abline(0,1)


## ----rkgls-bubble-diffs, fig.height=6, asp=1------------------------
meuse.sf$diff.gls.ols.resid <- residuals(m.gls) - residuals(m.lzn.ff.dist.i)
summary(meuse.sf$diff.gls.ols.resid)
ggplot(data = meuse.sf) +
    geom_sf(aes(size = abs(diff.gls.ols.resid),
                fill = ifelse(diff.gls.ols.resid > 0, "green", "red")),
            pch = 21) +
    labs(title = "GLS residuals - OLS residuals",
         size = "absolute difference", fill = "") +
    scale_fill_discrete(labels=c('+', '-'))


## ----rkgls-predict-gls----------------------------------------------
meuse.grid.sf$ols.pred <- predict(m.lzn.ff.dist.i, newdata=meuse.grid.sf)
meuse.grid.sf$gls.pred <- predict(m.gls, newdata=meuse.grid.sf)
summary(meuse.grid.sf$ols.pred)
summary(meuse.grid.sf$gls.pred)
meuse.grid.sf$diff.ols.gls.pred <- meuse.grid.sf$ols.pred - meuse.grid.sf$gls.pred
summary(meuse.grid.sf$diff.ols.gls.pred)


## ----rkgls-show-ols-gls-compare, tidy = FALSE-----------------------
zmax <- round(max(meuse.grid.sf$ols.pred,
                  meuse.grid.sf$gls.pred), 1) + 0.1
zmin <- round(min(meuse.grid.sf$ols.pred,
                  meuse.grid.sf$gls.pred), 1) - 0.1
ramp <- seq(from=zmin, to=zmax, by=.1)


## ----rkgls-spplot-ols-gls-diff, out.width='\\textwidth', fig.width=12, fig.height=4----
g.ols <- ggplot() + 
    geom_sf(data = meuse.grid.sf, aes(col = ols.pred)) +
    labs(title = "OLS", col = "log10(Zn)") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
g.gls <-ggplot() + 
    geom_sf(data = meuse.grid.sf, aes(col = gls.pred))  +
    labs(title = "GLS", col = "log10(Zn)") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
gridExtra::grid.arrange(g.ols, g.gls, nrow=1)


## ----rkgls-build-gls-resid-vgm--------------------------------------
meuse.sf$gls.resid <- residuals(m.gls)
(p.nugget <- intervals(m.gls)$corStruct["nugget","est."])
(t.sill <- var(meuse.sf$gls.resid))
(nugget <- t.sill * p.nugget)
(vmf.r.gls <- vgm(psill=t.sill-nugget,
                      model="Sph",
                      range=intervals(m.gls)$corStruct["range","est."],
                      nugget=nugget))
print(vrm2f)  


## ----rkgls-gls-resid-vgm, fig.width=8, fig.height=6, tidy = FALSE----
v.r.gls <- variogram(gls.resid ~ 1,
                     loc=meuse.sf, cutoff=1500, width=90)
# panel function to also show variogram and fitted model from OLS residuals
mypanel <- function(x, y, ...) {
    vgm.panel.xyplot(x, y, plot.numbers=TRUE, ...)
    panel.pointPairs(vr2$dist, vr2$gamma, col="red", pch=20)
    lattice::panel.lines(variogramLine(vrm2f, maxdist=1500), lty=2, col='red')
    }
plot(v.r.gls, pl=T, model=vmf.r.gls, 
     main="Variogram model fitted to GLS residuals",
     xlab = "separation",
     sub = "red: from OLS residuals; blue: from GLS residuals",
     panel = mypanel)


## ----rkgls-gls-resid-krige, tidy = FALSE----------------------------
k.gls.r <- krige(gls.resid ~ 1, loc=meuse.sf,
                 newdata=meuse.grid.sf, model=vmf.r.gls)
summary(k.gls.r)
k.gls.r$rk.gls.pred <- 
    meuse.grid.sf$gls.pred + k.gls.r$var1.pred
k.gls.r$diff.rk.gls.ked <- 
    k.gls.r$rk.gls.pred - kr240$var1.pred
summary(k.gls.r$diff.rk.gls.ked)
zmax <- round(max(k.gls.r$rk.gls.pred,
                  kr240$var1.pred), 1) + 0.1
zmin <- round(min(k.gls.r$rk.gls.pred,
                  kr240$var1.pred), 1) - 0.1
ramp <- seq(from=zmin, to=zmax, by=.1)


## ----rkgls-spplot-rkgls-ked-diff, out.width='\\textwidth', fig.width=10, fig.height=5----
g.glsrk <- ggplot() + 
    geom_sf(data = k.gls.r, aes(col = rk.gls.pred)) +
    labs(title = "RK/GLS prediction", col = "log10(Zn)") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
g.ked <-ggplot() + 
    geom_sf(data = kr40, aes(col = var1.pred))  +
    labs(title = "KED flood frequency", col = "log10(Zn)") +
    scale_color_gradientn(limits = range(ramp),
                         colors = sp::bpy.colors(length(ramp)))
gridExtra::grid.arrange(g.glsrk, g.ked, nrow=1)


## -------------------------------------------------------------------
g.glsrk.ked <- ggplot() + 
    geom_sf(data = k.gls.r, aes(col = diff.rk.gls.ked))  +
    labs(title = "Difference RK/GLS - KED", col = "log10(Zn)") +
    scale_color_gradientn(colors = topo.colors(64))
plot(g.glsrk.ked)


## ----child-gs_cv, child='gs_cv_sf.Rnw'------------------------------

## ----cv-set-parent, echo=FALSE, cache=FALSE-------------------------
set_parent('gs_intro_sf.Rnw')


## ----results='hide'-------------------------------------------------
kcv.ok <- krige.cv(logZn ~ 1, locations=meuse.sf, model=vmf)
kcv.rk <- krige.cv(logZn ~ ffreq, locations=meuse.sf, model=vrmf)
kcv.rk2 <- krige.cv(logZn ~ ffreq*dist, locations=meuse.sf, model=vrm2f)


## ----cv-result------------------------------------------------------
class(kcv.ok)
summary(kcv.ok)


## ----cv-resid-------------------------------------------------------
summary(kcv.ok$residual)
summary(kcv.rk$residual)
summary(kcv.rk2$residual)


## ----cv-rmse--------------------------------------------------------
sqrt(sum(kcv.ok$residual^2)/length(kcv.ok$residual))
sqrt(sum(kcv.rk$residual^2)/length(kcv.rk$residual))
sqrt(sum(kcv.rk2$residual^2)/length(kcv.rk$residual))


## -------------------------------------------------------------------
zmax <- round(max(abs(kcv.ok$residual), abs(kcv.rk$residual),
                   abs(kcv.rk2$residual)),2) + 0.01


## ----make-bubble----------------------------------------------------
cv.bubble.plot <- function(cv.result, kriging.type) {
    g <- ggplot(data = cv.result) +
        geom_sf(aes_(size = abs(cv.result$residual),
                    fill = ifelse(cv.result$residual > 0, "green", "red")),
                pch = 21) +
        labs(title = paste("X-validation,", kriging.type),
             size = "absolute residual", fill = "") +
        scale_fill_discrete(labels=c('+', '-')) +
        scale_size_continuous(limits = c(0, zmax))
    
    return(g)
}


## ----cv-plot-kcv, out.width='\\textwidth', fig.width=10, fig.height=5----
g.cv.ok <- cv.bubble.plot(kcv.ok, "OK")
g.cv.rk <- cv.bubble.plot(kcv.rk, "KED (univariate)")
g.cv.rk2 <- cv.bubble.plot(kcv.rk2, "KED (multivariate)")
gridExtra::grid.arrange(g.cv.ok, g.cv.rk, g.cv.rk2, nrow=1)


## -------------------------------------------------------------------
rm(zmax, g.cv.ok, g.cv.rk, g.cv.rk2, cv.bubble.plot)


## ----child-gs_gam, child='gs_gam_sf.Rnw'----------------------------

## ----gam-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('gs_intro_sf.Rnw')


## ----gam-load-data, echo=FALSE, results='hide'----------------------
require(sf)
require(gstat)
if (!exists("meuse")) {  
    data(meuse)
    meuse.sf <- meuse
    st_coordinates(meuse) <- c("x", "y")
}
if (!exists("meuse$logZn"))  {
    meuse$logZn <- log10(meuse$zinc)
    meuse.sf$logZn <- meuse$logZn
}


## ----gam-load-mgcv--------------------------------------------------
library(mgcv)


## ----gam-plot-Zn-elev-dist, out.width='\\textwidth', fig.width=10, fig.height=5----
library(ggplot2)
p1 <- ggplot(data = meuse, aes(dist, logZn)) +
    geom_point() +
    geom_smooth(method="loess") +
    geom_rug()  +
    labs(x = "normalized distance to river")
p2 <- ggplot(data = meuse, aes(elev, logZn)) +
    geom_point() +
    geom_smooth(method="loess")  +
    geom_rug() +
    labs(x = "elevation m.a.s.l.")
require(gridExtra)
grid.arrange(p1, p2, ncol=2)


## ----gam-smooth-one-add---------------------------------------------
library(mgcv)
m.g.dist <- gam(logZn ~ s(dist), data=meuse)
summary(m.g.dist)
summary(residuals(m.g.dist))
m.g.elev <- gam(logZn ~ s(elev), data=meuse)
summary(m.g.elev)
summary(residuals(m.g.elev))


## ----gam-plot-gam, out.width='\\textwidth', fig.width=10, fig.height=5----
par(mfrow=c(1,2))
plot.gam(m.g.dist, residuals=T, pch=20)
abline(h=0, lty=2)
plot.gam(m.g.elev, residuals=T, pch=20)
abline(h=0, lty=2)
par(mfrow=c(1,1))


## ----gam-scatter-elev-dist, fig.width=8, fig.height=6, out.width='0.65\\linewidth'----
ggplot(data = meuse, aes(dist, elev)) +
    geom_point() +
    geom_smooth() +
    geom_rug()
cor(meuse$elev, meuse$dist, method='pearson')
cor(meuse$elev, meuse$dist, method='spearman')


## ----build-gam------------------------------------------------------
m.g.dist.elev <- gam(logZn ~ s(dist) + s(elev),
                     data=meuse)
m.g.dist.elev.i <- gam(logZn ~ s(dist) + s(elev) + ti(dist,elev),
                       data=meuse)
summary(m.g.dist.elev)
summary(m.g.dist.elev.i)
summary(residuals(m.g.dist.elev))
summary(residuals(m.g.dist.elev.i))


## ----gam-m-g-dist-elev, out.width='\\textwidth', fig.width=10, fig.height=5----
plot.gam(m.g.dist.elev, pages=1)


## ----gam-m-g-dist-elev-i, out.width='\\textwidth', fig.width=10, fig.height=10----
plot.gam(m.g.dist.elev.i, pages=1)


## ----gam-vis-gam, tidy = FALSE, out.width="\\linewidth"-------------
vis.gam(m.g.dist.elev, theta=+60,
        plot.type="persp", color="terrain")


## ----gam-vis-gam-se, tidy = FALSE, out.width="\\linewidth"----------
vis.gam(m.g.dist.elev, theta=+60, color="terrain", se=1.96)


## ----gam-spplot-k-g-i, out.width='\\textwidth', fig.width=10, fig.height=5, tidy = FALSE----
tmp <- predict.gam(object=m.g.dist.elev.i, newdata=meuse.grid.sf, se.fit=TRUE)
names(tmp)
meuse.grid.sf$k.g.i <- tmp$fit
meuse.grid.sf$k.g.i.se <- tmp$se.fit

g.gam <- ggplot() + 
    geom_sf(data = meuse.grid.sf, aes(col = k.g.i)) +
    labs(title = "GAM prediction", col = "log10(Zn)") +
    scale_color_gradientn(colors = bpy.colors(64))
g.gam.se <- ggplot() +
    geom_sf(data = meuse.grid.sf, aes(col = k.g.i.se)) +
    labs(title = "GAM prediction standard error", col = "log10(Zn)") +
    scale_color_gradientn(colors = cm.colors(64))
grid.arrange(g.gam, g.gam.se, ncol=2)


## -------------------------------------------------------------------
summary(m.dist.elev.i <- lm(logZn ~ dist*elev, data=meuse))


## ----gam-m-dist-elev-i, out.width='\\textwidth', fig.width=12, fig.height=4----
par(mfrow=c(1,3))
plot(m.dist.elev.i, which=c(1,2,5))
par(mfrow=c(1,1))


## ----gam-plot-k-i, out.width='\\textwidth', fig.width=10, fig.height=5, tidy = FALSE----
tmp <- predict.lm(object=m.dist.elev.i,
                  newdata=meuse.grid.sf, se.fit=TRUE)
meuse.grid.sf$k.i <- tmp$fit
meuse.grid.sf$k.i.se <- tmp$se.fit

g.lm <- ggplot() + 
    geom_sf(data = meuse.grid.sf, aes(col = k.i)) +
    labs(title = "OLS prediction", col = "log10(Zn)") +
    scale_color_gradientn(colors = bpy.colors(64))
g.lm.se <- ggplot() +
    geom_sf(data = meuse.grid.sf, aes(col = k.i.se)) +
    labs(title = "OLS prediction standard error", col = "log10(Zn)") +
    scale_color_gradientn(colors = cm.colors(64))
gridExtra::grid.arrange(g.lm, g.lm.se, ncol=2)


## ----gam-highlev-hires, out.width='0.55\\textwidth', fig.width=4, fig.height=6----
ix <- which(row.names(meuse)=="76")
meuse[ix,c("zinc","elev","dist")]
log10(meuse[ix,"zinc"])
fitted(m.dist.elev.i)[ix]
fitted(m.g.dist.elev.i)[ix]
plot(st_coordinates(meuse.sf), asp=1, pch=21, cex=4*meuse$zinc/max(meuse$zinc),
     bg=ifelse(row.names(meuse)=="76","red","gray"))
data(meuse.riv)
lines(meuse.riv)
grid()


## ----gam-diff-lm-gam, tidy = FALSE----------------------------------
meuse.grid.sf$diff.gam.lm <- meuse.grid.sf$k.g.i - meuse.grid.sf$k.i
plot(meuse.grid.sf["diff.gam.lm"],
     pal = topo.colors, pch = 15,
     main="Difference, GAM prediction - linear prediction")


## ----gam-plot-gam-resid, tidy = FALSE-------------------------------
meuse.sf$resid.gam <- residuals(m.g.dist.elev.i)
ggplot(data = meuse.sf) +
    geom_sf(aes(size = abs(resid.gam),
                fill = ifelse(resid.gam > 0, "green", "red")),
            pch = 21) +
    labs(title = "GAM residuals",
         size = "absolute residual", fill = "") +
    scale_fill_discrete(labels=c('+', '-'))


## ----gam-plot-gam-resid-2, fig.width=8, fig.height=6, out.width='0.65\\textwidth', tidy = FALSE----
vr <- variogram(resid.gam ~ 1, locations=meuse.sf)
print(plot(vr, plot.numbers=T, main = "GAM residuals",
      pch = 20, xlab = "separation"))


## -------------------------------------------------------------------
max(vr$gamma)/max(variogram(logZn ~ 1, locations=meuse.sf)$gamma)


## ----echo=FALSE, results='hide'-------------------------------------
(vrmf <- fit.variogram(vr, model=vgm(0.15, "Sph", 800, 0.05)))
tmp <- krige(resid.gam ~ 1, locations=meuse.sf, newdata=meuse.grid.sf, model=vrmf, beta=0)
meuse.grid.sf$gam.r.sk <- tmp$var1.pred
meuse.grid.sf$gam.r.sk.var <- tmp$var1.var
meuse.grid.sf$rk.gam <- meuse.grid.sf$k.g.i + meuse.grid.sf$gam.r.sk


## ----gam-spplot-gam-r-sk, out.width='0.65\\textwidth', echo=FALSE----
plot(meuse.grid.sf["gam.r.sk"], pch = 15, nbreaks = 64,
     main = "GAM residuals", 
     pal = rainbow, reset = FALSE)
plot(meuse.sf["resid.gam"], pch=20, add = TRUE,
     cex=2*abs(meuse.sf$resid.gam)/max(abs(meuse.sf$resid.gam)),
              col=ifelse(meuse.sf$resid.gam > 0,"green", "red"))


## ----gam-pplot-meuse-rk-gam, out.width='0.65\\textwidth', echo=FALSE----
plot(meuse.grid.sf["rk.gam"], pch = 15, nbreaks = 64,
     main = "RK using GAM")



## ----child-gs_ans, child='gs_ans_sf.Rnw'----------------------------

## ----ans-set-parent, echo=FALSE, cache=FALSE------------------------
set_parent('gs_intro_sf.Rnw')


## ----compute-ffreq-Zn-mean------------------------------------------
round(tapply(meuse$logZn, meuse$ffreq, mean),3)


## ----m.ff.rpp-compute.success, echo=FALSE, results='hide'-----------
tmp.ff <- m.ff.rpp$frame[m.ff.rpp$frame$var=="<leaf>", c("n", "dev", "yval") ]
tmp.ff$correct <- tmp.ff$n - tmp.ff$dev
tmp.ff$prop <- (tmp.ff$correct)/tmp.ff$n
tmp.ff1 <- tmp.ff[tmp.ff$yval==1, ]
tmp.ff2 <- tmp.ff[tmp.ff$yval==2, ]
tmp.ff3 <- tmp.ff[tmp.ff$yval==3, ]


## ----show.lowest.pred.variance.distance-----------------------------
ix <- which.min(k.dist$var1.var)
k.dist[ix,]
meuse.grid[ix, "dist"]
mean(meuse$dist)