rm(list=ls())
library(lattice)
library(sp)
library(gstat)
# eine eigene Pausenfunktion:
pause <- function ()
{
    cat("Pause. Press <Enter> to continue...")
    readline()
    invisible()
}
#
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)
par.ylab.text <- trellis.par.get("par.ylab.text");par.ylab.text$cex <- 1.4;
trellis.par.set("par.ylab.text", par.ylab.text)
par.xlab.text <- trellis.par.get("par.xlab.text");par.xlab.text$cex <- 1.4;
trellis.par.set("par.xlab.text", par.xlab.text)
fontsize <- trellis.par.get("fontsize");fontsize$text =12;
trellis.par.set("fontsize", fontsize)
axis.text <- trellis.par.get("axis.text");axis.text$cex =1.2;
trellis.par.set("axis.text", axis.text)
plot.symbol <- trellis.par.get("plot.symbol");
plot.symbol$col ="blue"
plot.symbol$pch =16
plot.symbol$fill = "blue"
trellis.par.set("plot.symbol", plot.symbol)

#--------------------------------------------------------------------------
#   Einlesen der Daten
#--------------------------------------------------------------------------
filename <- file.choose()
dat <- read.table(filename,header=TRUE)
names(dat) <- c("X","Y","P","EC","pH","NitCo","NitFlx")
cc <- complete.cases(dat)
dat1 <- dat[cc,]
coordinates(dat1) <- c("X", "Y")
ZV <- dat1$pH
#--------------------------------------------------------------------------
#  variogram cloud
#--------------------------------------------------------------------------
par(mfrow=c(1,1))
g <- variogram(ZV~1,dat1,cloud=TRUE)
plot(g, plot.numbers = F, pch = 15, cex = 0.8
, xlim = range(0, 1.1*max(g$dist)), main = "Variogramm-Wolke")
pause()
#--------------------------------------------------------------------------
#  variogram cloud: select outliers
#--------------------------------------------------------------------------
g <- variogram(ZV~1,dat1,cloud=TRUE)
sel <- plot(g, plot.numbers = F, pch = 15, cex = 0.8
, xlim = range(0, 1.1*max(g$dist)), main = "Outlier selection"
, digitize = TRUE)
plot(sel,dat1)
cat("Selektierte Daten: sel$head: \n")
print(dat1[sel$head,])
cat("Selektierte Daten: sel$tail: \n")
print(dat1[sel$tail,])
#--------------------------------------------------------------------------
#  variogram estimation: standard method
#--------------------------------------------------------------------------
vg <- variogram(ZV~1,dat1)
# Ausgabe als plot
plot(vg, plot.numbers = TRUE, pch = 15, cex = 1.2
, xlim = range(0, 1.1*max(vg$dist))
, main = "Empirisches Variogramm\nStandardschätzung")
#--------------------------------------------------------------------------
#  variogram estimation: robust method (CRESSIE = TRUE)
#--------------------------------------------------------------------------
vg.cr <- variogram(ZV~1,dat1,cressie = TRUE)
# Ausgabe als plot
plot(vg.cr, plot.numbers = TRUE, pch = 15, cex = 1.2
, xlim = range(0, 1.1*max(vg$dist))
, main = "Empirisches Variogramm\nrobuste Schätzung")

#   Schluss