# Useful reference for these plots
#
# http://r-spatial.sourceforge.net/gallery/
#

library(gstat)
library(sp)

# Data points

data(meuse)
coordinates(meuse) = c('x', 'y')

# Grid of points where we will ultimately krige, say, zinc

data(meuse.grid) 
coordinates(meuse.grid) = c('x', 'y')

# River

data(meuse.riv)
meuse.sr = SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),"meuse.riv")))
rv = list("sp.polygons", meuse.sr, 
	fill = 'lightblue')

## coloured points plot with legend in plotting area and scales:
sp.theme(set=T)
spplot(meuse, 'zinc', key.space=list(x=0.2,y=0.9,corner=c(0,1)), scales=list(draw=T), sp.layout=list(rv), do.log=T)

# The variogram cloud

plot(variogram(log(zinc) ~1, meuse, cloud=T))

# Usual and robust estimators of the variogram

usual = variogram(log(zinc) ~1, meuse)
cressie = variogram(log(zinc) ~1, meuse, cressie=T)

plot(c(usual$dist, cressie$dist), c(usual$gamma, cressie$gamma), type='n', xlab='distance', ylab='semivariance')
points(usual$dist, usual$gamma, pch=23, bg='red')
points(cressie$dist, cressie$gamma, pch=22, bg='blue')

# Directional variation in the variogram

directional = variogram(log(zinc) ~1, meuse, alpha=c(0,45,90,135))
plot(directional)

# Different variogram models

show.vgms()
show.vgms(model='Mat', kappa.range=c(0.1,0.2,0.5,1,2,5,10), max=10)

sphere.usual = fit.variogram(usual, vgm(1, 'Sph', 800, 1))
plot(directional, sphere.usual)

# Using REML

sphere.reml = fit.variogram.reml(log(zinc) ~ 1, meuse, model=vgm(1, 'Sph', 800, 1))
plot(directional, sphere.reml)

# An exponential model

exp.usual = fit.variogram(usual, model=vgm(1, 'Exp', 800, 1))
plot(directional, exp.usual)

# A twice-differentiable Matern model

mat.2 = fit.variogram(usual, model=vgm(1, 'Mat', 800, 1, kappa=2))
plot(directional, mat.2)

# A very smooth Matern model

mat.10 = fit.variogram(usual, model=vgm(1, 'Mat', 800, 1, kappa=10))
plot(directional, mat.10)

# Kriging

# Using the spherical

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, sphere.reml)
spplot(zn.ok, 'var1.pred', sp.layout=list(rv), cuts=10)


# Using the exponential

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, exp.usual)
spplot(zn.ok, 'var1.pred', sp.layout=list(rv), cuts=10)


# Using the twice-differentiable Matern

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, mat.2)
spplot(zn.ok, 'var1.pred', sp.layout=list(rv), cuts=10)


# Using the smooth Matern

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, mat.10)
spplot(zn.ok, 'var1.pred', sp.layout=list(rv), cuts=10)

# A look at the kriging variance

# For the REML fit spherical

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, sphere.reml)
spplot(zn.ok, 'var1.var', sp.layout=list(rv), cuts=10, do.log=T)

# For the smooth Matern

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, mat.10)
spplot(zn.ok, 'var1.var', sp.layout=list(rv), cuts=10, do.log=T)


# Universal kriging, predictions

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ sqrt(dist), meuse, meuse.grid, mat.2)
spplot(zn.ok, 'var1.pred', sp.layout=list(rv), cuts=10, do.log=T)

# Universal kriging, variance

sp.theme(set=T)
zn.ok = krige(log(zinc) ~ sqrt(dist), meuse, meuse.grid, mat.2)
spplot(zn.ok, 'var1.var', sp.layout=list(rv), cuts=10, do.log=T)

# Co-kriging

g = gstat(NULL, "logCd", log(cadmium) ~ 1, meuse)
g = gstat(g, "logCu", log(copper) ~ 1, meuse)
g = gstat(g, "logPb", log(lead) ~ 1, meuse)
g = gstat(g, "logZn", log(zinc) ~ 1, meuse)

vm = variogram(g)

# The fitted model of the variogram
# enforces a symmetry condition on the
# variance-covariance matrix Sigma(h) of the multivariate
# increments Z(s+h)-Z(s). Namely Sigma(h)=Sigma(-h).

vm.fit = fit.lmc(vm, g, vgm(1, "Sph", 800, 1))
plot(vm, vm.fit)

cok.maps = predict(vm.fit, meuse.grid)
sp.theme(set=T)
spplot.vcov(cok.maps, cuts=10)
spplot(cok.maps, 'logZn.var', cuts=10)

X11()
sp.theme(set=T)
zn.ok = krige(log(zinc) ~ 1, meuse, meuse.grid, sphere.usual)
spplot(zn.ok, 'var1.var', sp.layout=list(rv), cuts=10)

# Compare the kriging variances, and predicted values

print(lm(zn.ok[['var1.var']] ~ cok.maps[['logZn.var']]))
plot(zn.ok[['var1.var']], cok.maps[['logZn.var']], pch=23, bg='red')
plot(zn.ok[['var1.pred']], cok.maps[['logZn.pred']], pch=23, bg='red')