library(sgeostat)
library(akima)
data(maas)

fivenum(maas$zinc)
zcl <- cut(maas$zinc,fivenum(maas$zinc),include.lowest=T)
plot(maas$x,maas$y, type = "n", ann = F, bty = "n", xaxt = "n", yaxt = "n")
axis(1); axis(2)
points(maas$x,maas$y, pch=codes(zcl))
legend(178700,333500,legend=levels(zcl), pch=c(1,2,3,4))

#ael: was missing:
maas.pts <- point(maas, x='x', y='y')  

maas.prs <- pair(maas.pts, num.lags=10, maxdist=2000)
spacebox(maas.pts, maas.prs, 'zinc')
maas.v <- est.variogram(maas.pts, maas.prs, 'zinc')
plot(maas.v)

#ael: made these plot.it calls FALSE - otherwise we get too many plots!
maas.vmod.spher <- fit.spherical(maas.v, c0=60000, cs=110000,as=800,iterations=25, plot.it=F)
maas.vmod.gauss <- fit.gaussian(maas.v, c0=60000, cg=110000, ag=800,
                      iterations=25, plot.it=F)
plot(maas.v, maas.vmod.spher)
plot(maas.v, maas.vmod.gauss)

#ael: use length=20, rather than by=100:
grid <- list(x=seq(min(maas$x), max(maas$x), length=20),
           y=seq(min(maas$y), max(maas$y), length=20))

grid$xy <- data.frame(cbind(
		c(matrix(grid$x, length(grid$x), length(grid$x))),
		c(matrix(grid$y, length(grid$y), length(grid$y), byrow=T))))
colnames(grid$xy) <- c("x", "y")

grid$pts <- point(grid$xy, x='x', y='y')
grid$krige.spher <- krige(grid$pts, maas.pts, "zinc", maas.vmod.spher,
                          maxdist=1000, extrap=F)
image(grid$x,grid$y,matrix(grid$krige.spher$zhat, length(grid$x), length(grid$y)))
contour(grid$x, grid$y,
          matrix(grid$krige.spher$zhat, length(grid$x), length(grid$y)),add=T)
points(maas$x,maas$y, pch=codes(zcl))
legend(178700,333500,legend=levels(zcl), pch=c(1,2,3,4))