####################################  INTERPOLATION OF TEMPERATURES  #######################################

############################  Script for manuscript analyses,tables and figures #######################################

#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                                                                      

#DATE: 08/15/2013            

#Version: 2

#PROJECT: Environmental Layers project                                     

#library used in the workflow production:

library(mgcv)                                # GAM package by Simon Wood

library(sp)                                  # Spatial pacakge with class definition by Bivand et al.

library(spdep)                               # Spatial pacakge with methods and spatial stat. by Bivand et al.

library(rgdal)                               # GDAL wrapper for R, spatial utilities

library(gstat)                               # Kriging and co-kriging by Pebesma et al.

library(fields)                              # NCAR Spatial Interpolation methods such as kriging, splines

library(raster)                              # Hijmans et al. package for raster processing

library(rasterVis)                           # Raster plotting functions

library(parallel)                            # Parallelization of processes with multiple cores

library(maptools)                            # Tools and functions for sp and other spatial objects e.g. spCbind

library(maps)                                # Tools and data for spatial/geographic objects

library(reshape)                             # Change shape of object, summarize results 

library(plotrix)                             # Additional plotting functions

library(plyr)                                # Various tools including rbind.fill

library(spgwr)                               # GWR method

library(automap)                             # Kriging automatic fitting of variogram using gstat

library(rgeos)                               # Geometric, topologic library of functions

#RPostgreSQL                                 # Interface R and Postgres, not used in this script

#Additional libraries not used in workflow

library(pgirmess)                            # Krusall Wallis test with mulitple options, Kruskalmc {pgirmess}  

#### FUNCTION USED IN SCRIPT

source(file.path(script_path,function_analyses_paper)) #source all functions used in this script.

in_dir3 <-"/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_kriging_day_lst_comb3_07112013"

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

in_dir7<- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gwr_daily_mults10_lst_comb3_08072013"

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

setwd(out_dir)

y_var_name <- "dailyTmax"

method_interpolation <- "gam_daily"

met_obj_file_1 <- "met_stations_outfiles_obj_gam_daily__365d_gam_day_lst_comb3_08132013.RData"

#raster_prediciton object for baseline 1 () s(lat,lon) + s(elev)) and baseline 2 (slat,lon))

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"

raster_obj_file_6 <- "raster_prediction_obj_kriging_daily_dailyTmax_365d_kriging_daily_mults10_lst_comb3_08062013.RData"

raster_obj_file_7 <- "raster_prediction_obj_gwr_daily_dailyTmax_365d_gwr_daily_mults10_lst_comb3_08072013.RData"

#raster_prediction_obj_gam_daily_dailyTmax_365d_gam_daily_mults10_lst_comb3_08082013.RData

infile_reg_outline <- covar_obj$infile_reg_outline

covar_names<- covar_obj$covar_names

#####

raster_prediction_obj_2 <-load_obj(file.path(in_dir2,raster_obj_file_2)) #comb4 (baseline 1)

raster_prediction_obj_3 <-load_obj(file.path(in_dir3,raster_obj_file_3)) #comb3/mod1 baseline 2, kriging

raster_prediction_obj_4 <-load_obj(file.path(in_dir4,raster_obj_file_4)) #comb3/mod1 baseline 2, gwr

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

##### 1) Generate: Table 3. Contribution of covariates using validation accuracy metrics

## This is a table of accuracy compared to baseline by difference 

#Check input covariates and model formula:

raster_prediction_obj_2$method_mod_obj[[1]]$formulas #models run for baseline 1

summary_metrics_v1<-raster_prediction_obj_1$summary_metrics_v

summary_metrics_v2<-raster_prediction_obj_2$summary_metrics_v

tb2 <-raster_prediction_obj_2$tb_diagnostic_v #365 days accuracy table for baseline 1

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

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

df3a$Model <-model_col

names(df3a)<- names_table_col

print(df3a) #show resulting table

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

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

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

df3b$Model <- model_col

names(df3b)<- names_table_col

print(df3b) #Part b of table 3

mod_compk2 <-kruskal.test(tb2$rmse ~ as.factor(tb2$pred_mod))

print(mod_compk1) #not significant

print(mod_compk2) #not significant

mod_compk1 <-kruskalmc(tb1$rmse ~ as.factor(tb1$pred_mod))

mod_compk2 <-kruskalmc(tb2$rmse ~ as.factor(tb2$pred_mod))

print(mod_compk1)

print(mod_compk2)

#Now write out table 3

file_name<-paste("table3b_baseline2_paper","_",out_prefix,".txt",sep="")

write.table(df3b,file=file_name,sep=",")

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

##### 2) Generate: Table 4. Comparison between interpolation methods using validation accuracy metrics

## This is a table of accuracy metric for the optimal model (baseline2) as identified in the previous step 

##Table 4: Interpolation methods comparison

#get sd for kriging, gam and gwr

#tb2 <-raster_prediction_obj_2$tb_diagnostic_v  #GAM baseline 2, loaded

sd1 <- calc_stat_from_raster_prediction_obj(raster_prediction_obj_1,"sd") # see function script

sd2 <- calc_stat_from_raster_prediction_obj(raster_prediction_obj_2,"sd") # standard deviation for baseline 2

sd3 <- calc_stat_from_raster_prediction_obj(raster_prediction_obj_3,"sd") # kriging

sd4 <- calc_stat_from_raster_prediction_obj(raster_prediction_obj_4,"sd") #gwr

table_sd <- do.call(rbind,list(sd1[1,],sd3[1,],sd4[1,])) #table containing the sd for the three mdethods for baseline 2

table_sd <- round(table_sd[,-1],digit=3) #round to three digits the differences

table_sd <- table_sd[,c("mae","rmse","me","r")] #column 5 contains m50, remove it

summary_metrics_v3<-raster_prediction_obj_3$summary_metrics_v  #Kriging methods

summary_metrics_v4<-raster_prediction_obj_4$summary_metrics_v  # GWR method

table_data3 <-summary_metrics_v3$avg[1,c("mae","rmse","me","r")] #first line mod1 (baseline)

table_data4 <-summary_metrics_v4$avg[1,c("mae","rmse","me","r")] #first line mod1 (baseline)

table_data1 <- table_data1[1,]

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

#combined tables with accuracy metrics and their standard deviations

table4_paper <-table_combined_symbol(table_ac,table_sd,"±")

#lapply(lapply(table_ac,FUN=paste,table_sd,FUN=paste,sep="±"),FUN=paste)

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

names_table_col<-c("MAE","RMSE","ME","R","Model")

names(table4_paper)<- names_table_col

write.table(as.data.frame(table4_paper),file=file_name,sep=",")

####### Now create figures #############

#figure 1: study area

#figure 2: methodological worklfow

#figure 3:Figure 3. MAE/RMSE and distance to closest fitting station.

#Figure 4. RMSE and MAE, mulitisampling and hold out for FSS and GAM.

#Figure 5. Overtraining tendency

# Workflow not generated in R

#Assign model's names

limit_val<-seq(0,150, by=10)

l1 <- distance_to_closest_fitting_station(raster_prediction_obj_1,names_mod,dist_classes=limit_val) #GAM method

l3 <- distance_to_closest_fitting_station(raster_prediction_obj_3,names_mod,dist_classes=limit_val) #Kriging method

l4 <- distance_to_closest_fitting_station(raster_prediction_obj_4,names_mod,dist_classes=limit_val) #GWR method

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

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

x_tick_labels <- limit_val<-seq(5,125, by=10)

metric_name <-"rmse_tb"

title_plot <- "RMSE and distance to closest station for baseline 2"

y_lab_text <- "RMSE (C)"

#quick test

list_param_plot<-list(list_dist_obj,col_t,pch_t,legend_text,mod_name,x_tick_labels,metric_name,title_plot,y_lab_text)

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

plot_dst_MAE(list_param_plot)

title_plot <- "MAE and distance to closest fitting station"

y_lab_text <- "MAE (C)"

#Now set up plotting device

res_pix<-480

col_mfrow<-2

row_mfrow<-1

png_file_name<- paste("Figure_3_accuracy_and_distance_to_closest_fitting_station_",out_prefix,".png", sep="")

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

#Figure 3a

list_param_plot<-list(list_dist_obj,col_t,pch_t,legend_text,mod_name,x_tick_labels,metric_name,title_plot,y_lab_text)

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

#Figure 3b: histogram

barplot(l1$n_tb$res_mod1,names.arg=limit_val,

        ylab="Number of observations",

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

title("Number of observation in term of distance bins")

box()

dev.off()

#Using baseline 2: lat,lon and elev

#Use run of 7 hold out proportions, 10 to 70% with 10 random samples and 12 dates...

#Use gam_day method

#Use comb3 i.e. using baseline s(lat,lon)+s(elev)

names_mod<-c("mod1")

prop_obj_kriging<-calc_stat_prop_tb(names_mod,raster_prediction_obj_6)

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

pch_t<- 1:length(col_t)

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

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

####

res_pix<-480

col_mfrow<-2

row_mfrow<-1

png_file_name<- paste("Figure_4_proportion_of_holdout_and_accuracy_",out_prefix,".png", sep="")

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)

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

plot_prop_metrics(list_param_plot)

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

      xlab="Hold out validation/testing proportion",

      ylab="MAE (C)")

#now figure 4b

metric_name<-"rmse"

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

plot_prop_metrics(list_param_plot)

      xlab="Hold out validation/testing proportion",

      ylab="RMSE (C)")

dev.off()

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

#########Figure 5. Overtraining tendency

#read in relevant data:

tb1_s<-extract_from_list_obj(raster_prediction_obj_1$validation_mod_obj,"metrics_s")

rownames(tb1_s)<-NULL #remove row names

tb1_s$method_interp <- "gam_daily" #add type of interpolation...out_prefix too??

tb3_s$method_interp <- "kriging_daily" #add type of interpolation...out_prefix too??

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

tb1 <-raster_prediction_obj_1$tb_diagnostic_v  #gam dailycontains the accuracy metrics for each run...

tb3 <-raster_prediction_obj_3$tb_diagnostic_v  #Kriging daily methods

tb4 <-raster_prediction_obj_4$tb_diagnostic_v  #gwr daily methods

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

title(main="Training and testing RMSE for hold out and interpolation 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

prop_obj_kriging_s<-calc_stat_prop_tb(names_mod,raster_prediction_obj_6,testing=FALSE)

prop_obj_gwr_s<-calc_stat_prop_tb(names_mod,raster_prediction_obj_7,testing=FALSE)

plot(prop_obj_gam_s$avg_tb$rmse ~ prop_obj_gam_s$avg_tb$prop, type="b",)

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

#### plot figure 5

#set up the output file to plot

res_pix<-480

col_mfrow<-2

row_mfrow<-1

png_file_name<- paste("Figure_5_overtraining_tendency_and_holdout_proportion_",out_prefix,".png", sep="")

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

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

pch_t<- 1:length(col_t)

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(prop_obj_gam$avg_tb$rmse ~ prop_obj_gam$avg_tb$prop, ylab="",xlab="",type="b",col=c("red"),pch=pch_t[1],ylim=y_range,lty=2)

lines(prop_obj_gam_s$avg_tb$rmse ~ prop_obj_gam_s$avg_tb$prop, ylab="",xlab="",type="b",pch=pch_t[1],ylim=y_range,col=c("red"))

lines(prop_obj_gwr_s$avg_tb$rmse ~ prop_obj_gwr_s$avg_tb$prop, ylab="",xlab="",type="b",ylim=y_range,pch=pch_t[3],col=c("black"))

lines(prop_obj_gwr$avg_tb$rmse ~ prop_obj_gam$avg_tb$prop, ylab="",xlab="",type="b",ylim=y_range,pch=pch_t[3],,col=c("black"),lty=2)

lines(prop_obj_kriging$avg_tb$rmse ~ prop_obj_kriging$avg_tb$prop,ylab="",xlab="", type="b",ylim=y_range,pch=pch_t[2],,col=c("blue"),lty=2)

lines(prop_obj_kriging_s$avg_tb$rmse ~ prop_obj_kriging_s$avg_tb$prop,ylab="",xlab="",type="b",ylim=y_range,pch=pch_t[2],col=c("blue"))

legend("topleft",legend=legend_text, 

       cex=0.9, pch=pch_t,col=col_t,lty=1,bty="n")

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

      xlab="Hold out validation/testing proportion",

      ylab="RMSE (C)")

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

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

      xlab="Interpolation method",

      ylab="MAE (C)")

#########Figure 6: Monthly RMSE averages for the three interpolation methods: GAM, GWR and Kriging.

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

                     gwr=tb4[tb4$pred_mod=="mod1",c("mae")])

lines(mae_tmp$kriging,col=c("blue"),type="b",pch=2)

lines(mae_tmp$gwr,col=c("black"),type="b",pch=2)

       cex=1.2, pch=pch_t,col=col_t,lty=1,bty="n")

arrange(x2,desc(mae))

#                     gwr=tb4[tb4$pred_mod=="mod1",c("mae")])

tb3_month<-raster_prediction_obj_3$summary_month_metrics_v[[1]]

tb4_month<- raster_prediction_obj_4$summary_month_metrics_v[[1]]

plot(1:12,tb1_month$mae,col=c("red"),type="b",ylim=y_range)

lines(1:12,tb3_month$mae,col=c("blue"),type="b")

lines(1:12,tb4_month$mae,col=c("black"),type="b")

  date<-strptime(tb$date, "%Y%m%d")   # interpolation date being processed

  month<-strftime(date, "%m")          # current month of the date being processed

}

tb1$month<-add_month_tag(tb1)

tb3$month<-add_month_tag(tb3)

tb4$month<-add_month_tag(tb4)

metric_name<-"mae"

month_data_list<-list(gam=tb1[tb1$pred_mod=="mod1",c(metric_name,"month")],

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

                      gwr=tb4[tb4$pred_mod=="mod1",c(metric_name,"month")])

y_range<-range(unlist(month_data_list))

#boxplot(test[[metric_ac]]~test[[c("month")]],outline=FALSE,horizontal=FALSE,cex=0.5,

#        ylab=paste(metric_ac,"in degree C",sep=" "),,axisnames=FALSE,axes=FALSE)

#boxplot(test[[metric_ac]]~test[[c("month")]],outline=FALSE,horizontal=FALSE,cex=0.5,

#        ylab=paste(metric_ac,"in degree C",sep=" "))

#axis(1, labels = FALSE)

## Create some text labels

#labels <- month.abb # abbreviated names for each month

## Plot x axis labels at default tick marks

#text(1:length(labels), par("usr")[3] - 0.25, srt = 45, adj = 1,

#     labels = labels, xpd = TRUE)

#axis(2)

#box()

#Now plot figure 6

res_pix<-480

col_mfrow<-2

row_mfrow<-1

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

xlab_tick <- unique(tb1$month)

xlab_text <-"Month"

plot(1:12,tb1_month$mae,col=c("red"),type="b",ylim=y_range,xlab=xlab_text,xaxt="n")

lines(1:12,tb3_month$mae,col=c("blue"),type="b")

lines(1:12,tb4_month$mae,col=c("black"),type="b")

axis(1,at=1:12,labels=xlab_tick)

title(main="Monthly average MAE")

xlab_text<-"Month"

#y_range<-range(month_data_list$gam$mae)

boxplot(mae~month,data=month_data_list$gam,main="GAM",ylab=ylab_text,outline=FALSE)

#boxplot(mae~month,data=month_data_list$kriging,ylim=y_range,main="Kriging",ylab=ylab_text,outline=FALSE)

#boxplot(mae~month,data=month_data_list$gwr,ylim=y_range,main="GWR",ylab=ylab_text,outline=FALSE)

dev.off()

length(tb1_month$mae)

names(tb1_month)

sd1 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_1,"sd") # see function script

sd2 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_2,"sd") # standard deviation for baseline 2

sd3 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_3,"sd") # kriging

sd4 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_4,"sd") #gwr

avg_v1 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_1,"mean") # see function script

avg_v2 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_2,"mean") # standard deviation for baseline 2

avg_v3 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_3,"mean") # kriging

avg_v4 <- calc_stat_month_from_raster_prediction_obj(raster_prediction_obj_4,"mean") #gwr

#Combined sd in one table for mod1 (baseline 2)

table_sd <- do.call(cbind,list(gam=sd1[sd1$pred_mod=="mod1",c("rmse")],

                               kriging=sd3[sd3$pred_mod=="mod1",c("rmse")],

                               gwr=sd4[sd4$pred_mod=="mod1",c("rmse")])) #table containing the sd for the three mdethods for baseline 2

table_sd <- as.data.frame(round(table_sd,digit=3)) #round to three digits the differences

#Combined mean in one table for mod1 (baseline 2)

table_avg <- do.call(cbind,list(gam=avg_v1[avg_v1$pred_mod=="mod1",c("rmse")],

                               kriging=avg_v3[sd3$pred_mod=="mod1",c("rmse")],

                               gwr=avg_v4[avg_v4$pred_mod=="mod1",c("rmse")])) #table containing the sd for the three mdethods for baseline 2

table_avg <- as.data.frame(round(table_avg,digit=3)) #round to three digits the differences

#combined tables with accuracy metrics and their standard deviations

table5_paper <-table_combined_symbol(table_avg,table_sd,"±")

table5_paper$month <- month.abb

file_name<-paste("table5_comparisons_monthly_averages_interpolation_methods_paper","_",out_prefix,".txt",sep="")

write.table(as.data.frame(table5_paper),file=file_name,sep=",")

y_var_name <-"dailyTmax"

index<-244 #index corresponding to January 1

lf4 <- raster_prediction_obj_4$method_mod_obj[[index]][[y_var_name]]

date_selected <- "20109101"

methods_names <-c("gam","kriging","gwr")

names_layers<-methods_names

lf <-list(lf1$mod1,lf3$mod1,lf4$mod1)

#lf <-lf[[1]]

pred_temp_s <-stack(lf)