rm(list=ls())
library(DAAG)
n.points <- 6
#
# das Kovarianzmodell
#
covmod <- function(h){
  covmod <- exp(-abs(h)/5)
}
#
par(mfrow=c(1,1))
plot(1,1,type="n",xlim=c(0,10),ylim=c(0,10),xlab="x",ylab="y",asp=1)
pause()
#
# Anklicken der Punkte des "Messgitters"
#
p.x <- 1:n.points
p.y <- 1:n.points
for (i in 1:n.points){
  p <- locator(1)
  p.x[i] <- p$x
  p.y[i] <- p$y
  points(p,cex=2.5,pch=20,col="darkgreen")
}
p.coords <- matrix(c(p.x,p.y),nrow=2,byrow=T)
#
# Aufbau der Krigematrix und der rechten Seite des LGS
#
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]
    D[i,j] <- covmod(sqrt(dist%*%dist))
  }
}
D.inv <- solve(D)
#
# Messwerte:
#
T.mean = 19
T.long <- c(15,24,17,25,18,23,19,22,16,25) 
# wenn n.points zu groß wird, müssen hier mehr Werte stehen
T <- T.long[1:n.points]
#===============================================================================
# Karte über das gesamte Gebiet
#-------------------------------------------------------------------------------
#   Erstellen eines Grids
#
ndataX <- 75
ndataY <- 75
x.min <- 0; x.max <- 10
y.min <- 0; y.max <- 10
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")
#stop()
#plot(dgr,cex=0.2,pch=16,asp=1)
#-------------------------------------------------------------------------------
# Berechnung der Schätzungen und der Schätzvarianz
#-------------------------------------------------------------------------------
n <- length(dgr[,1])
T.est <- 1:n
T.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(sqrt(dist%*%dist))
  }
  lamda <- D.inv%*%rs
  T.est[i] <- (T-T.mean)%*%lamda + T.mean
  T.var[i] <- covmod(0) - as.vector(lamda)%*%as.vector(rs)
}
results <- data.frame(w,T.est,T.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(T.est),max(T.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,], T, pmat = rotmat), col = 2, pch =16)
#
v <- results[,4]
datmat <- matrix(v,byrow=T,nrow=ndataX)
persp(x,y,datmat, zlim = c(min(T.var),max(T.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)

