#Analyses, figures, tables and data for the  paper are also produced in the script.

#AUTHOR: Benoit Parmentier 

#CREATED ON: 10/31/2013  

#Version: 3

#PROJECT: Environmental Layers project                                     

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

                                 "gam_fss","kriging_fss","gwr_fss")

y_var_name <- "dailyTmax"

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

out_dir <-paste(out_dir,"_",out_prefix,sep="")

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

####### PART 2: generate figures for paper #############

#figure 1: Study area OR

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

######### Figure 1: Oregon study area, add labeling names to  Willamette Valley an Mountain Ranges?

dev.off()

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

lst_md <- raster(ref_rast_name)

projection(lst_md)<-projection(s_raster)

lst_md<- lst_md - 273.16

lst_mm_01<-subset(s_raster,"mm_01")

par(mfrow=c(1,2))

plot(lst_md)

plot(interp_area,add=TRUE)

dev.off()

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

raster_prediction_obj_9 <-load_obj(file.path(in_dir9,raster_obj_file_9)) #comb5 gwr_fss

plot_names <- names(list_tb) #this is the default names for the plots

#now replace names for relevant figure used later on.

plot_names[7:12] <- c("RMSE for gam_FSS","RMSE for kriging_FSS","RMSE for gwr_FSS",

#Quick function to explore accuracy make this a function to create solo figure...and run only a subset...

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

dev.off()

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

#x<-tb_ms_list[[1]]

tb_ms_list_diff <- lapply(tb_ms_list,FUN=function(x){x[x$prop!=0,]})                           

## Now create boxplots...

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

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

#boxplot(diff_kriging_CAI$rmse,diff_gam_CAI$rmse,diff_gwr_CAI$rmse,names=c("kriging_CAI","gam_CAI","gwr_CAI"),

#        main="Difference between training and testing daily rmse")

#monthly CAI

boxplot(list_m_diff$kriging_CAI$rmse,list_m_diff$gam_CAI$rmse,list_m_diff$gwr_CAI$rmse,

        names=c("kriging_CAI","gam_CAI","gwr_CAI"),ylab="RMSE (°C)",xlab="Interpolation Method",

        main="Difference between training and testing for monhtly rmse")

legend("topleft",legend=c("a"),bty="n") #bty="n", don't put box around legend

#daily CAI

boxplot(list_diff$kriging_CAI$rmse,list_diff$gam_CAI$rmse,list_diff$gwr_CAI$rmse,

        names=c("kriging_CAI","gam_CAI","gwr_CAI"),ylab="RMSE (°C)",xlab="Interpolation Method",

        main="Difference between training and testing daily rmse")

legend("topleft",legend=c("b"),bty="n")

#daily fss

boxplot(list_diff$kriging_fss$rmse,list_diff$gam_fss$rmse,list_diff$gwr_fss$rmse,

        names=c("kriging_FSS","gam_FSS","gwr_FSS"),ylab="RMSE (°C)",xlab="Interpolation Method",

        main="Difference between training and testing daily rmse")

legend("topleft",legend=c("d"),bty="n")

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

y_var_name <-"dailyTmax"

index<-244 #index corresponding to Sept 1

names_layers <-c("mod1 = lat*long","mod2 = lat*long + LST","mod3 = lat*long + elev","mod4 = lat*long + N_w*E_w",

                 "mod5 = lat*long + elev + DISTOC","mod6 = lat*long + elev + LST","mod7 = lat*long + elev + LST*FOREST")

list_plots_spt <- vector("list",length=length(lf))

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

  pred_temp_s <-stack(lf[[i]])

  max_val <- 40

  min_val <- -10

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

  png(paste("Figure_",nb_fig[i],"_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

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

  p <- levelplot(pred_temp_s,main=methods_names[i], ylab=NULL,xlab=NULL,

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

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

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

  print(p)

  dev.off()

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

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

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

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

dev.off()

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

## Transects image location in OR  

layout_m<-c(1,1)

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

vx_text <- c(395770.1,395770.1,395770.1) #location of labels

vy_text  <-c(473000,297045.9,112611.5)

title("Transect Oregon")

dev.off()

nb_transect <- 3

list_transect2<-vector("list",nb_transect)

list_transect3<-vector("list",nb_transect)

layers_names3 <- names(rast_pred3)<-c("mod1","mod3","mod6","elev")

#pos<-c(1,2) # postions in the layer prection

#transect_list

                                           paste("mod1_2_6",collapse="_"),out_prefix,sep="_"))

                                           paste("mod1_2_6",collapse="_"),out_prefix,sep="_"))

                                           paste("mod1_2_6",collapse="_"),out_prefix,sep="_"))

                                           paste("mod1_3_6",collapse="_"),out_prefix,sep="_"))

                                           paste("mod1_3_6",collapse="_"),out_prefix,sep="_"))

                                           paste("mod1_3_6",collapse="_"),out_prefix,sep="_"))

names(list_transect2)<-c("Oregon Transect 1","Oregon Transect 2","Oregon Transect 3")

#title_plot2

#rast_pred2

#debug(plot_transect_m2)

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

#Do for january and September...?

#names_layers <-c("mod1 = lat*long","mod2 = lat*long + LST","mod3 = lat*long + elev","mod4 = lat*long + N_w*E_w",

r_stack_diff <-stack(c(diff_pred_date1_list,diff_pred_date2_list))

names(r_stack_diff) <- c("Jan_Daily","Jan_CAI","Jan_FSS","Sept_Daily","Sept_CAI","Sept_FSS")

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

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

dev.off()

###

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

#Figure 10: LST and Tmax at stations, elevation and land cover.