#different covariates using two baselines. Accuracy methods are added in the the function script to evaluate results.

#Figures, tables and data for the contribution of covariate paper are also produced in the script.

#AUTHOR: Benoit Parmentier                                                                      

#PROJECT: Environmental Layers project                                     

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

#### FUNCTION USED IN SCRIPT

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

#### Parameters and constants  

#gwr results:

in_dir4 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gwr_day_lst_comb3_part1_07122013"

#multisampling results (gam)

out_dir<-"/home/parmentier/Data/IPLANT_project/paper_analyses_tables_fig_08032013"

setwd(out_dir)

met_stations_outfiles_obj_file<-"/data/project/layers/commons/data_workflow/output_data_365d_gam_fus_lst_test_run_07172013/met_stations_outfiles_obj_gam_fusion__365d_gam_fus_lst_test_run_07172013.RData"

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

y_var_name <- "dailyTmax"

#method_interpolation <- "gam_daily"

covar_obj_file_1 <- "covar_obj__365d_gam_day_lst_comb3_08132013.RData"

raster_obj_file_3 <- "raster_prediction_obj_kriging_daily_dailyTmax_365d_kriging_day_lst_comb3_07112013.RData"

raster_obj_file_4 <- "raster_prediction_obj_gwr_daily_dailyTmax_365d_gwr_day_lst_comb3_part1_07122013.RData"

#multisampling using baseline lat,lon + elev

#Load objects containing training, testing, models objects 

raster_prediction_obj_6 <-load_obj(file.path(in_dir6,raster_obj_file_6)) #kriging daily multisampling 10 to 70% 

raster_prediction_obj_7 <-load_obj(file.path(in_dir7,raster_obj_file_7)) #gwr daily multisampling 10 to 70%

############### BEGIN SCRIPT #################

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

l1$mae_tb #contains

list_dist_obj <-list(l1,l3,l4)

col_t <- c("red","blue","black")

pch_t <- 1:length(col_t)

legend_text <- c("GAM","Kriging","GWR")

prop_obj_gwr<-calc_stat_prop_tb(names_mod,raster_prediction_obj_7)

list_prop_obj<-list(prop_obj_gam,prop_obj_kriging,prop_obj_gwr)

## plot setting for figure 4

png(filename=file.path(out_dir,png_file_name),

    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

par(mfrow=c(row_mfrow,col_mfrow))

metric_name<-"mae"

list_param_plot<-list(list_prop_obj,col_t,pch_t,legend_text,mod_name,metric_name,add_CI,CI_bar)

names(list_param_plot)<-c("list_prop_obj","col_t","pch_t","legend_text","mod_name","metric_name","add_CI","CI_bar")

rownames(tb4_s)<-NULL #remove row names

tb4_s$method_interp <- "gwr_daily" #add type of interpolation...out_prefix too??

#tb1_s <-raster_prediction_obj_1$tb_diagnostic_s  #gam dailycontains the accuracy metrics for each run...

#tb3_s <-raster_prediction_obj_3$tb_diagnostic_s  #Kriging daily methods

#tb4_s <-raster_prediction_obj_4$tb_diagnostic_s  #gwr daily methods

tb3 <-raster_prediction_obj_3$tb_diagnostic_v  #Kriging daily methods

tb4 <-raster_prediction_obj_4$tb_diagnostic_v  #gwr daily methods

#Calculate difference in MAE and RMSE for training and testing data: call diff_df function

diff_tb1 <-diff_df(tb1_s[tb1_s$pred_mod=="mod1",],tb1[tb1$pred_mod=="mod1",],c("mae","rmse")) #gam select differences for mod1

diff_tb3 <-diff_df(tb3_s[tb3_s$pred_mod=="mod1",],tb3[tb3$pred_mod=="mod1",],c("mae","rmse")) #kriging

diff_tb4 <-diff_df(tb4_s[tb4_s$pred_mod=="mod1",],tb4[tb4$pred_mod=="mod1",],c("mae","rmse")) #gwr

diff_mae_data <-data.frame(gam=diff_tb1$mae,kriging=diff_tb3$mae,gwr=diff_tb4$mae)

diff_rmse_data <-data.frame(gam=diff_tb1$rmse,kriging=diff_tb3$rmse,gwr=diff_tb4$rmse)

#Test the plot

boxplot(diff_mae_data)

boxplot(diff_rmse_data) #plot differences in training and testing accuracies for three methods

      xlab="Interpolation method",

      ylab="RMSE (C)")

tb5 <-raster_prediction_obj_5$tb_diagnostic_v  #gam dailycontains the accuracy metrics for each run...

tb6 <-raster_prediction_obj_6$tb_diagnostic_v  #Kriging daily methods

tb7 <-raster_prediction_obj_7$tb_diagnostic_v  #gwr daily methods

l_pch_t <- c(1,1,3,3,2,2)

l_lty_t <- c(2,1,1,2,2,1)

add_CI <- c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE)

y_range<-c(0.5,3)

plot_ac_holdout_prop(l_prop,l_col_t,l_pch_t,add_CI,y_range)

      ylab="RMSE (C)")

title(main="Difference between training and testing MAE",

      xlab="Interpolation method",

      ylab="MAE (C)")

dev.off()

############### STUDY TIME AND accuracy

mae_tmp<- data.frame(gam=tb1[tb1$pred_mod=="mod1",c("mae")],

                     kriging=tb3[tb3$pred_mod=="mod1",c("mae")],

# interpolation code.

#Figures and data for the contribution of covariate paper are also produced.

#AUTHOR: Benoit Parmentier                                                                      #

#Version: 2

#PROJECT: Environmental Layers project                                       #

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

    n_tb<-subset(l$n_tb,pred_mod==mod_name,select=metric_name) 

    sd_tb<-subset(l$sd_tb,pred_mod==mod_name,select=metric_name) #l$sd_abs_tb[,metric_name]

    no <- unlist(as.data.frame(n_tb))

    y <- unlist(as.data.frame(avg_tb))

    if(add_CI[i]==TRUE){

      if (i==1){

        lines(y~x, pch=pch_t[i],col=col_t[i],type="b")

      }else{

        lines(y~x, pch=pch_t[i],col=col_t[i],type="b")

        plotCI(y=y, x=x, uiw=ciw, col=col_t[i], barcol="blue", lwd=1,pch=pch_t[i],

      }

    }else{