#This script interpolates station values for the Oregon case study using a two stage regression.

#This program loads the station data from a shape file and perform 8 regressions using general additive model (GAM) followed by kriging on the residuals.

#It uses LST monthly averages as input variables.

#Note that this program:

#1)assumes that the shape file is in the current working 

#This scripts predicts tmax using ing GAM and LST derived from MOD11A1.

#Script created by Benoit Parmentier on May 6, 2012. 

library(gstat)

#path<-"/data/computer/parmentier/Data/IPLANT_project/data_Oregon_stations"

path<-"H:/Data/IPLANT_project/data_Oregon_stations"

infile2<-"dates_interpolation_03052012_2dates_test.txt"                                          #List of 10 dates for the regression

out_prefix<-"_05062012_Kr_LST"

CRS<-proj4string(ghcn)

mean_LST<- readGDAL("mean_day244_rescaled.rst")  #This reads the whole raster in memory and provide a grid for kriging

proj4string(mean_LST)<-CRS #Assigning coordinates information

results <- matrix(1,length(dates),14)            #This is a matrix containing the diagnostic measures from the GAM models.

results_AIC<- matrix(1,length(dates),length(models)+2)  

results_GCV<- matrix(1,length(dates),length(models)+2)

results_DEV<- matrix(1,length(dates),length(models)+2)

results_RMSE<- matrix(1,length(dates),length(models)+2)

results_RMSE_kr<- matrix(1,length(dates),length(models)+2)

cor_LST_LC1<-matrix(1,10,1)      #correlation LST-LC1

cor_LST_LC3<-matrix(1,10,1)      #correlation LST-LC3

cor_LST_tmax<-matrix(1,10,1)    #correlation LST-tmax

#coords<- ghcn[,c('x_OR83M','y_OR83M')]

  mod2<- gam(tmax~ s(lat,lon,ELEV_SRTM), data=data_s)

  #Automate this using a data frame of size??

  #Saving dataset in dataframes: residuals from RMSE

#   data_v$mod1<-y_mod1$fit

#   data_v$mod2<-y_mod2$fit

#   data_v$mod3<-y_mod3$fit

#   data_v$mod4<-y_mod4$fit

#   data_v$mod5<-y_mod5$fit

#   data_v$mod6<-y_mod6$fit

#   data_v$mod7<-y_mod7$fit

#   

#   data_s$mod1<-mod1$fit

#   data_s$mod2<-mod2$fit

#   data_s$mod3<-mod3$fit

#   data_s$mod4<-mod4$fit

#   data_s$mod5<-mod5$fit

#   data_s$mod6<-mod6$fit

#   data_s$mod7<-mod7$fit

#   

  data_v$res_mod1<-as.numeric(res_mod1)

  data_v$res_mod2<-as.numeric(res_mod2)

  data_v$res_mod3<-as.numeric(res_mod3)

  data_v$res_mod4<-as.numeric(res_mod4)

  data_v$res_mod5<-as.numeric(res_mod5)

  data_v$res_mod6<-as.numeric(res_mod6)

  data_v$res_mod7<-as.numeric(res_mod7)

  data_s$res_mod1<-as.numeric(mod1$residuals)

  data_s$res_mod2<-as.numeric(mod2$residuals)

  data_s$res_mod3<-as.numeric(mod3$residuals)

  data_s$res_mod4<-as.numeric(mod4$residuals)

  data_s$res_mod5<-as.numeric(mod5$residuals)

  data_s$res_mod6<-as.numeric(mod6$residuals)

  data_s$res_mod7<-as.numeric(mod7$residuals)

  ###BEFORE Kringing the data object must be transformed to SDF

  coords<- data_v[,c('x_OR83M','y_OR83M')]

  coordinates(data_v)<-coords

  proj4string(data_v)<-CRS  #Need to assign coordinates...

  coords<- data_s[,c('x_OR83M','y_OR83M')]

  coordinates(data_s)<-coords

  proj4string(data_s)<-CRS  #Need to assign coordinates..

  #Kriging residuals!!

  for (j in 1:length(models)){

    name<-paste("res_mod",j,sep="")

    data_s$residuals<-data_s[[name]]

    X11()

    hscat(residuals~1,data_s,(0:9)*20000) # 9 lag classes with 20,000m width

    v<-variogram(residuals~1, data_s)

    plot(v)

    v.fit<-fit.variogram(v,vgm(1,"Sph", 150000,1))

    res_krige<-krige(residuals~1, data_s,mean_LST, v.fit)#mean_LST provides the data grid/raster image for the kriging locations.

    # Kriging visualization of Residuals fit over space

    #spplot.vcov(co_kriged_surf)                           #Visualizing the covariance structure

    res_krig1_s <- overlay(res_krige,data_s)             #This overlays the kriged surface tmax and the location of weather stations

    res_krig1_v <- overlay(res_krige,data_v)             #This overlays the kriged surface tmax and the location of weather stations

    name2<-paste("pred_kr_mod",j,sep="")

    #Adding the results back into the original dataframes.

    data_s[[name2]]<-res_krig1_s$var1.pred

    data_v[[name2]]<-res_krig1_v$var1.pred  

    #Calculate RMSE and then krig the residuals....!

    res_mod_kr_s<- data_s$tmax - data_s[[name2]]           #Residuals from kriging.

    res_mod_kr_v<- data_v$tmax - data_v[[name2]]           #Residuals from cokriging.

    RMSE_mod_kr_s <- sqrt(sum(res_mod_kr_s^2,na.rm=TRUE)/(nv-sum(is.na(res_mod_kr_s))))                  #RMSE from kriged surface.

    RMSE_mod_kr_v <- sqrt(sum(res_mod_kr_v^2,na.rm=TRUE)/(nv-sum(is.na(res_mod_kr_v))))                  #RMSE from co-kriged surface.

    #(nv-sum(is.na(res_mod2)))

    #Writing out results

    results_RMSE_kr[i,1]<- dates[i]  #storing the interpolation dates in the first column

    results_RMSE_kr[i,2]<- ns        #number of stations used in the training stage

    results_RMSE_kr[i,j+2]<- RMSE_mod_kr_v

    #results_RMSE_kr[i,3]<- res_mod_kr_v

    name3<-paste("res_kr_mod",j,sep="")

    data_s[[name3]]<-res_mod_kr_s

    data_v[[name3]]<-res_mod_kr_v #Writing residuals from kriging

    }

  data_name<-paste("ghcn_v_",out_prefix,"_",dates[[i]],sep="")

  write.table(data_v, file= paste(path,"/",data_name,".txt",sep=""), sep=" ")

  #write out a new shapefile (including .prj component)

  outfile<-sub(".shp","",data_name)   #Removing extension if it is present

  writeOGR(data_v,".", outfile, driver ="ESRI Shapefile")

  data_name<-paste("ghcn_s_",out_prefix,"_",dates[[i]],sep="")

  assign(data_name,data_s)

  write.table(data_s, file= paste(path,"/",data_name,".txt",sep=""), sep=" ")

  outfile<-sub(".shp","",data_name)   #Removing extension if it is present

  writeOGR(data_s,".", outfile, driver ="ESRI Shapefile")

  # end of the for loop #1

mode(results_RMSEnum)<- "numeric"         # Make it numeric first

mode(results_AICnum)<- "numeric"          # Now turn it into a data.frame...

colnames(results_table_RMSE)<-c("dates","ns","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7")

colnames(results_table_AIC)<-c("dates","ns","mod1", "mod2","mod3", "mod4", "mod5", "mod6", "mod7")

results_RMSE_kr_num <-results_RMSE_kr

mode(results_RMSE_kr_num)<- "numeric"         # Make it numeric first

results_table_RMSE_kr<-as.data.frame(results_RMSE_kr_num)

colnames(results_table_RMSE_kr)<-c("dates","ns","mod1k", "mod2k","mod3k", "mod4k", "mod5k", "mod6k", "mod7k")

write.table(results_table_RMSE_kr, file= paste(path,"/","results_GAM_Assessment_kr",out_prefix,".txt",sep=""), sep=",")

##Visualization of results##

for(i in 1:length(dates)){

  X11()

  RMSE_kr<-results_table_RMSE_kr[i,]

  RMSE_ga<-results_table_RMSE[i,]

  RMSE_kr<-RMSE_kr[,1:length(models)+2]

  RMSE_ga<-RMSE_ga[,1:length(models)+2]

  colnames(RMSE_kr)<-names(RMSE_ga)

  height<-rbind(RMSE_ga,RMSE_kr)

  rownames(height)<-c("GAM","GAM_KR")

  height<-as.matrix(height)

  barplot(height,ylim=c(0,105),ylab="RMSE in tenth deg C",beside=TRUE,

          legend.text=rownames(height),

          args.legend=list(x="topright"),

          main=paste("RMSE for date ",dates[i], sep=""))

  savePlot(paste("Barplot_results_RMSE_GAM_KR_",dates[i],out_prefix,".png", sep=""), type="png")

  }

# End of script##########