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

#

#AUTHOR: Benoit Parmentier                                                                       

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

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

# Modify code for stage 2, make it a function and fully automated (distoc var)

# 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   

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

clim_script <- file.path(script_path,"climatology_05142013.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,"Database_stations_covariates_processing_function_05212013.R"))

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

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

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

var<-"TMAX" # variable being interpolated

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

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

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

if (!file.exists(out_path)){

  dir.create(out_path)

  #} else{

#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

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

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

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

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

start_year = "2001"

end_year = "2010"

#hdfdir =  '/home/layers/commons/modis/MOD11A1_tiles' #destination file where hdf files are stored locally after download.

hdfdir <- "/data/project/layers/commons/modis/MOD11A1_tiles"

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

                                            var,grass_setting_script,modis_download_script, clim_script,

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

if (stages_to_run[1]==1){

}

#Collect LST climatology list as output???

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

#list of 17 parameters

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

                                  out_suffix,out_suffix_modis,ref_rast_name,hdfdir,covar_names) 

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

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

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

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

names(list_param_prep)<-cnames

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

#names(outfiles_obj)<- c("loc_stations","loc_stations_ghcn","daily_covar_ghcn_data","monthly_covar_ghcn_data")

#Set additional parameters

#Input for sampling function...

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

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

               "y_var ~ s(LST)",

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

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

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

                                interpolation_method)

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

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

                                "interpolation_method")

raster_prediction_obj <-raster_prediction_fun(list_param_raster_prediction)

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

                                  infile_covariates,covar_names,date_selected_results,var,out_prefix)

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

  #Call as function...

}