library(gstat)
library(maptools)
library(rgdal)
library(lattice)
# trellis.par.set(sp.theme())

# load data:
data(meuse)
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")
# load grids:
data(meuse.grid)
meuse.grid$X <- meuse.grid$x
meuse.grid$Y <- meuse.grid$y
meuse.grid$soil <- as.factor(meuse.grid$soil)
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE
fullgrid(meuse.grid) <- TRUE
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
# get values from grids to points:
meuse.ov <- overlay(meuse.grid, meuse)
meuse$X <- meuse.ov$X
meuse$Y <- meuse.ov$Y
meuse$soil <- meuse.ov$soil
names(meuse@data)

#------------------------------------------------------------
# Generate predictions
#------------------------------------------------------------

# inverse distances:
zinc.id <- krige(zinc~1, meuse, meuse.grid)
# regression on coordinates (trend):
zinc.tr <- krige(zinc~1, meuse, meuse.grid, degree=2, nmax=20)
zinc.tr$zinc.pred <- ifelse(zinc.tr$var1.pred<0, 0, zinc.tr$var1.pred)
# linear regression:
zinc.lm <- krige(log1p(zinc)~sqrt(dist), meuse, meuse.grid)
zinc.lm$zinc.pred <- expm1(zinc.lm$var1.pred)
# ordinary kriging:
vt.fit <- fit.variogram(variogram(log1p(zinc)~1, meuse), vgm(1, "Exp", 300, 1))
plot(variogram(log1p(zinc)~1, meuse), vt.fit)
zinc.ok <- krige(log1p(zinc)~1, meuse, meuse.grid, vt.fit)
zinc.ok$zinc.pred <- expm1(zinc.ok$var1.pred)
zinc.ok$svar <- zinc.ok$var1.var/var(log1p(meuse$zinc))
# GLM prediction (all possible predictors):
m.glm <- glm(zinc~dist+soil+ffreq+X+Y, meuse, family=poisson())
summary(m.glm)
p.glm <- predict(m.glm, newdata=meuse.grid, type="response", se.fit=T)
zinc.glm <- as(meuse.grid["soil"], "SpatialPointsDataFrame")
zinc.glm$var1.pred <- p.glm$fit
zinc.glm$var1.var <- p.glm$se.fit
zinc.glm$svar <- p.glm$se.fit^2/(m.glm$null.deviance/m.glm$df.null)
gridded(zinc.glm) <- TRUE
fullgrid(zinc.glm) <- TRUE
# regression-kriging:
vr.fit <- fit.variogram(variogram(log1p(zinc)~sqrt(dist), meuse), vgm(1, "Exp", 300, 1))
plot(variogram(log1p(zinc)~sqrt(dist), meuse), vr.fit)
zinc.rk <- krige(log1p(zinc)~sqrt(dist), meuse, meuse.grid, vr.fit)
zinc.rk$zinc.pred <- expm1(zinc.rk$var1.pred)

# plot all predictions next to each other:
meuse.grid$zinc.id <- zinc.id$var1.pred
meuse.grid$zinc.tr <- zinc.tr$zinc.pred
meuse.grid$zinc.lm <- zinc.lm$zinc.pred
meuse.grid$zinc.ok <- zinc.ok$zinc.pred
meuse.grid$zinc.glm <- zinc.glm$var1.pred
meuse.grid$zinc.rk <- zinc.rk$zinc.pred

spplot(meuse.grid[c("zinc.id", "zinc.tr", "zinc.lm", "zinc.ok")], col.regions=grey(rev(seq(0,1,0.025))))