soil {RandomFields} | R Documentation |
Soil physical and chemical data collected on a field in the Weissenstaedter Becken, Germany
data(soil)
This data frame contains the following columns:
For technical reasons some of the data were obtained as differences of two measurements (which are not available anymore). Therefore, some of the data have negative values.
The data were collected by Wolfgang Falk, Soil Physics Group, http://www.geo.uni-bayreuth.de/bodenphysik/Welcome.html, University of Bayreuth, Germany.
Falk, W. (2000) Kleinskalige raeumliche Variabilitaet von Lachgas und bodenchemischen Parameters [Small Scale Spatial Variability of Nitrous Oxide and Pedo-Chemical Parameters], Master thesis, University of Bayreuth, Germany.
################################################################ ## ## ## a geostatistical analysis that demonstrates ## ## features of the package `RandomFields' ## ## ## ## *** THIS EXAMPLE IS NOT RUN AUTOMATICALLY *** ## ## *** SINCE IT IS TIME CONSUMING *** ## ## ## ################################################################ ## Not run: 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(bg="white",cex=1, cex.lab=1.4, cex.axis=1.4, 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.4) ## empirical variogram 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, Mean=mean(d)) ## fit parameters of the whittlematern model by MLE mle <- fitvario(x=pts, data=d, model="whittle", par=rep(NA,5)) ## 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=by.eye)) lines(gx, Variogram(gx, model=mle$mle), col=2) lines(gx, Variogram(gx, model=mle$lsq), col=3) lines(gx, Variogram(gx, model=mle$nlsq), col=4) lines(gx, Variogram(gx, model=mle$slsq), col=5) legend(120, 4, c("empirical", "by eye", "MLE", "lsq", "sqrt(n) lsq", "sd^-0.5 lsq"), lty=c(-1, rep(1, 5)), pch=c(1, rep(-1, 5)), col=c(1, 1, 2, 3, 4, 5), cex=1.4) ## map of expected values k <- Kriging("O", x=x, y=y, grid=TRUE, model=mle$mle, given=pts, data=d) image(x, y, k, col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]") par(bg="white") my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1.4) ## 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=mle$mle, given=pts, data=d, n=10) # better n=100 or n=1000 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.4) mean(apply(cs<=24, 3, all)) ## about 40 percent ... ## End(Not run)