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

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

#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

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

#

# TO DO:

  ##Function to predict temperature interpolation with 21 input parameters

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

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

  #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

  #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(spgwr)   #GWR method

  #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

  #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

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

    y_var_name<-"dailyTmax"                                       

  }

    y_var_name<-"dailyTmin"                                       

  }

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

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

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

  if (dates_selected==""){

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

  }

  if (dates_selected!=""){

    dates<-dates_selected #dates to be predicted 

  }

  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

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

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

  #First predict at the monthly time scale: climatology

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

      file=log_fname,sep="\n")

  j=12

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

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

    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

    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

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

    #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

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

  #rast_clim_yearlist<-list_tmp

      file=log_fname,sep="\n")

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

    #test<-runGAMFusion(1,list_param=list_param_runGAMFusion)

    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 = 9) #This is the end bracket from mclapply(...) statement

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

    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:11,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

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

  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)

  }

      file=log_fname,sep="\n")

  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

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

  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

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

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

  #close log_file connection and add meta data

  #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=="gam_daily" | interpolation_method=="kriging_daily" | interpolation_method=="gwr_daily"){

  }

  return(raster_prediction_obj)

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