rm(list=ls())
library(DAAG)
#===============================================================================
# das Variogrammodell
#===============================================================================
svg.sph <- function(nugg,sill,rnge,x){
  if (x <= rnge) {
    xhlp <- nugg + sill*(1.5*x/rnge -0.5*(x/rnge)^3)
  } else {
    xhlp <- nugg + sill
  }
  if (x == 0) {
      xhlp <- 0
  }
#      
  svg.sph <- xhlp
}
#
# und noch das zugehörige Kovarianzmodell:
#
covmod <- function(nugg,sill,rnge,x) {
  covmod <- nugg + sill - svg.sph(nugg,sill,rnge,x)
}
#===============================================================================
# Einlesen und bearbeiten der Messdaten
#===============================================================================
filename <- file.choose()
dat <- read.table(filename,header=F)
names(dat) <- c("X","Y","P","EC","pH","NitCo","NitFlx")
dat$NitFlx <- dat$P*dat$NitCo                             
cln <- 7                  # Spaltennummer
dat <- dat[,c(1,2,cln)]   # Auswahl der Spalte cln
# Bem: dat enthält jetzt noch genau 3 Spalten: X,Y,und 1 Spalte mit Messwerten 
cc <- complete.cases(dat) # Ermittlung fehlender Werte
dat1 <- dat[cc,]          # Ausschluss fehlender Werte
vartext <- names(dat1)[3] # Spaltenüberschrift für Grafik
ZV <- dat1[,3]            # Verwenden der Spalte 3 für das Variogramm
X <- dat1[,1]
Y <- dat1[,2]
#===============================================================================
# Berechnung des empirischen Variogramms
#===============================================================================
# Variogram-Cloud
#--------------------------------------------------------------------------
n <- length(X)
n.count <- 0
n.pairs <- (n-1)*n/2         # Anzahl der Paarkombinationen
h.var   <- rep(0,n.pairs)    # Array für Punktabstände
semvar  <- rep(0,n.pairs)    # Array für Semivarianzen
for (i in 1:(n-1)) {
  for (j in (i+1):n) {
    n.count <- n.count + 1
    h.var[n.count]  <- sqrt((X[j]-X[i])^2 + (Y[j]-Y[i])^2)
    semvar[n.count] <- 0.5*(ZV[i]-ZV[j])^2 
}}
#--------------------------------------------------------------------------
# Variogram-empirisch-isotrop
#--------------------------------------------------------------------------
n.hclass <- 20                   # Anzahl der lag-Klassen
count.class <- rep(0,n.hclass)   
sem.class   <- rep(0,n.hclass)
dh <- 1.01*max(h.var)/n.hclass   # Inkrement der lag-Klassen
#
# Klassenmitten der lag-Klassen
#
h.class.mean <- c(-dh/2 + dh*1:(n.hclass)) # immer nachprüfen, ob das stimmt!!
#
# Zählen der Wertepaare in den Klassen:
#
for (i in 1:length(semvar)) {
  nr <- trunc(h.var[i]/dh) + 1
  count.class[nr] <- count.class[nr] + 1
  sem.class[nr]   <- sem.class[nr] + semvar[i]
} 
semvar.class <- sem.class/count.class
#
# Plots:
#
par(mfrow=c(1,2),mar=c(5 , 4, 4, 2) + 0.1) # mar macht Ränder um die Graphik
#
# counts in den Klassen:
#
plot(h.class.mean,count.class,main=paste("classcounts: ",vartext))
nvar.max <- trunc(0.5*n.hclass)
#
# Semivarianzen in den Klassen:
#
plot(h.class.mean[1:nvar.max],semvar.class[1:nvar.max],
ylim=c(0,max(semvar.class[1:nvar.max])),pch=15,col="blue",
xlab="lag",ylab="semivariance",main=paste("variogram: ",vartext))
#===============================================================================
# Anpassung des Variogrammmodells mit nls()
#===============================================================================
# Wichtungsfaktoren:
#wft <- count.class/(h.class.mean^2)
wft <- count.class/(semvar.class^2)
sum.wft <- sum(wft)
wft <- 100*wft/sum.wft  # Wichtungsfaktoren in %
#
n.max <- nvar.max #: siehe einige Zeilen weiter oben
x <- h.class.mean[1:n.max]
s <- semvar.class[1:n.max]
wts <- wft[1:n.max]
#
# Anfangsschätzungen optisch: hier muss man von Hand eingreifen!
#
ng <- 100
sl <- 300
rg <- 500
model <- nls(s~svg.sph(nugg,sill,rnge,x),model=TRUE,weigths=wts,
             start=list(nugg=ng,sill=sl,rnge=rg))
print(summary(model))
par(mfrow=c(1,1))
lag <- seq(0,max(x),20)
#
# Plot Modell gegen Messung
#
plot(x,s,ylim=c(0,1.2*max(s)),pch=15,col="blue",cex=1.5,
main="Fit: Sphärisches Modell")
#
# Abgreifen der optimierten Parameter aus dem angepassten Modell
#
par <- model$m$getAllPars()
ngo <- par[1]
slo <- par[2]
rgo <- par[3]
n.est <- length(lag)
svest <- 1:n.est
for(i in (1:n.est)) {
  svest[i] <- svg.sph(ngo,slo,rgo,lag[i]) 
  # hab ich ohne Schleife nicht hingekriegt
} 
lines(lag,svest,lty=1,lwd=1,col="gray4")
#===============================================================================
# Aufbau der Krigematrix
#===============================================================================
p.coords <- t(as.matrix(dat1[,c(1,2)]))
n.points <- length(ZV)
D <- matrix(1:n.points^2,ncol=n.points)
rs <- 1:n.points
for(i in 1:n.points) {
  for(j in 1:n.points) {
    dist <- p.coords[,i]-p.coords[,j]
    cat(i,j,"\n")
    D[i,j] <- covmod(ngo,slo,rgo,sqrt(dist%*%dist))
  }
}
D.inv <- solve(D)
#===============================================================================
# Karte über das gesamte Gebiet
#===============================================================================
#   Erstellen eines Grids
#
ndataX <- 75
ndataY <- 75
x.min <- min(X); x.max <- max(X)
y.min <- min(Y); y.max <- max(Y)
x <- seq(x.min,x.max,(x.max-x.min)/(ndataX-1))
y <- seq(y.min,y.max,(y.max-y.min)/(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)
gridmat <- as.matrix(dgr)
plot(w$X1,w$X2,cex=0.2,pch=16,asp=1,xlab="X",ylab="Y")
#plot(dgr,cex=0.2,pch=16,asp=1)
#-------------------------------------------------------------------------------
# Berechnung der Schätzungen und der Schätzvarianz
#-------------------------------------------------------------------------------
n <- length(dgr[,1])
ZV.est <- 1:n
ZV.var <- 1:n
#
# Berechnung der rechten Seiten und der Lösung des Gleichungssystems 
#
for (i in (1:n)) {
  for(j in 1:n.points) {
    dist <- gridmat[i,]-p.coords[,j]
    rs[j] <- covmod(ngo,slo,rgo,sqrt(dist%*%dist))
  }
  lamda <- D.inv%*%rs
#
# Hier kommen zum ersten Mal die Messwerte ins Spiel:
#
ZV.mean = mean(ZV)
  ZV.est[i] <- (ZV-ZV.mean)%*%lamda + ZV.mean
  ZV.var[i] <- covmod(ngo,slo,rgo,0) - as.vector(lamda)%*%as.vector(rs)
}
results <- data.frame(w,ZV.est,ZV.var)
#===============================================================================
#  Contourplots
#===============================================================================
v <- results[,3]
datmat <- matrix(v,byrow=T,nrow=ndataX)
filled.contour(x,y,datmat,asp=1,color = topo.colors,nlevels=50,
, main="Simple Kriging: Schätzwerte"
,plot.axes={axis(1);axis(2)
; contour(x,y,datmat,nlevels=20,add=T,lwd=1,col="brown", axes=F)
; points(p.coords[1,],p.coords[2,],pch=16,col="white",cex=1,lwd=1)
})
v <- results[,4]
datmat <- matrix(v,byrow=T,nrow=ndataX)
filled.contour(x,y,datmat,asp=1,color = topo.colors,nlevels=50,
, main="Simple Kriging: Schätzvarianz"
,plot.axes={axis(1);axis(2)
; contour(x,y,datmat,nlevels=20,add=T,lwd=1,col="brown", axes=F)
; points(p.coords[1,],p.coords[2,],pch=16,col="white",cex=1,lwd=1)
})
#===============================================================================
#  3D-Plots
#===============================================================================
v <- results[,3]
datmat <- matrix(v,byrow=T,nrow=ndataX)
persp(x,y,datmat, zlim = c(min(ZV.est),max(ZV.est)),theta = 30, phi = 20, expand = 0.7, col = "lightblue",
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "X-Value", ylab = "\nY-Value", zlab = "Temp", main="Simple Kriging: Schätzwerte",
     r=4,border="lightblue") -> rotmat
#   Hinzufügen gemessener Datenpunkte
points(trans3d(p.coords[1,],p.coords[2,], ZV, pmat = rotmat), col = 2, pch =16)
#
v <- results[,4]
datmat <- matrix(v,byrow=T,nrow=ndataX)
persp(x,y,datmat, zlim = c(0,max(ZV.var)),theta = 30, phi = 20, expand = 0.7, col = "lightblue",
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "X-Value", ylab = "\nY-Value", zlab = "Temp", main="Simple Kriging: Schätzvarianz",
     r=4,border="lightblue") -> rotmat
points(trans3d(p.coords[1,],p.coords[2,], rep(0,n.points), pmat = rotmat), col = 2, pch =16)

