################################################################
## ##
## a geostatistical analysis that demonstrates ##
## features of the package `RandomFields' ##
## ##
################################################################
data(soil)
names(soil)
pts <- soil[,c(1,2)]
d <- soil$moisture
## define some graphical parameters first
close.screen(close.screen())
maxbin <- max(dist(pts)) / 2
(zlim <- range(d))
cn <- 100
colour <- rainbow(cn)
par(cex=1, cex.lab=1.3, cex.axis=1.3, mar=c(4.3,4.3,0.8,0.8))
lu.x <- min(soil$x)
lu.y <- max(soil$y)
y <- x <- seq(min(soil$x), max(soil$x), l=121)
## ... and a certain appearance of the legend
my.legend <- function(lu.x, lu.y, zlim, col, cex=1) {
## uses already the legend code of R-1.3.0
cn <- length(col)
filler <- vector("character", length=(cn-3)/2)
legend(lu.x, lu.y, y.i=0.03, x.i=0.1,
legend=c(format(zlim[2], dig=2), filler,
format(mean(zlim), dig=2), filler,
format(zlim[1], dig=2)),
lty=1, col=col[length(col):1],cex=cex)
}
## plot the data first
plot(pts, col=colour[1+99*((d-min(d))/diff(zlim))], pch=16,
xlab="x [cm]", ylab="y [cm]")
my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1.3)
ev <- EmpiricalVariogram(pts, data=d, grid=FALSE,
bin=c(-1,seq(0,maxbin,l=30)))
## show all models,
by.eye <- ShowModels(x=x, y=y, emp=ev, col=colour, zlim=zlim)
## fit paramters of the whittlematern model by MLE
p <- mleRF(pts, d, "whittle", par=rep(NA,5), lower.k=0.01,
upper.k=30)
## plot the fitted model and the empirical variogram
plot(ev$c, ev$emp.var, ylim=c(0,11), ylab="variogram", xlab="lag")
gx <- seq(0.001, max(ev$c), l=100)
if(!is.null(by.eye))
lines(gx, Variogram(gx, model="whittle", par=by.eye$p))
lines(gx, Variogram(gx, model="whittle", par=p), col=2)
legend(120, 4, c("empirical", "by eye", "MLE"),
lty=c(-1, 1, 1), pch=c(1, -1, -1), col=c(1, 1, 2), cex=1.3)
## map of expected values
k <- Kriging("O", x=x, y=y, grid=TRUE, model="whittle", par=p,
given=pts, data=d)
image(x, y, k, col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]")
my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1.3)
## what is the probability that at no point of the
## grid given by x and y the moisture is greater than 24 percent?
cs <- CondSimu("O", x=x, y=y, grid=TRUE, model="whittle",
param=p, given=pts, data=d, n=100)
image(x, y, cs[,,1], col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]")
my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1.3)sum(apply(cs<=24, 3, all)) ## about 40 percent ...
Run the code above in your browser using DataLab