data3<-read.table(paste(path,"/","ghcn_data_TMAXy1980_2010_OR_0602012.txt",sep=""),sep=",", header=TRUE)

prop<-0.3                                                                           #Proportion of testing retained for validation   

seed_number<- 100                                                                   #Seed number for random sampling

setwd(path)

############ START OF THE SCRIPT ##################

mean_LST<- readGDAL(infile5)                 #Reading the whole raster in memory. This provides a grid for kriging

proj4string(mean_LST)<-CRS                   #Assigning coordinate information to prediction grid.

set.seed(seed_number)                        #Using a seed number allow results based on random number to be compared...

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

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

ghcn_all<-ghcn

ghcn_test2<-subset(ghcn_test,ghcn_test$ELEV_SRTM>0)

ghcn<-ghcn_test2

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

  date2<-as.POSIXlt(as.Date(date1))

  data3$date<-date2

  d<-subset(data3,year>=2000 & mflag=="0" ) #Selecting dataset 2000-2010 with good quality

  d1<-aggregate(value~station+month, data=d, mean)  #Calculate monthly mean for every station in OR

  id<-as.data.frame(unique(d1$station))     #Unique station in OR for year 2000-2010      

  #This allows to change only one name of the data.frame

  dst$TMax<-dst$TMax/10                #TMax is the average max temp for months

  #dstjan=dst[dst$month==9,]  #dst contains the monthly averages for tmax for every station over 2000-2010

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

  #########

  sta_bias=sta_tmax_from_lst-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean

  bias_xy=project(as.matrix(sta_lola),proj_str)

  # windows()

  ##########################Commented out by Benoit

  #added by Benoit 

  d<-data_s

  d$dailyTmax=(as.numeric(d$dailyTmax))/10 #stored as 1/10 degree C to allow integer storage

  names(modst)[1]<-c("id")       #modst contains the average tmax per month for every stations...

  names(dmoday)[4]<-c("lat")

  ###

  #dmoday contains the daily tmax values with TMax being the monthly station tmax mean

  # windows()

  #### Added by Benoit on 06/19

  data_s<-dmoday #put the 

  data_s$daily_delta<-daily_delta

  #data_s$y_var<-daily_delta  #y_var is the variable currently being modeled, may be better with BIAS!!

  #Model and response variable can be changed without affecting the script

  mod1<- gam(y_var~ s(lat) + s (lon) + s (ELEV_SRTM), data=data_s)

  mod6<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC1), data=data_s)

  mod7<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST)+s(LC3), data=data_s)

  mod8<- gam(y_var~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST) + s(LC1), data=data_s)

  #### Added by Benoit ends

  #########

  plot(daily_delta_rast,main="Raster Daily Delta")

  #tmax_predicted=themolst+daily_delta_rast+bias_rast #Added by Benoit, why is it -bias_rast

  plot(tmax_predicted,main="Predicted daily")

  sta_pred=lookup(tmax_predicted,data_v$lat,data_v$lon)

  #sta_pred=lookup(tmax_predicted,daily_sta_lola$lat,daily_sta_lola$lon)

  #rmse=RMSE(sta_pred,dmoday$dailyTmax)

  rmse=RMSE(sta_pred,tmax)

  #plot(sta_pred~dmoday$dailyTmax,xlab=paste("Actual daily for",datelabel),ylab="Pred daily",main=paste("RMSE=",rmse))

  X11()

  results_RMSE_f[i,3]<- "RMSE"

  results_RMSE_f[i,j+3]<- rmse_fit  #Storing RMSE for the model j

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

    results_DEV[i,3]<- "DEV"

    results_DEV[i,j+3]<- mod$deviance

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

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

    results_RMSE_f[i,3]<- "RSME"

    ##Model assessment: general diagnostic/metrics

    ##validation: using the testing data

    #This was modified on 06192012

    y_mod<- predict(mod, newdata=data_v, se.fit = TRUE) #Using the coeff to predict new values.

    #tmax_predicted=sta_LST+sta_bias-y_mod$fit

    #data_v$tmax<-(data_v$tmax)/10

    RMSE_mod <- sqrt(sum(res_mod^2)/nv)                 #RMSE FOR REGRESSION STEP 1: GAM     

    MAE_mod<- sum(abs(res_mod))/nv                     #MAE, Mean abs. Error FOR REGRESSION STEP 1: GAM   

    ME_mod<- sum(res_mod)/nv                            #ME, Mean Error or bias FOR REGRESSION STEP 1: GAM

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

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

    #Saving residuals and prediction in the dataframes: tmax predicted from GAM

    pred<-paste("pred_mod",j,sep="")

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

    data_v[[name2]]<-as.numeric(res_mod)

    #end of loop calculating RMSE

  }