#########################    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 GAM+FUSION or GAM+CAI and other options added

#The interpolation is done first at the monthly time scale then delta surfaces are added.

#AUTHOR: Benoit Parmentier                                                                        

#DATE: 07/16/2013                                                                                 

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

#

# TO DO:

# 1) modify to make it general for any method i.e. make call to method e.g. gam_fus, gam_cai etc.

# 2) simplify and bundle validation steps, make it general--method_mod_validation

# 3) solve issues with log file recordings

# 4) output location folder on the fly???

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

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:

  #3)infile_daily

  #4)infile_locs:

  #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)dates_selected

  #

  #6 additional parameters for monthly climatology and more

  #16)list_models: model formulas in character vector

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

  #18)n_path

  #19)out_path

  #20)script_path: path to script

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

  ###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

  #6 additional parameters for monthly climatology and more

  list_models<-list_param_raster_prediction$list_models

  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

  setwd(out_path)

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

  if (var=="TMAX"){

    y_var_name<-"dailyTmax"                                       

  }

  if (var=="TMIN"){

    y_var_name<-"dailyTmin"                                       

  }

  ################# 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 OTEHR DATASETS  #################

  #infile_daily = gsub("/data/project/layers/commons/","/nobackup/aguzman4/climate/interp/testdata/",infile_daily)

  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

  if (dates_selected==""){

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

  }

  if (dates_selected!=""){

    dates<-dates_selected #dates to be predicted 

  }

  #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

  #Getting this name from command line

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

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

  #Input for sampling function...

  #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)

  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!!

  sampling_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=12

  if (interpolation_method=="gam_fusion"){

    list_param_runClim_KGFusion<-list(j,s_raster,covar_names,lst_avg,list_models,dst,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","var","y_var_name","out_prefix","out_path")

    print("Monthly")

    clim_method_mod_obj<-mclapply(1:12, list_param=list_param_runClim_KGFusion, runClim_KGFusion,mc.preschedule=FALSE,mc.cores = 6) #This is the end bracket from mclapply(...) statement

    print("done")

    file2<-file.path(out_path,paste(interpolation_method,"_mod_",y_var_name,out_prefix,".RData",sep=""))

    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

  }

  if (interpolation_method=="gam_CAI"){

    list_param_runClim_KGCAI<-list(j,s_raster,covar_names,lst_avg,list_models,dst,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","var","y_var_name","out_prefix","out_path")

    clim_method_mod_obj<-mclapply(1:12, list_param=list_param_runClim_KGCAI, runClim_KGCAI,mc.preschedule=FALSE,mc.cores = 6) #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="")))

    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)

  ################## 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

  print("daily")

  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...

  if (interpolation_method=="gam_CAI" | interpolation_method=="gam_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,sampling_obj,dst,var,y_var_name, interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_daily_deviation)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","interpolation_method","out_prefix","out_path")

    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 6) #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")

    method_mod_obj<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 12) #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_day),list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 6) #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_day),list_param=list_param_run_prediction_gwr_daily,run_interp_day_fun,mc.preschedule=FALSE,mc.cores = 12) #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

  print("validation")

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

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

    tmp<-method_mod_obj[[i]][[y_var_name]]  #y_var_name is the variable predicted (dailyTmax or dailyTmin)

    list_tmp[[i]]<-tmp

  }

  rast_day_yearlist<-list_tmp #list of predicted images over full year... 

  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")

  list_param_validation<-list(i,rast_day_yearlist,method_mod_obj,y_var_name, out_prefix, out_path)

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

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

  #test_val<-calculate_accuracy_metrics(1,list_param_validation)

  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)

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

  ##Create data.frame with valiation metrics for a full year

  tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v")

  rownames(tb_diagnostic_v)<-NULL #remove row names

  #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)

  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)

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

  #Will add more information to be returned

  if (interpolation_method=="gam_CAI" | interpolation_method=="gam_fusion"){

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

                                summary_metrics_v,summary_month_metrics_v)

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

                                    "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,

                                summary_metrics_v,summary_month_metrics_v)

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

                                    "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="")))

  }

  return(raster_prediction_obj)

}

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

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

  rasterOptions(maxmemory=max_mem,timer=TRUE)

  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

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

  j=6 #12

  #quit()

  print("stage 4 DONE")

##################  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: 05/07/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) #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)

      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,k=5)) #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) 

}

####

#TODO:

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

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

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

  prop_month<-0 #proportion retained for validation

  run_samp<-1

  #### 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]]

  #Create formulas object from list of characters...

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

  #TMax to model...

  if (var=="TMAX"){   

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

  }

  if (var=="TMIN"){   

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

  }

  #Fit gam models using data and list of formulas

  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?

  names(mod_list)<-cname

  #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(mod_list))

  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_",j,"_",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

  }

  #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

  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

  #Saving kriged surface in raster images

  #data_name<-paste("clim_month_",j,"_",model_name,"_",prop_month,

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

  #raster_name_clim<-paste("fusion_",data_name,out_prefix,".tif", sep="")

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

  #now climatology layer

  #clim_rast<-LST-bias_rast

  data_name<-paste(var,"_clim_month_",j,"_",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,mod_list,list_formulas)

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

  save(clim_obj,file= file.path(out_path,paste("clim_obj_CAI_month_",j,"_",var,"_",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)

  options(error=function()traceback(2))

  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

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

  prop_month<-0 #proportion retained for validation

  run_samp<-1 #This option can be added later on if/when neeeded

  #### 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

  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

  }

  #### STEP3:  NOW FIT AND PREDICT  MODEL

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

  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?

  names(mod_list)<-cname

  #Now generate file names for the predictions...

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

  names(list_out_filename)<-cname  

  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_",j,"_",cname[k],"_",prop_month,

                     "_",run_samp,sep="")

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

    list_out_filename[[k]]<-raster_name

  }

  print("predict raster model")

  #now 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)

  print("done")

  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_",j,"_",names(rast_clim_list)[k],"_",prop_month,

                     "_",run_samp,sep="")

    raster_name<-file.path(out_path,paste("fusion_",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_",j,"_",model_name,"_",prop_month,

                     "_",run_samp,sep="")

    raster_name_bias<-file.path(out_path,paste("fusion_",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_",j,"_",model_name,"_",prop_month,

                     "_",run_samp,sep="")

    raster_name_clim<-file.path(out_path,paste("fusion_",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")){

    print("Error with FitBias")

    #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

  clim_obj<-list(rast_bias_list,rast_clim_list,data_month,mod_list,list_formulas)

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

  save(clim_obj,file= file.path(out_path,paste("clim_obj_month_",j,"_",var,"_",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

  #

  #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_day

  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

  ##########

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

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

  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

  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[[i]]

  mod_LST <- ghcn.subsets[[i]][,match(LST_month, names(ghcn.subsets[[i]]))]  #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[[i]]

  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[i]

  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"))

  ##########

  # STEP 2 - JOIN DAILY AND MONTHLY STATION INFORMATION

  ##########

  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

    }

  }

  d[1] <- NULL

  x[1] <- NULL

  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...

  mergeTry<-try(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

  #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

  ##########

  # STEP 3 - interpolate daily delta across space

  ##########

  #Change to take into account TMin and TMax

  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_xy<-as.matrix(cbind(dmoday$x,dmoday$y))

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

  mod_krtmp2<-fitdelta

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

  data_s<-dmoday #put the 

  data_s$daily_delta<-daily_delta

  #########

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

  #########

  rast_clim_list<-rast_clim_yearlist[[mo]]  #select relevant month

  rast_clim_month<-raster(rast_clim_list[[1]])

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

  #Saving kriged surface in raster images

  data_name<-paste("daily_delta_",y_var_name,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],

                   "_",sampling_dat$run_samp[i],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)

  #Now predict daily after having selected the relevant month

  temp_list<-vector("list",length(rast_clim_list))  

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

    rast_clim_month<-raster(rast_clim_list[[j]])

    temp_predicted<-rast_clim_month+daily_delta_rast

    data_name<-paste(y_var_name,"_predicted_",names(rast_clim_list)[j],"_",

                     sampling_dat$date[i],"_",sampling_dat$prop[i],

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

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

    writeRaster(temp_predicted, filename=raster_name,overwrite=TRUE) 

    temp_list[[j]]<-raster_name

  }

  ##########

  # STEP 5 - Prepare output object to return

  ##########

  mod_krtmp2<-fitdelta

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

  names(temp_list)<-names(rast_clim_list)

  delta_obj<-list(temp_list,rast_clim_list,raster_name_delta,data_s,

                  data_v,sampling_dat[i,],mod_krtmp2)

  obj_names<-c(y_var_name,"clim","delta","data_s","data_v",

               "sampling_dat",model_name)

  names(delta_obj)<-obj_names 

  save(delta_obj,file= file.path(out_path,paste("delta_obj_",var,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],

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

  return(delta_obj)

}

  if (interpolation_method=="gam_fusion"){

  }

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

  parser.add_option("-d","--daynight",dest="dn",default=0,

      help="Use 0 for both(default),1 for day only,2 for night only")