###################################################
### chunk number 1: 
###################################################
rm(list=ls())
library(digest)
.owidth <- getOption("width")
options("width"=70)
owarn <- options("warn")$warn
options(warn=1)
.epsNo <- 0


###################################################
### chunk number 2: figreset eval=FALSE
###################################################
## .iwidth <- 5
## .iheight <- 6
## .ipointsize <- 12


###################################################
### chunk number 3: 
###################################################
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 4: afig eval=FALSE
###################################################
## .epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
## postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)


###################################################
### chunk number 5: zfig eval=FALSE
###################################################
## dev.null <- dev.off()
## system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
## cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")


###################################################
### chunk number 6: afig_png eval=FALSE
###################################################
## .epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".png", sep="")
## png(filename=file, width = .iwidth, height = .iheight, pointsize = .ipointsize)


###################################################
### chunk number 7: zfig_png eval=FALSE
###################################################
## dev.null <- dev.off()
## cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")


###################################################
### chunk number 8:  eval=FALSE
###################################################
## library(lattice)
## library(sp)
## data(meuse)
## coordinates(meuse) <- c("x", "y")
## spplot(meuse, "zinc", do.log = T)
## bubble(meuse, "zinc", do.log = T, key.space = "bottom")


###################################################
### chunk number 9: 
###################################################
library(lattice)
trellis.par.set(canonical.theme(color = FALSE))
library(sp)
data(meuse)
coordinates(meuse) <- c("x", "y")
#print(spplot(meuse, "zinc", do.log = T), split = c(1,1,2,1), more = TRUE)
#print(bubble(meuse, "zinc", do.log = T, key.space = "bottom"), split = c(2,1,2,1))


###################################################
### chunk number 10:  eval=FALSE
###################################################
## xyplot(log(zinc)~sqrt(dist), as.data.frame(meuse))
## zn.lm <- lm(log(zinc)~sqrt(dist), meuse)
## meuse$fitted.s <- predict(zn.lm, meuse) - mean(predict(zn.lm, meuse))
## meuse$residuals <- residuals(zn.lm)
## spplot(meuse, c("fitted.s", "residuals"))


###################################################
### chunk number 11: 
###################################################
.iwidth <- 6
.iheight <- 4
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
print(xyplot(log(zinc)~sqrt(dist), as.data.frame(meuse), asp = 1), split = c(1,1,2,1), more = TRUE)
zn.lm <- lm(log(zinc)~sqrt(dist), meuse)
meuse$fitted.s <- predict(zn.lm, meuse) - mean(predict(zn.lm, meuse))
meuse$residuals <- residuals(zn.lm)
# RSB quietening greys
print(spplot(meuse, c("fitted.s", "residuals"), cuts=4, col.regions=grey.colors(5, 0.95, 0.55, 2.2)), split = c(2,1,2,1))
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 12: 
###################################################
data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")


###################################################
### chunk number 13: 
###################################################
library(gstat)
idw.out <- idw(zinc~1, meuse, meuse.grid, idp = 2.5)
as.data.frame(idw.out)[1:5,]


###################################################
### chunk number 14: 
###################################################
zn.lm <- lm(log(zinc)~sqrt(dist), meuse)
meuse.grid$pred <- predict(zn.lm, meuse.grid)
meuse.grid$se.fit <- predict(zn.lm, meuse.grid, se.fit=TRUE)$se.fit
#$


###################################################
### chunk number 15: 
###################################################
meuse.lm <- krige(log(zinc)~sqrt(dist), meuse, meuse.grid)


###################################################
### chunk number 16: 
###################################################
meuse.tr2 <- krige(log(zinc)~1, meuse, meuse.grid, degree = 2)


###################################################
### chunk number 17:  eval=FALSE
###################################################
## lm(log(zinc)~I(x^2)+I(y^2)+I(x*y) + x + y, meuse)


###################################################
### chunk number 18:  eval=FALSE
###################################################
## lm(log(zinc)~poly(x,y,2), meuse)


###################################################
### chunk number 19:  eval=FALSE
###################################################
## hscat(log(zinc)~1,meuse,(0:9)*100)


###################################################
### chunk number 20: 
###################################################
.iwidth <- 6
.iheight <- 6
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
hscat(log(zinc)~1,meuse,(0:9)*100,pch=3,cex=.6)
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 21:  eval=FALSE
###################################################
## library(gstat)
## variogram(log(zinc) ~ 1, meuse, cloud = TRUE)


###################################################
### chunk number 22: 
###################################################
.iwidth <- 6
.iheight <- 6
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
library(gstat)
cld <- variogram(log(zinc) ~ 1, meuse, cloud = TRUE)
svgm <- variogram(log(zinc) ~ 1, meuse)
d <- data.frame(gamma = c(cld$gamma, svgm$gamma),
	dist = c(cld$dist, svgm$dist),
	id = c(rep("cloud", nrow(cld)), rep("sample variogram", nrow(svgm)))
	)
xyplot(gamma ~ dist | id, d,
	scales = list(y = list(relation = "free", ylim = list(NULL, c(-.005,0.7)))),
	layout = c(1, 2), as.table = TRUE,
	panel = function(x,y, ...) {
		if (panel.number() == 2)
			ltext(x+10, y, svgm$np, adj = c(0,0.5)) #$
		panel.xyplot(x,y,...)
	}, 
	xlim = c(0, 1590),
	cex = .5, pch = 3
)
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 23:  eval=FALSE
###################################################
## plot(variogram(log(zinc) ~ 1, meuse))


###################################################
### chunk number 24:  eval=FALSE
###################################################
## sel <- plot(variogram(zinc ~ 1, meuse, cloud = TRUE), digitize = TRUE)
## plot(sel, meuse)


###################################################
### chunk number 25: 
###################################################
.iwidth <- 6
.iheight <- 4
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
sel <-
structure(list(x = c(145.291968730077, 266.107479142605, 320.156523274526, 
339.232656497557, 323.335878811698, 212.058435010685, 135.753902118561, 
46.7319470777507, 78.5255024494688, 142.112613192905), y = c(574649.690841889, 
581256.265954825, 627502.29174538, 822396.257577002, 1053626.38652977, 
1278249.94036961, 1255126.92747433, 792666.669568789, 634108.866858316, 
577952.978398357)), .Names = c("x", "y"))
v <- variogram(zinc ~ 1, meuse, cloud = TRUE)
v$gamma <- v$gamma/1e6
sel$y <- sel$y/1e6
p1 <- xyplot(gamma~dist, v, 
	panel = function(x, y, ...) {
		panel.xyplot(x, y, ...)
		llines(sel$x, sel$y)
	},
	pch=3, cex = .5, asp = 1, ylab = "gamma (x 1e6)")
x <-
structure(list(head = c(40, 40, 40, 54, 55, 54, 47, 80, 55, 55, 
54, 53, 54, 55, 59, 59), tail = c(41, 42, 43, 57, 57, 58, 59, 
99, 121, 122, 123, 125, 125, 125, 125, 132)), .Names = c("head", 
"tail"), row.names = as.integer(c(NA, 16)), class = c("pointPairs", 
"data.frame"))
p2 = plot(x, meuse, scales=list(draw=F), col.line = 1)
print(p1, split = c(1,1,2,1), more = TRUE)
print(p2, split = c(2,1,2,1))
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 26: 
###################################################
.iwidth <- 6
.iheight <- 4
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
v <- variogram(log(zinc) ~ 1, meuse)
# INTERACTIVE mode out-commented:
#plot(v, type = 'b', pch = 3)
#fn = function(n = 100) {
#        for (i in 1:n) {
#        	meuse$random = sample(meuse$zinc)
#        	v = variogram(log(random) ~ 1, meuse)
#        	trellis.focus("panel", 1, 1, highlight = FALSE)
#        	llines(v$dist, v$gamma, col = 'grey')
#        	trellis.unfocus()
#        }
#}
#fn()
#trellis.focus("panel", 1, 1, highlight = FALSE)
#lpoints(v$dist, v$gamma, col = 'black', type = 'b', lwd = 2, pch=3)
#trellis.unfocus()
print(xyplot(gamma ~ dist, v, pch = 3, type = 'b', lwd = 2,
	panel = function(x, y, ...) {
        for (i in 1:100) {
        	meuse$random = sample(meuse$zinc)
        	v = variogram(log(random) ~ 1, meuse)
        	llines(v$dist, v$gamma, col = 'grey')
		}
		panel.xyplot(x, y, ...)
	},
	ylim = c(0, 0.75), xlab = 'distance', ylab = 'semivariance'
))
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 27:  eval=FALSE
###################################################
## plot(variogram(log(zinc) ~ 1, meuse))


###################################################
### chunk number 28:  eval=FALSE
###################################################
## plot(variogram(log(zinc) ~ 1, meuse, alpha = c(0, 45, 90, 135)))


###################################################
### chunk number 29:  eval=FALSE
###################################################
## plot(variogram(log(zinc) ~ 1, meuse, cutoff = 1000, width = 50))


###################################################
### chunk number 30:  eval=FALSE
###################################################
## variogram(log(zinc) ~ 1, meuse, boundaries = c(0,50,100,seq(250,1500,250)))


###################################################
### chunk number 31:  eval=FALSE
###################################################
## show.vgms()
## show.vgms(model = "Mat", kappa.range = c(.1, .2, .5, 1, 2, 5, 10), max = 10)


###################################################
### chunk number 32: 
###################################################
vgm(1, "Sph", 300)
vgm(1, "Sph", 300, 0.5)
v1 <- vgm(1, "Sph", 300, 0.5)
v2 <- vgm(0.8, "Sph", 800, add.to = v1)
v2
vgm(0.5, "Nug", 0)


###################################################
### chunk number 33: 
###################################################
vgm()


###################################################
### chunk number 34:  eval=FALSE
###################################################
## v <- variogram(log(zinc) ~ 1, meuse)
## plot(v)


###################################################
### chunk number 35: 
###################################################
fit.variogram(v, vgm(1, "Sph", 800, 1))


###################################################
### chunk number 36: 
###################################################
fit.variogram(v, vgm(1, "Sph", 10, 1))


###################################################
### chunk number 37:  eval=FALSE
###################################################
## v.fit <- fit.variogram(v, vgm(1, "Sph", 10, 1))
## if (attr(v.fit, "singular"))
## 	stop("singular fit")


###################################################
### chunk number 38: 
###################################################
.iwidth <- 6
.iheight <- 3.5
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
v <- variogram(log(zinc) ~ 1, meuse)
v.fit <- fit.variogram(v, vgm(1, "Sph", 800, 1))
plot(v, v.fit, pch = 3)
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 39: 
###################################################
attr(v.fit, "SSErr")


###################################################
### chunk number 40:  eval=FALSE
###################################################
## library(geoR)
## v.eye <- eyefit(variog(as.geodata(meuse["zinc"]),max.dist=1500))
## ve.fit <- as.vgm.variomodel(v.eye[[1]])


###################################################
### chunk number 41: 
###################################################
fit.variogram(v, vgm(1, "Sph", 800, 0.06), fit.sills = c(FALSE, TRUE))


###################################################
### chunk number 42: 
###################################################
fit.variogram.reml(log(zinc)~1, meuse, model=vgm(0.6, "Sph", 800, 0.06))


###################################################
### chunk number 43: 
###################################################
v.dir <- variogram(log(zinc)~1,meuse,alpha=(0:3)*45) 
v.anis <- vgm(.6, "Sph", 1600, .05, anis=c(45, 0.3))


###################################################
### chunk number 44:  eval=FALSE
###################################################
## plot(v.dir, v.anis)


###################################################
### chunk number 45: 
###################################################
.iwidth <- 6
.iheight <- 3
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
v.dir <- variogram(log(zinc)~1,meuse,alpha=(0:3)*45) 
v.anis <- vgm(.6, "Sph", 1600, .05, anis=c(45,.3))
print(plot(v.dir, v.anis, pch=3))
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 46:  eval=FALSE
###################################################
## plot(variogram(log(zinc)~1,meuse, map=TRUE, cutoff=1000, width=100))


###################################################
### chunk number 47: 
###################################################
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)
g 
vm <- variogram(g)
vm.fit <- fit.lmc(vm, g, vgm(1, "Sph", 800, 1))


###################################################
### chunk number 48:  eval=FALSE
###################################################
## plot(vm, vm.fit)


###################################################
### chunk number 49: 
###################################################
cor(as.data.frame(meuse)[c("cadmium", "copper", "lead", "zinc")])


###################################################
### chunk number 50: 
###################################################
.iwidth <- 7
.iheight <- 7
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
print(plot(vm, vm.fit))
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 51:  eval=FALSE
###################################################
## variogram(log(zinc) ~ sqrt(dist), meuse)


###################################################
### chunk number 52: 
###################################################
f <- log(zinc) ~ sqrt(dist)
vt <- variogram(f, meuse)
vt.fit <- fit.variogram(vt, vgm(1, "Exp", 300, 1))
vt.fit
g.wls <- gstat(NULL, "log-zinc", f, meuse, model=vt.fit, set = list(gls=1))
(variogram(g.wls)$gamma - vt$gamma)/ mean(vt$gamma)
#$


###################################################
### chunk number 53: 
###################################################
lz.sk <- krige(log(zinc)~1, meuse, meuse.grid, v.fit, beta = 5.9)
lz.ok <- krige(log(zinc)~1, meuse, meuse.grid, v.fit)
lz.uk <- krige(log(zinc)~sqrt(dist), meuse, meuse.grid, vt.fit)


###################################################
### chunk number 54: 
###################################################
cok.maps <- predict(vm.fit, meuse.grid)
names(cok.maps)


###################################################
### chunk number 55:  eval=FALSE
###################################################
## spplot.vcov(cok.maps)


###################################################
### chunk number 56: 
###################################################
.iwidth <- 8
.iheight <- 8
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
print(spplot.vcov(cok.maps, cuts=8,
 col.regions=grey.colors(9, 0.9, 0.45, 2.2)))
# RSB quietening greys
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 57:  eval=FALSE
###################################################
## vm2.fit <- vm.fit
## vm2.fit$model[[3]]$range  = c(0, 900)
## predict(vm2.fit, meuse.grid)


###################################################
### chunk number 58: 
###################################################
vm2.fit <- vm.fit
vm2.fit$model[[3]]$range  = c(0, 900)


###################################################
### chunk number 59: 
###################################################
vm2.fit$set <- list(nocheck=1) #$
x <- predict(vm2.fit, meuse.grid)
names(as.data.frame(x))
any(as.data.frame(x)[c(2,4,6,8)] < 0)


###################################################
### chunk number 60: 
###################################################
g.cc <- gstat(NULL, "log.zinc", log(zinc)~1, meuse, model = v.fit)
meuse.grid$distn <- meuse.grid$dist - mean(meuse.grid$dist) + mean(log(meuse$zinc))
vd.fit <- v.fit
vov <- var(meuse.grid$distn) / var(log(meuse$zinc))
vd.fit$psill <- v.fit$psill * vov
g.cc <- gstat(g.cc, "distn", distn ~ 1, meuse.grid, nmax = 1, model=vd.fit, 
	merge = c("log.zinc","distn"))
vx.fit <- v.fit
vx.fit$psill <- sqrt(v.fit$psill * vd.fit$psill) * cor(meuse$dist, log(meuse$zinc)) #$
g.cc <- gstat(g.cc, c("log.zinc", "distn"), model = vx.fit) 
x <- predict(g.cc, meuse.grid)


###################################################
### chunk number 61: 
###################################################
.iwidth <- 6
.iheight <- 6
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
x$lz.uk <- lz.uk$var1.pred
x$lz.ok <- lz.ok$var1.pred
print(spplot(x, c("log.zinc.pred", "lz.ok", "lz.uk"),
	names.attr = c("collocated", "ordinary", "universal"),
	cuts=7, col.regions=grey.colors(8, 0.55, 0.95, 2.2)
))
# RSB quietening greys
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 62:  eval=FALSE
###################################################
## demo(pcb)


###################################################
### chunk number 63: 
###################################################
lz.ok <- krige(log(zinc)~1, meuse, meuse.grid, v.fit, block = c(40, 40))


###################################################
### chunk number 64: 
###################################################
xy <- expand.grid(x = seq(-20, 20, 4), y = seq(-20, 20, 4))
xy <- xy[(xy$x^2 + xy$y^2) <= 20^2, ]
lz.ok <- krige(log(zinc)~1, meuse, meuse.grid, v.fit, block = xy)


###################################################
### chunk number 65:  eval=FALSE
###################################################
## lz.pols <- krige(log(zinc)~1, meuse, meuse.polygons, v.fit)


###################################################
### chunk number 66:  eval=FALSE
###################################################
## spsample(polygon, n = 500, type = "regular", offset = c(.5, .5))


###################################################
### chunk number 67: 
###################################################
options("width"=55)


###################################################
### chunk number 68: 
###################################################
meuse$part.a <- idw(part.a~1, meuse.grid, meuse, nmax=1)$var1.pred


###################################################
### chunk number 69: 
###################################################
options("width"=70)


###################################################
### chunk number 70: 
###################################################
meuse$part.a <- meuse.grid$part.a[overlay(meuse.grid, meuse)]


###################################################
### chunk number 71: 
###################################################
x1 <- krige(log(zinc)~1, meuse[meuse$part.a == 0,],
meuse.grid[meuse.grid$part.a == 0,], model = vgm(.548, "Sph", 900, .0654),
nmin = 20, nmax = 40, maxdist = 1000)
x2 <- krige(log(zinc)~1, meuse[meuse$part.a == 1,],
meuse.grid[meuse.grid$part.a == 1,], model = vgm(.716, "Sph", 900),
nmin = 20, nmax = 40, maxdist = 1000)
lz.stk <- rbind(as.data.frame(x1), as.data.frame(x2))
coordinates(lz.stk) <- c("x", "y")
lz.stk <- as(x, "SpatialPixelsDataFrame")


###################################################
### chunk number 72:  eval=FALSE
###################################################
## spplot(lz.stk["var1.pred"], main = "stratified kriging predictions")


###################################################
### chunk number 73: 
###################################################
g.tr <- gstat(formula = log(zinc) ~ sqrt(dist), data = meuse, model = v.fit)
predict(g.tr, meuse[1,])
predict(g.tr, meuse[1,], BLUE = TRUE)


###################################################
### chunk number 74:  eval=FALSE
###################################################
## predict(g, meuse[1,], BLUE = TRUE, debug = 32)


###################################################
### chunk number 75: 
###################################################
meuse$Int <- rep(1, 155) #$
g.tr <- gstat(formula = log(zinc) ~ -1+Int+sqrt(dist), data = meuse, model = v.fit)
rn <- c("Intercept", "beta1")
df <- data.frame(Int = c(0,1), dist = c(1,0), row.names=rn)
spdf <- SpatialPointsDataFrame(SpatialPoints(matrix(0, 2, 2)), df)
spdf
predict(g.tr, spdf, BLUE = TRUE)


###################################################
### chunk number 76:  eval=FALSE
###################################################
## library(MASS)
## boxcox(zinc~sqrt(dist), data=as.data.frame(meuse))


###################################################
### chunk number 77: 
###################################################
meuse$zinc.ns <- qqnorm(meuse$zinc, plot.it = FALSE)$x
#$


###################################################
### chunk number 78: 
###################################################
ind.f <- I(zinc < 500) ~ 1
ind.fit <- fit.variogram(variogram(ind.f, meuse), vgm(1, "Sph", 800, 1))
ind.kr <- krige(ind.f, meuse, meuse.grid, ind.fit)
summary(ind.kr$var1.pred)
#$


###################################################
### chunk number 79:  eval=FALSE
###################################################
## meuse.dup <- rbind(as.data.frame(meuse)[1,], as.data.frame(meuse))
## coordinates(meuse.dup)=~x+y
## krige(log(zinc)~1, meuse.dup, meuse[1,], v.fit)


###################################################
### chunk number 80: 
###################################################
meuse.dup <- rbind(as.data.frame(meuse)[1,], as.data.frame(meuse))
coordinates(meuse.dup)=~x+y


###################################################
### chunk number 81: 
###################################################
zd <- zerodist(meuse.dup)
zd
krige(log(zinc)~1, meuse.dup[-zd[,1], ], meuse[1,], v.fit)


###################################################
### chunk number 82: 
###################################################
setL <- list(cn_max=1e10)
krige(log(zinc)~1, meuse.dup, meuse[1,], v.fit, set = setL)


###################################################
### chunk number 83: 
###################################################
set.seed(1357531)


###################################################
### chunk number 84: 
###################################################
sel100 <- sample(1:155, 100)
m.model <- meuse[sel100,]
m.valid <- meuse[-sel100,]
v100.fit <- fit.variogram(variogram(log(zinc)~1, m.model), vgm(1, "Sph", 800, 1))
m.valid.pr <- krige(log(zinc)~1, m.model, m.valid, v100.fit)
resid.kr <- log(m.valid$zinc) - m.valid.pr$var1.pred
summary(resid.kr)
resid.mean <- log(m.valid$zinc) - mean(log(m.valid$zinc))
R2 <- 1 - sum(resid.kr^2)/sum(resid.mean^2)
R2


###################################################
### chunk number 85: 
###################################################
m.valid.pr$res <- resid.kr #$


###################################################
### chunk number 86:  eval=FALSE
###################################################
## bubble(m.valid.pr, "res")


###################################################
### chunk number 87: 
###################################################
nfold <- 3
part <- sample(1:nfold, 155, replace = TRUE)
sel <- (part != 1)
m.model <- meuse[sel,]
m.valid <- meuse[-sel,]


###################################################
### chunk number 88: 
###################################################
v.fit <- vgm(.59, "Sph", 874, .04)
cv155 <- krige.cv(log(zinc)~1, meuse, v.fit, nfold=5)


###################################################
### chunk number 89:  eval=FALSE
###################################################
## bubble(cv155, "residual", main = "log(zinc): 5-fold CV residuals")


###################################################
### chunk number 90: 
###################################################
.iwidth <- 6
.iheight <- 6
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
print(bubble(cv155, "residual", main = "log(zinc): 5-fold CV residuals",
	maxsize = 1.5, col = c(grey(.15), grey(.7))))
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 91: 
###################################################
summary(cv155)


###################################################
### chunk number 92:  eval=FALSE
###################################################
## g.cv <- gstat.cv(g, nmax=40)


###################################################
### chunk number 93: 
###################################################
v1.fit <- vgm(0.591, "Sph", 897, .0507)
v2.fit <- vgm(0.591, "Sph", 897, add.to = vgm(0.0507, "Sph", 40))


###################################################
### chunk number 94: 
###################################################
set.seed(13331)
cv155.1 <- krige.cv(log(zinc)~1, meuse, v1.fit, nfold=5)
set.seed(13331)
cv155.2 <- krige.cv(log(zinc)~1, meuse, v2.fit, nfold=5)
summary(cv155.1$residual-cv155.2$residual)


###################################################
### chunk number 95: 
###################################################
b1 <- krige(log(zinc)~1, meuse, meuse.grid, v1.fit, block = c(40,40))$var1.var
b2 <- krige(log(zinc)~1, meuse, meuse.grid, v2.fit, block = c(40,40))$var1.var
summary((b1-b2)/b1)


###################################################
### chunk number 96: 
###################################################
lzn.sim <- krige(log(zinc)~1, meuse, meuse.grid, v.fit, nsim = 6, nmax=40)


###################################################
### chunk number 97:  eval=FALSE
###################################################
## spplot(lzn.sim)


###################################################
### chunk number 98: 
###################################################
.iwidth <- 6
.iheight <- 5
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
#print(spplot(lzn.sim, col.regions=grey(c(.1,.3,.4,.5,.7,.8,.9,.95)),cuts=6))
print(spplot(lzn.sim, col.regions=grey.colors(7, 0.55, 0.95, 2.2),cuts=6))
# RSB quietening greys
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 99: 
###################################################
`area` <-
structure(c(181130.059662363, 180917.131281033, 180769.007189672, 
180676.429632572, 180611.625342602, 180463.501251241, 180389.439205561, 
180556.078808342, 180694.945143992, 180796.780456802, 180824.553723932, 
180898.615769613, 181000.451082423, 181130.059662363, 181241.152730884, 
181268.925998014, 181194.863952334, 181130.059662363, 333718.312421053, 
333491.307789474, 333245.386105263, 332990.005894737, 332781.918315789, 
332554.913684211, 332498.162526316, 332413.035789474, 332592.747789474, 
332810.293894737, 333008.922947368, 333188.634947368, 333358.888421053, 
333491.307789474, 333566.976, 333614.268631579, 333793.980631579, 
333718.312421053), .Dim = as.integer(c(18, 2)))
area.sp = SpatialPolygons(list(Polygons(list(Polygon(area)), "area")))


###################################################
### chunk number 100: 
###################################################
nsim <- 1000
cutoff <- 500
grd <- overlay(meuse.grid, area.sp)
sel.grid <- meuse.grid[!is.na(grd),]
lzn.sim <- krige(log(zinc)~1, meuse, sel.grid, v.fit, nsim = nsim, nmax=40)
res <- apply(as.data.frame(lzn.sim)[1:nsim], 2, function(x) mean(x > log(cutoff)))


###################################################
### chunk number 101:  eval=FALSE
###################################################
## hist(res, main = paste("fraction above", cutoff), xlab = NULL, ylab = NULL)


###################################################
### chunk number 102: 
###################################################
bkr <- krige(log(zinc)~1, meuse, area.sp, v.fit)
1 - pnorm(log(cutoff), bkr$var1.pred, sqrt(bkr$var1.var))


###################################################
### chunk number 103: 
###################################################
.iwidth <- 6
.iheight <- 3
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
meuse.grid$part <- meuse.grid$part.a + meuse.grid$part.b
layout(matrix(1:2, 1, 2))
omar <- par("mar")
par(mar = rep(0,4))
image(meuse.grid["part"], col = 'gray') #$
lines(area)
par(mar = c(2,2,0.5,0)+.1)
hist(res, main = NULL, xlab = NULL, ylab = NULL)
par(mar = omar)
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 104:  eval=FALSE
###################################################
## cok.sims <- predict(vm.fit, meuse.grid, nsim = 1000)


###################################################
### chunk number 105: 
###################################################
table(meuse$soil) #$
s1.fit <- fit.variogram(variogram(I(soil==1)~1,meuse), vgm(1, "Sph", 800, 1))
s1.sim <- krige(I(soil==1)~1, meuse, meuse.grid, s1.fit, nsim = 6, indicators = TRUE, nmax = 40)


###################################################
### chunk number 106:  eval=FALSE
###################################################
## spplot(s1.sim)


###################################################
### chunk number 107: 
###################################################
.iwidth <- 6
.iheight <- 5
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
print(spplot(s1.sim, cuts=1, col.regions=grey.colors(2, 0.6, 0.90, 2.2)))
# RSB quietening greys
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 108:  eval=FALSE
###################################################
## m1 <- sapply(1:155, function(x) mean(krige(log(zinc)~1, meuse[-x,], 
## 	meuse.grid, v.fit)$var1.var))
## which(m1 == min(m1)) #$


###################################################
### chunk number 109:  eval=FALSE
###################################################
## plot(sort(m1))


###################################################
### chunk number 110:  eval=FALSE
###################################################
## cutoff <- 1000
## f <- function(x) {
## 	kr = krige(log(zinc)~1, meuse[-x,], meuse.grid, v.fit)
## 	mean(abs(pnorm((kr$var1.pred - log(cutoff))/sqrt(kr$var1.var)) - 0.5))
## }
## m2 <- sapply(1:155, f)
## which(m2 == max(m2))


###################################################
### chunk number 111: 
###################################################
.iwidth <- 6
.iheight <- 4
.epsNo <- .epsNo + 1; file <- paste("Fig-geos-", .epsNo, ".eps", sep="")
postscript(file=file, onefile = TRUE, paper="special", width = .iwidth, height = .iheight, pointsize = .ipointsize, horizontal=FALSE)
source("m1m2.R")
layout(matrix(1:2, 1, 2))
omar <- par("mar")
par(mar = rep(0,4))
plot(meuse)
points(meuse[m1 < quantile(m1,.1),], pch=1)
points(meuse[m1 > quantile(m1,.9),], pch=16)
plot(meuse)
points(meuse[m2 < quantile(m2,.1),], pch=16)
points(meuse[m2 > quantile(m2,.9),], pch=1)
par(mar = omar)
dev.null <- dev.off()
system(paste("./mungeEps.sh", file, "3 > outfile ; mv outfile", file, sep=" "))
cat("\\includegraphics[width=0.95\\textwidth]{", file, "}\n\n", sep="")
.iwidth <- 5
.iheight <- 6
.ipointsize <- 12


###################################################
### chunk number 112: 
###################################################
library(RandomFields)
ver <- packageDescription("RandomFields")$Version


###################################################
### chunk number 113: 
###################################################
PrintModelList()


###################################################
### chunk number 114:  eval=FALSE
###################################################
## library(geoR)
## plot(variog(as.geodata(meuse["zinc"]), max.dist=1500))


###################################################
### chunk number 115: 
###################################################
sI <- capture.output(print(sessionInfo()))
cat(paste("%", sI[1:(length(sI)/2)], sep=" "), sep="\n")
sT <- capture.output(print(Sys.time()))
cat("%", sT, "\n")