rm(list=ls())
library(lattice)
library(sp)
library(gstat)
# eine eigene Pausenfunktion:
pause <- function ()
{
    cat("Pause. Press <Enter> to continue...")
    readline()
    invisible()
}
#--------------------------------------------------------------------------
#   Einlesen der Daten
#--------------------------------------------------------------------------
setwd("C:\\Dokumente und Einstellungen\\Huwe\\Eigene Dateien\\Lehre\\3_Geoökologie_Bsc\\Lehrveranstaltungen\\VÜ2_Geostatistik\\Übungen_GS\\Datensätze\\Bruchsal")
filename <- file.choose()
dat <- read.table(filename,header=F)
names(dat) <- c("X","Y","P","EC","pH","NitCo","NitFlx")
cc <- complete.cases(dat)
dat1 <- dat[cc,]
ZV <- dat1$pH
#--------------------------------------------------------------------------
# Lineare und Quadratische Regressionen
#--------------------------------------------------------------------------
par(mfrow=c(2,2))
dat1.lm1 <- lm(ZV~X + Y,dat1)
print(summary(dat1.lm1))
#res.lin <- as.vector(dat1.lm1$residuals)
res.lin <- dat1.lm1$residuals
#
dat1.lm2 <- lm(ZV~ I(X^2) + I(Y^2) + I(X*Y) + X + Y,dat1)
print(summary(dat1.lm2))
res.quad <- dat1.lm2$residuals
#
# Erzeugen der Data-Frames
#
dfr.res.lin <- data.frame(dat1[,c(1,2)],res.lin)
names(dfr.res.lin) <- c("X","Y","ZVres")
coordinates(dfr.res.lin) <- c("X", "Y")
dfr.res.quad <- data.frame(dat1[,c(1,2)],res.quad)
names(dfr.res.quad) <- c("X","Y","ZVres")
coordinates(dfr.res.quad) <- c("X", "Y")
#
#   Residuenvariogramme
#
vg <- variogram(ZVres~1,dfr.res.lin)
plot(vg, plot.numbers = TRUE, pch = 15, cex = 1.2
, xlim = range(0, 1.1*max(vg$dist)), main = "Residuals\nlinear trend")
pause()
vg <- variogram(ZVres~1,dfr.res.quad)
plot(vg, plot.numbers = TRUE, pch = 15, cex = 1.2
, xlim = range(0, 1.1*max(vg$dist)), main = "Residuals\nquadratic trend")
vg <- variogram(ZVres~1,dfr.res.lin,cressie = TRUE)
plot(vg, plot.numbers = TRUE, pch = 15, cex = 1.2
, xlim = range(0, 1.1*max(vg$dist)), main = "Residuals robust\nlinear trend")
vg <- variogram(ZVres~1,dfr.res.quad,cressie = TRUE)
plot(vg, plot.numbers = TRUE, pch = 15, cex = 1.2
, xlim = range(0, 1.1*max(vg$dist)), main = "Residuals robust\nquadratic trend")