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

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

}

##################    Master script for temperature predictions  #######################################

############################ TMIN AND TMAX predictions ##########################################

#                           

##This script produces intperpolated surface of TMIN and TMAX for specified processing region(s) given sets 

#of inputs and parameters.

#STAGE 1: LST climatology downloading and/or calculation

#STAGE 2: Covariates preparation for study/processing area: calculation of covariates (spect,land cover,etc.) and reprojection

#STAGE 3: Data preparation: meteorological station database query and extraction of covariates values from raster brick

#STAGE 4: Raster prediction: run interpolation method (-- gam fusion, gam CAI, ...) and perform validation 

#STAGE 5: Output analyses: assessment of results for specific dates...

#

#AUTHOR: Benoit Parmentier                                                                       

#DATE: 07/18/2013                                                                                 

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363, TASK$568--   

## TODO:

# Modify code for stage 1 and call python script from R in parallel

# Add options to run only specific stage + additional out_suffix?

# Make master script a function?

# Add log file for master script,add function to collect inputs and outputs

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

###Loading R library and packages   

library(RPostgreSQL)

library(maps)

library(maptools)

library(parallel)

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(spgwr)

library(reshape)

library(plotrix)

######## PARAMETERS FOR WORK FLOW #########################

### Need to add documentation ###

#Adding command line arguments to use mpiexec

args<-commandArgs(TRUE)

script_path<-"/nobackup/aguzman4/climateLayers/climateCode/environmental-layers/climate/research/oregon/interpolation"

dataHome<-"/nobackup/aguzman4/climateLayers/interp/testdata/"

script_path2<-"/nobackup/aguzman4/climateLayers/climateCode/environmental-layers/climate/research/world/interpolation"

#CALLED FROM MASTER SCRIPT:

modis_download_script <- file.path(script_path,"modis_download_05142013.py") # LST modis download python script

clim_script <- file.path(script_path,"climatology_05312013.py") # LST climatology python script

grass_setting_script <- file.path(script_path,"grass-setup.R") #Set up system shell environment for python+GRASS

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

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

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

source(file.path(script_path2,"GAM_fusion_analysis_raster_prediction_multisampling.R"))

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

#source(file.path(script_path,"results_covariates_database_stations_output_analyses_04012013.R")) #to be completed

#FUNCTIONS CALLED FROM GAM ANALYSIS RASTER PREDICTION ARE FOUND IN...

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

source(file.path(script_path2,"GAM_fusion_function_multisampling.R")) #Include GAM_CAI

source(file.path(script_path,"interpolation_method_day_function_multisampling.R")) #Include GAM_day

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

#stages_to_run<-c(1,2,3,4,5) #May decide on antoher strategy later on...

stages_to_run<-c(0,0,0,4,0) #MRun only raster fitting, prediction and assessemnt (providing lst averages, covar brick and met stations)

#If stage 2 is skipped then use previous covar object

covar_obj_file<-args[6]

#If stage 3 is skipped then use previous met_stations object

met_stations_outfiles_obj_file<-args[7]

var<-"TMAX"    

out_prefix<-args[1] #User defined output prefix

out_suffix<-args[2] #Regional suffix

out_suffix_modis<-args[3] #pattern to find tiles produced previously

interpolation_method<-c("gam_fusion")#,"gam_CAI") #Testing mem usage

out_path<-"/nobackup/aguzman4/climateLayers/output/"

out_path<-paste(out_path,out_prefix,sep="")

if(!file.exists(out_path)){

   dir.create(out_path)

}

lc_path<-"/nobackup/aguzman4/climateLayers/interp/testdata/data_workflow/inputs/lc-consensus-global"

infile_modis_grid<-"/nobackup/aguzman4/climateLayers/interp/testdata/data_workflow/inputs/modis_grid/modis_sinusoidal_grid_world.shp" #modis grid tiling system, global

infile_elev<-"/nobackup/aguzman4/climateLayers/interp/testdata/data_workflow/inputs/dem-cgiar-srtm-1km-tif/srtm_1km.tif"  #elevation at 1km, global extent to be replaced by the new fused product 

infile_canheight<-"/nobackup/aguzman4/climateLayers/interp/testdata/data_workflow/inputs/treeheight-simard2011/Simard_Pinto_3DGlobalVeg_JGR.tif"         #Canopy height, global extent

infile_distoc <- "/nobackup/aguzman4/climateLayers/interp/testdata/data_workflow/inputs/distance_to_coast/GMT_intermediate_coast_distance_01d_rev.tif" #distance to coast, global extent at 0.01 deg

#infile_reg_outline <- "/nobackup/aguzman4/climateLayers/subsetWGS84/shapefiles/shp_40.00_-115.00_1929.shp"

infile_reg_outline <- args[4]

ref_rast_name<-args[5]

buffer_dist<-0 #not in use yet, must change climatology step to make sure additional tiles are downloaded and LST averages

               #must also be calculated for neighbouring tiles.

list_tiles_modis <- c("h08v04,h09v04") #tiles for Oregon

CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84

CRS_interp <-CRS_locs_WGS84 

out_region_name<-""

#The names of covariates can be changed...these names should be output/input from covar script!!!

rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev_s","slope","aspect","CANHGHT","DISTOC")

lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")

lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12",

               "nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",

               "nobs_09","nobs_10","nobs_11","nobs_12")

covar_names<-c(rnames,lc_names,lst_names)

list_val_range <-c("lon,-180,180","lat,-90,90","N,-1,1","E,-1,1","N_w,-1,1","E_w,-1,1","elev_s,0,6000","slope,0,90",

                   "aspect,0,360","DISTOC,-0,10000000","CANHGHT,0,255","LC1,0,100","LC5,0,100","mm_01,-15,50",

                   "mm_02,-15,50","mm_03,-15,50","mm_04,-15,50","mm_05,-15,50","mm_06,-15,50","mm_07,-15,50",

                   "mm_08,-15,50","mm_09,-15,50","mm_10,-15,50","mm_11,-15,50","mm_12,-15,50")

############ STAGE 1: LST Climatology ###############

#print("Starting stage 1")

#Parameters,Inputs from R to Python??

start_year = "2001"

end_year = "2010"

hdfdir <- "/nobackup/aguzman4/climate/interp/testdata/data_workflow/inputs/MOD11A1_tiles" #path directory to MODIS data

download=0 #download MODIS product if 1

clim_calc=1 #calculate lst averages/climatology if 1

list_param_download_clim_LST_script <- list(list_tiles_modis,start_year,end_year,hdfdir,

                                            var,grass_setting_script,modis_download_script, clim_script,

                                            download,clim_calc,out_suffix_modis)

names(list_param_download_clim_LST_script)<-c("list_tiles_modis","start_year","end_year","hdfdir",

                                              "var","grass_setting_script","modis_download_script","clim_script",

                                              "download","clim_calc","out_suffix_modis")

no_tiles <- length(unlist(strsplit(list_tiles_modis,",")))  # transform string into separate element in char vector

if (stages_to_run[1]==1){

  clim_production_obj <-lapply(1:no_tiles, list_param=list_param_download_clim_LST_script, download_calculate_MODIS_LST_climatology)

}

#Collect LST climatology list as output???

############ STAGE 2: Covariate production ################

#list of 18 parameters

list_param_covar_production<-list(var,out_path,lc_path,infile_modis_grid,infile_elev,infile_canheight,

                                  infile_distoc,list_tiles_modis,infile_reg_outline,CRS_interp,CRS_locs_WGS84,out_region_name,

                                  buffer_dist,list_val_range,out_suffix,out_suffix_modis,ref_rast_name,hdfdir,covar_names) 

names(list_param_covar_production)<-c("var","out_path","lc_path","infile_modis_grid","infile_elev","infile_canheight",

                                      "infile_distoc","list_tiles_modis","infile_reg_outline","CRS_interp","CRS_locs_WGS84","out_region_name",

                                      "buffer_dist","list_val_range","out_suffix","out_suffix_modis","ref_rast_name","hdfdir","covar_names") 

## Modify to store infile_covar_brick in output folder!!!

if (stages_to_run[2]==2){

  covar_obj <- covariates_production_temperature(list_param_covar_production)

  infile_covariates <- covar_obj$infile_covariates

  infile_reg_outline <- covar_obj$infile_reg_outline

  covar_names<- covar_obj$covar_names

}else{

  covar_obj <-load_obj(covar_obj_file)

  infile_covariates <- covar_obj$infile_covariates

  #Region passed as input from command line

  #infile_reg_outline <- covar_obj$infile_reg_outline

  covar_names<- covar_obj$covar_names

}

#Note that if stages_to_run[2]!=2, then use values defined at the beginning of the script for infile_covariates and infile_reg_outline

############# STAGE 3: Data preparation ###############

#specific to this stage

db.name <- "ghcn" #       # name of the Postgres database

range_years<-c("2010","2011") #right bound not included in the range!!

range_years_clim<-c("2001","2011") #right bound not included in the range!!

infile_ghncd_data <-"/nobackup/aguzman4/climateLayers/stations/ghcnd-stations.txt"                              #This is the textfile of station locations from GHCND

#qc_flags_stations<-c("0","S")    #flags allowed for screening after the query from the GHCND??

qc_flags_stations<-c("0")    #flags allowed for screening after the query from the GHCND??

#infile_covariates and infile_reg_outline defined in stage 2 or at the start of script...

#list of 12 parameters for input in the function...

list_param_prep<-list(db.name,var,range_years,range_years_clim,infile_reg_outline,infile_ghncd_data,infile_covariates,CRS_locs_WGS84,out_path,covar_names,qc_flags_stations,out_prefix)

cnames<-c("db.name","var","range_years","range_years_clim","infile_reg_outline","infile_ghncd_data","infile_covariates","CRS_locs_WGS84","out_path","covar_names","qc_flags_stations","out_prefix")

names(list_param_prep)<-cnames

##### RUN SCRIPT TO GET STATION DATA WITH COVARIATES #####

if (stages_to_run[3]==3){

  list_outfiles<-database_covariates_preparation(list_param_prep)

}else{

  list_outfiles <-load_obj(met_stations_outfiles_obj_file) 

}

############### STAGE 4: RASTER PREDICTION #################

#Collect parameters from the previous stage: data preparation stage

#3 parameters from output

infile_monthly<-list_outfiles$monthly_covar_ghcn_data #outile4 from database_covar script

infile_daily<-list_outfiles$daily_covar_ghcn_data  #outfile3 from database_covar script

infile_locs<- list_outfiles$loc_stations_ghcn #outfile2? from database covar script

list_param_data_prep <- list(infile_monthly,infile_daily,infile_locs,infile_covariates,covar_names,var,out_prefix,CRS_locs_WGS84)

names(list_param_data_prep) <- c("infile_monthly","infile_daily","infile_locs","infile_covariates","covar_names","var","out_prefix","CRS_locs_WGS84")

#Set additional parameters

#Input for sampling function...

seed_number<- 100  #if seed zero then no seed?     

nb_sample<-1           #number of time random sampling must be repeated for every hold out proportion

step<-0         

constant<-0             #if value 1 then use the same samples as date one for the all set of dates

prop_minmax<-c(0.3,0.3)  #if prop_min=prop_max and step=0 then predicitons are done for the number of dates...

#dates_selected<-c("20100101","20100102","20100103","20100901") # Note that the dates set must have a specific format: yyymmdd

dates_selected<-"" # if empty string then predict for the full year specified earlier

screen_data_training<-FALSE #screen training data for NA and use same input training for all models fitted

#Models to run...this can be changed for each run

#LC1: Evergreen/deciduous needleleaf trees

#Combination for test run:

#list_models<-c("y_var ~ s(elev_s)",

#                "y_var ~ s(LST)",

#                "y_var ~ s(lat,lon)+ s(elev_s)",

#                "y_var ~ te(lat,lon,elev_s)",

#                "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST)", 

#                "y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC1)") 

list_models<-c("y_var ~ s(lat,lon) + s(elev_s)",

               "y_var ~ s(lat,lon) + s(elev_s) + s(LST)")#,

               #"y_var ~ s(lat,lon) + s(elev_s) + s(N_w,E_w) + s(LST) + s(LC1)")

#Default name of LST avg to be matched               

lst_avg<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12")  

#Collect all parameters in a list

list_param_raster_prediction<-list(list_param_data_prep,screen_data_training,

                                seed_number,nb_sample,step,constant,prop_minmax,dates_selected,

                                list_models,lst_avg,out_path,script_path,

                                interpolation_method)

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

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

                                "list_models","lst_avg","out_path","script_path",

                                "interpolation_method")

raster_prediction_obj <-raster_prediction_fun(list_param_raster_prediction)

############## STAGE 5: OUTPUT ANALYSES ##################

date_selected_results<-c("20100101") 

list_param_results_analyses<-list(out_path,script_path,raster_prediction_obj,interpolation_method,

                                  infile_covariates,covar_names,date_selected_results,var,out_prefix)

names(list_param_results_analyses)<-c("out_path","script_path","raster_prediction_obj","interpolation_method",

                     "infile_covariates","covar_names","date_selected_results","var","out_prefix")

#plots_assessment_by_date<-function(j,list_param){

if (stages_to_run[5]==5){

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

  #Use lapply or mclapply

  summary_v_day <-plots_assessment_by_date(1,list_param_results_analyses)

  #Call as function...

}

###############   END OF SCRIPT   ###################

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

# #LAND COVER INFORMATION

# LC1: Evergreen/deciduous needleleaf trees

# LC2: Evergreen broadleaf trees

# LC3: Deciduous broadleaf trees

# LC4: Mixed/other trees

# LC5: Shrubs

# LC6: Herbaceous vegetation

# LC7: Cultivated and managed vegetation

# LC8: Regularly flooded shrub/herbaceous vegetation

# LC9: Urban/built-up

# LC10: Snow/ice

# LC11: Barren lands/sparse vegetation

# LC12: Open water