#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                                                                      

#Version: 5

#PROJECT: Environmental Layers project                                     

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

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"

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

#setwd(out_dir)

interp_area_WGS84 <-spTransform(interp_area,CRS_locs_WGS84)         # Project from WGS84 to new coord. system

usa_map <- getData('GADM', country='USA', level=1) #Get US map

usa_map_2 <- usa_map[usa_map$NAME_1!="Alaska",] #remove Alaska

usa_map_2 <- usa_map_2[usa_map_2$NAME_1!="Hawaii",] #remove Hawai 

usa_map_OR <- usa_map_2[usa_map_2$NAME_1=="Oregon",] #get OR

#set up the output file to plot

res_pix<-960

col_mfrow<-1

row_mfrow<-1

png(filename=paste("Figure1_contribution_covariates_study_area_",out_prefix,".png",sep=""),

#opar <- par(fig=c(0.85,0.95,0.8, 0.9), new=TRUE)

#opar <- par(fig=c(0.8,0.95,0.75, 0.9), new=TRUE)

#opar <- par(fig=c(0.75,0.95,0.7, 0.9), new=TRUE)

plot(usa_map_2,border="black") #border and lwd are options of graphics package polygon object

plot(usa_map_OR,col="dark grey",add=T)

#temp.colors <- matlab.like(no_brks)

temp.colors <- colorRampPalette(c('blue', 'khaki', 'red'))

par(mfrow=c(1,2))

plot(lst_md,col=temp.colors(25),axes=F) #use axes=F to remove lat and lon or coordinates

plot(interp_area,add=TRUE)

title("(b) Mean for month of January")

dev.off()

par(mfrow=c(1,2))

p <-levelplot(lst_md,col.regions=temp.colors(25),axes=F,margin=F,scales = list(draw = FALSE), #use list(draw=F) so as  not to dispaly coordinates!!!

y_lab_text <- "MAE (°C)"

add_CI <- c(TRUE,TRUE,TRUE)

#Now set up plotting device

col_mfrow<-2

row_mfrow<-1

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

##### plot figure 5 for paper

####

col_mfrow<-2

row_mfrow<-1

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

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

#### plot figure 6

#set up the output file to plot

col_mfrow<-2

row_mfrow<-1

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

y_range<-range(unlist(month_data_list))

#Now plot figure 14

col_mfrow<-2

#row_mfrow<-2

row_mfrow<-1

layout_m<-c(4,3) #one row two columns

png(paste("Figure_7a_spatial_pattern_tmax_prediction_models_baseline1_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

p<- levelplot(pred_temp_s,main="Interpolated Surfaces Model Comparison baseline 1", ylab=NULL,xlab=NULL,

          par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m,

layout_m<-c(4,3) #one row two columns

png(paste("Figure_7b_spatial_correlogram_tmax_prediction_models_baseline1_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

p<-xyplot(data ~ lag | which, dd,type="b",

          par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=layout_m,

layout_m<-c(1,2) #one row two columns

png(paste("Figure_8a_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

#X11(height=7*layout_m[1],width=7*layout_m[2])

names_layers <-c("mod1 = s(lat,long)+s(elev)","mod4 = s(lat,long)+s(LST)")

#col.regions=temp.colors(25))

dev.off()

diff<-raster(lf1$mod1)-raster(lf1$mod4)

names_layers <- c("mod1-mod4")

diff<-stack(diff)

png(paste("Figure_8b_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

#plot(diff,col=temp.colors(100),main=names_layers)

          margin=F,

          par.settings = list(axis.text = list(font = 2, cex = 1.3),layout=c(1,1),

                              par.main.text=list(font=2,cex=2),strip.background=list(col="white")),par.strip.text=list(font=2,cex=1.5),

          names.attr=names_layers,col.regions=temp.colors)

dev.off()

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

names(tb_mod4_LST_rec3) <- c("month","p_val_rec3")

tb_mod4_LST_rec3$month <- 1:12

layout_m <- c(1,3) # works if set to this?? ok set the resolution...      

#date_selected <- c("20100101","20100901")

png(paste("Figure_15_paper_","LST_elev_",out_prefix,".png", sep=""),

layout_m <- c(2,2) # works if set to this?? ok set the resolution...      

#date_selected <- c("20100101","20100901")

png(paste("Figure10_paper_","_variogram_",out_prefix,".png", sep=""),

    #height=480*6,width=480*4)

p_hist <-histogram(tb_variogram$model,

layout_m<-c(4,3) # works if set to this?? ok set the resolution...

png(paste("Figure_11_paper_","_residuals_",date_selected,"_",out_prefix,".png", sep=""),

grid.arrange(list_p[[1]],list_p[[2]],list_p[[3]],

             list_p[[4]],list_p[[5]],list_p[[6]],

layout_m<-c(1,3) # works if set to this?? ok set the resolution...

png(paste("Figure_11a_paper_","_residuals_",date_selected,"_",out_prefix,".png", sep=""),

grid.arrange(list_p[[1]],list_p[[2]],list_p[[3]],

             ncol=3)  

dev.off() 

layout_m<-c(1,3) # works if set to this?? ok set the resolution...

png(paste("Figure_11b_paper_","_residuals_",date_selected,"_",out_prefix,".png", sep=""),

grid.arrange(list_p[[4]],list_p[[5]],list_p[[6]],

             ncol=3)  

dev.off() 

layout_m<-c(1,3) # works if set to this?? ok set the resolution...

png(paste("Figure_11c_paper_","_residuals_",date_selected,"_",out_prefix,".png", sep=""),

grid.arrange(list_p[[7]],list_p[[8]],list_p[[9]],

             ncol=3)  

dev.off() 

layout_m<-c(1,3) # works if set to this?? ok set the resolution...

png(paste("Figure_11d_paper_","_residuals_",date_selected,"_",out_prefix,".png", sep=""),

grid.arrange(list_p[[10]],

             ncol=3)  

layout_m<-c(1,3) # works if set to this?? ok set the resolution...

png(paste("Figure12_paper_","average_MAE_","_",out_prefix,".png", sep=""),

grid.arrange(list_p_mae[[1]],list_p_mae[[2]],list_p_mae[[3]],ncol=3)

dev.off()   

layout_m <- c(1,3) # works if set to this?? ok set the resolution...      

png(paste("Figure_13_paper_","residuals_MAE_",out_prefix,".png", sep=""),

    #height=480*6,width=480*4)

p_bw1<-bwplot(data_v_ag$res_mod1~elev_rcstat,do.out=F,ylim=c(-15,15),

#max_val<- -10

#min_val <- 0

layout_m<-c(1,3) #one row, three columns

png(paste("Figure_0a_graphic_abstact_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

    height=res_pix*layout_m[1],width=res_pix*layout_m[2])

dd$lag <- mydata$lag 

layout_m<-c(1,3) #one row three columns

png(paste("Figure_0b_graphic_abstract_spatial_correlogram_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

    height=res_pix*layout_m[1],width=res_pix*layout_m[2])