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=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
#--------------------------------------------------------------------------
# verfügbare Variogrammmodelle
#--------------------------------------------------------------------------
show.vgms()
pause()
#Modelle aus der Matern-Klasse:
show.vgms(model="Mat", kappa.range = c(0.1,0.2,0.5,1,2,5,10),max=10)
pause()
#--------------------------------------------------------------------------
# Überlagerung von Variogrammmodellen
#--------------------------------------------------------------------------
vgm1 <- vgm(1, "Sph", 300)
print(vgm1)
#
vgm2 <- vgm(1, "Sph", 300, 0.5)
print(vgm2)
#
vgm3 <- vgm(0.8, "Sph", 800, 0.5, add.to = vgm2)
print(vgm3)
#
vgm4 <- vgm(0.5, "Nug", 0)
print(vgm4)
# Liste der verfügbaren Modelle:
vgm()
#--------------------------------------------------------------------------
#  Emp. Variogramm
#--------------------------------------------------------------------------
cut.val = 2000
wdth.val = cut.val/10
vg1 <- variogram(ZV~1, 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()
#--------------------------------------------------------------------------
# Anpassung eines Gauss'schen Variogrammmodells mit nugget
#--------------------------------------------------------------------------
mod1 <- vgm(0, "Gau", 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\nGauss'sches Modell")
print(v.fit)
attr(v.fit,"SSErr")