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

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

#                           

#of inputs and parameters.

#STAGE 1: LST climatology calculation

#STAGE 2: Covariates preparation: aspect, land cover, distance to coast etc.

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

#

#AUTHOR: Benoit Parmentier                                                                       

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

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

library(plotrix)

script_path<-"/home/parmentier/Data/IPLANT_project/env_layers_scripts/"

#script_path<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/"

#list_script_files<-

#CALLED FROM MASTER SCRIPT:

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

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

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

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

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

#Change input path?? for LST stage ??

in_path  <- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/input_data/"

#in_path <- out_path

out_path <- paste("/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data",

                  out_prefix,"/",sep="")

if (!file.exists(out_path)){

  dir.create(out_path)

lc_path<-"/home/layers/data/land-cover/lc-consensus-global"

infile_modis_grid<-"modis_sinusoidal_grid_world.shp" #Give path!!! NEED TO CHANGE THIS...

infile_elev<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif/srtm_1km.tif"  #this is the global file: replace later with the input produced by the DEM team

infile_canheight<-"/home/layers/data/land-cover/treeheight-simard2011/Simard_Pinto_3DGlobalVeg_JGR.tif"              #Canopy height

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

infile_reg_outline=""  #input region outline defined by polygon: none for Venezuel

#infile_reg_outline <- "OR83M_state_outline.shp"  #input region outline defined by polygon: Oregon

CRS_interp<-"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs";

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

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

out_suffix<-"_VE_05092013"

ref_rast_name<-""  #local raster name defining resolution, exent, local projection--. set on the fly??

#ref_rast_name<-"mean_day244_rescaled.rst"  #local raster name defining resolution, exent: oregon

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

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

#lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10") #use older version for continuity check to be changed

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

covar_names<-c(rnames,lc_names,lst_names)

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

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

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

if (stages_to_run[2]==2){

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

infile1<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/outline_venezuela_region__VE_01292013.shp"      #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"  #input region outline defined by polygon: Oregon

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

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

infile2<-"/home/layers/data/climate/ghcn/v2.92-upd-2012052822/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??

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

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

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

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

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

                                list_models,lst_avg,in_path,out_path,script_path,

                                interpolation_method)

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

raster_prediction_obj <-raster_prediction_fun(list_param_raster_prediction)

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

date_selected_results<-c("20100101") ##This is for year 2000!!!

#raster_prediciton_obj<-load_obj(paste("raster_prediction_obj","_","interpolation_method,

#                                y_var_name,out_prefix,sep="")

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