rm(list=ls())
#--------------------------------------------------------------------------
#   Laden der packages
#--------------------------------------------------------------------------
library(lattice)
library(sp)
library(gstat)
# eine eigene Pausenfunktion:
pause <- function ()
{
    cat("Pause. Press <Enter> to continue...")
    readline()
    invisible()
}
#--------------------------------------------------------------------------
#   Voreinstellungen
#--------------------------------------------------------------------------
# Festsetzen von Graphikparametern für das Trellis-package
trellis.par.set(theme = col.whitebg())
axis.line <- trellis.par.get("axis.line");
axis.line$col <- "black";axis.line$lwd <- 4;
trellis.par.set("axis.line", axis.line)
#
fontsize <- trellis.par.get("fontsize");fontsize$text =10;
trellis.par.set("fontsize", fontsize)
par.ylab.text <- trellis.par.get("par.ylab.text");par.ylab.text$cex <- 1.8;
trellis.par.set("par.ylab.text", par.ylab.text)
par.xlab.text <- trellis.par.get("par.xlab.text");par.xlab.text$cex <- 1.8;
trellis.par.set("par.xlab.text", par.xlab.text)
axis.text <- trellis.par.get("axis.text");axis.text$cex =1.4;
trellis.par.set("axis.text", axis.text)
par.main.text <- trellis.par.get("par.main.text");par.main.text$cex =2.0;
trellis.par.set("par.main.text", par.main.text)
#
plot.symbol <- trellis.par.get("plot.symbol");
plot.symbol$col ="blue"
plot.symbol$pch =16
plot.symbol$fill = "blue"
plot.symbol$cex = 0.8
trellis.par.set("plot.symbol", plot.symbol)
#--------------------------------------------------------------------------
#   Einlesen der Daten
#--------------------------------------------------------------------------
filename <- file.choose()
dat <- read.table(filename,header=FALSE)
names(dat) <- c("X","Y","P","EC","pH","NitCo","NitFlx")
dat <- dat[,c(1,2,5)]
cc <- complete.cases(dat)
dat1 <- dat[cc,]
ZV <- dat1$pH
#--------------------------------------------------------------------------
#  Emp. Variogramm
#--------------------------------------------------------------------------
cut.val = 2000
wdth.val = cut.val/10
vg1 <- variogram(ZV~1, loc = ~X+Y, dat1, cutoff = cut.val, width = wdth.val)
a <- plot(vg1, plot.numbers = TRUE, pch = 15, cex = 1.0, col="blue"
, xlim = range(0, 1.1*max(vg1$dist)), main = "Empirisches Variogramm")
print(a)
pause()
#--------------------------------------------------------------------------
# Anpassung eines sphärischen Variogrammmodells mit nugget
#--------------------------------------------------------------------------
mod1 <- vgm(0, "Sph", 500,1)
v.fit <- fit.variogram(vg1, model=mod1, fit.method=7, fit.sills=TRUE, fit.ranges=TRUE)
plot(vg1, model=v.fit, lwd=4,pch=15,cex=1.8, col="black", xlim = range(0, 1.1*max(vg1$dist)),
     main = "Bruchsal T1 pH\nSphärisches Modell")
print(v.fit)
attr(v.fit,"SSErr")
pause()
#--------------------------------------------------------------------------
#   Erstellen eines Grids
#--------------------------------------------------------------------------
ndataX <- 100
ndataY <- 200
x <- seq(min(dat1$X),max(dat1$X),(max(dat1$X)-min(dat1$X))/(ndataX-1))
y <- seq(min(dat1$Y),max(dat1$Y),(max(dat1$Y)-min(dat1$Y))/(ndataY-1))
rpv <- rep(ndataY,ndataX)
u<-rep(1:ndataX,rpv);v<-rep(1:ndataY,ndataX)
w <- data.frame(cbind(x[u],y[v]))
dgr <- data.frame(X = w$X1, Y=w$X2)
plot(w$X1,w$X2,cex=0.2,pch=16,asp=1)
#--------------------------------------------------------------------------
#   Ordinary Kriging
#--------------------------------------------------------------------------
ZV.okr <- krige(ZV ~ 1, ~ X + Y, dat1, dgr, model = v.fit)
#
# Contourplot
#cl <- colors()[seq(280,360,1)]
#cl <- topo.colors(100)
title1 <- paste(names(dat)[3],"Ordinary kriging\n predictions")
title2 <- paste(names(dat)[3],"Ordinary kriging\n estimation variance")
cl <- rainbow(100,start=.7,end=.25)
cts <- 40
pl1 <- levelplot(var1.pred~X+Y, ZV.okr, aspect = "iso",
        main = title1,contour=FALSE,region=TRUE,pretty=TRUE,labels=FALSE,col.regions=cl,cuts=cts)
print(pl1)
#--------------------------------------------------------------------------
#    Schätzvarianz
#--------------------------------------------------------------------------
pl2 <- contourplot(var1.var~X+Y, ZV.okr, aspect = "iso",
        main = title2,contour=FALSE,region=TRUE,
        pretty=TRUE,labels=FALSE,col.regions=cl,cuts=cts)
print(pl2)
#--------------------------------------------------------------------------
#   alternative Plots
#--------------------------------------------------------------------------
# Schätzwerte: Contour
v <- ZV.okr$var1.pred
datmat <- matrix(v,byrow=T,nrow=ndataX)
filled.contour(x,y,datmat,asp=1,color = topo.colors,nlevels=50,
, main=title1,plot.axes={axis(1);axis(2)
; contour(x,y,datmat,nlevels=20,add=T,lwd=1,col="brown", axes=F)
; points(dat1$X,dat1$Y,pch=16,col="white",cex=1,lwd=1)
})
# Schätzwerte: 3D perspektivisch
title <- paste(names(dat)[3],"Ordinary kriging\n Estimation variance")
persp(x,y,datmat, zlim = c(min(v),1.2*max(v)),theta = 30, phi = 20, expand = 0.7, col = "lightblue",
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "X-Value", ylab = "\nY-Value", zlab = "ZV-Value", main=title1,
     r=4,border="lightblue4") -> rotmat
#   Hinzufügen gemessener Datenpunkte
points(trans3d(dat1$X, dat1$Y, ZV, pmat = rotmat), col = 2, pch =16)

# Schätzvarianz: Contour
v <- ZV.okr$var1.var
datmat <- matrix(v,byrow=T,nrow=ndataX)
filled.contour(x,y,datmat,asp=1,color = topo.colors,nlevels=50,
, main=title2,plot.axes={axis(1);axis(2)
; contour(x,y,datmat,nlevels=20,add=T,lwd=1,col="brown", axes=F)
; points(dat1$X,dat1$Y,pch=16,col="white",cex=1,lwd=1)
})
# Schätzvarianz: 3D perspektivisch
persp(x,y,datmat, zlim = c(0,1.2*max(v)),theta = 30, phi = 20, expand = 0.7, col = "lightblue",
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "X-Value", ylab = "\nY-Value", zlab = "ZV-Value", main=title2,
     r=4,border="lightblue4") -> rotmat
#   Hinzufügen gemessener Datenpunkte
#ZV.var <- ZV
#ZV.var[] <- 0.
#points(trans3d(dat1$X, dat1$Y, ZV.var, pmat = rotmat), col = 2, pch =16)
