#########################    Raster prediction    ####################################

############################ Interpolation of temperature for given processing region ##########################################

#This script interpolates temperature values using MODIS LST, covariates and GHCND station data.                      

#It requires the text file of stations and a shape file of the study area.           

#Note that the projection for both GHCND and study area is lonlat WGS84.       

#Options to run this program are:

#1) Multisampling: vary the porportions of hold out and use random samples for each run

#2)Constant sampling: use the same sample over the runs

#3)over dates: run over for example 365 dates without mulitsampling

#4)use seed number: use seed if random samples must be repeatable

#5)possibilty of running single and multiple time scale methods:

   # gam_daily, kriging_daily,gwr_daily,gam_CAI,gam_fusion,kriging_fusion,gwr_fusion and other options added.

#For multiple time scale methods, the interpolation is done first at the monthly time scale then delta surfaces are added.

#AUTHOR: Benoit Parmentier                                                                        

#DATE: 11/03/2013                                                                                 

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#568--     

#

# TO DO:

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

raster_prediction_fun <-function(list_param_raster_prediction){

  ##Function to predict temperature interpolation with 21 input parameters

  #9 parameters used in the data preparation stage and input in the current script

  #1)list_param_data_prep: used in earlier code for the query from the database and extraction for raster brick

  #2)infile_monthly: monthly averages with covariates for GHCND stations obtained after query

  #3)infile_daily: daily GHCND stations with covariates, obtained after query

  #4)infile_locs: vector file with station locations for the processing/study area (ESRI shapefile)

  #5)infile_covariates: raster covariate brick, tif file

  #6)covar_names: covar_names #remove at a later stage...

  #7)var: variable being interpolated-TMIN or TMAX

  #8)out_prefix

  #9)CRS_locs_WGS84

  #10)screen_data_training

  #

  #6 parameters for sampling function

  #10)seed_number

  #11)nb_sample

  #12)step

  #13)constant

  #14)prop_minmax

  #15)seed_number_month

  #16)nb_sample_month

  #17)step_month

  #18)constant_month

  #19)prop_minmax_month

  #20)dates_selected

  #

  #6 additional parameters for monthly climatology and more

  #21)list_models: model formulas in character vector

  #22)lst_avg: LST climatology name in the brick of covariate--change later

  #23)n_path

  #24)out_path

  #25)script_path: path to script

  #26)interpolation_method: c("gam_fusion","gam_CAI") #other otpions to be added later

  #27) use_clim_image

  #28) join_daily

  #29)list_models2: models' formulas as string vector for daily devation

  #30)interp_method2: intepolation method for daily devation step

  #31)num_cores: How many cores to use

  #32)max_mem: Max memory to use for predict step

  ###Loading R library and packages     

  library(gtools)                                         # loading some useful tools 

  library(mgcv)                                           # GAM package by Simon Wood

  library(sp)                                             # Spatial pacakge with class definition by Bivand et al.

  library(spdep)                               # Spatial pacakge with methods and spatial stat. by Bivand et al.

  library(rgdal)                               # GDAL wrapper for R, spatial utilities

  library(gstat)                               # Kriging and co-kriging by Pebesma et al.

  library(fields)                             # NCAR Spatial Interpolation methods such as kriging, splines

  library(raster)                              # Hijmans et al. package for raster processing

  library(rasterVis)

  library(parallel)                            # Urbanek S. and Ripley B., package for multi cores & parralel processing

  library(reshape)

  library(plotrix)

  library(maptools)

  library(gdata) #Nesssary to use cbindX

  library(automap)  #autokrige function

  library(spgwr)   #GWR method

  ### Parameters and arguments

  #PARSING INPUTS/ARGUMENTS

#   

#   names(list_param_raster_prediction)<-c("list_param_data_prep",

#                                          "seed_number","nb_sample","step","constant","prop_minmax","dates_selected",

#                                          "list_models","lst_avg","in_path","out_path","script_path",

#                                          "interpolation_method")

  #9 parameters used in the data preparation stage and input in the current script

  list_param_data_prep<-list_param_raster_prediction$list_param_data_prep

  infile_monthly<-list_param_data_prep$infile_monthly

  infile_daily<-list_param_data_prep$infile_daily

  infile_locs<-list_param_data_prep$infile_locs

  infile_covariates<-list_param_data_prep$infile_covariates #raster covariate brick, tif file

  covar_names<- list_param_data_prep$covar_names #remove at a later stage...

  var<-list_param_data_prep$var

  out_prefix<-list_param_data_prep$out_prefix

  CRS_locs_WGS84<-list_param_data_prep$CRS_locs_WGS84

  #6 parameters for sampling function

  seed_number<-list_param_raster_prediction$seed_number

  nb_sample<-list_param_raster_prediction$nb_sample

  step<-list_param_raster_prediction$step

  constant<-list_param_raster_prediction$constant

  prop_minmax<-list_param_raster_prediction$prop_minmax

  dates_selected<-list_param_raster_prediction$dates_selected

  seed_number_month <-list_param_raster_prediction$seed_number_month

  nb_sample_month <-list_param_raster_prediction$nb_sample_month

  step_month <-list_param_raster_prediction$step_month

  constant_month <-list_param_raster_prediction$constant_month

  prop_minmax_month <-list_param_raster_prediction$prop_minmax_month

  #6 additional parameters for monthly climatology and more

  list_models<-list_param_raster_prediction$list_models

  list_models2<-list_param_raster_prediction$list_models2

  interp_method2 <- list_param_raster_prediction$interp_method2

  lst_avg<-list_param_raster_prediction$lst_avg

  out_path<-list_param_raster_prediction$out_path

  script_path<-list_param_raster_prediction$script_path

  interpolation_method<-list_param_raster_prediction$interpolation_method

  screen_data_training <-list_param_raster_prediction$screen_data_training

  use_clim_image <- list_param_raster_prediction$use_clim_image # use predicted image as a base...rather than average Tmin at the station for delta

  join_daily <- list_param_raster_prediction$join_daily # join monthly and daily station before calucating delta

  #cores and memory usage options

  num_cores <- list_param_raster_prediction$num_cores

  max_mem<- as.numeric(list_param_raster_prediction$max_mem)

  rasterOptions(maxmemory=max_mem,timer=TRUE)

  setwd(out_path)

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

  #This should not be set here...? master script, modify for precip

  if (var=="TMAX"){

    y_var_name<-"dailyTmax"

    y_var_month<-"TMax"

  }

  if (var=="TMIN"){

    y_var_name<-"dailyTmin"

    y_var_month <-"TMin"

  }

  ################# CREATE LOG FILE #####################

  #create log file to keep track of details such as processing times and parameters.

  #log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")

  log_fname<-paste("R_log_raster_prediction",out_prefix, ".log",sep="")

  #sink(log_fname) #create new log file

  file.create(file.path(out_path,log_fname)) #create new log file

  time1<-proc.time()    #Start stop watch

  cat(paste("Starting script at this local Date and Time: ",as.character(Sys.time()),sep=""),

             file=log_fname,sep="\n")

  cat("Starting script process time:",file=log_fname,sep="\n",append=TRUE)

  cat(as.character(time1),file=log_fname,sep="\n",append=TRUE)    

  ############### READING INPUTS: DAILY STATION DATA AND OTHER DATASETS  #################

  ghcn<-readOGR(dsn=dirname(infile_daily),layer=sub(".shp","",basename(infile_daily)))

  CRS_interp<-proj4string(ghcn)                       #Storing projection information (ellipsoid, datum,etc.)

  stat_loc<-readOGR(dsn=dirname(infile_locs),layer=sub(".shp","",basename(infile_locs)))

  #dates2 <-readLines(file.path(in_path,dates_selected)) #dates to be predicted, now read directly from the file

  #Should clean this up, reduce the number of if

  if (dates_selected==""){

    dates<-as.character(sort(unique(ghcn$date))) #dates to be predicted 

  }

  if (dates_selected!=""){

    dates<-dates_selected #dates to be predicted 

  }

  if(class(dates_selected)=="numeric"){ #select n every  observation, may change this later.

    dates<-as.character(sort(unique(ghcn$date))) #dates to be predicted 

    dates <- dates[seq(1, length(dates), dates_selected)]

  }

  #Reading in covariate brickcan be changed...

  s_raster<-brick(infile_covariates)                   #read in the data brck

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

  #Reading monthly data

  dst<-readOGR(dsn=dirname(infile_monthly),layer=sub(".shp","",basename(infile_monthly)))

  #construct date based on input end_year !!!

  day_tmp <- rep("15",length=nrow(dst))

  year_tmp <- rep(as.character(end_year),length=nrow(dst))

  #dates_month <-do.call(paste,c(list(day_tmp,sprintf( "%02d", dst$month ),year_tmp),sep="")) #reformat integer using formatC or sprintf

  dates_month <-do.call(paste,c(list(year_tmp,sprintf( "%02d", dst$month ),day_tmp),sep="")) #reformat integer using formatC or sprintf

  dst$date <- dates_month

  ########### CREATE SAMPLING -TRAINING AND TESTING STATIONS ###########

  #dates #list of dates for prediction

  #ghcn_name<-"ghcn" #infile daily data 

  list_param_sampling<-list(seed_number,nb_sample,step,constant,prop_minmax,dates,ghcn)

  #list_param_sampling<-list(seed_number,nb_sample,step,constant,prop_minmax,dates,ghcn_name)

  names(list_param_sampling)<-c("seed_number","nb_sample","step","constant","prop_minmax","dates","ghcn")

  #run function, note that dates must be a character vector!! Daily sampling

  sampling_obj<-sampling_training_testing(list_param_sampling)

  #Now run monthly sampling

  dates_month<-as.character(sort(unique(dst$date)))

  list_param_sampling<-list(seed_number_month,nb_sample_month,step_month,constant_month,prop_minmax_month,dates_month,dst)

  #list_param_sampling<-list(seed_number,nb_sample,step,constant,prop_minmax,dates,ghcn_name)

  names(list_param_sampling)<-c("seed_number","nb_sample","step","constant","prop_minmax","dates","ghcn")

  sampling_month_obj <- sampling_training_testing(list_param_sampling)

  ########### PREDICT FOR MONTHLY SCALE  ##################

  #First predict at the monthly time scale: climatology

  cat("Predictions at monthly scale:",file=log_fname,sep="\n", append=TRUE)

  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),

      file=log_fname,sep="\n")

  t1<-proc.time()

  j=6 #12

  #browser() #Missing out_path for gam_fusion list param!!!

  #if (interpolation_method=="gam_fusion"){

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

    list_param_runClim_KGFusion<-list(j,s_raster,covar_names,lst_avg,list_models,dst,sampling_month_obj,var,y_var_name, out_prefix,out_path)

    names(list_param_runClim_KGFusion)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","sampling_month_obj","var","y_var_name","out_prefix","out_path")

    #debug(runClim_KGFusion)

    #test<-runClim_KGFusion(1,list_param=list_param_runClim_KGFusion)

    clim_method_mod_obj<-mclapply(1:length(sampling_month_obj$ghcn_data), list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))

    #Use function to extract list

    list_tmp<-vector("list",length(clim_method_mod_obj))

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

      tmp<-clim_method_mod_obj[[i]]$clim

      list_tmp[[i]]<-tmp

    }

    clim_yearlist<-list_tmp

  }

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

    list_param_runClim_KGCAI<-list(j,s_raster,covar_names,lst_avg,list_models,dst,sampling_month_obj,var,y_var_name, out_prefix,out_path)

    names(list_param_runClim_KGCAI)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","sampling_month_obj","var","y_var_name","out_prefix","out_path")

    clim_method_mod_obj<-mclapply(1:length(sampling_month_obj$ghcn_data), list_param=list_param_runClim_KGCAI, runClim_KGCAI,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    #test<-runClim_KGCAI(1,list_param=list_param_runClim_KGCAI)

    save(clim_method_mod_obj,file= file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep="")))

    list_tmp<-vector("list",length(clim_method_mod_obj))

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

      tmp<-clim_method_mod_obj[[i]]$clim

      list_tmp[[i]]<-tmp

    }

    clim_yearlist<-list_tmp

  }

  t2<-proc.time()-t1

  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)

  #Getting rid of raster temp files

  removeTmpFiles(h=0)

  #quit()

  ################## PREDICT AT DAILY TIME SCALE #################

  #Predict at daily time scale from single time scale or multiple time scale methods: 2 methods availabe now

  #put together list of clim models per month...

  #rast_clim_yearlist<-list_tmp

  #Second predict at the daily time scale: delta

  #method_mod_obj<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement

  cat("Predictions at the daily scale:",file=log_fname,sep="\n", append=TRUE)

  t1<-proc.time()

  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),

      file=log_fname,sep="\n")

  #TODO : Same call for all functions!!! Replace by one "if" for all multi time scale methods...

  #The methods could be defined earlier as constant??

  #Create data.frame combining sampling at daily and monthly time scales:

  daily_dev_sampling_dat <- combine_sampling_daily_monthly_for_daily_deviation_pred(sampling_obj,sampling_month_obj)

  #use_clim_image<- TRUE

  #use_clim_image<-FALSE

  #join_daily <- FALSE

  #join_daily <- TRUE

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

    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn

    i<-1

    list_param_run_prediction_daily_deviation <-list(i,clim_yearlist,daily_dev_sampling_dat,sampling_month_obj,sampling_obj,dst,list_models2,interp_method2,

                                                     s_raster,use_clim_image,join_daily,var,y_var_name, interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","daily_dev_sampling_dat","sampling_month_obj","sampling_obj","dst","list_models2","interp_method2",

                                                        "s_raster","use_clim_image","join_daily","var","y_var_name","interpolation_method","out_prefix","out_path")

    #method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data),list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

    #debug(run_prediction_daily_deviation)

    #test <- run_prediction_daily_deviation(1,list_param=list_param_run_prediction_daily_deviation) #This is the end bracket from mclapply(...) statement

    #test <- mclapply(1:9,list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

    method_mod_obj<-mclapply(1:nrow(daily_dev_sampling_dat),list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

  }

  #TODO : Same call for all functions!!! Replace by one "if" for all daily single time scale methods...

  if (interpolation_method=="gam_daily"){

    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn

    i<-1

    list_param_run_prediction_gam_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,screen_data_training,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_gam_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","screen_data_training","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")

    #test <- runGAM_day_fun(1,list_param_run_prediction_gam_daily)

    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data),list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    #method_mod_obj<-mclapply(1:22,list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

  }

  if (interpolation_method=="kriging_daily"){

    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn

    i<-1

    list_param_run_prediction_kriging_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_kriging_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")

    #test <- runKriging_day_fun(1,list_param_run_prediction_kriging_daily)

    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data),list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

  }

  if (interpolation_method=="gwr_daily"){

    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn

    i<-1

    list_param_run_prediction_gwr_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_gwr_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")

    #test <- run_interp_day_fun(1,list_param_run_prediction_gwr_daily)

    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data),list_param=list_param_run_prediction_gwr_daily,run_interp_day_fun,mc.preschedule=FALSE,mc.cores = num_cores) #This is the end bracket from mclapply(...) statement

    #method_mod_obj<-mclapply(1:22,list_param=list_param_run_prediction_gwr_daily,run_interp_day_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

  }

  t2<-proc.time()-t1

  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)

  #browser()

  ############### NOW RUN VALIDATION #########################

  #SIMPLIFY THIS PART: one call

  cat("Validation step:",file=log_fname,sep="\n", append=TRUE)

  t1<-proc.time()

  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),

      file=log_fname,sep="\n")

  if (interpolation_method=="gam_daily" | interpolation_method=="kriging_daily" | interpolation_method=="gwr_daily"){

    multi_time_scale <- FALSE

    list_data_v <- extract_list_from_list_obj(method_mod_obj,"data_v")

    list_data_s <- extract_list_from_list_obj(method_mod_obj,"data_s")

    rast_day_yearlist <- extract_list_from_list_obj(method_mod_obj,y_var_name) #list_tmp #list of predicted images over full year...

    list_sampling_dat <- extract_list_from_list_obj(method_mod_obj,"sampling_dat")

    list_param_validation<-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_name, multi_time_scale,out_prefix, out_path)

    names(list_param_validation)<-c("list_index","rast_day_year_list",

                                  "list_data_v","list_data_s","list_sampling_dat","y_ref","multi_time_scale","out_prefix", "out_path") #same names for any method

    #debug(calculate_accuracy_metrics)

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation)

    validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = num_cores) 

    save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep="")))

    t2<-proc.time()-t1

    cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)

  }

  ### Run monthly validation if multi-time scale methods and add information to daily...

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

    multi_time_scale <- TRUE

    i<-1

    ## daily time scale

    list_data_v <- extract_list_from_list_obj(method_mod_obj,"data_v")

    list_data_s <- extract_list_from_list_obj(method_mod_obj,"data_s")

    rast_day_yearlist <- extract_list_from_list_obj(method_mod_obj,y_var_name) #list_tmp #list of predicted images over full year...

    list_sampling_dat <- extract_list_from_list_obj(method_mod_obj,"daily_dev_sampling_dat")

    list_param_validation<-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_name, multi_time_scale,out_prefix, out_path)

    names(list_param_validation)<-c("list_index","rast_day_year_list",

                                    "list_data_v","list_data_s","list_sampling_dat","y_ref","multi_time_scale","out_prefix", "out_path") #same names for any method

    #debug(calculate_accuracy_metrics)

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation)

    validation_mod_obj <-mclapply(1:length(method_mod_obj), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = num_cores) 

    save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep="")))

    ### monthly time scale

    list_data_v <- extract_list_from_list_obj(clim_method_mod_obj,"data_month_v") #extract monthly testing/validation dataset

    list_data_s <- extract_list_from_list_obj(clim_method_mod_obj,"data_month") #extract monthly training/fitting dataset

    rast_day_yearlist <- extract_list_from_list_obj(clim_method_mod_obj,"clim") #list_tmp #list of predicted images over full year at monthly time scale

    list_sampling_dat <- extract_list_from_list_obj(clim_method_mod_obj,"sampling_month_dat")

    #list_param_validation_month <-list(i,clim_yearlist,clim_method_mod_obj,y_var_name, multi_time_scale ,out_prefix, out_path)

    #names(list_param_validation_month)<-c("list_index","rast_day_year_list","method_mod_obj","y_var_name","multi_time_scale","out_prefix", "out_path") #same names for any method

    list_param_validation_month <-list(i,rast_day_yearlist,list_data_v,list_data_s,list_sampling_dat,y_var_month, multi_time_scale,out_prefix, out_path)

    names(list_param_validation_month)<-c("list_index","rast_day_year_list",

                                    "list_data_v","list_data_s","list_sampling_dat","y_ref","multi_time_scale","out_prefix", "out_path") #same names for any method

    #debug(calculate_accuracy_metrics)    

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation_month)

    validation_mod_month_obj <- mclapply(1:length(clim_method_mod_obj), list_param=list_param_validation_month, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = num_cores) 

    #test_val<-calculate_accuracy_metrics(1,list_param_validation)

    save(validation_mod_month_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_month_obj_",y_var_name,out_prefix,".RData",sep="")))

    ##Create data.frame with validation and fit metrics for a full year/full numbe of runs

    tb_month_diagnostic_v <- extract_from_list_obj(validation_mod_month_obj,"metrics_v") 

    #tb_diagnostic_v contains accuracy metrics for models sample and proportion for every run...if full year then 365 rows maximum

    rownames(tb_month_diagnostic_v)<-NULL #remove row names

    tb_month_diagnostic_v$method_interp <- interpolation_method

    tb_month_diagnostic_s<-extract_from_list_obj(validation_mod_month_obj,"metrics_s")

    rownames(tb_month_diagnostic_s)<-NULL #remove row names

    tb_month_diagnostic_s$method_interp <- interpolation_method #add type of interpolation...out_prefix too??

  }

  #Cleaning raster temp files

  removeTmpFiles(h=0)

  #################### ASSESSMENT OF PREDICTIONS: PLOTS OF ACCURACY METRICS ###########

  ##Create data.frame with validation and fit metrics for a full year/full numbe of runs

  tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v") 

  #tb_diagnostic_v contains accuracy metrics for models sample and proportion for every run...if full year then 365 rows maximum

  rownames(tb_diagnostic_v)<-NULL #remove row names

  tb_diagnostic_v$method_interp <- interpolation_method

  tb_diagnostic_s<-extract_from_list_obj(validation_mod_obj,"metrics_s")

  rownames(tb_diagnostic_s)<-NULL #remove row names

  tb_diagnostic_s$method_interp <- interpolation_method #add type of interpolation...out_prefix too??

  #Call functions to create plots of metrics for validation dataset

  metric_names<-c("rmse","mae","me","r","m50")

  summary_metrics_v<- boxplot_from_tb(tb_diagnostic_v,metric_names,out_prefix,out_path) #if adding for fit need to change outprefix

  names(summary_metrics_v)<-c("avg","median")

  summary_month_metrics_v<- boxplot_month_from_tb(tb_diagnostic_v,metric_names,out_prefix,out_path)

  #################### CLOSE LOG FILE  ####################

  #close log_file connection and add meta data

  cat("Finished script process time:",file=log_fname,sep="\n", append=TRUE)

  time2<-proc.time()-time1

  cat(as.character(time2),file=log_fname,sep="\n", append=TRUE)

  #later on add all the parameters used in the script...

  cat(paste("Finished script at this local Date and Time: ",as.character(Sys.time()),sep=""),

             file=log_fname,sep="\n", append=TRUE)

  cat("End of script",file=log_fname,sep="\n", append=TRUE)

  #close(log_fname)

  ################### PREPARE RETURN OBJECT ###############

  #Will add more information to be returned

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

    raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,validation_mod_month_obj, tb_diagnostic_v,

                                tb_diagnostic_s,tb_month_diagnostic_v,tb_month_diagnostic_s,summary_metrics_v,summary_month_metrics_v)

    names(raster_prediction_obj)<-c("clim_method_mod_obj","method_mod_obj","validation_mod_obj","validation_mod_month_obj","tb_diagnostic_v",

                                    "tb_diagnostic_s","tb_month_diagnostic_v","tb_month_diagnostic_s","summary_metrics_v","summary_month_metrics_v")  

    save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))

  }

  #use %in% instead of "|" operator

  if (interpolation_method=="gam_daily" | interpolation_method=="kriging_daily" | interpolation_method=="gwr_daily"){

    raster_prediction_obj<-list(method_mod_obj,validation_mod_obj,tb_diagnostic_v,

                                tb_diagnostic_s,summary_metrics_v,summary_month_metrics_v)

    names(raster_prediction_obj)<-c("method_mod_obj","validation_mod_obj","tb_diagnostic_v",

                                    "tb_diagnostic_s","summary_metrics_v","summary_month_metrics_v")  

    save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))

  }

  print("stage 4 DONE")

  return(raster_prediction_obj)

}

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

######################## END OF SCRIPT/FUNCTION #####################

##################  Functions for use in the raster prediction stage   #######################################

############################ Interpolation in a given tile/region ##########################################

#This script contains 5 functions used in the interpolation of temperature in the specfied study/processing area:                             

# 1)predict_raster_model<-function(in_models,r_stack,out_filename)                                                             

# 2)fit_models<-function(list_formulas,data_training)           

# 3)runClim_KGCAI<-function(j,list_param) : function that peforms GAM CAI method

# 4)runClim_KGFusion<-function(j,list_param) function for monthly step (climatology) in the fusion method

# 5)runGAMFusion <- function(i,list_param) : daily step for fusion method, perform daily prediction

#

#AUTHOR: Benoit Parmentier                                                                       

#DATE: 10/03/2013                                                                                 

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363--   

##Comments and TODO:

#This script is meant to be for general processing tile by tile or region by region.

# Note that the functions are called from GAM_fusion_analysis_raster_prediction_mutlisampling.R.

# This will be expanded to other methods.

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

predict_raster_model<-function(in_models,r_stack,out_filename){

  #This functions performs predictions on a raster grid given input models.

  #Arguments: list of fitted models, raster stack of covariates

  #Output: spatial grid data frame of the subset of tiles

  list_rast_pred<-vector("list",length(in_models))

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

    mod <-in_models[[i]] #accessing GAM model ojbect "j"

    raster_name<-out_filename[[i]]

    if (inherits(mod,"gam")) {           #change to c("gam","autoKrige")

      raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE,block.size=1000) #Using the coeff to predict new values.

      raster_pred<- predict(object=r_stack,model=mod,na.rm=FALSE) #Using the coeff to predict new values.

      names(raster_pred)<-"y_pred"  

      writeRaster(raster_pred, filename=raster_name,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

      #print(paste("Interpolation:","mod", j ,sep=" "))

      list_rast_pred[[i]]<-raster_name

    }

  }

  if (inherits(mod,"try-error")) {

    print(paste("no gam model fitted:",mod[1],sep=" ")) #change message for any model type...

  }

  return(list_rast_pred)

}

fit_models<-function(list_formulas,data_training){

  #This functions several models and returns model objects.

  #Arguments: - list of formulas for GAM models

  #           - fitting data in a data.frame or SpatialPointDataFrame

  #Output: list of model objects 

  list_fitted_models<-vector("list",length(list_formulas))

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

    formula<-list_formulas[[k]]

    mod<- try(gam(formula, data=data_training)) #change to any model!!

    #mod<- try(autoKrige(formula, input_data=data_s,new_data=s_sgdf,data_variogram=data_s))

    model_name<-paste("mod",k,sep="")

    assign(model_name,mod) 

    list_fitted_models[[k]]<-mod

  }

  return(list_fitted_models) 

}

#Function to glue all methods together...still need to separate fit and training for gwr and kriging, ok for now

interpolate_area_fun <- function(method_interp,list_models,s_raster,list_out_filename,data_df){

  ##Function to fit and predict an interpolation surface

  ##Author: Benoit Parmentier

  ##Function depends on other functions!!!

  #inpputs:

  #method_interp: interpolation method with value "gam","gwr","kriging"

  #list_models: models to fit and predict as string (i.e.vector char)

  #s_raster: stack with covariate variables, must match in name the data.frame input

  #data_df: spatial point data.frame with covariates, must be projected match names of covariates

  #list_out_filename: list of char containing output names for models

  #Conver to formula object

  list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

  cname<-paste("mod",1:length(list_formulas),sep="") #change to more meaningful name?

  names(list_out_filename)<-cname  

  ##Now carry out prediction

  if(method_interp=="gam"){    

    #First fitting

    mod_list<-fit_models(list_formulas,data_df) #only gam at this stage

    names(mod_list)<-cname

    #if raster provided then predict surface

    if(!is.null(s_raster)){

      #Second predict values for raster image...by providing fitted model list, raster brick and list of output file names

      rast_pred_list<-predict_raster_model(mod_list,s_raster,list_out_filename)

      names(rast_pred_list)<-cname

    }

  }

  if(method_interp%in%c("gwr","kriging")){

    #Call funciton to fit and predict gwr and/or kriging

    #month_prediction_obj<-predict_auto_krige_raster_model(list_formulas,s_raster,data_month,list_out_filename)

    rast_prediction_obj<-predict_autokrige_gwr_raster_model(method_interp,list_formulas,s_raster,data_df,list_out_filename)

    mod_list <-rast_prediction_obj$list_fitted_models

    rast_pred_list <-rast_prediction_obj$list_rast_pred

    names(rast_pred_list)<-cname

  }

  #Now prepare to return object

  interp_area_obj <-list(mod_list,list_formulas,rast_pred_list)

  names(interp_area_obj) <- c("mod_list","list_formulas","rast_pred_list")

  return(interp_area_obj)

} 

####

#TODO:

#Add log file and calculate time and sizes for processes-outputs

#Can combine runClim_KGFusion and runClim_KGCAI

runClim_KGCAI <-function(j,list_param){

  #Make this a function with multiple argument that can be used by mcmapply??

  #Arguments: 

  #1)list_index: j 

  #2)covar_rast: covariates raster images used in the modeling

  #3)covar_names: names of input variables 

  #4)lst_avg: list of LST climatogy names, may be removed later on

  #5)list_models: list input models for bias calculation

  #6)dst: data at the monthly time scale

  #7)var: TMAX or TMIN, variable being interpolated

  #8)y_var_name: output name, not used at this stage

  #9)out_prefix

  #10) out_path

  #The output is a list of four shapefile names produced by the function:

  #1) clim: list of output names for raster climatogies 

  #2) data_month: monthly training data for bias surface modeling

  #3) mod: list of model objects fitted 

  #4) formulas: list of formulas used in bias modeling

  ### PARSING INPUT ARGUMENTS

  #list_param_runGAMFusion<-list(i,clim_yearlist,sampling_obj,var,y_var_name, out_prefix)

  index<-list_param$j

  s_raster<-list_param$covar_rast

  covar_names<-list_param$covar_names

  lst_avg<-list_param$lst_avg

  list_models<-list_param$list_models

  dst<-list_param$dst #monthly station dataset

  var<-list_param$var

  y_var_name<-list_param$y_var_name

  out_prefix<-list_param$out_prefix

  out_path<-list_param$out_path

  #inserted #

  sampling_month_obj<-list_param$sampling_month_obj

  ghcn.month.subsets<-sampling_month_obj$ghcn_data

  sampling_month_dat <- sampling_month_obj$sampling_dat

  sampling_month_index <- sampling_month_obj$sampling_index

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

  #prop_month<-0 #proportion retained for validation...

  #run_samp<-1 #sample number, can be introduced later...

  prop_month <- sampling_month_dat$prop[j] #proportion retained for validation...

  run_samp <- sampling_month_dat$run_samp[j] #sample number if multisampling...will need create mulitple prediction at daily!!! could be complicated

                                       #possibility is to average per proportion !!!

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

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

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

  LST_name <-LST_month  

  #### STEP 2: PREPARE DATA

  #change here...use training data...

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

  #LST_name <-lst_avg[j] # name of LST month to be matched

  #data_month$LST<-data_month[[LST_name]]

  dataset_month <-ghcn.month.subsets[[j]]

  mod_LST <- ghcn.month.subsets[[j]][,match(LST_month, names(ghcn.month.subsets[[j]]))]  #Match interpolation date and monthly LST average

  dataset_month$LST <- as.data.frame(mod_LST)[,1] #Add the variable LST to the dataset

  #change here...

  dst$LST<-dst[[LST_month]] #Add the variable LST to the monthly dataset

  proj_str<-proj4string(dst) #get the local projection information from monthly data

  #TMax to model..., add precip later

  if (var=="TMAX"){   

    dataset_month$y_var<-dataset_month$TMax #Adding TMax as the variable modeled

  }

  if (var=="TMIN"){   

    dataset_month$y_var<-dataset_month$TMin #Adding TMin as the variable modeled

  }

  ind.training <- sampling_month_index[[j]]

  ind.testing  <- setdiff(1:nrow(dataset_month), ind.training)

  data_month_s <- dataset_month[ind.training, ]   #Training dataset currently used in the modeling

  data_month_v <- dataset_month[ind.testing, ]    #Testing/validation dataset using input sampling

  data_month <- data_month_s #training data for  monthhly predictions...

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

  #mo<-as.integer(strftime(date_proc, "%m"))          # current month of the date being processed

  #day<-as.integer(strftime(date_proc, "%d"))

  #year<-as.integer(strftime(date_proc, "%Y"))

  ## end of pasted

  #end of insert...

  #Fit gam models using data and list of formulas

  list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

  cname<-paste("mod",1:length(list_formulas),sep="") #change to more meaningful name?

  #mod_list<-fit_models(list_formulas,data_month) #only gam at this stage

  #cname<-paste("mod",1:length(mod_list),sep="") #change to more meaningful name?

  #Adding layer LST to the raster stack  

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

  s_raster<-dropLayer(s_raster,pos)      # If it exists drop layer

  LST<-subset(s_raster,LST_name)

  names(LST)<-"LST"

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

  #Now generate file names for the predictions...

  list_out_filename<-vector("list",length(list_formulas))

  names(list_out_filename)<-cname  

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

    #j indicate which month is predicted

    data_name<-paste(var,"_clim_month_",as.integer(month_no),"_",cname[k],"_",prop_month,

                     "_",run_samp,sep="")

    raster_name<-file.path(out_path,paste("CAI_",data_name,out_prefix,".tif", sep=""))

    list_out_filename[[k]]<-raster_name

  }

  ## Select the relevant method...

  if (interpolation_method=="gam_CAI"){

    #First fitting

    mod_list<-fit_models(list_formulas,data_month) #only gam at this stage

    names(mod_list)<-cname

    #Second predict values for raster image...by providing fitted model list, raster brick and list of output file names

    #now predict values for raster image...

    rast_clim_list<-predict_raster_model(mod_list,s_raster,list_out_filename)

    names(rast_clim_list)<-cname

    #Some models will not be predicted because of the lack of training data...remove empty string from list of models

  }

  if (interpolation_method %in% c("kriging_CAI","gwr_CAI")){

    if(interpolation_method=="kriging_CAI"){

      method_interp <- "kriging"

    }else{

      method_interp <- "gwr"

    }

    #Call function to fit and predict gwr and/or kriging

    #month_prediction_obj<-predict_auto_krige_raster_model(list_formulas,s_raster,data_month,list_out_filename)

    month_prediction_obj<-predict_autokrige_gwr_raster_model(method_interp,list_formulas,s_raster,data_month,list_out_filename)

    mod_list <-month_prediction_obj$list_fitted_models

    rast_clim_list <-month_prediction_obj$list_rast_pred

    names(rast_clim_list)<-cname

  }

  rast_clim_list<-rast_clim_list[!sapply(rast_clim_list,is.null)] #remove NULL elements in list

  #Adding Kriging for Climatology options 

  clim_xy<-coordinates(data_month)

  fitclim<-Krig(clim_xy,data_month$y_var,theta=1e5) #use TPS or krige 

  #fitclim<-Krig(clim_xy,data_month$TMax,theta=1e5) #use TPS or krige 

  mod_krtmp1<-fitclim

  model_name<-"mod_kr"

  clim_rast<-interpolate(LST,fitclim) #interpolation using function from raster package

  #Write out modeled layers

  data_name<-paste(var,"_clim_month_",as.integer(month_no),"_",model_name,"_",prop_month,

                   "_",run_samp,sep="")

  raster_name_clim<-file.path(out_path,paste("CAI_",data_name,out_prefix,".tif", sep=""))

  writeRaster(clim_rast, filename=raster_name_clim,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

  #Adding to current objects

  mod_list[[model_name]]<-mod_krtmp1

  #rast_bias_list[[model_name]]<-raster_name_bias

  rast_clim_list[[model_name]]<-raster_name_clim

  #Prepare object to return

  clim_obj<-list(rast_clim_list,data_month,data_month_v,sampling_month_dat[j,],mod_list,list_formulas)

  names(clim_obj)<-c("clim","data_month","data_month_v","sampling_month_dat","mod","formulas")

  save(clim_obj,file= file.path(out_path,paste("clim_obj_CAI_month_",as.integer(month_no),"_",var,"_",prop_month,

                                               "_",run_samp,"_",out_prefix,".RData",sep="")))

  return(clim_obj) 

}

#

runClim_KGFusion<-function(j,list_param){

  #Make this a function with multiple argument that can be used by mcmapply??

  #Arguments: 

  #1)list_index: j 

  #2)covar_rast: covariates raster images used in the modeling

  #3)covar_names: names of input variables 

  #4)lst_avg: list of LST climatogy names, may be removed later on

  #5)list_models: list input models for bias calculation

  #6)dst: data at the monthly time scale

  #7)var: TMAX or TMIN, variable being interpolated

  #8)y_var_name: output name, not used at this stage

  #9)out_prefix

  #

  #The output is a list of four shapefile names produced by the function:

  #1) clim: list of output names for raster climatogies 

  #2) data_month: monthly training data for bias surface modeling

  #3) mod: list of model objects fitted 

  #4) formulas: list of formulas used in bias modeling

  ### PARSING INPUT ARGUMENTS

  #list_param_runGAMFusion<-list(i,clim_yearlist,sampling_obj,var,y_var_name, out_prefix)

  index<-list_param$j

  s_raster<-list_param$covar_rast

  covar_names<-list_param$covar_names

  lst_avg<-list_param$lst_avg

  list_models<-list_param$list_models

  dst<-list_param$dst #monthly station dataset

  var<-list_param$var

  y_var_name<-list_param$y_var_name

  out_prefix<-list_param$out_prefix

  out_path<-list_param$out_path

  #inserted #

  sampling_month_obj<-list_param$sampling_month_obj

  ghcn.month.subsets<-sampling_month_obj$ghcn_data

  sampling_month_dat <- sampling_month_obj$sampling_dat

  sampling_month_index <- sampling_month_obj$sampling_index

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

  #prop_month<-0 #proportion retained for validation...

  #run_samp<-1 #sample number, can be introduced later...

  prop_month <- sampling_month_dat$prop[j] #proportion retained for validation...

  run_samp <- sampling_month_dat$run_samp[j] #sample number if multisampling...

  #will need create mulitple prediction at daily!!! could be complicated

  #possibility is to average per proportion !!!

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

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

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

  LST_name <-LST_month  

  #### STEP 2: PREPARE DATA

  #change here...use training data...

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

  #LST_name <-lst_avg[j] # name of LST month to be matched

  #data_month$LST<-data_month[[LST_name]]

  dataset_month <-ghcn.month.subsets[[j]]

  mod_LST <- ghcn.month.subsets[[j]][,match(LST_month, names(ghcn.month.subsets[[j]]))]  #Match interpolation date and monthly LST average

  dataset_month$LST <- as.data.frame(mod_LST)[,1] #Add the variable LST to the dataset

  #change here...

  dst$LST<-dst[[LST_month]] #Add the variable LST to the monthly dataset

  proj_str<-proj4string(dst) #get the local projection information from monthly data

  ind.training <- sampling_month_index[[j]]

  ind.testing  <- setdiff(1:nrow(dataset_month), ind.training)

  data_month_s <- dataset_month[ind.training, ]   #Training dataset currently used in the modeling

  data_month_v <- dataset_month[ind.testing, ]    #Testing/validation dataset using input sampling

  data_month <- data_month_s #training data for  monthhly predictions...

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

  #mo<-as.integer(strftime(date_proc, "%m"))          # current month of the date being processed

  #day<-as.integer(strftime(date_proc, "%d"))

  #year<-as.integer(strftime(date_proc, "%Y"))

  ## end of pasted

  #end of insert...09/04

  #### STEP 2: PREPARE DATA

  #data_month<-dst[dst$month==j,] #Subsetting dataset for the relevant month of the date being processed

  #LST_name<-lst_avg[j] # name of LST month to be matched

  #data_month$LST<-data_month[[LST_name]]

  #Adding layer LST to the raster stack  

  covar_rast<-s_raster

  #names(s_raster)<-covar_names

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

  s_raster<-dropLayer(s_raster,pos)      # If it exists drop layer

  LST<-subset(s_raster,LST_name)

  names(LST)<-"LST"

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

  #LST bias to model...

  if (var=="TMAX"){

    data_month$LSTD_bias<-data_month$LST-data_month$TMax

    data_month$y_var<-data_month$LSTD_bias #Adding bias as the variable modeled

  }

  if (var=="TMIN"){

    data_month$LSTD_bias<-data_month$LST-data_month$TMin

    data_month$y_var<-data_month$LSTD_bias #Adding bias as the variable modeled

  }

  #If CAI model then...

  #TMax to model..., add precip later

  #if (var=="TMAX"){   

  #  dataset_month$y_var<-dataset_month$TMax #Adding TMax as the variable modeled

  #}

  #if (var=="TMIN"){   

  #  dataset_month$y_var<-dataset_month$TMin #Adding TMin as the variable modeled

  #}

  #### STEP3:  NOW FIT AND PREDICT  MODEL

  list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

  cname<-paste("mod",1:length(list_formulas),sep="") #change to more meaningful name?

  #Now generate file names for the predictions...

  list_out_filename<-vector("list",length(list_formulas))

  names(list_out_filename)<-cname  

  ##Change name...

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

    #j indicate which month is predicted, var indicates TMIN or TMAX

    data_name<-paste(var,"_bias_LST_month_",as.integer(month_no),"_",cname[k],"_",prop_month,

                     "_",run_samp,sep="")

    raster_name<-file.path(out_path,paste("fusion_",interpolation_method,"_",data_name,out_prefix,".tif", sep=""))

    list_out_filename[[k]]<-raster_name

  }

  #for (k in 1:length(list_out_filename)){

  #  #j indicate which month is predicted

  #  data_name<-paste(var,"_clim_month_",as.integer(month_no),"_",cname[k],"_",prop_month,

  #                   "_",run_samp,sep="")

  #  raster_name<-file.path(out_path,paste("CAI_",data_name,out_prefix,".tif", sep=""))

  #  list_out_filename[[k]]<-raster_name

  #}

  ## Select the relevant method...

  if (interpolation_method=="gam_fusion"){

    #First fitting

    mod_list<-fit_models(list_formulas,data_month) #only gam at this stage

    names(mod_list)<-cname

    #Second predict values for raster image...by providing fitted model list, raster brick and list of output file names

    rast_bias_list<-predict_raster_model(mod_list,s_raster,list_out_filename)

    names(rast_bias_list)<-cname

  }

  if (interpolation_method %in% c("kriging_fusion","gwr_fusion")){

    if(interpolation_method=="kriging_fusion"){

      method_interp <- "kriging"

    }else{

      method_interp <- "gwr"

    }

    #Call funciton to fit and predict gwr and/or kriging

    #month_prediction_obj<-predict_auto_krige_raster_model(list_formulas,s_raster,data_month,list_out_filename)

    month_prediction_obj<-predict_autokrige_gwr_raster_model(method_interp,list_formulas,s_raster,data_month,list_out_filename)

    mod_list <-month_prediction_obj$list_fitted_models

    rast_bias_list <-month_prediction_obj$list_rast_pred

    names(rast_bias_list)<-cname

  }

  #Some modles will not be predicted...remove them

  rast_bias_list<-rast_bias_list[!sapply(rast_bias_list,is.null)] #remove NULL elements in list

  mod_rast<-stack(rast_bias_list)  #stack of bias raster images from models

  rast_clim_list<-vector("list",nlayers(mod_rast))

  names(rast_clim_list)<-names(rast_bias_list)

  for (k in 1:nlayers(mod_rast)){

    clim_fus_rast<-LST-subset(mod_rast,k)

    data_name<-paste(var,"_clim_LST_month_",as.integer(month_no),"_",names(rast_clim_list)[k],"_",prop_month,

                     "_",run_samp,sep="")

    raster_name<-file.path(out_path,paste("fusion_",interpolation_method,"_",data_name,out_prefix,".tif", sep=""))

    rast_clim_list[[k]]<-raster_name

    writeRaster(clim_fus_rast, filename=raster_name,overwrite=TRUE)  #Wri

  }

  #### STEP 4:Adding Kriging for Climatology options

  bias_xy<-coordinates(data_month)

  #fitbias<-Krig(bias_xy,data_month$LSTD_bias,theta=1e5) #use TPS or krige 

  fitbias<-try(Krig(bias_xy,data_month$LSTD_bias,theta=1e5)) #use TPS or krige 

  model_name<-"mod_kr"

  if (inherits(fitbias,"Krig")){

    #Saving kriged surface in raster images

    bias_rast<-bias_rast<-interpolate(LST,fitbias) #interpolation using function from raster package

    data_name<-paste(var,"_bias_LST_month_",as.integer(month_no),"_",model_name,"_",prop_month,

                     "_",run_samp,sep="")

    raster_name_bias<-file.path(out_path,paste("fusion_",interpolation_method,"_",data_name,out_prefix,".tif", sep=""))

    writeRaster(bias_rast, filename=raster_name_bias,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

    #now climatology layer

    clim_rast<-LST-bias_rast

    data_name<-paste(var,"_clim_LST_month_",as.integer(month_no),"_",model_name,"_",prop_month,

                     "_",run_samp,sep="")

    raster_name_clim<-file.path(out_path,paste("fusion_",interpolation_method,"_",data_name,out_prefix,".tif", sep=""))

    writeRaster(clim_rast, filename=raster_name_clim,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

    #Adding to current objects

    mod_list[[model_name]]<-fitbias

    rast_bias_list[[model_name]]<-raster_name_bias

    rast_clim_list[[model_name]]<-raster_name_clim

  }

  if (inherits(fitbias,"try-error")){

    #NEED TO DEAL WITH THIS!!!

    #Adding to current objects

    mod_list[[model_name]]<-NULL

    rast_bias_list[[model_name]]<-NULL

    rast_clim_list[[model_name]]<-NULL

  }

  #### STEP 5: Prepare object and return

  #Prepare object to return

  clim_obj<-list(rast_bias_list,rast_clim_list,data_month,data_month_v,sampling_month_dat[j,],mod_list,list_formulas)

  names(clim_obj)<-c("bias","clim","data_month","data_month_v","sampling_month_dat","mod","formulas")

  save(clim_obj,file= file.path(out_path,paste("clim_obj_fusion_month_",as.integer(month_no),"_",var,"_",prop_month,

                                               "_",run_samp,"_",out_prefix,".RData",sep=""))) 

  return(clim_obj)

}

## Run function for kriging...?

#runGAMFusion <- function(i,list_param) {            # loop over dates

run_prediction_daily_deviation <- function(i,list_param) {            # loop over dates

  #This function produce daily prediction using monthly predicted clim surface.

  #The output is both daily prediction and daily deviation from monthly steps.

  #### Change this to allow explicitly arguments...

  #Arguments: 

  #1)index: loop list index for individual run/fit

  #2)clim_year_list: list of climatology files for all models...(12*nb of models)

  #3)sampling_obj: contains, data per date/fit, sampling information

  #4)dst: data at the monthly time scale

  #5)var: variable predicted -TMAX or TMIN

  #6)y_var_name: name of the variable predicted - dailyTMax, dailyTMin

  #7)out_prefix

  #8)out_path

  #9)list_models2 : interpolation model's formulas as string

  #10)interp_methods2: "gam","gwr","kriging"

  #11)s_raster: stack for covariates and toher variables

  #The output is a list of four shapefile names produced by the function:

  #1) list_temp: y_var_name

  #2) rast_clim_list: list of files for temperature climatology predictions

  #3) delta: list of files for temperature delta predictions

  #4) data_s: training data

  #5) data_v: testing data

  #6) sampling_dat: sampling information for the current prediction (date,proportion of holdout and sample number)

  #7) mod_kr: kriging delta fit, field package model object

  ### PARSING INPUT ARGUMENTS

  #list_param_runGAMFusion<-list(i,clim_yearlist,sampling_obj,var,y_var_name, out_prefix)

  rast_clim_yearlist<-list_param$clim_yearlist

  sampling_obj<-list_param$sampling_obj

  ghcn.subsets<-sampling_obj$ghcn_data

  sampling_dat <- sampling_obj$sampling_dat

  sampling <- sampling_obj$sampling_index

  var<-list_param$var

  y_var_name<-list_param$y_var_name

  out_prefix<-list_param$out_prefix

  dst<-list_param$dst #monthly station dataset

  out_path <-list_param$out_path

  list_models2 <-list_param$list_models2

  interp_method2 <- list_param$interp_method2

  s_raster <- list_param$s_raster

  sampling_month_obj <- list_param$sampling_month_obj

  daily_dev_sampling_dat <- list_param$daily_dev_sampling_dat

  index_d <- daily_dev_sampling_dat$index_d[i]

  index_m <- daily_dev_sampling_dat$index_m[i]

  use_clim_image <- list_param$use_clim_image # use predicted image as a base...rather than average Tmin at the station for delta

  join_daily <- list_param$join_daily # join monthly and daily station before calucating delta

  #use_clim_image

  ##########

  # STEP 1 - Read in information and get traing and testing stations

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

  #use index_d and index_m

  date<-strptime(daily_dev_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

  proj_str<-proj4string(dst) #get the local projection information from monthly data

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

  data_day<-ghcn.subsets[[index_d]]

  mod_LST <- ghcn.subsets[[index_d]][,match(LST_month, names(ghcn.subsets[[index_d]]))]  #Match interpolation date and monthly LST average

  data_day$LST <- as.data.frame(mod_LST)[,1] #Add the variable LST to the dataset

  dst$LST<-dst[[LST_month]] #Add the variable LST to the monthly dataset

  ind.training<-sampling[[index_d]]

  ind.testing <- setdiff(1:nrow(data_day), ind.training)

  data_s <- data_day[ind.training, ]   #Training dataset currently used in the modeling

  data_v <- data_day[ind.testing, ]    #Testing/validation dataset using input sampling

  ns<-nrow(data_s)

  nv<-nrow(data_v)

  #i=1

  date_proc<-sampling_dat$date[index_d]

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

  mo<-as.integer(strftime(date_proc, "%m"))          # current month of the date being processed

  day<-as.integer(strftime(date_proc, "%d"))

  year<-as.integer(strftime(date_proc, "%Y"))

  #Adding layer LST to the raster stack  

  #names(s_raster)<-covar_names

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

  s_raster<-dropLayer(s_raster,pos)      # If it exists drop layer

  LST<-subset(s_raster,LST_month)

  names(LST)<-"LST"

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

  #Now get monthly data...

  ghcn.month.subsets<-sampling_month_obj$ghcn_data

  sampling_month_dat <- sampling_month_obj$sampling_dat

  sampling_month_index <- sampling_month_obj$sampling_index

  dataset_month <-ghcn.month.subsets[[index_m]]

  mod_LST <- ghcn.month.subsets[[index_m]][,match(LST_month, names(ghcn.month.subsets[[index_m]]))]  #Match interpolation date and monthly LST average

  dataset_month$LST <- as.data.frame(mod_LST)[,1] #Add the variable LST to the dataset

  #change here...

  dst$LST<-dst[[LST_month]] #Add the variable LST to the monthly dataset

  proj_str<-proj4string(dst) #get the local projection information from monthly data

  ind.training_month <- sampling_month_index[[index_m]]

  ind.testing_month  <- setdiff(1:nrow(dataset_month), ind.training_month)

  data_month_s <- dataset_month[ind.training_month, ]   #Training dataset currently used in the modeling

  data_month_v <- dataset_month[ind.testing_month, ]    #Testing/validation dataset using input sampling

  modst <- data_month_s #training data for  monthhly predictions...

  ##########

  # STEP 2 - CLEAN DATA AND JOIN DAILY TO MONTHLY STATION INFORMATION

  ##########

  #if use join

  #modst<-dst[dst$month==mo,] #Subsetting dataset for the relevant month of the date being processed

  if (var=="TMIN"){

    modst$LSTD_bias <- modst$LST-modst$TMin; #That is the difference between the monthly LST mean and monthly station mean

  }

  if (var=="TMAX"){

    modst$LSTD_bias <- modst$LST-modst$TMax; #That is the difference between the monthly LST mean and monthly station mean    

  }

  #This may be unnecessary since LSTD_bias is already in dst?? check the info

  #Some loss of observations: LSTD_bias for January has only 56 out of 66 possible TMIN!!! We may need to look into this issue

  #to avoid some losses of station data...

  #Clearn out this part: make this a function call

  x<-as.data.frame(data_v)

  d<-as.data.frame(data_s)

  for (j in 1:nrow(x)){

    if (x$value[j]== -999.9){

      x$value[j]<-NA

    }

  }

  for (j in 1:nrow(d)){

    if (d$value[j]== -999.9){

      d$value[j]<-NA

    }

  }

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

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

  names(d)[pos]<-y_var_name

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

  names(x)[pos]<-y_var_name

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

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

  pos<-match("station",names(x)) #Find column with name station ID

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

  pos<-match("station",names(modst)) #Find column with name station ID

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

  ##########

  # STEP 3 - interpolate daily delta across space

  ##########

  #if used images

  # extract from image 

  #Change to take into account TMin and TMax

  if(use_clim_image==FALSE){

    #must join daily and monthly data first...

    dmoday <-merge(modst,d,by="id",suffixes=c("",".y2"))  

    xmoday <-merge(modst,x,by="id",suffixes=c("",".y2"))  

    mod_pat<-glob2rx("*.y2")   #remove duplicate columns that have ".y2" in their names

    var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names

    dmoday<-dmoday[,-var_pat] #dropping relevant columns

    mod_pat<-glob2rx("*.y2")   

    var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names

    xmoday<-xmoday[,-var_pat] #Removing duplicate columns

    data_v<-xmoday

    #coords <-dmoday[,c("coords.x1","coords.x2")]

    coords <-dmoday[,c("x","y")]

    coordinates(dmoday)<-coords

    proj4string(dmoday)<-proj_str

    #dmoday contains the daily tmax values for training with TMax/TMin being the monthly station tmax/tmin mean

    #xmoday contains the daily tmax values for validation with TMax/TMin being the monthly station tmax/tmin mean

    if (var=="TMIN"){

      daily_delta <-dmoday$dailyTmin-dmoday$TMin #daily detl is the difference between monthly and daily temperatures

    }

    if (var=="TMAX"){

      daily_delta <- dmoday$dailyTmax-dmoday$TMax

    }  

    #daily_delta <- dmoday[[y_var_name]] - 

    #only one delta in this case!!!

    #list(mod)

    if(is.null(list_models2)){ #change here...

      list_daily_delta_rast <- vector("list",length=1) #only one delta surface in this case!!

      list_mod_krtmp2 <- vector("list",length=1) #only one delta model in this case!!

      model_name<-paste("mod_stat_kr",sep="_")

      daily_delta_xy<-as.matrix(cbind(dmoday$x,dmoday$y))

      fitdelta<-Krig(daily_delta_xy,daily_delta,theta=1e5) #use TPS or krige

      mod_krtmp2 <- fitdelta

      #names(mod_krtmp2)[k] <- model_name

      #data_s$daily_delta<-daily_delta

      #rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

      rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

      rast_clim_mod <- stack(rast_clim_list)

      names(rast_clim_mod) <- names(rast_clim_list)

      rast_clim_month <- subset(rast_clim_mod,1) #example layer to interpolate to

      daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the of the daily devation

      #there is only one daily devation (delta) sruface in this case

      #To many I/O out of swap memory on atlas

      #Saving kriged surface in raster images

      data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                       sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

      raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

      writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

      list_daily_delta_rast[[1]] <- raster_name_delta

      list_mod_krtmp2[[1]] <- mod_krtmp2

    }

    if(!is.null(list_models2)){ #change here...

      list_daily_delta_rast <- vector("list",length=1) #several delta surfaces in this case but stored as one list!!

      list_mod_krtmp2 <- vector("list",length=1) #several delta model in this case but stored as one list!!

      dev_mod_name<-paste("dev_mod",1:length(list_models2),sep="") #change to more meaningful name?

      model_name<-paste("mod_stat_",sep="_")

      #Now generate file names for the predictions...

      list_out_filename<-vector("list",length(list_models2))

      names(list_out_filename)<- dev_mod_name  

      ##Change name...

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

        #j indicate which month is predicted, var indicates TMIN or TMAX

        data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                         sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],

                         "_",interp_method2,"_",dev_mod_name[j],sep="")

        raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

        list_out_filename[[j]]<-raster_name_delta

      }

      #Now call function

      #for (j in 1:length(list_models2)){

      dmoday$y_var <- daily_delta

      #coordinates(data_s)<-cbind(data_s$x,data_s$y)

      #proj4string(data_s)<-proj_str

      #coordinates(data_v)<-cbind(data_v$x,data_v$y)

      #proj4string(data_v)<-proj_str

      interp_area_obj <-interpolate_area_fun(interp_method2,list_models2,s_raster,list_out_filename,dmoday)

      rast_pred_list <- interp_area_obj$rast_pred_list

      rast_pred_list <-rast_pred_list[!sapply(rast_pred_list,is.null)] #remove NULL elements in list

      list_daily_delta_rast[[1]] <-rast_pred_list

      #names(list_daily_delta_rast) <- names(daily_delta_df)

      list_mod_krtmp2[[1]] <-interp_area_obj$mod_list      

    }

  }

  if(use_clim_image==TRUE){

    # User can choose to join daily and monthly station before interpolation:

    #-this ensures that the delta difference is "more" exact since its starting point is basesd on average value but there is risk to loose some stations

    #may need to change this option later!!

    #if jion_Daily is true then daily station used as training will match monthly station used as training

    if(join_daily==TRUE){

      dmoday <-merge(modst,d,by="id",suffixes=c("",".y2"))  

      xmoday <-merge(modst,x,by="id",suffixes=c("",".y2"))  

      mod_pat<-glob2rx("*.y2")   #remove duplicate columns that have ".y2" in their names

      var_pat<-grep(mod_pat,names(dmoday),value=FALSE) # using grep with "value" extracts the matching names

      dmoday<-dmoday[,-var_pat] #dropping relevant columns

      mod_pat<-glob2rx("*.y2")   

      var_pat<-grep(mod_pat,names(xmoday),value=FALSE) # using grep with "value" extracts the matching names

      xmoday<-xmoday[,-var_pat] #Removing duplicate columns

      data_v<-xmoday

    }else{

      dmoday<-d

      data_v<-x

    }

    #dmoday contains the daily tmax values for training with TMax/TMin being the monthly station tmax/tmin mean

    #xmoday contains the daily tmax values for validation with TMax/TMin being the monthly station tmax/tmin mean

    #coords <-dmoday[,c("coords.x1","coords.x2")]

    coords <-dmoday[,c("x","y")]

    coordinates(dmoday)<-coords

    proj4string(dmoday)<-proj_str

    #Now compute daily delta deviation from climatology layer:

    rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

    rast_clim_mod <- stack(rast_clim_list)

    names(rast_clim_mod) <- names(rast_clim_list)

    extract_data_s <-extract(rast_clim_mod,dmoday,df=TRUE)

    #list_daily_delta

    daily_delta_df <- dmoday[[y_var_name]] - extract_data_s  

    daily_delta_df <- daily_delta_df[,-1]

    names(daily_delta_df) <- paste(names(daily_delta_df),"_del",sep="")

    names(extract_data_s) <- paste(names(extract_data_s),"_m",sep="") # "m" for monthly predictions...

    dmoday <-spCbind(dmoday,extract_data_s) #contains the predicted clim at locations

    dmoday <-spCbind(dmoday,daily_delta_df) #contains the predicted clim at locations

    #Now krige  forevery model !! loop

    list_mod_krtmp2 <- vector("list",length=nlayers(rast_clim_mod)) 

    list_daily_delta_rast <- vector("list",length=nlayers(rast_clim_mod)) 

    names(list_daily_delta_rast) <- names(daily_delta_df)

    names(list_mod_krtmp2) <- names(daily_delta_df)

    for(k in 1:nlayers(rast_clim_mod)){

      daily_delta <- daily_delta_df[[k]] #Current daily deviation being process: the reference monthly prediction varies...

      #model_name<-paste("mod_kr","day",sep="_")

      model_name<- names(daily_delta_df)[k]

      if(is.null(list_models2)){

        daily_delta_xy<-as.matrix(cbind(dmoday$x,dmoday$y))

        fitdelta<-Krig(daily_delta_xy,daily_delta,theta=1e5) #use TPS or krige

        list_mod_krtmp2[[k]] <-fitdelta

        names(list_mod_krtmp2)[k] <- model_name

        #data_s$daily_delta<-daily_delta

        #rast_clim_list<-rast_clim_yearlist[[index_m]]  #select relevant monthly climatology image ...

        rast_clim_month <- subset(rast_clim_mod,1) #example layer to interpolate to

        daily_delta_rast<-interpolate(rast_clim_month,fitdelta) #Interpolation of the bias surface...

        #list_daily_delta_rast[[k]] <- raster_name_delta

        data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                         sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

        raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

        writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

        #writeRaster(r_spat, NAflag=NA_flag_val,filename=raster_name,bylayer=TRUE,bandorder="BSQ",overwrite=TRUE)   

        #raster_name_delta <- list_daily_delta_rast

        #mod_krtmp2 <- list_mod_krtmp2

        list_daily_delta_rast[[k]] <- raster_name_delta

      }

      if (!is.null(list_models2)){

        #list_formulas<-lapply(list_models,as.formula,env=.GlobalEnv) #mulitple arguments passed to lapply!!

        dev_mod_name<-paste("dev_mod",1:length(list_models2),sep="") #change to more meaningful name?

        #Now generate file names for the predictions...

        list_out_filename<-vector("list",length(list_models2))

        names(list_out_filename)<- dev_mod_name  

        ##Change name...

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

          #j indicate which month is predicted, var indicates TMIN or TMAX

          data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                           sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],

                           "_",interp_method2,"_",dev_mod_name[j],sep="")

          raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

          list_out_filename[[j]]<-raster_name_delta

        }

        #Now call function

        #for (j in 1:length(list_models2)){

        dmoday$y_var <- daily_delta

        #coordinates(data_s)<-cbind(data_s$x,data_s$y)

        #proj4string(data_s)<-proj_str

        #coordinates(data_v)<-cbind(data_v$x,data_v$y)

        #proj4string(data_v)<-proj_str

        interp_area_obj <-interpolate_area_fun(interp_method2,list_models2,s_raster,list_out_filename,dmoday)

        rast_pred_list <- interp_area_obj$rast_pred_list

        names(rast_pred_list) <- dev_mod_name

        rast_pred_list <-rast_pred_list[!sapply(rast_pred_list,is.null)] #remove NULL elements in list

        list_daily_delta_rast[[k]] <-rast_pred_list

        names(list_daily_delta_rast) <- names(daily_delta_df)

        mod_list <-interp_area_obj$mod_list

        names(mod_list) <- dev_mod_name

        list_mod_krtmp2[[k]] <-interp_area_obj$mod_list

      }

    }

    #Too many I/O out of swap memory on atlas

    #Saving kriged surface in raster images

    #delta_rast_s <-stack(list_daily_delta_rast)

    #names(delta_rast_s) <- names(daily_delta_df)

    #Should check that all delta images have been created for every model!!! remove from list empty elements!!

    #data_name<-paste("daily_delta_",y_var_name,"_",model_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

    #                 sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

    #raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

    #writeRaster(daily_delta_rast, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

    #data_name<-paste("daily_delta_",y_var_name,"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

    #                 sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

    #raster_name_delta<-file.path(out_path,paste(interpolation_method,"_",var,"_",data_name,out_prefix,".tif", sep=""))

    #writeRaster(delta_rast_s, filename=raster_name_delta,overwrite=TRUE)  #Writing the data in a raster file format...(IDRISI)

    #writeRaster(r_spat, NAflag=NA_flag_val,filename=raster_name,bylayer=TRUE,bandorder="BSQ",overwrite=TRUE)   

    #raster_name_delta <- list_daily_delta_rast

    #mod_krtmp2 <- list_mod_krtmp2

  }

  #########

  # STEP 4 - Calculate daily predictions - T(day) = clim(month) + delta(day)

  #########

  #if(use_clim_image==FALSE){

  #  list_daily_delta_rast <- rep(raster_name_delta,length=nlayers(rast_clim_mod))

  #}

  #Now predict daily after having selected the relevant month

  temp_list<-vector("list",nlayers(rast_clim_mod))  

  for (k in 1:nlayers(rast_clim_mod)){

    if(use_clim_image==TRUE){

      if (is.null(list_models2)){ 

        daily_delta_rast <- raster(list_daily_delta_rast[[k]]) #There is only one image of deviation per model if list_models2 is NULL

      }

      if (!is.null(list_models2)){ #then possible multiple daily dev predictions

        daily_delta_rast <- stack(unlist(list_daily_delta_rast[[k]]))

      }

      #daily_delta_rast <- subset(delta_rast_s,k)

    }

    #if use_clim_image==FALSE then daily__delta_rast already defined earlier...

    if(use_clim_image==FALSE){

      if (is.null(list_models2)){ 

        daily_delta_rast <- raster(list_daily_delta_rast[[1]]) #There is only one image of deviation per model if list_models2 is NULL

      }

      if (!is.null(list_models2)){ #then possible multiple daily dev predictions hence use stack

        daily_delta_rast <- stack(unlist(list_daily_delta_rast[[1]]))

      }

      #daily_delta_rast <- subset(delta_rast_s,k)

    }

    #rast_clim_month<-raster(rast_clim_list[[k]])

    rast_clim_month <- subset(rast_clim_mod,k) #long term monthly prediction

    if (is.null(list_models2)){ 

      temp_predicted<-rast_clim_month + daily_delta_rast

      data_name<-paste(y_var_name,"_predicted_",names(rast_clim_mod)[k],"_",sampling_month_dat$prop[index_m],"_",sampling_month_dat$run_samp[index_m],"_",

                       sampling_dat$date[index_d],"_",sampling_dat$prop[index_d],"_",sampling_dat$run_samp[index_d],sep="")

      raster_name<-file.path(out_path,paste(interpolation_method,"_",data_name,out_prefix,".tif", sep=""))