#MODIFIED ON: 12/23/2013            

#Version: 2

function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_12232013.R"

## Results from monthly holdout 0 to 70%

#List of object with 0% monthly hold out

list_raster_obj_files  <- list(file.path(in_dir1,raster_obj_file_1),file.path(in_dir2,raster_obj_file_2),

                               file.path(in_dir3a,raster_obj_file_3a),file.path(in_dir3b,raster_obj_file_3b),

                               file.path(in_dir4,raster_obj_file_4),file.path(in_dir5,raster_obj_file_5),

                               file.path(in_dir6,raster_obj_file_6),file.path(in_dir7,raster_obj_file_7),

                               file.path(in_dir8,raster_obj_file_8),file.path(in_dir9,raster_obj_file_9))

names(list_raster_obj_files)<- c("gam_daily","kriging_daily","gwr_daily_a","gwr_daily_b",

                                 "gam_CAI","kriging_CAI","gwr_CAI",

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

y_var_name <- "dailyTmax"

out_prefix<-"analyses_12232013"

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

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

if (!file.exists(out_dir)){

  dir.create(out_dir)

  #} else{

  #  out_path <-paste(out_path..)

}

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

############### BEGIN SCRIPT ##################################

#############################PART I: Generate tables for paper ##########

#Table 1: Covariates used (not produced in R)

#Table 2: Covariates models and functional form (not produced in R)

#Table 3: Combinations of models tested (not produced in R)

#Table 4. Average RMSE for all models acroos year 2010 wiht 30% daily hold out

#Table 5: Average monhtly RMSE for CAI model 7, climatology surfaces,  CAI and single time scale methods 

##### Table 4: Contribution of covariates using validation accuracy metrics

table_gwr <- summary_metrics_v_list[grep("gwr",names(summary_metrics_v_list))] #select models related to gwr

table_gwr$gwr_daily <- c(table_gwr[[1]],table_gwr[[2]]) #combine both gwr object results

table_gwr <- table_gwr[c(5,3,4)] #reorder list and drop first two objects

table4_paper<-table4_paper[(c(5,4,1,2,3))]  #reordering columns                           

################ Table 5 ################

#### Monthly RMSE for GAM CAI Climatology, GAM CAI daily and Single  time scale

## This is a table of accuracy  

#raster_prediction_obj_4 <-load_obj(file.path(in_dir4,raster_obj_file_4)) #comb5 gam_dir

tb_m_gam_CAI <- load_obj(file.path(in_dir13,raster_obj_file_13))$tb_month_diagnostic_v #getting objet

tb_gam_CAI <- load_obj(file.path(in_dir13,raster_obj_file_13))$tb_diagnostic_v #getting objet

tb_gam_daily <- load_obj(file.path(in_dir1,raster_obj_file_1))$tb_diagnostic_v #getting objet

names_mod <- paste("mod",1:7,sep="")

list_stat_tb_gam_CAI <-calc_stat_by_month_tb(names_mod,tb_gam_CAI,month_holdout=T)

list_stat_tb_gam_daily <-calc_stat_by_month_tb(names_mod,tb_gam_daily,month_holdout=F)

list_stat_tb_gam_CAI$avg

clim_rmse <- subset(tb_m_gam_CAI,prop==30 & pred_mod== "mod7",select="rmse")

CAI_rmse <- subset(list_stat_tb_gam_CAI$avg,prop_month==30 & pred_mod== "mod7",select="rmse")

daily_rmse <- subset(list_stat_tb_gam_daily$avg, pred_mod== "mod7",select="rmse")

table5 <- data.frame(clim=clim_rmse,CAI_month=CAI_rmse,daily_month=daily_rmse)

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

table5$month <- month.abb

table5 <- table5[,c(4,1,2,3)]

###Table 5, writeout

names_table_col <-c("Month","Long term Climatology RMSE",

                    "Average Monthly RMSE from Daily pred (CAI)"

                    ,"Average Monthly RMSE from Daily pred (Single time scale)") # 

names(table5)<- names_table_col

print(table5) #show resulting table

#Now write out table 5

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

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

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

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

#Figure 11: Spatial lag profiles for predicted surfaces 

######### Figure 2: LST averaging: daily mean compared to monthly mean

lst_md <- raster(ref_rast_name)

png(filename=paste("Figure2_paper_Comparison_daily_monthly_mean_lst",out_prefix,".png",sep=""),width=960,height=480)

title("Mean for January 1")

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

######### Figure 3: Prediction procedures: direct, CAI and FSS 

#(figure created outside R)

#plot_names <- c("RMSE for gam_FSS","RMSE for kriging_FSS","RMSE for gwr_FSS",

#               "RMSE for gam_CAI","RMSE for kriging_CAI","RMSE for gwr_CAI")

names_mod <- paste("mod",1:7,sep="")

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

               "RMSE for gam_CAI","RMSE for kriging_CAI","RMSE for gwr_CAI")

list_plots <- plot_accuracy_by_holdout_fun(list_tb,ac_metric,plot_names,names_mod)

##For paper...combine figures... tb_v for GWR, Kriging and GAM for both FSS and CAI

png(paste("Figure4_paper_accuracy_",ac_metric,"_prop_month","_",out_prefix,".png", sep=""),

png(paste("Figure5_paper_boxplot_overtraining_",out_prefix,".png", sep=""),

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

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

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

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

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

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

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

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

nb_fig<- c("6a","6b","6c")

png(paste("Figure6_paper_","_spatial_pattern_tmax_prediction_models_gam_levelplot_",date_selected,out_prefix,".png", sep=""),

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

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

list_transect2[[1]]<-c(transect_list[1],paste("Figure8a_paper_tmax_elevation_transect1_OR_",date_selected,

list_transect2[[2]]<-c(transect_list[2],paste("Figure8b_tmax_elevation_transect2_OR_",date_selected,

list_transect2[[3]]<-c(transect_list[3],paste("Figure8c_paper_tmax_elevation_transect3_OR_",date_selected,

list_transect3[[1]]<-c(transect_list[1],paste("Figure8a_tmax_elevation_transect1_OR_",date_selected,

list_transect3[[2]]<-c(transect_list[2],paste("Figure8b_tmax_elevation_transect2_OR_",date_selected,

list_transect3[[3]]<-c(transect_list[3],paste("Figure8c_tmax_elevation_transect3_OR_",date_selected,

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

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

title_plot2<-paste(names(list_transect2),"on",date_selected,sep=" ")

#title_plot2<-paste(names(list_transect2),"on",date_selected,sep=" ")

#title_plot2<-paste(rep("Oregon transect on ",3), date_selected,sep="")

#title_plot3<-paste(names(list_transect3),date_selected,sep=" ")

#title_plot3<-paste(rep("Oregon transect on ",3), date_selected,sep="")

trans_data3 <-plot_transect_m2(list_transect3,rast_pred3,title_plot2,disp=FALSE,m_layers_sc)

#Figure 9: Spatial pattern: Image differencing  

png(paste("Figure9_paper_difference_image_",out_prefix,".png", sep=""),

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

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

png(paste("Figure10_paper_diff_prediction_tmax_difference_land cover,ac_metric","_",out_prefix,".png", sep=""),

  #make it slanted...

  axis(2,cex.axis=1.4,font=2)

        cex=1.4, col=c("black","red","green","blue","darkgreen","purple"),bty="n",

  title(paste(title_plots_list[i],sep=""),cex=1.6, font=2)

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

#### Figure 11: Spatial lag profiles  

#This figure is generated to show the spatial Moran'I for 10 spatial 

#for Jan 1 and Sept 1 in 2010 for all models (1 to 7) and methods

index<-1 #index corresponding to Jan 1 #For now create Moran's I for only one date...

lf_moran_list_date1 <-lapply(list_raster_obj_files[c("gam_daily","gam_CAI","gam_fss")],

                               FUN=function(x){x<-load_obj(x);x$method_mod_obj[[index]][[y_var_name]]})                           

index<-244 #index corresponding to Sept 1 #For now create Moran's I for only one date...

lf_moran_list_date2 <-lapply(list_raster_obj_files[c("gam_daily","gam_CAI","gam_fss")],

                               FUN=function(x){x<-load_obj(x);x$method_mod_obj[[index]][[y_var_name]]})                           

#date_selected <- "20100901"

#date_selected <- "20100101"

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

methods_names <-c("gam_daily","gam_CAI","gam_FSS")

names_layers<-methods_names

#Subset images to eliminate mod_kr

lf1 <- list(lf_moran_list_date1[[1]],lf_moran_list_date1[[2]][1:7],lf_moran_list_date1[[3]][1:7])

lf2 <- list(lf_moran_list_date2[[1]],lf_moran_list_date2[[2]][1:7],lf_moran_list_date2[[3]][1:7])

names(lf1)<-c("gam_daily","gam_CAI","gam_FSS")

names(lf2)<-c("gam_daily","gam_CAI","gam_FSS")

### Now extract Moran's I for a range of lag using a list of images

#set maximum lag range

nb_lag <-10

#Provide list of raster images:

list_lf <- list(lf1,lf2)

list_moran_df1 <- calculate_moranI_profile(list_lf[[1]],nb_lag) #for January 1

list_moran_df2 <- calculate_moranI_profile(list_lf[[2]],nb_lag) #for September 1

names(list_moran_df1)<-c("gam_daily","gam_CAI","gam_FSS")

names(list_moran_df2)<-c("gam_daily","gam_CAI","gam_FSS")

#Run accross two dates...

#list_moran lapply(list_lf,FUN=calculate_moranI_profile,nb_lag=nb_lag)

### Prepare to plot lag Moran's I profiles

#generate automatic title for exploration if necessary!!

list_title_plot<- list(c("Spatial lag profile on January 1, 2010"),

                  c("Spatial lag profile on September 1, 2010"))

#name used in the panel!!!

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

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

list_moran_df <- list(list_moran_df1,list_moran_df2)

list_param_plot_moranI_profile_fun <- list(list_moran_df,list_title_plot,names_panel_plot,layout_m)

names(list_param_plot_moranI_profile_fun) <- c("list_moran_df","list_title_plot","names_panel_plot","layout_m")

#debug(plot_moranI_profile_fun)

#p<- plot_moranI_profile_fun(1,list_param=list_param_plot_moranI_profile_fun)

list_moranI_plots <- lapply(1:2,FUN=plot_moranI_profile_fun,list_param=list_param_plot_moranI_profile_fun)

#This function uses list moran_df object from calculate_moranI_profile function!!

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

png(paste("Figure11_paper_spatial_correlogram_prediction_models_levelplot_",out_prefix,".png", sep=""),

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

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

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

#png(paste("Figure11_paper_spatial_correlogram_prediction_models_levelplot_",out_prefix,".png", sep=""),

#    height=480,width=480)

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

p1 <- list_moranI_plots[[1]]

p2 <- list_moranI_plots[[2]]

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

dev.off()

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

#Figure 12: Monthly climatology, Daily deviation and bias

lf_delta_gam_fss <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_fss"]])$method_mod_obj,"delta") #getting objet

test01 <-stack(lf_delta_gam_fss[[287]]) #Ot.14

test0 <-stack(lf_delta_gam_fss[[288]]) #Ot.15

test1 <-stack(lf_delta_gam_fss[[289]]) #Ot.16

test2<-stack(lf_delta_gam_fss[[290]]) #Ot.17

plot(test01)

plot(test0)

plot(test1)

plot(test2)

#Figure 13: Tmax and LST averagees

y_range_nobs <- c(0,385)

y_range_avg <- c(-15,60)

var_name <- "LST"

#It works for any variable with stack of monthly values...!!

LST_stat <- plot_mean_nobs_r_stack_by_month(var_s=LST_s,var_nobs=LST_nobs,y_range_nobs,y_range_avg,var_name,out_prefix)

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

####### Figure 13 : disussion  

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

####### Figure 13a: spatial variation -MoranI and  std  

#Use data from accuraccy with no monthly hold out

#Extract monthly data frame used in fitting 

list_data_month_s <-extract_list_from_list_obj(raster_obj$validation_mod_month_obj,"data_s") 

year_nbs <- sapply(list_data_month_s,FUN=length) #number of observations per month

#Convert sp data.frame and combined them in one unique df, see function define earlier

data_month_all <-convert_spdf_to_df_from_list(list_data_month_s) #long rownames

#LSTD_bias_avgm<-aggregate(LSTD_bias~month,data=data_month_all,mean)

#LSTD_bias_sdm<-aggregate(LSTD_bias~month,data=data_month_all,sd)

LST_avgm<-aggregate(LST~month,data=data_month_all,mean)

TMax_avgm<-aggregate(TMax~month,data=data_month_all,mean)

#plot(LST_avgm,type="b")

#lines(TMax_avgm,col="blue",add=T)

plot(data_month_all$month,data_month_all$LST, xlab="month",ylab="LST at station")

title(paste("Monthly LST at stations in Oregon", "2010",sep=" "))

png(paste("LST_at_stations_per_month_",out_prefix,".png", sep=""), type="png")

statistics_LST_s<- LST_stat$avg #extracted earlier!!!

png(paste("Monthly_mean_TMax_LST_at_Stations_",out_prefix,".png", sep=""), type="png")

plot(TMax_avgm,type="b",ylim=c(0,35),col="red",xlab="month",ylab="tmax (degree C)")

lines(1:12,LST_avgm$LST,type="b",col="blue")

#lines(1:12,statistics_LST_s$mean,type="b",col="darkgreen")

text(TMax_avgm[,1],TMax_avgm[,2],rownames(statistics_LST_s),cex=0.8,pos=2)

legend("topleft",legend=c("TMax","LST"), cex=1, col=c("red","blue"),

              lty=1,bty="n")

title(paste("Monthly mean tmax and LST at stations in Oregon", "2010",sep=" "))           

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

####### Figure 13b: spatial variation -MoranI and  std  

#LST for area with FOrest 50> and grass >50

#Account for forest

data_month_all_forest<-data_month_all[data_month_all$LC1>=50,]

data_month_all_grass<-data_month_all[data_month_all$LC6>=50,]

#data_month_all_grass<-data_month_all[data_month_all$LC6>=10,]

data_month_all_urban<-data_month_all[data_month_all$LC9>=50,]

LST_avgm_forest<-aggregate(LST~month,data=data_month_all_forest,mean)

LST_avgm_grass<-aggregate(LST~month,data=data_month_all_grass,mean)

LST_avgm_urban<-aggregate(LST~month,data=data_month_all_urban,mean)

plot(TMax_avgm,type="b",ylim=c(-7,42))

lines(LST_avgm$LST,col="blue",add=T)

lines(LST_avgm_forest,col="green",add=T)

lines(LST_avgm_grass,col="red",add=T)

lines(LST_avgm_urban,col="pink",add=T)

legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass","LST_urban"), cex=0.8, col=c("black","blue","green","red","pink"),

       lty=1,bty="n")

title(paste("Monthly average tmax for stations in Oregon ", "2010",sep=""))

savePlot(paste("Temporal_profile_res",id[i],out_prefix,".png", sep=""), type="png")  

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

####### Figure 13c: spatial variation -MoranI and  std  

#Takes 1 hour to get the average moran's I for the whole year so load moran_std_tt_fss2.RData

#tt_fss <- mclapply(1:365, list_param=list_param_stat_moran, FUN=stat_moran_std_raster_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

#save(tt_fss,file=file.path(out_dir,"moran_std_tt_fss.RData"))

#tx<- do.call(rbind,tt_fss)

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

#mo<-as.integer(strftime(dates, "%m"))          # current month of the date being processed

dates<-unique(raster_obj$tb_diagnostic_v$date)

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

####### Figure 13d: correlation between elevation and LST

LST_s <- subset(s_raster,c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12"))

LST1 <- subset(LST_s,1)

test<-stack(LST_s,elev)

names(test) <- c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12","elev")

x_cor <-layerStats(test,stat="pearson",na.rm=T)

corr(values(LST1),values(elev))

plot(1:12,x_cor[[1]][1:12,13],type="b",ylab="corelation",xlab="month",main="correlation between elevation and LST")

#### Now elev?

#LC1<-mask(LC1,mask_ELEV_SRTM)

#  cellStats(LC1,"countNA")        #Check that NA have been assigned to water and areas below 0 m

#LC1_50_m<- LC1>50

#LC1_100_m<- LC1>=100

#LC1_50_m[LC1_50_m==0]<-NA

#LC1_100_m[LC1_100_m==0]<-NA

#LC1_50<-LC1_50_m*LC1

#LC1_100<-LC1_100_m*LC1

#avl<-c(0,500,1,500,1000,2,1000,1500,3,1500,2000,4,2000,4000,5)

#rclmat<-matrix(avl,ncol=3,byrow=TRUE)

#elev_rec<-reclass(ELEV_SRTM,rclmat)  #Loss of layer names when using reclass

#elev_rec_forest<-elev_rec*LC1_100_m

#avg_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="mean",na.rm=TRUE)

#std_elev_rec<-zonal(rast_diff,zones=elev_rec,stat="sd",na.rm=TRUE)

#avg_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="mean",na.rm=TRUE)

#std_elev_rec_forest<-zonal(rast_diff,zones=elev_rec_forest,stat="sd",na.rm=TRUE)

## CREATE plots

#X11()

#plot(avg_elev_rec[,1],avg_elev_rec[,2],type="b",ylim=c(-10,1),

#       ylab="",xlab="",axes=FALSE)

#mtext("tmax difference between FSS and CAI (degree C)",side=2,cex=1.2,line=3,font=2)

#mtext("elevation classes (m)",side=1,cex=1.2,line=3,font=2)

#lines(avg_elev_rec_forest[,1],avg_elev_rec_forest[,2],col="green",type="b") #Elevation and 100% forest...

#breaks_lab<-avg_elev_rec[,1]

#elev_lab<-c("0-500","500-1000","1000-1500","1500-2000","2000-4000")

#axis(side=1,las=1,

#       at=breaks_lab,labels=elev_lab, cex=1.5,font=2) #reduce number of labels to Jan and June

#axis(2,cex.axis=1.2,font=2)

#legend("bottomleft",legend=c("Elevation", "elev_forest"), 

#         cex=1, lwd=1.3,col=c("black","green"),bty="n",

#         lty=1)

#box()

#title(paste("Prediction tmax difference (",mf_selected,"-",mc_selected,") and elevation ",sep=""),cex=1.4,font=2)

#savePlot(paste("fig7_diff_prediction_tmax_difference_elevation",mf_selected,mc_selected,date_selected,out_prefix,".png", sep="_"), type="png")

#dev.off()