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

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,"covariates_production_temperatures_05132013.R"))

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

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

stages_to_run<-c(0,0,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

#Change input path?? for LST stage ??

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

list_tiles_modis <- c("h11v08,h11v07,h12v07,h12v08,h10v07,h10v08") #tile for Venezuela and surrounding area

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

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

hdfdir =  '/home/parmentier/Data/IPLANT_project/MOD11A1_tiles'

download=1

clim_calc=0

#end_month= "12"

#start_month= "1"

}else{

  #Provide brick of covariates if stage 2 is not run

  infile_covariates<-"/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script and used in stage 3 and stage 4

  #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

}

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

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

#covar_names see stage 2

#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("2000","2011") #right bound not included in the range!!

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

names(list_param_prep)<-cnames

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

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

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

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