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

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

plot(lst_mm_01,col=temp.colors(25),axes=F)

plot(interp_area,add=TRUE)

dev.off()

## Calucate the proprotion of missing pixel for January 1 mean climatotology image

plot_dst_MAE(list_param_plot)

metric_name <-"mae_tb"

y_lab_text <- "MAE (°C)"

add_CI <- c(TRUE,TRUE,TRUE)

#Now set up plotting device

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

        ylab="Number of observations",

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

box()

dev.off()

#debug(plot_prop_metrics)

plot_prop_metrics(list_param_plot)

      xlab="Hold out validation/testing proportion",

      ylab="MAE (°C)")

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

plot_prop_metrics(list_param_plot)

      xlab="Hold out validation/testing proportion",

      ylab="RMSE (°C)")

legend("top",legend=legend_text_data, 

       cex=0.9, lty=c(1,2),bty="n")

      xlab="Hold out validation/testing proportion",

      ylab="RMSE (°C)")

boxplot(diff_mae_data_mult[-4]) #plot differences in training and testing accuracies for three methods

      xlab="Interpolation method",

      ylab="MAE (°C)")

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)

legend("topleft",legend=legend_text, 

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

#Second plot

ylab_text<-"MAE (°C)"

xlab_text<-"Month"

dev.off()

#Now generate table 5

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

          strip=strip.custom(factor.levels=names_layers),

          xlab=list(label="Spatial lag neighbor", cex=2,font=2),

print(p)

dev.off()

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

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

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

          names.attr=names_layers,col.regions=temp.colors,at=seq(max_val,min_val,by=0.01))

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

print(p)

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

dev.off()

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

names(diff) <- names_layers

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

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

########## Prepare table 6

plot(tb_sig_p_val_rec2)

#Now prepare 

s_table_LST_mod4 <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(LST)")

plot(tb_mod4_LST_rec3,type="l",ylim=c(0.2,1))

s_table_elev_mod4 <- subset(s.table_term_tb,mod_name=="mod4" & term_name=="s(elev_s)")

#tb_mod4_elev_rec3 <- aggregate(tb_mod4_elev_rec3$p_val_rec2 ~ tb_mod4_elev_rec3$month,FUN=mean)

plot(tb_mod4_elev_rec3)

lines(tb_mod4_elev_rec3)

test1 <- subset(s.table_term_tb,mod_name=="mod1" & term_name=="s(elev_s)")

          type="b",

          ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5),

          xlab=list(label="Month",cex=1.5),

p_prop <- update(p_prop,par.settings = list(axis.text = list(font = 2, cex = 1.3),

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

          type="b",

          ylab=list(label="\u0394RMSE between mod1 and mod4",cex=1.5),

          xlab=list(label="Month",cex=1.5),

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

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

p_dif <- update(p_dif, panel = function(...) {

          ylim=c(-1,1),

          ylab=list(label="Pearson Correlation",cex=1.5),

          xlab=list(label="Month",cex=1.5),

          auto.key = list("topright", corner = c(0,1),# col=c("black","red"),

                     border = FALSE, lines = TRUE,cex=1.2))

p_cor <- update(p_cor, panel = function(...) {

p1<-plot(raster_prediction_obj_3$method_mod_obj[[1]]$mod[[1]]$exp_var,raster_prediction_obj_3$method_mod_obj[[1]]$mod[[1]]$var_model,

         ylim=c(0,9),

         ylab=list(label="Semivariance",cex=1.5),

         par.settings = list(axis.text = list(font = 2, cex = 1.3), #control the font size!!

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

                par.strip.text=list(font=2,cex=1.5)

p241<-plot(raster_prediction_obj_3$method_mod_obj[[241]]$mod[[1]]$exp_var,raster_prediction_obj_3$method_mod_obj[[241]]$mod[[1]]$var_model,

         ylim=c(0,9),

         ylab=list(label="Semivariance",cex=1.5),

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

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

                par.strip.text=list(font=2,cex=1.5)

          col=c("grey"),

          ylab=list(label="Percent of total",cex=1.5),

          xlab=list(label="Variogram model",cex=1.5),

p_bw1<- bwplot(tb_variogram$range~tb_variogram$month,do.out=F,ylim=c(0,250000),

         ylab=list(label="Range of variograms",cex=1.5),

         xlab=list(label="Month",cex=1.5),

p_bw2<-bwplot(tb_variogram$Nug~tb_variogram$month,do.out=F,ylim=c(0,12),

         ylab=list(label="Nugget of variograms",cex=1.5),

         xlab=list(label="Month",cex=1.5),

p_bw3<-bwplot(tb_variogram$psill~tb_variogram$month,do.out=F,ylim=c(0,30),

         ylab=list(label="Sill of variograms",cex=1.5),

         xlab=list(label="Month",cex=1.5),

#grid.arrange(p1,p2,p3,ncol=1)

grid.arrange(p_hist,p_bw1,p_bw2,p_bw3,ncol=2)

day<-as.integer(strftime(date_proc, "%d"))

year<-as.integer(strftime(date_proc, "%Y"))

datelabel=format(ISOdate(year,mo,day),"%b %d, %Y")

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

list_p <- vector("list", length(names_mod))

for (k in 1:length(names_mod)){

grid.arrange(list_p[[10]],

             ncol=3)  

dev.off() 

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

### Figure 12: map of MAE by stations over 365 days to summarize residuals information

elev <- subset(s_raster,"elev_s")

list_p_mae <- vector("list", 3)

names_var <- c("mod1","mod2","mod5")

for (k in 1:length(names_var)){

  model_name <- names_var[k]

  res_model_name <- paste("res",model_name,sep="_")

  #res_model_name <- "res_mod1"

  p1 <- levelplot(elev,scales = list(draw = FALSE), colorkey = FALSE,par.settings = GrTheme)

  #tt <- na.omit(data_v_mae)

  #df_tmp=subset(data_v_mae,data_v_mae$res_mod1!="NaN")

  df_tmp=subset(data_v_mae,data_v_mae[[res_model_name]]!="NaN")

               col=matlab.like(5),na.rm=TRUE)

  p3 <- p2 + p1 + p2 #to force legend...

  #print(p3)

}

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

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

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

dev.off()   

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

         ylab=list(label="Residuals (deg C)",cex=1.5),

         xlab=list(label="Elevation classes (meter)",cex=1.5),

         scales = list(x = list(at = c(1, 2, 3, 4), #provide tick location and labels

                               labels = c("0-500","500-1000","1000-1500","1500-2000"))),

                       par.settings = list(axis.text = list(font = 2, cex = 1.3), #control the font size!!

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

         ylab=list(label="Residuals (deg C)",cex=1.5),

         xlab=list(label="Elevation classes (meter)",cex=1.5),

         scales = list(x = list(at = c(1, 2, 3, 4), 

                               labels = c("0-500","500-1000","1000-1500","1500-2000"))),

         par.settings = list(axis.text = list(font = 2, cex = 1.3), #control the font size!!

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

         ylab=list(label="Residuals (deg C)",cex=1.5),

         xlab=list(label="Elevation classes (meter)",cex=1.5),

         scales = list(x = list(at = c(1, 2, 3, 4), 

                               labels = c("0-500","500-1000","1000-1500","1500-2000"))),

         par.settings = list(axis.text = list(font = 2, cex = 1.3), #control the font size!!