library(leaps)
library(MPV)
library(MASS)
library(prodlim)
library(akima)
library(fields)
source("GLM_fit.txt")

######################
#### Read in Data ####
######################

## Read data as table from local directory 
precip = read.table("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/climatol-ann.txt")      # ann avg precip 246 x 4
grid = read.table("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/india-grid-topo.txt")     # high res DEM 3384 x 3
rajeevan = read.table("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/Rajeevan-grid.txt")   # rajeevan grid 357 x 2

## Organize data into a data frame with names for columns/variables:
precip_dt = data.frame(precip)                          # ann avg precip at 246 locations
colnames(precip_dt) = c("lon", "lat", "elev", "pmm")
india_grid = data.frame(grid)                           # high resolution DEM (more than just India)
colnames(india_grid) = c("lon", "lat", "elev")
raj_grid = data.frame(rajeevan[,2], rajeevan[,1])       # lat/lon for just India
colnames(raj_grid) = c("lon","lat")

########## Add this the the list of functions!!! ######## WJR

######################################
#### i. Fit Best Regression Model ####
######################################

# Define variables Y and X
pmm = precip_dt[,4]     # indepedent variable
X = precip_dt[,1:3]     # all the predictor set
X = as.data.frame(X)
family = "gamma"

# Find best GLM model 
bestmod = GLM_fit(pmm, X, family)

gamshape = gamma.shape(bestmod)[1]


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

 X = as.data.frame(rajeevgridel)
 names(X)=c("lon","lat","elev")
 ypred <- predict(bestmod,newdata=X,se.fit=TRUE, response=TRUE)
predscale = ypred$fit/as.numeric(gamshape)
ypredfit = qgamma(0.5,scale = predscale,shape=gamshape)
ypredfit[ypredfit == "NaN"]=0

quilt.plot(X[,1],X[,2],ypredfit,xlab="Longitude (m)",ylab="Latitude (m)",main='Predicted Precipitation on DEM Grid (mm)')
 world(add=T,lwd=1)

quilt.plot(X[,1],X[,2],ypred$se.fit,xlab="Longitude (m)",ylab="Latitude (m)",main='Predicted Standard Error on DEM Grid (mm)')
 world(add=T,lwd=1)