library(fields)
library(MASS)
library(prodlim)

### 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..

data=read.table(file='http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/climatol-ann.txt')

Y=data$V4      #ann avg precip at 246 locations
Y = Y/10		#convert to cm
X=data[,1:3]   #Parameters (log,lat,elev)
names(X)=c("lon","lat","elev")

lat=data$V2
lon=data$V1
elev=data$V3
precip=data$V4/10   #convert to cm

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=TRUE)
 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 # 
 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=5000,
			 minFactor=1e-12) )
 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 = rho*(1 - exp(-1*rdist(xs,xs)/theta))
 K1 = K + diag( sigma^2, nobs)
 #Kstar<- Matern(rdist( xgrid, x)/.3, smoothness=1.5)
 Kstar = K
 nobs=length(yHW1)
ghat<- Kstar %*% solve( K1) %*% 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..

#### Create the Rajeevan Grid with elevation..

## read te India DEM and estimate on it
### Create Rajeevan grid with elevation
## read the Rajeevan grid (lat and lon)

### India and Central Asia topo
test1=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/india-grid-topo.txt"),ncol=3,byrow=T)

### Rajeevan grid
test2=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/Rajeevan-grid.txt"),ncol=2,byrow=T)
### put longitude first
test3=cbind(test2[,2],test2[,1])

xv1=data.frame(test1[,1:2])
xv2=data.frame(test3[,1:2])

## find common points in the top grid that contains
## Rajeevan grid points and the missing points
##

zzint=row.match(xv2,xv1)
indexnotna = which(!is.na(zzint))
indexna = which( is.na(zzint) )

## create the new rajeevan grid
rajeevgridel = cbind(test3[indexnotna,1:2],test1[zzint[indexnotna],3])  

lat = rajeevgridel[,2]
lon = rajeevgridel[,1]
elev = rajeevgridel[,3]

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


#Kstar = sigma^2 + rho*(1 - exp(-1*rdist(xps,xs)/theta))
Kstar = rho*(1 - exp(-1*rdist(xps,xs)/theta))
ghat<- Kstar %*% solve( K1) %*% yHW1
kse = sqrt(var(yHW1) - diag((Kstar) %*% solve( K1) %*% 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..

quilt.plot(lon,lat,ghat)   #spatial estimates
quilt.plot(lon,lat,kse)		#standard error

##############################################################
### Model 2  - Lecture 4
#### Spatial trend with a linear model + Kriging residuals #####
### Estimated together..
## equivalent to fitting a lm on precip and lat and lon and
### Kriging the residuals

## xs and xse already defined..
yHW1 = precip
 nobs=length(yHW1)
xx = cbind(rep(1,nobs),xs)

nobs = length(yHW1)

sigma=1


 look<- vgram(xs, yHW1, N=15, lon.lat=TRUE)
 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))
K = 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
chat = solve( rho*K + diag( sigma^2, nobs)) %*% (yHW1-xx%*%dhat)

ghat<- xx%*%dhat + rho*Kstar %*% chat

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

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

#### Create the Rajeevan Grid with elevation..

## read te India DEM and estimate on it
### Create Rajeevan grid with elevation
## read the Rajeevan grid (lat and lon)

### India and Central Asia topo
test1=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/india-grid-topo.txt"),ncol=3,byrow=T)

### Rajeevan grid
test2=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/Rajeevan-grid.txt"),ncol=2,byrow=T)
### put longitude first
test3=cbind(test2[,2],test2[,1])

xv1=data.frame(test1[,1:2])
xv2=data.frame(test3[,1:2])

## find common points in the top grid that contains
## Rajeevan grid points and the missing points
##

zzint=row.match(xv2,xv1)
indexnotna = which(!is.na(zzint))
indexna = which( is.na(zzint) )

## create the new rajeevan grid
rajeevgridel = cbind(test3[indexnotna,1:2],test1[zzint[indexnotna],3])  

lat = rajeevgridel[,2]
lon = rajeevgridel[,1]
elev = rajeevgridel[,3]

nlocs = length(lat)


xps=rajeevgridel[,1:2]
xpse = cbind(lon,lat,elev)

Kstar = rho*(1 - exp(-1*rdist(xps,xs)/theta))

xx = cbind(rep(1,nlocs),xps)

ghat<- xx%*%dhat + rho*Kstar %*% chat

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=2)
yp = predict.Krig(zz,x=xps,drop.Z=TRUE)


### spatial map of estimates and errors - modify to add beautify..

quilt.plot(lon,lat,ghat)   #spatial estimates
quilt.plot(lon,lat,kse)		#standard error

#############################################
#### Model 3
############## GLM on elevation - elev
### to get the mean function + Kriging the residuals
#Explicitly

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


lat=data$V2
lon=data$V1
elev=data$V3
precip = data$V4/10   #convert to cm

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

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


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


sigma=1
 look<- vgram(xs, yHW1, N=15, lon.lat=TRUE)
 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 = rho*(1 - exp(-1*rdist(xs,xs)/theta))
 K1 = K + diag( sigma^2, nobs)
 
 #Kstar<- Matern(rdist( xgrid, x)/.3, smoothness=1.5)
 Kstar = K
 nobs=length(yHW1)
ghat1 <- Kstar %*% solve( K1) %*% yHW1
ghat = predict(zglm) + ghat1

## 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..


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

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

#### Create the Rajeevan Grid with elevation..

## read te India DEM and estimate on it
### Create Rajeevan grid with elevation
## read the Rajeevan grid (lat and lon)

### India and Central Asia topo
test1=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/india-grid-topo.txt"),ncol=3,byrow=T)

### Rajeevan grid
test2=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/Rajeevan-grid.txt"),ncol=2,byrow=T)
### put longitude first
test3=cbind(test2[,2],test2[,1])

xv1=data.frame(test1[,1:2])
xv2=data.frame(test3[,1:2])

## find common points in the top grid that contains
## Rajeevan grid points and the missing points
##

zzint=row.match(xv2,xv1)
indexnotna = which(!is.na(zzint))
indexna = which( is.na(zzint) )

## create the new rajeevan grid
rajeevgridel = cbind(test3[indexnotna,1:2],test1[zzint[indexnotna],3])  

lat = rajeevgridel[,2]
lon = rajeevgridel[,1]
elev = rajeevgridel[,3]

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

Kstar = rho*(1 - exp(-1*rdist(xps,xs)/theta))
ghat<- Kstar %*% solve( K1) %*% yHW1
ghat = ghat + predict(zglm,newdata=as.data.frame(elev))

kse1 = predict(zglm,newdata=as.data.frame(elev),se.fit=TRUE)$se.fit

kse2 = sqrt(var(yHW1) - diag((Kstar) %*% solve( K1) %*% t(Kstar)))

kse = sqrt(kse1^2 + kse2^2)



#Kstar = sigma^2 + rho*(1 - exp(-1*rdist(xpse,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 
## Might have to ply with rho, theta..
zz=Krig(xs,precip,Z=xse[,3],rho=rho,theta=theta,sigma2=sigma,m=2)
yp = predict.Krig(zz,x=xps,Z=elev)


### spatial map of estimates and errors - modify to add beautify..

quilt.plot(lon,lat,ghat)   #spatial estimates
quilt.plot(lon,lat,kse)		#standard error

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