date<-strptime(sampling_dat$date[i], "%Y%m%d")   # interpolation date being processed

  month<-strftime(date, "%m")          # current month of the date being processed

  LST_month<-paste("mm_",month,sep="") # name of LST month to be matched

  #Adding layer LST to the raster stack

  pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", if not present pos=NA

  #r1[r1 < (min_val)]<-NA

  s_raster<-addLayer(s_raster,r1)            #Adding current month

  ###Regression part 1: Creating a validation dataset by creating training and testing datasets

  ind.training<-sampling[[i]]

  ns<-nrow(data_s)

  nv<-nrow(data_v)

  #min_val<-(-15)     #Screening for extreme values

  #themolst[themolst < (min_val)]<-NA

  ###########

  # STEP 2 - Weather station means across same days: Monthly mean calculation

  ###########

  #sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month..

  #proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";

  #sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon) #Extracted values of LST for the stations

  sta_tmax_from_lst<-modst$LST

  #########

  #Added by Benoit

  modst$LSTD_bias<-sta_bias  #Adding bias to data frame modst containning the monthly average for 10 years

  png(paste("LST_TMax_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i],"_",sampling_dat$run_samp[i], out_prefix,".png", sep=""))

  plot(modst$TMax,sta_tmax_from_lst,xlab="Station mo Tmax",ylab="LST mo Tmax",main=paste("LST vs TMax for",datelabel,sep=" "))

  abline(0,1)

  dev.off()

  #added by Benoit 

  #x<-ghcn.subsets[[i]]  #Holds both training and testing for instance 161 rows for Jan 1

  pos<-match("value",names(d)) #Find column with name "value"

  #names(d)[pos]<-c("dailyTmax")

  names(d)[pos]<-y_var_name

  names(x)[pos]<-y_var_name

  #names(x)[pos]<-c("dailyTmax")

  pos<-match("station",names(d)) #Find column with name "value"

  names(d)[pos]<-c("id")

  names(x)[pos]<-c("id")

  #fitbias<-Tps(bias_xy,sta_bias) #use TPS or krige

  #Adding options to use only training stations : 07/11/2012

  #bias_xy2=project(as.matrix(c(dmoday$lon,dmoday$lat),proj_str)

  if(bias_val==1){

    sta_bias<-dmoday$LSTD_bias

    bias_xy<-cbind(dmoday$x_OR83M,dmoday$y_OR83M)

  daily_sta_lola=dmoday[,c("lon","lat")] #could be same as before but why assume merge does this - assume not

  daily_sta_xy=project(as.matrix(daily_sta_lola),proj_str)

  daily_delta=dmoday$dailyTmax-dmoday$TMax

  #fitdelta<-Tps(daily_sta_xy,daily_delta) #use TPS or krige

  ########

  # check: assessment of results: validation

  ########

  #ME_fun<-function(x,y){return(mean(abs(y)))}

  #FIT ASSESSMENT

  sta_pred_data_s=lookup(tmax_predicted,data_s$lat,data_s$lon)

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

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

  #names(data_v)[pos]<-c("dailyTmax")

  tmax<-data_v$dailyTmax

  #data_v$dailyTmax<-tmax

  r2<-cor(sta_pred,tmax)^2              #R2, coef. of var

  png(paste("Predicted_tmax_versus_observed_scatterplot_",sampling_dat$date[i],"_",sampling_dat$prop[i],

            "_",sampling_dat$run_samp[i],out_prefix,".png", sep=""))

  ###BEFORE GAM prediction the data object must be transformed to SDF

  coordinates(data_v)<-coords

  coordinates(data_s)<-coords

  ns<-nrow(data_s) #This is added to because some loss of data might have happened because of the averaging...

  nv<-nrow(data_v)

  list_formulas<-vector("list",nmodels)

  list_formulas[[2]] <- as.formula("y_var ~ s(LST)", env=.GlobalEnv)

  #list_formulas[[1]] <- as.formula("y_var ~ s(ELEV_SRTM)", env=.GlobalEnv)

  #list_formulas[[2]] <- as.formula("y_var ~ s(lat,lon)", env=.GlobalEnv)