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

#

#AUTHOR: Benoit Parmentier                                                                       

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

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

#CALLED FROM MASTER SCRIPT:

#/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/climatology_03192013.py

source(file.path(script_path,"covariates_production_temperatures_03212013.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"))

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

##Paths to inputs and output

var<-"TMAX"

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

#in_path <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/input"

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('h11v08','h11v07','h12v07','h12v08','h10v07','h10v08') #tile for Venezuel and surrounding area

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_region_name<-"_venezuela_region" #generated on the fly

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

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

                                      "out_suffix","ref_rast_name","covar_names") 

if (stages_to_run[2]==2){

  #Transform into function...

  #/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates_production_temperatures_03212013.R

  infile_covar_brick<-covariates_production_temperature(list_param_covar_production)

}

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

#set up earlier

#var <- "TMIN"           # name of the variables to keep: TMIN, TMAX or PRCP --already set up earlier

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

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

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

#covar_names see stage 2

#specific to this stage

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

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

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

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

names(list_param_prep)<-cnames

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

               "y_var ~ s(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(LC2)",

#Choose interpolation method...

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

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

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

raster_prediction_obj <-raster_prediction_fun(list_param_raster_prediction)

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

#raster_prediction_obj<-load_obj("raster_prediction_obj_dailyTmin_365d_GAM_fus5_all_lstd_03292013.RData")

#raster_prediction_obj<-load_obj("raster_prediction_obj_gam_CAI_dailyTmax_365d_GAM_CAI_all_lst_04162013.RData")

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

#                                y_var_name,out_prefix,sep="")

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

#plots_assessment_by_date<-function(j,list_param){

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

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