###Parameters and arguments

infile1<-"ghcn_or_tmax_b_04142012_OR83M.shp"

setwd(path) 

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

prop<-0.3                                                                            #Proportion of testing retained for validation   

infile3<-"LST_dates_var_names.txt"

infile4<-"models_interpolation_04032012b.txt"

filename<-sub(".shp","",infile1)              #Removing the extension from file.

ghcn<-readOGR(".", filename)                  #reading shapefile 

ghcn = transform(ghcn,Northness = cos(ASPECT)) #Adding a variable to the dataframe

ghcn = transform(ghcn,Eastness = sin(ASPECT))  #adding variable to the dataframe.

LST_dates <-readLines(paste(path,"/",infile3, sep=""))

models <-readLines(paste(path,"/",infile4, sep=""))

#Screening for bad values

ghcn_all<-ghcn

ghcn_test<-subset(ghcn,ghcn$tmax>-150 & ghcn$tmax<400)

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

ghcn<-ghcn_test2

month_var<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09", "mm_10", "mm_11", "mm_12")

ghcn.subsets <-lapply(dates, function(d) subset(ghcn, date==d)) #this creates a list of 10 subsets data

  ####Regression part 2: GAM models

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

  mod3<- gam(tmax~ s(lat) + s (lon) + s (ELEV_SRTM) +  s (Northness)+ s (Eastness) + s(DISTOC), data=data_s)

  mod4<- gam(tmax~ s(lat) + s (lon) + s(ELEV_SRTM) + s(Northness) + s (Eastness) + s(DISTOC) + s(LST), data=data_s)

  mod5<- gam(tmax~ s(lat,lon) +s(ELEV_SRTM) + s(Northness,Eastness) + s(DISTOC) + s(LST), data=data_s)

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

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

  ####Regression part 3: Calculating and storing diagnostic measures

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

  results_AIC[i,7]<- AIC (mod5)

  results_AIC[i,8]<- AIC (mod6)

  results_AIC[i,9]<- AIC (mod7)

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

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

  results_GCV[i,7]<- mod5$gcv.ubre

  results_GCV[i,8]<- mod6$gcv.ubre

  results_GCV[i,9]<- mod7$gcv.ubre

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

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

  results_DEV[i,6]<- mod4$deviance

  results_DEV[i,7]<- mod5$deviance

  results_DEV[i,8]<- mod6$deviance

  #####VALIDATION: Prediction checking the results using the testing data########

  y_mod5<- predict(mod5, newdata=data_v, se.fit = TRUE) 

  y_mod6<- predict(mod6, newdata=data_v, se.fit = TRUE)

  y_mod7<- predict(mod7, newdata=data_v, se.fit = TRUE)

  res_mod1<- data_v$tmax - y_mod1$fit #Residuals for GMA model that resembles the ANUSPLIN interpolation

  res_mod2<- data_v$tmax - y_mod2$fit   #Residuals for GAM model that resembles the PRISM interpolation                               

  res_mod5<- data_v$tmax - y_mod5$fit

  res_mod6<- data_v$tmax - y_mod6$fit

  res_mod7<- data_v$tmax - y_mod7$fit

  RMSE_mod1 <- sqrt(sum(res_mod1^2)/nv)          

  RMSE_mod2 <- sqrt(sum(res_mod2^2)/nv)

  RMSE_mod5 <- sqrt(sum(res_mod5^2)/nv)

  RMSE_mod6 <- sqrt(sum(res_mod6^2)/nv)

  RMSE_mod7 <- sqrt(sum(res_mod7^2)/nv)

  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

  results_RMSE[i,7]<- RMSE_mod5

  results_RMSE[i,8]<- RMSE_mod6

  results_RMSE[i,9]<- RMSE_mod7

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

  assign(data_name,data_v)

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

  assign(data_name,data_s)

  }

## Plotting and saving diagnostic measures

f<-file(paste(path,"/","results_GAM_Assessment",out_prefix,"2.txt",sep=""),"w")

write(results_table_RMSE, file=f)

close(f)

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

#Selecting dates and files based on names

#   cor_LST_LC3[i]<-cor(ghcn.subsets[[i]]$LST,ghcn.subsets[[i]]$LC3)

# }