#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                                                                        

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

#

# TO DO:

  library(plotrix)

  library(maptools)

  library(gdata) #Nesssary to use cbindX

  ### Parameters and arguments

  #PARSING INPUTS/ARGUMENTS

#   

  interpolation_method<-list_param_raster_prediction$interpolation_method

  setwd(out_path)

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

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

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

    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

    #gamclim_fus_mod<-mclapply(1:6, list_param=list_param_runClim_KGFusion, runClim_KGFusion,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))

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

      file=log_fname,sep="\n")

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

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

    list_param_run_prediction_daily_deviation <-list(i,clim_yearlist,sampling_obj,dst,var,y_var_name, interpolation_method,out_prefix,out_path)

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

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

    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)

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

  }

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

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

#

#AUTHOR: Benoit Parmentier                                                                       

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

##SCRIPT USED FOR THE PREDICTIONS: Source or list all scripts here to avoid confusion on versions being run!!!!

#CALLED FROM MASTER SCRIPT:

modis_download_script <- file.path(script_path,"modis_download_05142013.py") # LST modis download 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_05302013.R"))

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

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

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

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

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

source(file.path(script_path,"GAM_fusion_function_multisampling_05212013.R")) #Include GAM_CAI

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

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

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

var<-"TMAX" # variable being interpolated

#interpolation_method<-c("gam_fusion","gam_CAI","gam_daily") #other otpions to be added later

#interpolation_method<-c("gam_fusion") #other otpions to be added later

#interpolation_method<-c("gam_daily") #other otpions to be added later

#interpolation_method<-c("kriging_daily") #other otpions to be added later

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

if (!file.exists(out_path)){

#infile_reg_outline<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/outline_venezuela_region__VE_01292013.shp" 

#infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region_TMIN__VE_03192013.tif" #covariates stack for TMIN

#infile_covariates<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates_Oregon_region_TMAX__OR_04052013.tif" #Oregon covar TMAX from earlier codes...for continuity

#This is the shape file of outline of the study area                                                      #It is an input/output of the covariate script

#infile_reg_outline <-"OR83M_state_outline.shp" #remove this parameter!!!

#this may be redundant with infile_reg_outline

#covar_names see stage 2

#"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80"

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

#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","CANHEIGHT","DISTOC")

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

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

infile_ghncd_data <-"/home/layers/data/climate/ghcn/v2.92-upd-2012052822/ghcnd-stations.txt"                              #This is the textfile of station locations from GHCND

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

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

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