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

###################################################
### code chunk number 1: exA1.Rnw:20-24
###################################################
options(prompt="> ", continue="+ ", digits=5, width=70, show.signif.stars=T)
par(mfrow=c(1,1))
rm(list=ls())
setwd("~/data/edu/dgeostats/ex")


###################################################
### code chunk number 2: exA1.Rnw:28-32
###################################################
require(sp)
require(gstat)
require(lattice)
require(rgdal)


###################################################
### code chunk number 3: exA1.Rnw:56-58
###################################################
grid.128 <- SpatialGrid(GridTopology(cellcentre.offset=c(0.5, 0.5), cellsize=c(1, 1), cells.dim=c(128, 128)))
str(grid.128)


###################################################
### code chunk number 4: exA1.Rnw:69-70
###################################################
vm <- vgm(psill=0.8,,model="Exp",range=8,nugget=0.2)


###################################################
### code chunk number 5: exA1.Rnw:87-96 (eval = FALSE)
###################################################
## set.seed(621)
## system.time(
##             k.e8.128 <- krige(z~1, loc=NULL,
##                               newdata=grid.128,
##                               model=vm, beta=0,
##                               dummy=T, nsim=1,
##                               maxdist=24)
##             )
## names(k.e8.128) <- "z"


###################################################
### code chunk number 6: exA1.Rnw:115-116 (eval = FALSE)
###################################################
## save(k.e8.128, file="./ds/sim/exp24.RData")


###################################################
### code chunk number 7: exA1.Rnw:120-121
###################################################
load("./ds/sim/exp24.RData")


###################################################
### code chunk number 8: exA1.Rnw:130-133
###################################################
summary(k.e8.128)
sim.plot.128 <- spplot(k.e8.128, col.regions=bpy.colors(64), main='Simulated field, vgm(0.8,"Exp",8,0.2), 128x=128')
print(sim.plot.128)


###################################################
### code chunk number 9: exA1.Rnw:147-148
###################################################
var(k.e8.128$z)


###################################################
### code chunk number 10: exA1.Rnw:166-168
###################################################
k.e8.128.pts <- SpatialPointsDataFrame(SpatialPoints(k.e8.128), data=k.e8.128@data)
str(k.e8.128.pts)


###################################################
### code chunk number 11: exA1.Rnw:193-194
###################################################
(block <- subset(k.e8.128.pts,(coordinates(k.e8.128.pts)[,1] %in% c(0.5,1.5)) & (coordinates(k.e8.128.pts)[,2] %in% c(0.5,1.5))))


###################################################
### code chunk number 12: exA1.Rnw:199-200
###################################################
(vc <- variogram(z ~ 1, loc=block, cutoff=sqrt(2), cloud=T))


###################################################
### code chunk number 13: exA1.Rnw:217-218
###################################################
mean(vc$gamma)


###################################################
### code chunk number 14: exA1.Rnw:223-225
###################################################
(v <- variogram(z ~ 1, loc=block, cutoff=2*sqrt(2)))
sum(v$np*v$gamma)/sum(v$np)


###################################################
### code chunk number 15: exA1.Rnw:230-240
###################################################
bbox(k.e8.128.pts)
dv.2 <- vector(mode="numeric",length=12^2); i <- 1
for (x in seq(0.5,23.5,by=2)) {
  for (y in seq(0.5,23.5,by=2)) {
    block <- subset(k.e8.128.pts,(coordinates(k.e8.128.pts)[,1] %in% c(x,x+1)) & (coordinates(k.e8.128.pts)[,2] %in% c(y, y+1)));
    v <- variogram(z ~ 1, loc=block, cutoff=2*sqrt(2));
    dv.2[i] <- sum(v$np*v$gamma)/sum(v$np); i <- i + 1
  }
}
length(dv.2);head(dv.2)


###################################################
### code chunk number 16: exA1.Rnw:252-257
###################################################
summary(dv.2)
hist(
    dv.2,
    main="Dispersion variances of 2 x 2 blocks, upper-left 24 x 24 block",
    xlab="Dispersion variance", breaks=seq(0,1.6,by=.05))


###################################################
### code chunk number 17: exA1.Rnw:298-306
###################################################
dv.3 <- vector(mode="numeric",length=12^2); i <- 1
for (x in seq(0.5,34.5,by=3)) {
  for (y in seq(0.5,34.5,by=3)) {
    block <- subset(k.e8.128.pts,(coordinates(k.e8.128.pts)[,1] %in% c(x, x+1, x+2)) & (coordinates(k.e8.128.pts)[,2] %in% c(y, y+1, y+2)));
    v <- variogram(z ~ 1, loc=block, cutoff=2*sqrt(2));
    dv.3[i] <- sum(v$np*v$gamma)/sum(v$np); i <- i + 1
  }
}


###################################################
### code chunk number 18: exA1.Rnw:311-316
###################################################
summary(dv.3)
hist(
    dv.3,
    main="Dispersion variances of 3 x 3 blocks, upper-left 36x36 block",
    xlab="Dispersion variance", breaks=seq(0,1.4,by=.05))


###################################################
### code chunk number 19: exA1.Rnw:354-363
###################################################
dv.4 <- vector(mode="numeric",length=12^2); i <- 1
for (x in seq(0.5,45.5,by=4)) {
  for (y in seq(0.5,45.5,by=4)) {
    block <- subset(k.e8.128.pts,(coordinates(k.e8.128.pts)[,1] %in% c(x,x+1, x+2, x+3)) & (coordinates(k.e8.128.pts)[,2] %in% c(y, y+1, y+2, y +3)));
    v <- variogram(z ~ 1, loc=block, cutoff=4*sqrt(2), width=0.1);
    dv.4[i] <- sum(v$np*v$gamma)/sum(v$np); i <- i + 1
  }
}
summary(dv.4)


###################################################
### code chunk number 20: exA1.Rnw:366-370
###################################################
hist(
    dv.4,
    main="Dispersion variances of 4 x 4 blocks, upper-left 48x48 block",
    xlab="Dispersion variance", breaks=seq(0,1.4,by=.05))


###################################################
### code chunk number 21: exA1.Rnw:391-400
###################################################
dv.8 <- vector(mode="numeric",length=12^2); i <- 1
for (x in seq(0.5,88.5,by=8)) {
  for (y in seq(0.5,88.5,by=8)) {
    block <- subset(k.e8.128.pts,(coordinates(k.e8.128.pts)[,1] %in% c(x,x+1, x+2, x+3, x+4, x+5, x+6, x+7)) & (coordinates(k.e8.128.pts)[,2] %in% c(y, y+1, y+2, y +3, y+4, y+5, y+6, y+7)));
    v <- variogram(z ~ 1, loc=block, cutoff=8*sqrt(2), width=0.1);
    dv.8[i] <- sum(v$np*v$gamma)/sum(v$np); i <- i + 1
  }
}
summary(dv.8)


###################################################
### code chunk number 22: exA1.Rnw:403-407
###################################################
hist(
    dv.8,
    main="Dispersion variances of 8 x 8 blocks, upper-left 96 x 96 block",
    xlab="Dispersion variance", breaks=seq(0,2,by=.05))


###################################################
### code chunk number 23: exA1.Rnw:428-432
###################################################
dv.means <- rbind("2x2"=mean(dv.2), "3x3"=mean(dv.3), "4x4"=mean(dv.4), "8x8"=mean(dv.8))
dv.means.cells <- c(2^2, 3^2, 4^2, 8^2)
dv.means.std <- dv.means/dv.means.cells
data.frame("means"=dv.means, "std means"=dv.means.std)


###################################################
### code chunk number 24: exA1.Rnw:436-440
###################################################
par(mfrow=c(1,2))
plot(dv.means ~ dv.means.cells, type="b", main="", xlab="Number of cells", ylab="Mean dispersion variance")
plot(dv.means.std ~ dv.means.cells, type="b", main="", xlab="Number of cells", ylab="Standardized dispersion variance")
par(mfrow=c(1,2))


###################################################
### code chunk number 25: exA1.Rnw:456-466
###################################################
par(mfrow=c(2,2))
hist(dv.2, freq=F, main="Dispersion variances of 2 x 2 blocks", xlab="Dispersion variance", breaks=seq(0,1.9,by=.05), ylim=c(0,5))
abline(v=mean(dv.2), col="red")
hist(dv.3, freq=F, main="Dispersion variances of 3 x 3 blocks", xlab="Dispersion variance", breaks=seq(0,1.9,by=.05), ylim=c(0,5))
abline(v=mean(dv.3), col="red")
hist(dv.4, freq=F, main="Dispersion variances of 4 x 4 blocks", xlab="Dispersion variance", breaks=seq(0,1.9,by=.05), ylim=c(0,5))
abline(v=mean(dv.4), col="red")
hist(dv.8, freq=F, main="Dispersion variances of 8 x 8 blocks", xlab="Dispersion variance", breaks=seq(0,1.9,by=.05), ylim=c(0,5))
abline(v=mean(dv.8), col="red")
par(mfrow=c(1,1))


###################################################
### code chunk number 26: exA1.Rnw:491-495
###################################################
grid.64  <- SpatialGrid(GridTopology(cellcentre.offset=c(2, 2), cellsize=c(4, 4), cells.dim=c(32, 32)))
grid.32 <- SpatialGrid(GridTopology(cellcentre.offset=c(4, 4), cellsize=c(8, 8), cells.dim=c(16, 16)))
grid.16  <- SpatialGrid(GridTopology(cellcentre.offset=c(8, 8), cellsize=c(16, 16), cells.dim=c(8,8)))
summary(grid.128); summary(grid.64); summary(grid.32); summary(grid.16)


###################################################
### code chunk number 27: exA1.Rnw:510-516
###################################################
k.e8.64 <- idw(z~1, loc=k.e8.128, newdata=grid.64, idp=0, nmax=4)
names(k.e8.64) <- c("z","nr")
k.e8.32 <- idw(z~1, loc=k.e8.128, newdata=grid.32, idp=0, nmax=16)
names(k.e8.32) <- c("z","nr")
k.e8.16 <- idw(z~1, loc=k.e8.128, newdata=grid.16, idp=0, nmax=32)
names(k.e8.16) <- c("z","nr")


###################################################
### code chunk number 28: exA1.Rnw:526-529
###################################################
tmp <- idw(z~1, loc=k.e8.128, newdata=grid.64, idp=0, nmax=4)
summary(tmp$var1.pred); summary(k.e8.64$z)
rm(tmp)


###################################################
### code chunk number 29: exA1.Rnw:539-544
###################################################
z.at <- seq(-4.4,4.4,by=.2)
sim.plot.128 <- spplot(k.e8.128, zcol="z", col.regions=bpy.colors(64), main='Simulated field, vgm(1,"Exp",8,0), 128x128', at=z.at)
sim.plot.64 <- spplot(k.e8.64, zcol="z", col.regions=bpy.colors(64), main='Simulated field, vgm(1,"Exp",8,0), 64x64', at=z.at)
sim.plot.32 <- spplot(k.e8.32, zcol="z", col.regions=bpy.colors(64), main='Simulated field, vgm(1,"Exp",8,0), 32x32', at=z.at)
sim.plot.16 <- spplot(k.e8.16, zcol="z", col.regions=bpy.colors(64), main='Simulated field, vgm(1,"Exp",8,0), 16x16', at=z.at)


###################################################
### code chunk number 30: exA1.Rnw:549-553
###################################################
print(sim.plot.128, split=c(1,1,2,2), more=T)
print(sim.plot.64, split=c(2,1,2,2), more=T)
print(sim.plot.32, split=c(1,2,2,2), more=T)
print(sim.plot.16, split=c(2,2,2,2), more=F)


###################################################
### code chunk number 31: exA1.Rnw:572-581
###################################################
cbind(rbind(summary(k.e8.128$z),
            summary(k.e8.64$z),
            summary(k.e8.32$z),
            summary(k.e8.16$z)),
      "Var."=c(
          var(k.e8.128$z),
          var(k.e8.64$z),
          var(k.e8.32$z),
          var(k.e8.16$z)))


###################################################
### code chunk number 32: exA1.Rnw:596-602
###################################################
par(mfrow=c(2,2))
hist(k.e8.128$z, xlim=c(-4.5, 4.5), main="128 x 128")
hist(k.e8.64$z, xlim=c(-4.5, 4.5), main="64 x 64")
hist(k.e8.32$z, xlim=c(-4.5, 4.5), main="32 x 32")
hist(k.e8.16$z, xlim=c(-4.5, 4.5), main="16 x 16")
par(mfrow=c(1,1))


###################################################
### code chunk number 33: exA1.Rnw:609-615
###################################################
boxplot(data.frame(res.128=k.e8.128$z,
                   res.64=k.e8.64$z,
                   res.32=k.e8.32$z,
                   res.16=k.e8.16$z),
        main="Four aggregation levels",
        horizontal=T, boxwex=.5, cex=.6)


###################################################
### code chunk number 34: exA1.Rnw:661-670
###################################################
v.128 <- variogram(z ~ 1, loc=k.e8.128, cutoff=32)
v.64 <- variogram(z ~ 1, loc=k.e8.64, cutoff=32)
v.32 <- variogram(z ~ 1, loc=k.e8.32, cutoff=32)
v.16 <- variogram(z ~ 1, loc=k.e8.16, cutoff=32)
myTheme <- simpleTheme(col = "slateblue4", pch=20, cex=0.9)
v.128.pl <- plot(v.128, ylim=c(0, 1.0), pl=F, model=vm, main="128x128")
v.64.pl <- plot(v.64, ylim=c(0, 1.0), pl=F, model=vm, main="64x64")
v.32.pl <- plot(v.32, ylim=c(0, 1.0), pl=F, model=vm, main="32x32")
v.16.pl <- plot(v.16, ylim=c(0, 1.0), pl=F, model=vm, main="16x16")


###################################################
### code chunk number 35: printv (eval = FALSE)
###################################################
## trellis.par.set(myTheme)
## print(v.128.pl, split=c(1,1,2,2), more=T)
## print(v.64.pl, split=c(2,1,2,2), more=T)
## print(v.32.pl, split=c(1,2,2,2), more=T)
## print(v.16.pl, split=c(2,2,2,2), more=F)


###################################################
### code chunk number 36: exA1.Rnw:684-685
###################################################
trellis.par.set(myTheme)
print(v.128.pl, split=c(1,1,2,2), more=T)
print(v.64.pl, split=c(2,1,2,2), more=T)
print(v.32.pl, split=c(1,2,2,2), more=T)
print(v.16.pl, split=c(2,2,2,2), more=F)


###################################################
### code chunk number 37: exA1.Rnw:702-710
###################################################
(vmf.128 <- fit.variogram(v.128, vm))
(vmf.64 <- fit.variogram(v.64, vm))
(vmf.32 <- fit.variogram(v.32, vm))
(vmf.16 <- fit.variogram(v.16, vm))
v.128.pl <- plot(v.128, ylim=c(0, 1.0), pl=F, model=vmf.128, main="128x128")
v.64.pl <- plot(v.64, ylim=c(0, 1.0), pl=F, model=vmf.64, main="64x64")
v.32.pl <- plot(v.32, ylim=c(0, 1.0), pl=F, model=vmf.32, main="32x32")
v.16.pl <- plot(v.16, ylim=c(0, 1.0), pl=F, model=vmf.16, main="16x16")


###################################################
### code chunk number 38: exA1.Rnw:716-717
###################################################
trellis.par.set(myTheme)
print(v.128.pl, split=c(1,1,2,2), more=T)
print(v.64.pl, split=c(2,1,2,2), more=T)
print(v.32.pl, split=c(1,2,2,2), more=T)
print(v.16.pl, split=c(2,2,2,2), more=F)


###################################################
### code chunk number 39: exA1.Rnw:755-756
###################################################
save.image(file="exA1.RData")