## ----setup, include=FALSE, cache=FALSE-----------------------------------
## ----fileShow, eval=FALSE------------------------------------------------
## file.show("sandford.txt")

## ----readTable1, eval=FALSE----------------------------------------------
## ds <- read.table("sandford.txt", skip=2)
## names(ds) <- c("seq","clay1","silt1","clay2","silt2","clay3","silt3")

## ----readTable2, echo=FALSE----------------------------------------------
ds <- read.table("./ds/sandford/sandford.txt", skip=2)
names(ds) <- c("seq","clay1","silt1","clay2","silt2","clay3","silt3")
str(ds)

## ----showTransect, fig.width=10, fig.height=5, out.width="\\linewidth"----
plot(ds$clay2, type="l", ylim=c(0,100), xlim=c(0,320), main="Layer 2", xlab="Station on transect", ylab="weight %")
lines(ds$silt2, lty=2, col="blue")
legend("topright", c("clay","silt"), lty=1:2, col=c("black","blue"))
grid()

## ----sBasis--------------------------------------------------------------
require(splines)
(n <- length(ds$clay2))
basis.10 <- ns(ds$seq, df=floor(n/10))
str(basis.10)

## ------------------------------------------------------------------------
for (i in 148:152) {
    print(paste("Station",i))
    print(round(basis.10[i,basis.10[i,]!=0],4))
}

## ----sLm-----------------------------------------------------------------
m.10 <- lm(ds$clay2 ~ basis.10)
summary(m.10)$adj.r.squared

## ------------------------------------------------------------------------
summary(pred.10 <- fitted(m.10))
summary(ds$clay2)

## ----sPredict, fig.width=10, fig.height=5, out.width='\\linewidth'-------
plot(pred.10, type="l",ylab="clay, weight %", xlab="Position on transect")
lines(ds$clay2, col="darkgray")
legend("bottomright", lty=1,
       col=c("black","darkgray"),
       c("spline fit","actual"))

## ----sDouble-------------------------------------------------------------
basis.20 <- ns(ds$seq, df=floor(n/20))
basis.5 <- ns(ds$seq, df=floor(n/5))
m.20 <- lm(ds$clay2 ~ basis.20)
summary(m.20)$adj.r.squared
m.5 <- lm(ds$clay2 ~ basis.5)
summary(m.5)$adj.r.squared

## ----sPlot3, fig.width=10, fig.height=5, out.width='\\linewidth'---------
plot(pred.10, type="l", ylab="clay, weight %", xlab="Position on transect",
     ylim=c(-10,100))
grid()
abline(h=0, col="darkgray"); abline(h=100, col="darkgray")
lines(predict(m.20), col="blue")
lines(predict(m.5), col="darkgreen")
lines(ds$clay2, col="darkgray")
legend("bottomright", lty=1, col=c("black","blue","darkgreen","darkgray"),
       c("spline (10)","spline (20)", "spline (5)", "actual"))

## ------------------------------------------------------------------------
pts <- c(145.5,150.5,155.5)
str(ispl <- interpSpline(ds$seq, ds$clay2, bSpline=TRUE))
(pred3 <- predict(ispl, x=pts)$y)
ds$clay2[145:156]

## ----sPlot3Interp, out.width='0.7\\linewidth'----------------------------
plot(pred.10, type="l", ylab="clay, weight %",
     xlab="Position on transect", xlim=c(140,160))
grid()
lines(predict(m.20), col="blue", type="l")
lines(predict(m.5), col="darkgreen", type="l")
lines(ds$clay2, col="darkgray", type="b")
legend("topleft", lty=1, col=c("black","blue","darkgreen","darkgray","red"),
       c("spline (10)","spline (20)", "spline (5)", "actual", "interpolation"))
points(x=pts, y=pred3, col="red")
rm(pred3)

## ----fitSmooth-----------------------------------------------------------
(spl.fit <- smooth.spline(ds$clay2))
round(length(ds$seq)/spl.fit$df,1)

## ----plotSmooth, fig.width=10, fig.height=5, out.width="\\linewidth"-----
plot(ds$clay2, type="l", col="darkgray", ylim=c(0,100), xlim=c(0,320),
     xlab="Station on transect", ylab="clay, weight %")
lines(spl.fit$y)
lines(predict(m.5), col="darkgreen", type="l")
points(spl.fit$fit$knot*length(ds$clay2), rep(0, length(spl.fit$fit$knot)),
       spl.fit$nk, pch="|", col="blue")
legend("topright", c("clay","smooth spline fit", "natural spline fit, df=5"),
       lty=1, col=c("darkgray","black", "darkgreen"))
grid()

## ----plotSmoothZoom, fig.width=10, fig.height=5, out.width="\\linewidth"----
plot(ds$clay2, type="b", col="darkgray", ylim=c(0,100), xlim=c(0,100), xlab="Station on transect", ylab="clay, weight %")
lines(spl.fit$y)
lines(predict(m.5), col="darkgreen")
points(spl.fit$fit$knot*length(ds$clay2), rep(0, length(spl.fit$fit$knot)),
       spl.fit$nk, pch="|", col="blue")
legend("topright", c("clay","smoothing spline", "natural spline, df=5"),
       lty=1, col=c("darkgray","black", "darkgreen"))
grid()

## ----loadFields----------------------------------------------------------
require(fields)

## ----loadData------------------------------------------------------------
require(gstat)
data(jura)
ls(pattern="jura")

## ----strJura-------------------------------------------------------------
str(jura.pred)
str(jura.grid)

## ----makeMatrix----------------------------------------------------------
jura.pred$coords <- matrix(c(jura.pred$Xloc, jura.pred$Yloc), byrow=F, ncol=2)
str(jura.pred$coords)

## ----tps-----------------------------------------------------------------
surf.1 <-Tps(jura.pred$coords, jura.pred$Co)
class(surf.1)

## ----tps-summary---------------------------------------------------------
summary(surf.1)

## ----lambda--------------------------------------------------------------
print(surf.1$lambda)

## ----grid----------------------------------------------------------------
jura.grid$coords <- matrix(c(jura.grid$Xloc, jura.grid$Yloc), byrow=F, ncol=2)
surf.1.pred <- predict.Krig(surf.1, jura.grid$coords)
summary(as.vector(surf.1.pred))

## ----plotPredBw, fig.width=7, fig.height=7, out.width="0.8\\linewidth", tidy=FALSE----
surf.1.pred.m <- matrix(surf.1.pred,
                        length(jura.grid$coords[,1]),
                        length(jura.grid$coords[,2]),
                        byrow=F)
str(surf.1.pred.m)

require(lattice)
# plot(surf.1.pred.m, wireframe=TRUE, col.regions=bpy.colors(96))

str(surf.1.pred.m)
getwd()
pdf("TPS_wireframe.pdf", width=12, height=12)
wireframe(surf.1.pred.m,
      shade = TRUE, aspect = c(61/87, 0.4),
     light.source = c(10, 0, 10), zoom=1.1, box=F,
     scales=list(draw=F), xlab="", ylab="", zlab="")
dev.off()


plot(jura.grid$coords, pch=20, asp=1,
     cex=0.8*surf.1.pred/max(surf.1.pred),
     xlab="E", ylab="N",
     main="relative Co concentration")
grid(col="gray")

## ----plotPredCol, fig.width=7, fig.height=7, out.width="0.8\\linewidth", tidy=FALSE----
require(sp)
plot(jura.grid$coords, pch=20, asp=1, cex=.6,
     col=bpy.colors(256)[cut(surf.1.pred, 256)],
     xlab="E", ylab="N",
     main="relative Co concentration")
grid(col="gray")

## ----valTps, tidy=FALSE--------------------------------------------------
jura.val$coords <-
    matrix(c(jura.val$Xloc, jura.val$Yloc), byrow=F, ncol=2)
surf.1.val <- 
    predict.Krig(surf.1, jura.val$coords)
surf.1.val.res <- 
    (jura.val$Co - surf.1.val)
summary(as.vector(surf.1.val.res))
require(sp)
surf.1.val.res.sp <-
    SpatialPointsDataFrame(coords = jura.val$coords,
                           data=data.frame(res=surf.1.val.res))
(rmse.tps <- with(surf.1.val.res.sp@data, sum(res^2)/length(res)))
rmse.tps/median(jura.val$Co)
(me.tps <- with(surf.1.val.res.sp@data, sum(res)/length(res)))
bubble(surf.1.val.res.sp, zcol="res", pch=1, 
       main="Residuals from thin-plate spline validation",
       sub="Co concentration, mg kg-1")

## ----vgm-----------------------------------------------------------------
jura.pred.sp <- jura.pred
coordinates(jura.pred.sp) <- ~Xloc + Yloc
v <- variogram(Co ~ 1, loc=jura.pred.sp)
(vmf <- fit.variogram(v, vgm(12.5, "Pen", 1.2, 1.5)))

## ----vgmplot, tidy=FALSE, out.width='0.45\\linewidth'--------------------
plot(v, plot.numbers=T, model=vmf,
     main="Fitted variogram model, Jura Co",
     xlab="Separation (km)")

## ----vaOK----------------------------------------------------------------
jura.val.sp <- jura.val
coordinates(jura.val.sp) <- ~Xloc + Yloc
k <- krige(Co ~ 1, loc=jura.pred.sp, newdata=jura.val.sp, model=vmf)

## ----vaOKstats-----------------------------------------------------------
summary(k$var1.pred)
summary(as.vector(surf.1.val))
summary(k$var1.pred - as.vector(surf.1.val))

## ----vaOKresidstats------------------------------------------------------
summary(k$res <- (jura.val.sp$Co - k$var1.pred))
summary(as.vector(surf.1.val.res))
summary(k$res - as.vector(surf.1.val.res))

## ------------------------------------------------------------------------
(rmse.k <- with(k@data, sum(res^2)/length(res)))
print(rmse.tps)
(me.k <- with(k@data, sum(res)/length(res)))
print(me.tps)

## ----grid-ok-------------------------------------------------------------
jura.grid.sp <- jura.grid
coordinates(jura.grid.sp) <- ~Xloc + Yloc
gridded(jura.grid.sp) <- TRUE
k <- krige(Co ~ 1, loc=jura.pred.sp, newdata=jura.grid.sp, model=vmf)

## ----tps-as-spdf---------------------------------------------------------
class(k)
k$tps.pred <- surf.1.pred
str(k)

## ----vaOKgrid, out.width='0.48\\linewidth', tidy=FALSE-------------------
zlim <- c(floor(min(k$var1.pred, k$tps.pred)),
    ceiling(max(k$var1.pred, k$tps.pred)))
ats <- seq(zlim[1], zlim[2], length=256)
spplot(k, zcol="var1.pred", col.regions=bpy.colors(256), at=ats,
       main="Co concentration, mg kg-1", sub="OK")
spplot(k, zcol="tps.pred", col.regions=bpy.colors(256), at=ats,
       main="Co concentration, mg kg-1", sub="TPS")

## ----compare-interp, out.width='0.7\\linewidth', tidy=FALSE--------------
summary(k$diff <- k$tps.pred - k$var1.pred)
spplot(k, zcol="diff", col.regions=terrain.colors(64),
       main="Difference, TPS-OK", sub="Co, mg kg-1")