#DATE MODIFIED: 05/21/2014                                                                                 

  #1) in_path_tile: location of files, if NULL then code is run on NEX node or as stage 5

  #date_selected_results <- c("x")

  in_path_tile <- list_param$in_path_tile

  if(!is.null(in_path_tile)){

    covar_obj <- load_obj(list_param$covar_obj)

    infile_covariates <- file.path(in_path_tile,basename(covar_obj$infile_covariates))

    covar_names <- covar_obj$covar_names

  }else{ #we are on the node or running as stage 5

    infile_covariates <- list_param$covar_obj$infile_covariates #already loaded in memory

    covar_names<-list_param$covar_obj$covar_names

  }

  #if raster_obj has not been loaded in memory then we have

  #the name of the RData object for a specific tile

  if(class(raster_prediction_obj)=="character"){

    raster_prediction_obj <- load_obj(raster_prediction_obj)

  }

  method_mod_obj <- raster_prediction_obj$method_mod_obj

  validation_mod_obj <-raster_prediction_obj$validation_mod_obj

  if(interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){

    multi_timescale <- TRUE

  }

  #add precip option later...

  s_raster <- brick(infile_covariates) #stack produced for specific tile

  names(s_raster)<- covar_names               #Assigning names to the raster layers: making sure it is included in the extraction

  LC_mask <- subset(s_raster,"LC12")

  LC_mask[LC_mask==100]<- NA

  LC_mask_rec[is.na(LC_mask_rec)]<- 0

  date <- strptime(date_selected[j], "%Y%m%d")   # interpolation date being processed

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

  index <- i_dates[[j]]

  ##The path of production is not the same if input_path_tile is not NULL

  if(!is.null(in_path_tile)){

    #infile_covariates <- file.path(in_path_tile,basename(list_param$covar_obj$infile_covariates))

    pred_temp <- basename(as.character(method_mod_obj[[index]][[y_var_name]])) #list of daily prediction files with path included

    pred_temp <- file.path(in_path_tile,pred_temp)

  }else{

    pred_temp <- as.character(method_mod_obj[[index]][[y_var_name]]) #list of daily prediction files with path included

  }

  rast_pred_temp_s <- stack(pred_temp) #stack of temperature predictions from models (daily)

  rast_pred_temp <- mask(rast_pred_temp_s,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE)

  data_v<-validation_mod_obj[[index]]$data_v #testing with residuals

  #no formula if multi-timescale method

  if(multi_timescale==TRUE){

    formulas<- raster_prediction_obj$clim_method_mod_obj[[as.integer(month)]]$formulas #models ran

  }else{

    formulas<- method_mod_obj[[index]]$formulas #models ran

  }

  #if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){

  if (multi_timescale==TRUE){

    data_month <-clim_method_mod_obj[[mo]]$data_month

    plot(data_month[[y_var_month]],data_month$LST,xlab=paste("Station mo ",y_var_month,sep=""),ylab=paste("LST month ",mo," ",sep=""))

    title(paste("LST vs ", y_var_month,"for month ",mo,sep=" "))

    #This is not stored in data_s$TMax?

    #png(file.path(out_path,paste("Month_day_scatterplot_",y_var_name,"_",y_var_month,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,

    #                             out_prefix,".png", sep="")))

    #plot(data_s[[y_var_name]]~data_s[[y_var_month]],xlab=paste("Month") ,ylab=paste("Daily for",datelabel),main="across stations in VE")

    #nb_point<-paste("ns=",length(data_s[[y_var_month]]),sep="")

    #nb_point2<-paste("ns_obs=",length(data_s[[y_var_month]])-sum(is.na(data_s[[y_var_name]])),sep="")

    #nb_point3<-paste("n_month=",length(data_month[[y_var_month]]),sep="")

    #legend("topleft",legend=c(nb_point,nb_point2,nb_point3),bty="n",cex=0.8)

    #dev.off()

    plot(raster(rast_pred_temp,layer=1))

    title("Monthly ghcn station in tile for January")

  } #End of if multi_timescale=TRUE

  png(file.path(out_path,paste("Raster_prediction_levelplot_",y_var_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,

  png(file.path(out_path,paste("Raster_prediction_plot_",y_var_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",sampling_dat$run_samp,

  ## Figure 6: training and testing daily stations used

  plot(raster(rast_pred_temp,layer=1))

  title(paste("Daily stations ", datelabel,sep=""))

  nb_point1<-paste("ns_obs=",nrow(data_s),sep="")

  nb_point2<-paste("nv_obs=",nrow(data_v),sep="")

  legend("bottomright",legend=c(nb_point1,nb_point2),bty="n",cex=0.8)

      ##The path of production is not the same if input_path_tile is not NULL

    if(!is.null(in_path_tile)){

      #infile_covariates <- file.path(in_path_tile,basename(list_param$covar_obj$infile_covariates))

      bias_lf <- basename(as.character(clim_method_mod_obj[[index]]$bias)) #list of daily prediction files with path included

      bias_lf <- file.path(in_path_tile,bias_lf)

      delta_lf <- basename(unlist(method_mod_obj[[index]]$delta))

      delta_lf <- file.path(in_path,delta_lf)

    }else{

      bias_lf <- clim_method_mod_obj[[index]]$bias #list of daily prediction files with path included

      delta_lf <- method_mod_obj[[index]]$delta

    }

    bias_rast<-stack(bias_lf)

    delta_rast<-raster(delta_lf) #only one delta image!!!

  #if CAI

      ##The path of production is not the same if input_path_tile is not NULL

    if(!is.null(in_path_tile)){

      #infile_covariates <- file.path(in_path_tile,basename(list_param$covar_obj$infile_covariates))

      clim_lf <- basename(as.character(clim_method_mod_obj[[mo]]$clim)) #list of daily prediction files with path included

      clim_lf <- file.path(in_path_tile,clim_lf)

      delta_lf <- basename(unlist(method_mod_obj[[index]]$delta))

      delta_lf <- file.path(in_path,delta_lf)

    }else{

      clim_lf <- clim_method_mod_obj[[index]]$clim #list of monthly prediction files with path included

      delta_lf <- method_mod_obj[[index]]$delta

    }    

    clim_rast<-stack(clim_lf)

    delta_rast<-stack(delta_lf) #this is a stack now... delta images!!!

    names(delta_rast)<- paste(names_mod,"_delta",sep="")

    names(clim_rast) <- paste(names_mod,"_month",mo,sep="")

    #layers_names <- c(names(clim_rast),"delta")

    #rast_temp_date<-stack(clim_rast,delta_rast)

    #rast_temp_date<-mask(rast_temp_date,LC_mask,file=file.path(out_path,"test.tif"),overwrite=TRUE) #loosing names here

    #names(rast_temp_date) <-layers_names

    #plot(rast_temp_date)

    plot(clim_rast)

    #title("Climatology for month ", mo, sep="")

    dev.off()

    png(file.path(out_path,paste("delta_surface_",y_var_name,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],

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

    plot(delta_rast) 

    dev.off()

  }

  ### Figure 9: map of residuals...

  for (k in 1:length(names_mod)){

    model_name <- names_mod[k]

    #fig_name <- file.path(out_path,paste("Figure_residuals_map_",y_var_name,"_",model_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",

    #                             sampling_dat$run_samp,out_prefix,".png", sep=""))

    png(file.path(out_path,paste("Figure_residuals_map_",y_var_name,"_",model_name,"_",sampling_dat$date,"_",sampling_dat$prop,"_",

                                 sampling_dat$run_samp,out_prefix,".png", sep="")))

    res_model_name <- paste("res",model_name,sep="_")

    elev <- subset(s_raster,"elev_s")

    p1 <- levelplot(elev,scales = list(draw = FALSE), colorkey = FALSE,col.regions=rev(terrain.colors(255)))

    #add legend..

    cx <- ((data_v[[res_model_name]])^2)/10

    p2 <- spplot(data_v,res_model_name, cex=1 * cx,main=paste("Residuals for ",res_model_name," ",datelabel,sep=""))

    p3 <-p2+p1+p2 #to force legend...

    #p2

    print(p3)

  #Figure 9: histogram for all images/residuals...