library(fields)
library(MASS)

### Doug's Lecture Notes 1 through 4
### Chapter 8 of Springer online book
#### Applied Spatial Data Analysis with R
#Authors:
#    Roger S. Bivand,
#    Edzer Pebesma,
#    Virgilio Gómez-Rubio 

### Three Models 
#### Fit a Kriging (i.e. spatial) Model to the precipitation data
####  Predict on the DEM grid

### Model 2  ### Spatial Additive Model - Linear term for the mean function
### plus the spatial model on the residual

#### Model 3 - Fit a GLM to lat, lon and elevation and fit a spatial model to
#### the residuals

### Model 2 estimates the mean component and residual spatial model parameters
### simultaneously, while Model 3 does it in two steps - i.e. Hierarchy!.

#### Spatial trend, mean function all refer to the same..

#### Kriging is an exact estimator - i.e., at the observation locations Yhat will
### be exactly equal to Y without nugget effect, but with one it will be close.

### Therefore you have to do cross-validation to evaluate the performance.

## You can write a loop to perform cross validation for all these models

### Problem 8 - Model 1

### Problem 9 - Models 2 and 3

##### Model 1 - Lectures 2B and 3 
## Kriging on the precipitation..

colopcp=read.table("http://cires.colorado.edu/~aslater/CVEN_6833/colo_precip.dat")

lat=colopcp$V1
lon=colopcp$V2
elev=colopcp$V3
precip=colopcp$V4
vars = cbind(lon,lat,elev)
vars=as.data.frame(vars)
varsdata = vars

fitall<-glm(precip~., data=vars)

xs = cbind(lon,lat)
xse = cbind(lon,lat,elev)

yHW1 = precip

### Also If Kriging with elevation then replace xps with xpse below

## Compute empirical variogram and plot ...
par(mar=c(4,4,1,1))
 look<- vgram(xs, yHW1, N=15, lon.lat=FALSE)
 bplot.xy( look$d, look$vgram, breaks= look$breaks,
          outline=FALSE,
          xlab="distance (degrees)", ylab="variogram")
 points( look$centers, look$stats[2,], col="blue")
 
# fit of exponential by nonlinear least squares
 xd<- look$centers
 ym<- look$stats[2,]
 sigma<- 1.0 # 
 nls.fit<- nls( ym ~ sigma^2 + rho*( 1- exp(-xd/theta)),
               start= list(rho=max(look$stats[2,]), theta=quantile(xd,0.1)), control=list( maxiter=200) )
 pars<- coefficients(nls.fit)
 rho<-pars[1]
 theta<-pars[2]
 xr = round(max(xd)+sd(xd))
 dgrid<- seq(0,xr,,400)
 lines(dgrid, sigma^2 + rho*(1 - exp(-1*dgrid/theta)), col="blue", lwd=3)
 
 # same fit using maximum likelihood - if it works..
 #MLE.fit<- MLE.Matern( xs, yHW1, smoothness=.5, m=1)
 #pars<- MLE.fit$pars
 #rho<-pars[1]
 #theta<-pars[2]
 #sigma<-pars[3]
 #lines(dgrid, sigma^2 + rho*(1 - exp(-1*dgrid/theta)), col="green", lwd=3)
 
 ### Predict at observation points.. is sigma = 0 then it will be exact.
 K = sigma^2 + rho*(1 - exp(-1*rdist(xs,xs)/theta))
 #Kstar<- Matern(rdist( xgrid, x)/.3, smoothness=1.5)
 Kstar = K
 nobs=length(yHW1)
ghat<- Kstar %*% solve( K + diag( sigma^2, nobs)) %*% yHW1

## Check with Krig and predict.krig
zz=Krig(xs,yHW1,rho=rho,theta=theta,m=1,sigma2=sigma)
yp = predict.Krig(zz)

## Predict on the grid.. and standard error..
predpts=read.table("http://cires.colorado.edu/~aslater/CVEN_6833/colo_dem.dat")

lat=predpts$V1
lon=predpts$V2
elev=predpts$V3

xps=cbind(lon,lat)
xpse = cbind(lon,lat,elev)


Kstar = sigma^2 + rho*(1 - exp(-1*rdist(xps,xs)/theta))
ghat<- Kstar %*% solve( K + diag( sigma^2, nobs)) %*% yHW1
kse = sqrt(var(yHW1) - diag((Kstar) %*% solve( K + diag( sigma^2, nobs)) %*% t(Kstar)))

## Check with Krig and predict.Krig
zz=Krig(xs,yHW1,rho=rho,theta=theta,sigma2=sigma,m=1)
yp = predict.Krig(zz,x=xps,drop.Z=TRUE)

### spatial map of estimates and errors - modify to add beautify..
xlon = sort(unique(lon))
nrs = length(xlon)
ylat = sort(unique(lat))
ncs = length(ylat)
zmat = matrix(ghat, nrow=ncs, ncol= nrs) 
image.plot(xlon,ylat,t(zmat),zlim=range(precip))
contour(xlon,ylat,t(zmat),add=T)
zmat = matrix(kse, nrow=ncs, ncol= nrs) 
image.plot(xlon,ylat,t(zmat))
contour(xlon,ylat,t(zmat),add=T)

##############################################################
### Model 2  - Lecture 3
#### Spatial trend with a linear model + Kriging residuals #####
### Estimated together..

xx = cbind(rep(1,nobs),colopcp$V2,colopcp$V1)
xx = cbind(rep(1,nobs),colopcp$V2,colopcp$V1,colopcp$V3)
xe = cbind(colopcp$V2,colopcp$V1,colopcp$V3)

## select the desired xx
xs = cbind(colopcp$V2,colopcp$V1)

nobs = length(yHW1)
sigma=1

 look<- vgram(xs, yHW1, N=15, lon.lat=FALSE)
 bplot.xy(look$d, look$vgram, breaks= look$breaks,
          outline=FALSE,
          xlab="distance (degrees)", ylab="variogram")
 points( look$centers, look$stats[2,], col="blue")
 
# fit of exponential by nonlinear least squares
 xd<- look$centers
 ym<- look$stats[2,]
 sigma<- 1.0 # 
 nls.fit<- nls( ym ~ sigma^2 + rho*( 1- exp(-xd/theta)),
    start= list(rho=max(look$stats[2,]), theta=quantile(xd,0.1)), control=list( maxiter=200) )
 pars<- coefficients(nls.fit)
 rho<-pars[1]
 theta<-pars[2]
 xr=round(max(xd)+sd(xd))
 dgrid<- seq(0,xr,,400)
 lines(dgrid, sigma^2 + rho*(1 - exp(-1*dgrid/theta)), col="blue", lwd=3)

 
 K = sigma^2 + rho*(1 - exp(-1*rdist(xs,xs)/theta))
MM = K + diag(sigma^2, nobs)
dhat = solve(t(xx) %*% solve(MM) %*% xx) %*% t(xx) %*% solve(MM) %*% yHW1

 Kstar = K
 nobs=length(yHW1)
ghat<- xx%*%dhat + Kstar %*% solve( K + diag( sigma^2, nobs)) %*% (yHW1-xx%*%dhat)


## Predict on the grid.. and standard error..
predpts=read.table("http://cires.colorado.edu/~aslater/CVEN_6833/colo_dem.dat")

lat=predpts$V1
lon=predpts$V2
elev=predpts$V3
varpreds = cbind(lon,lat,elev)
npred = length(elev)

xps = cbind(lon,lat)
xp=cbind(rep(1,npred),lon,lat,elev)
Kstar = sigma^2 + rho*(1 - exp(-1*rdist(xps,xs)/theta))
ghat<- xp%*%dhat + Kstar %*% solve( K + diag( sigma^2, nobs)) %*% (yHW1-xx%*%dhat)
kse = sqrt(var(yHW1) - diag((Kstar) %*% solve( K + diag( sigma^2, nobs)) %*% t(Kstar)))


## Check with Krig and predict.Krig
zz=Krig(xe[,1:2],yHW1,Z=xe[,3],rho=rho,theta=theta,sigma2=sigma,m=2)
yp = predict.Krig(zz,x=varpreds[,1:2],Z=varpreds[,3],drop.Z=FALSE)

### spatial map of estimates and errors - modify to add beautify..
xlon = sort(unique(lon))
nrs = length(xlon)
ylat = sort(unique(lat))
ncs = length(ylat)
zmat = matrix(ghat, nrow=ncs, ncol= nrs) 
image.plot(xlon,ylat,t(zmat),zlim=range(precip))
contour(xlon,ylat,t(zmat),add=T)
zmat = matrix(kse, nrow=ncs, ncol= nrs) 
image.plot(xlon,ylat,t(zmat))
contour(xlon,ylat,t(zmat),add=T)

#############################################
#### Model 3
############## GLM to get the mean function + Kriging the residuals
#######################


lat=colopcp$V1
lon=colopcp$V2
elev=colopcp$V3
precip = colopcp$V4

xe = cbind(lon,lat,elev)
yHW1 = precip

### fit GLM and get residuals..
zglm = glm(yHW1 ~ ., data=as.data.frame(xe))
yres = residuals(zglm)


xs = cbind(colopcp$V2,colopcp$V1)
yHW1 = yres
nobs = length(yHW1)


sigma=1
 look<- vgram(xs, yHW1, N=15, lon.lat=FALSE)
 bplot.xy(look$d, look$vgram, breaks= look$breaks,
          outline=FALSE,
          xlab="distance (degrees)", ylab="variogram")
 points( look$centers, look$stats[2,], col="blue")
 
# fit of exponential by nonlinear least squares
 xd<- look$centers
 ym<- look$stats[2,]
 sigma<- 1.0 # 
 nls.fit<- nls( ym ~ sigma^2 + rho*( 1- exp(-xd/theta)),
    start= list(rho=max(look$stats[2,]), theta=quantile(xd,0.1)), control=list( maxiter=200) )
 pars<- coefficients(nls.fit)
 rho<-pars[1]
 theta<-pars[2]
 xr=round(max(xd)+sd(xd))
 dgrid<- seq(0,xr,,400)
 lines(dgrid, sigma^2 + rho*(1 - exp(-1*dgrid/theta)), col="blue", lwd=3)


 ### Predict at observation points.. is sigma = 0 then it will be exact.
 K = sigma^2 + rho*(1 - exp(-1*rdist(xs,xs)/theta))
 #Kstar<- Matern(rdist( xgrid, x)/.3, smoothness=1.5)
 Kstar = K
 nobs=length(yHW1)
ghat<- Kstar %*% solve( K + diag( sigma^2, nobs)) %*% yHW1
ghat = predict(zglm) + ghat

## Check with Krig and predict.krig
zz=Krig(xs,yHW1,rho=rho,theta=theta,m=1,sigma2=sigma)
yp = predict(zglm) + predict.Krig(zz)

## Predict on the grid.. and standard error..

predpts=read.table("http://cires.colorado.edu/~aslater/CVEN_6833/colo_dem.dat")

lat=predpts$V1
lon=predpts$V2
elev=predpts$V3
xe = cbind(lon,lat,elev)

xps=cbind(lon,lat) 
Kstar = sigma^2 + rho*(1 - exp(-1*rdist(xps,xs)/theta))
ghat<- Kstar %*% solve( K + diag( sigma^2, nobs)) %*% yHW1
kse = sqrt(var(precip) - diag((Kstar) %*% solve( K + diag( sigma^2, nobs)) %*% t(Kstar)))
ghat = ghat + predict(zglm,newdata=as.data.frame(xe))

## Check with Krig and predict.Krig
zz=Krig(xs,yHW1,rho=rho,theta=theta,sigma2=sigma,m=1)
yp = predict(zglm,newdata=as.data.frame(xe))+ predict.Krig(zz,x=xps,drop.Z=TRUE)


### spatial map of estimates and errors - modify to add beautify..
xlon = sort(unique(lon))
nrs = length(xlon)
ylat = sort(unique(lat))
ncs = length(ylat)
zmat = matrix(ghat, nrow=ncs, ncol= nrs) 
image.plot(xlon,ylat,t(zmat),zlim=range(precip))
contour(xlon,ylat,t(zmat),add=T)
zmat = matrix(kse, nrow=ncs, ncol= nrs) 
image.plot(xlon,ylat,t(zmat))
contour(xlon,ylat,t(zmat),add=T)

######################