### R code from vignette source 'exA2.Rnw'

###################################################
### code chunk number 1: exA2.Rnw:18-19
###################################################
require(sp); require(gstat)


###################################################
### code chunk number 2: exA2.Rnw:22-23
###################################################
load("./exA1.RData")


###################################################
### code chunk number 3: exA2.Rnw:75-76
###################################################
vm <- vgm(psill = 0.8, , model = "Exp", range = 8, nugget = 0.2)


###################################################
### code chunk number 4: exA2.Rnw:114-121
###################################################
set.seed(316318)
scheme <- data.frame(round(kmeans(x=coordinates(grid.128),centers=2^4,
                                  iter.max=1000,nstart=2^4,
                                  algorithm="Hartigan-Wong")$centers))
names(scheme) <- c("x","y"); coordinates(scheme) <- ~x + y;
pts.z <- overlay(k.e8.128, scheme)
summary(pts.z)


###################################################
### code chunk number 5: exA2.Rnw:129-133
###################################################
layout.2 <- list("sp.points", pts.z, pch = 21, cex = 1.5, col = "black", fill = "white")
print(spplot(k.e8.128, col.regions=bpy.colors(64),
               sp.layout=list(layout.2)
              ))


###################################################
### code chunk number 6: exA2.Rnw:146-154
###################################################
vc <- variogram(z ~ 1, data=pts.z, cloud=T,
                cutoff=128*sqrt(2))
mean(vc$gamma)
trellis.par.set(myTheme)
print(plot(
    vc,
    main="Variogram cloud, 16 grid points in 128 x 128 grid",
    model=vm))


###################################################
### code chunk number 7: exA2.Rnw:174-189
###################################################
set.seed(1768)
dv <- vector(mode="numeric", length=128)
for (i in 1:128) {
  pts <- spsample(k.e8.128, 4^2, type="regular")
  pts.z <- overlay(k.e8.128, pts)
  dv[i] <- mean(variogram(z ~ 1, data=pts.z,
                          cutoff=128*sqrt(2),
                          cloud=T)$gamma)
}
summary(dv)
hist(
    dv,
    main="Dispersion variances for different 4 x 4 discretizations")
print(paste("Naive variance:",
            round(var(k.e8.128$z),4)))


###################################################
### code chunk number 8: exA2.Rnw:204-206
###################################################
dv <- matrix(0, nrow=128, ncol=6)
colnames(dv) <- c("4x4", "6x6", "8x8", "10x10", "12x12", "16x16")


###################################################
### code chunk number 9: exA2.Rnw:212-224
###################################################
set.seed(1814)
n <- c(4^2, 6^2, 8^2, 10^2, 12^2, 16^2)
for (p in 1:6) {
  for (i in 1:128) {
    pts <- spsample(k.e8.128, n[p], type="regular")
    pts.z <- overlay(k.e8.128, pts);
    v <- variogram(z ~ 1, data=pts.z, width=1, cutoff=128*sqrt(2));
    dv[i,p] <- sum(v$np*v$gamma)/sum(v$np);
}}
summary(dv)
apply(dv, 2, sd)
apply(dv, 2, function(x) 100*sd(x)/mean(x))


###################################################
### code chunk number 10: exA2.Rnw:251-254
###################################################
bbox(k.e8.128)
k.e8.64 <- k.e8.128[1:64,1:64]
bbox(k.e8.64)


###################################################
### code chunk number 11: exA2.Rnw:261-268
###################################################
set.seed(250)
repeat {
  try(scheme <- spsample(k.e8.64, n=128, type="random"), silent=T)
  if (class(.Last.value) != "try-error") break
}
plot(coordinates(scheme))
grid()


###################################################
### code chunk number 12: exA2.Rnw:273-275
###################################################
pts.z <- overlay(k.e8.128, scheme)
summary(pts.z)


###################################################
### code chunk number 13: exA2.Rnw:280-282
###################################################
v.e <- variogram(z ~ 1, pts.z)
print(plot(v.e, pl=T, model=vm))


###################################################
### code chunk number 14: exA2.Rnw:294-297
###################################################
print(vm)
(v.emf <- fit.variogram(v.e, model=vgm(1, "Exp", 15/3, .2)))
print(plot(v.e, pl=T, model=v.emf))


###################################################
### code chunk number 15: exA2.Rnw:314-315
###################################################
(dv <- mean(variogram(z ~ 1, data=pts.z, cutoff=8*sqrt(2), cloud=T)$gamma))


###################################################
### code chunk number 16: exA2.Rnw:324-329
###################################################
(v.emf.reg <- v.emf)
v.emf.reg[1,"psill"] <- 0
v.emf.reg[2,"psill"] <- v.emf[2,"psill"] - dv + v.emf[1,"psill"]
v.emf.reg[2,"range"] <- v.emf.reg[2,"range"] + (8/3)
print(v.emf.reg)


###################################################
### code chunk number 17: exA2.Rnw:343-344
###################################################
print(vmf.16)