library(fields); 
library(geoR);  
library(MASS)
#install.packages("akima")
library(akima)
# The rdist function is used to find the distance between all the grid points.

n=20
kombi = expand.grid(seq(1,n, by=1), seq(1,n, by=1)) #Make n x n grid
dist = rdist(kombi, kombi)

#Sample explanatoryvariabl elevation h(s) 
#Matern Correlation function
set.seed(100)
theta1E=1
mynuE=10
sigma1E=70000
sigmaElev=sigma1E*matrix(Matern(dist, range = theta1E,nu=mynuE),ncol=n*n, nrow=n*n)
muE  = rep(500,n*n) #Expected value

OurElev = mvrnorm(n = 1, muE, sigmaElev) #Sample one realisation.
# Display the altitude
Result=interp(x=kombi$Var1,y=kombi$Var2,z=OurElev)
image.plot(Result, asp=1)
title("Elevation (m)")

#Sample true temperature
beta0=0
beta1=-0.5/100 #0.5 deg C less per 100 m elevation
muTemp=0+beta1*OurElev

set.seed(200)
theta1T=1
mynuT=5
sigma1T=2
sigmaT=sigma1T*matrix(Matern(dist, range = theta1T,nu=mynuT),ncol=n*n, nrow=n*n)
delta=mvrnorm(n = 1, rep(0,n*n), sigmaT)
OurTrueTemp=muTemp + delta
# Display the altitude
Result=interp(x=kombi$Var1,y=kombi$Var2,z=OurTrueTemp)
image.plot(Result, asp=1, title("Temperature"))
title("Temperature (deg C)")

#Sample observations at locations (5,3), (20,15), (16,18)
s1=45
s2=300
s3=356

muZ=c(OurTrueTemp[s1], OurTrueTemp[s2], OurTrueTemp[s3])
sigmaEpsilon=0.5
ObsTemp=mvrnorm(n=1,mu=muZ,Sigma=sigmaEpsilon*diag(3))
locWithObs=matrix(data=0,nrow=3,ncol=2)
locWithObs[1,]=c(kombi$Var1[s1], kombi$Var2[s1])
locWithObs[2,]=c(kombi$Var1[s2], kombi$Var2[s2])
locWithObs[3,]=c(kombi$Var1[s3], kombi$Var2[s3])

muZ
ObsTemp

locToPredict=c(10,10)
nlocWithObs=3

distLocWithObs=rdist(locWithObs, locWithObs)
Sigma22=sigma1T*matrix(Matern(distLocWithObs, range = theta1T,nu=mynuT), ncol=nlocWithObs, nrow=nlocWithObs)+sigmaEpsilon*diag(3) 
Sigma11=sigma1T
distLocObsPred=rdist(locWithObs, locToPredict)
Sigma12=sigma1T*matrix(Matern(distLocWithObs, range = theta1T,nu=mynuT), ncol=nlocWithObs, nrow=1)

muPred=beta0+beta1*OurElev[10*n+10-1]
condMu=muPred-Sigma12%*%solve(Sigma22)%*%(muZ-ObsTemp  )
condMu
#Conditional variance
condVar=Sigma11 - Sigma12%*%solve(Sigma22)%*%c(Sigma12)
sqrt(condVar)
OurTrueTemp[(10*n+10-1)]