#########################    Raster prediction GAM FUSION    ####################################

############################ 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)GAM fusion: 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: 03/12/2013                                                                                 

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

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

raster_prediction_gam_fusion<-function(list_param_raster_prediction){

  ##Function to predict temperature interpolation with 21 input parameters

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

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

  in_path<-list_param_raster_prediction$in_path

  out_path<-list_param_raster_prediction$out_path

  script_path<-list_param_raster_prediction$script_path

  interpolation_method<-list_param_raster_prediction$interpolation_method

  setwd(in_path)

  source(file.path(script_path,"sampling_script_functions_03122013.R"))

  source(file.path(script_path,"GAM_fusion_function_multisampling_03122013.R"))

  source(file.path(script_path,"GAM_fusion_function_multisampling_validation_metrics_03122013.R"))

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

  if (file.exists(log_fname)){  #Stop the script???

    file.remove(log_fname)

    log_file<-file(log_fname,"w")

  }

  if (!file.exists(log_fname)){

    log_file<-file(log_fname,"w")

  }

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

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

             con=log_file,sep="\n")

  writeLines("Starting script process time:",con=log_file,sep="\n")

  writeLines(as.character(time1),con=log_file,sep="\n")    

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

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

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

  stat_loc<-readOGR(dsn=in_path,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 of covariate brick covariates can 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

  pos<-match("elev",names(s_raster))

  names(s_raster)[pos]<-"elev_1"

  #Screen for extreme values": this needs more thought, min and max val vary with regions

  #min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets)

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

  #Reading monthly data

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

  dst_all<-data3

  dst<-data3

  ### TO DO -important ###

  #Cleaning/sceerniging functions for daily stations, monthly stations and covariates?? do this during the preparation stage!!!??

  ###

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

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

  sampling_obj<-sampling_training_testing(list_param_sampling)

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

  #First predict at the monthly time scale: climatology

  writeLines("Predictions at monthly scale:",con=log_file,sep="\n")

  t1<-proc.time()

  j=12

  #browser()

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

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

  #source(file.path(script_path,"GAM_fusion_function_multisampling_03122013.R"))

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

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

  save(gamclim_fus_mod,file= paste("gamclim_fus_mod",out_prefix,".RData",sep=""))

  t2<-proc.time()-t1

  writeLines(as.character(t2),con=log_file,sep="\n")

  #now get list of raster clim layers

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

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

    tmp<-gamclim_fus_mod[[i]]$clim

    list_tmp[[i]]<-tmp

  }

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

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

  #rast_clim_yearlist<-list_tmp

  clim_yearlist<-list_tmp

  #Second predict at the daily time scale: delta

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

  writeLines("Predictions at the daily scale:",con=log_file,sep="\n")

  t1<-proc.time()

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

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

  names(list_param_runGAMFusion)<-c("list_index","clim_yearlist","sampling_obj","dst","var","y_var_name","out_prefix")

  #test<-mclapply(1:18, runGAMFusion,list_param=list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9)

  gam_fus_mod<-mclapply(1:length(sampling_obj$ghcn_data_day),list_param=list_param_runGAMFusion,runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

  #gam_fus_mod<-mclapply(1:length(sampling_obj$ghcn_data_day),runGAMFusion,list_param_runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

  #gam_fus_mod<-mclapply(1:length(ghcn.subsets), runGAMFusion,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

  save(gam_fus_mod,file= paste("gam_fus_mod",out_prefix,".RData",sep=""))

  t2<-proc.time()-t1

  writeLines(as.character(t2),con=log_file,sep="\n")

  #browser()

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

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

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

    tmp<-gam_fus_mod[[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

  writeLines("Validation step:",con=log_file,sep="\n")

  t1<-proc.time()

  #calculate_accuary_metrics<-function(i)

  list_param_validation<-list(i,rast_day_yearlist,gam_fus_mod,y_var_name, out_prefix)

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

  #gam_fus_validation_mod<-mclapply(1:length(gam_fus_mod), calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

  gam_fus_validation_mod<-mclapply(1:length(gam_fus_mod), list_param=list_param_validation, calculate_accuracy_metrics,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

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

  save(gam_fus_validation_mod,file= paste("gam_fus_validation_mod",out_prefix,".RData",sep=""))

  t2<-proc.time()-t1

  writeLines(as.character(t2),con=log_file,sep="\n")

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

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

  tb_diagnostic_v<-extract_from_list_obj(gam_fus_validation_mod,"metrics_v")

  rownames(tb_diagnostic_v)<-NULL #remove row names

  #Call function to create plots of metrics for validation dataset

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

  summary_metrics<-boxplot_from_tb(tb_diagnostic_v,metric_names,out_prefix)

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

  ##Write out information concerning accuracy and predictions

  outfile<-file.path(in_path,paste("assessment_measures_",out_prefix,".txt",sep=""))

  write.table(tb_diagnostic_v,file= outfile,row.names=FALSE,sep=",")

  write.table(summary_metrics[[1]], file= outfile, append=TRUE,sep=",") #write out avg

  write.table(summary_metrics[[2]], file= outfile, append=TRUE,sep=",") #write out median

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

  #close log_file connection and add meta data

  writeLines("Finished script process time:",con=log_file,sep="\n")

  time2<-proc.time()-time1

  writeLines(as.character(time2),con=log_file,sep="\n")

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

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

             con=log_file,sep="\n")

  writeLines("End of script",con=log_file,sep="\n")

  close(log_file)

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

  #Will add more information to be returned

  raster_prediction_obj<-list(gamclim_fus_mod,gam_fus_mod,gam_fus_validation_mod,tb_diagnostic_v)

  names(raster_prediction_obj)<-c("gamclim_fus_mod","gam_fus_mod","gam_fus_validation_mod","tb_diagnostic_v")  

  save(raster_prediction_obj,file= paste("raster_prediction_obj_",out_prefix,".RData",sep=""))

  return(raster_prediction_obj)

}

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

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