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

library(sp)
library(lattice)
library(gstat)



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

###  - Model 1

### 

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

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

#ann avg precip at 246 locations

data1[,4]=data1[,4]/10 #convert to cm

names(data1)=c("lon","lat","elev","rain")

lat=data1$lat
lon=data1$lon
elev=data1$elev
precip=data1$rain   

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



X=data1[,1:3]   #Parameters (log,lat,elev)
Y = data1[,4]

coordinates(data1) = ~lon+lat

v=variogram(rain ~ lon+lat, data1)
v=variogram(rain ~ 1, data1)
vv = fit.variogram(v,vgm("Sph"))
plot(v,vv)

plot(v)	#plot the variogram


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

yHW1 = precip


 ### Predict at observation points.. is sigma = 0 then it will be exact.
 
 sigma=1 #add a small nugget
 rho = vv[2,2]
 theta=vv[2,3]
 sigma=vv[1,2]
 sigma2=sigma^2
 
plot(yHW1, ghat)

## Check with Krig and predict.krig
zk = krige(rain ~ lon+lat, data1, data1, model=vv)
plot(yHW1, zk$var1.pred)

zz=mKrig(xs,yHW1,aRange=theta, lambda=sigma2/rho,m=1)
surface(zz)
plot(yHW1, zz$fitted.values)


## 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-2019/HW-1/india-grid-topo.txt"),ncol=3,byrow=T)

### Rajeevan grid
test2=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs-2019/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)


## Check with Krig and predict.Krig
zz=mKrig(xs,yHW1, aRange=theta, lambda=sigma2/rho,m=1)
yp = predict(zz,xnew=xps)
kse = predictSE(zz, xnew=xps)



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

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

##############################################################
### Model 2  - Lecture 4