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"

###Table 3a, baseline 1: s(lat,lon) 

#Change here !!! need  to reorder rows based on  mod first

model_col<-c("Baseline1","Elevation","Northing","Easting","LST","DISTOC","Forest","CANHEIGHT","LST*Forest","LST*CANHEIGHT") # 

df3a<- as.data.frame(sapply(table_data2,FUN=function(x) x-x[1]))

df3a<- round(df3a,digit=3) #roundto three digits teh differences

df3a$Model <-model_col

###Table 3b, baseline 2: s(lat,lon) + s(elev)

model_col<-c("Baseline2","Northness","Eastness","LST","DISTOC","Forest","CANHEIGHT","LST*Forest","LST*CANHEIGHT")

df3b <- as.data.frame(sapply(table_data1,FUN=function(x) x-x[1])) #difference between baseline (line 1) and other models

df3b <- round(df3b,digit=3) #roundto three digits the differences

### Figure 1: Oregon study area

#3 parameters from output

ghcn_dat <- readOGR(dsn=dirname(met_stations_obj$monthly_covar_ghcn_data),

        sub(".shp","",basename(met_stations_obj$monthly_covar_ghcn_data)))

                      title_plot,y_lab_text,add_CI)

names(list_param_plot)<-c("list_dist_obj","col_t","pch_t","legend_text","mod_name","x_tick_labels","metric_name",

                          "title_plot","y_lab_text","add_CI")

plot_dst_MAE(list_param_plot)

title(xlab="Distance to closest fitting station (km)")

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

mod_name<-c("mod1","mod1","mod1")#selected models

add_CI <- c(TRUE,TRUE,TRUE)

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

##### plot figure 4 for paper

par(mfrow=c(row_mfrow,col_mfrow))

metric_name<-"mae"

#debug(plot_prop_metrics)

plot_prop_metrics(list_param_plot)

#now figure 4b

metric_name<-"rmse"

plot_prop_metrics(list_param_plot)

title(main="RMSE for hold out and methods",

      xlab="Hold out validation/testing proportion",

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

pch_t<- 1:length(col_t)

##Make this a function???

y_range<-range(prop_obj_kriging$avg_tb$rmse,prop_obj_kriging_s$avg_tb$rmse)

#y_range<-range(prop_obj_gam$avg_tb$rmse,prop_obj_gam_s$avg_tb$rmse)

plot_ac_holdout_prop<- function(l_prop,l_col_t,l_pch_t,add_CI,y_range){

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

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

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

y_range<-c(0.5,3)

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

legend("topleft",legend=legend_text, 

title(main="Training and testing RMSE for hold out and methods",

      xlab="Hold out validation/testing proportion",

      ylab="RMSE (C)")

  list_mod_name<-list_param$mod_name

  metric_name<-list_param$metric_name

  add_CI <- list_param$add_CI

  for (i in 1:length(list_obj)){

    x<- l$avg_tb$prop

    y_sd <- unlist(as.data.frame(sd_tb)) #sd_tb

    #ciw2   <- qt(0.975, n) * y_sd2 / sqrt(n)

    #plotCI(y=y, x=x, uiw=ciw, col="red", main=paste(" MAE for ",mod_name,sep=""), barcol="blue", lwd=1,

    #       ylab="RMSE (C)", xlab=xlab_text)

  axis(1,at=1:length(stat_tb[,1]),labels=stat_tb[,1])

}

#Calculate the difference between training and testing in two different data.frames. Columns to substract are provided.

diff_df<-function(tb_s,tb_v,list_metric_names){

  tb_diff<-vector("list", length(list_metric_names))