library(fields)
library(maps)
library(prodlim)

#### Example spatial modeling on India Rainfall data.
test=read.table("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/climatol-ann.txt")
X = test[,1:2]
Y=test[,4]
Z1=test[,3]
fit1=spatialProcess(X,Y,Z=Z1)	#with elevation
set.panel(2,2)
plot(fit1)	#diagnostics plot

### Plot the surface quickly..
surface(fit1,ylim=range(0,35))
map('world2',wrap=c(-180,180), add=TRUE)
grid()
###

####  Predict on Rajeevan deg grid over India with Elevation

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

### Rajeevan grid over India is 1 x 1 degree grid that does not contain elevation
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 elevation
## 
zzint=row.match(xv2,xv1)
indexnotna = which(!is.na(zzint))
indexna = which( is.na(zzint) )

## create the new 0.5 degree grid over India
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)

## Predict using the spatial process model fitted
yp = predict(fit1,x=xps,Z=elev)

## quilt plot
par(mfrow=c(1,1))
quilt.plot(lon, lat, yp)
map('world2',wrap=c(-180,180), add=TRUE)
grid()


#################### Simple Kriging

Yscale = scale(Y)
yy=Krig(X,Yscale,m=1)
plot(yy)
map('world2',wrap=c(-180,180), add=TRUE)
#grid(72,24,lwd=1,lty=2)
grid()

######