#MODIFIED ON: 03/06/2014            

function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_12252013.R"

out_prefix<-"analyses_03062013"

#Figure annex: LST climatology production: daily mean compared to monthly mean: moved to appendix

#figure 1: Study area OR

#Figure 2. Accuracy and monthly hold out for FSS and CAI procedures and GWR, Kriging and GAM methods.

#Figure 5: Transect location (transect map)

#Figure 6: Transect profiles for one day 

#Figure 7: Image differencing and land cover: spatial patterns   

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

#Figure 9: Spatial lag profiles for predicted surfaces 

#Figure 10: Daily deviation

#Figure 11: Spatial correlogram at stations, LST and elevation every 5 lags

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

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

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

#This is now moved in Annex

lst_md <- raster(ref_rast_name)

projection(lst_md)<-projection(s_raster)

lst_mm_09<-subset(s_raster,"mm_09")

lst_md<-raster(ref_rast_d001)

lst_md<- lst_md - 273.16

lst_mm_01<-subset(s_raster,"mm_01")

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

par(mfrow=c(1,2))

plot(lst_md)

plot(interp_area,add=TRUE)

title("Mean for January 1")

plot(lst_mm_01)

plot(interp_area,add=TRUE)

title("Mean for month of January")

dev.off()

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

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

            index.cond=list(c(1,5,3,2,6,4)),    #this provides the order of the panels)

            pch=1:length(avg_tb$pred_mod),

            pch=1:length(avg_tb$pred_mod))),

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

        ylim=c(-9,1), outline=TRUE,

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

        ylim=c(-9,1),outline=TRUE,

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

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

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

        ylim=c(-12,1.5),outline=TRUE,

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

        ylim=c(-12,1.5),outline=TRUE,

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

names_layers<-methods_names

lf<-list(lf_list[[1]],lf_list[[2]][1:7],lf_list[[3]][1:7])

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

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

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

p3 <- list_plots_spt[[3]]

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

dev.off()

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

#Figure 5 and Figure 6: Spatial transects for one day (maps)

#######Figure 5: Map of transects

## Transects image location in OR  

layout_m<-c(1,1)

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

    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)

plot(elev)

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

  filename<-sub(".shp","",transect_list[i])             #Removing the extension from file.

  transect<-readOGR(dirname(filename), basename(filename))                 #reading shapefile 

  #transect2 <- elide(transect, shift=c(0, -24000))

  #plot(transect2,add=TRUE)

  #writeOGR(transect2,dsn=".",layer="transect_OR_1_12282013",driver="ESRI Shapefile")

  plot(transect,add=TRUE)

}

text(x=vx_text,y=vy_text,labels=c("Transect 1","Transect 2","Transect 3"))

title("Transect Oregon")

dev.off()

######Figure 6: TRANSECT PROFILES

nb_transect <- 3

list_transect2<-vector("list",nb_transect)

list_transect3<-vector("list",nb_transect)

list_transect4<-vector("list",nb_transect)

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

rast_pred<-stack(lf[[2]]) #GAM_CAI

rast_pred_selected2<-subset(rast_pred,c(1,2,6)) #3 is referring to FSS, plot it first because it has the

rast_pred_selected3<-subset(rast_pred,c(1,3,6)) #3 is referring to FSS, plot it first because it has the

                                          # the largest range.

rast_pred2 <- stack(rast_pred_selected2,subset(s_raster,"elev_s"))

rast_pred3 <- stack(rast_pred_selected3,subset(s_raster,"elev_s"))

#layers_names<-layerNames(rast_pred2)<-c("lat*lon","lat*lon + elev + LST","elev")

layers_names2 <- names(rast_pred2)<-c("mod1","mod2","mod6","elev")

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

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

#transect_list

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

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

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

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

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

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

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

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

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

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

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

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

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

names(rast_pred2)<-layers_names2

names(rast_pred3)<-layers_names3

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

#r_stack<-rast_pred

#m_layers_sc<-c(3) #elevation in the third layer

m_layers_sc<-c(4) #elevation in the third layer

#title_plot2

#rast_pred2

#debug(plot_transect_m2)

trans_data2 <- plot_transect_m2(list_transect2,rast_pred2,title_plot2,disp=FALSE,m_layers_sc)

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

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

#Figure7: Spatial pattern: Image differencing  

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

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

methods_name <-c("gam_daily","gam_CAI","gam_fss")

index<-244 #index corresponding to Sept 1

y_var_name <-"dailyTmax"

ref_mod <- 3 #mod3

alt_mod <- 6 #mod6

file_format <- ".rst"

NA_flag_val <- -9999

list_param_diff <- list(index,list_raster_obj_files,methods_name,y_var_name,ref_mod,alt_mod,NA_flag_val,file_format,out_dir,out_prefix)

names(list_param_diff) <- c("index","list_raster_obj_files","methods_name","y_var_name","ref_mod","alt_mod","NA_flag_val","file_format","out_dir","out_prefix")

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

diff_pred_date1_list<- diff_date_rast_pred_fun(1,list_param_diff)

diff_pred_date2_list<- diff_date_rast_pred_fun(244,list_param_diff)

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

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

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

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

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

#plot(r_stack_diff,col=temp.colors(25))

#r_diff1<-subset(r_stack_diff,1:3)

#r_diff2<-subset(r_stack_diff,4:6)

#p1<-levelplot(r_diff1, col.regions=matlab.like(100)) #January

#p2<-levelplot(r_diff2,col.regions=matlab.like(100)) #January

#p1<-levelplot(r_stack_diff,layers=1:3, col.regions=matlab.like(100)) #January

#p2<-levelplot(r_stack_diff,layers=4:6, col.regions=matlab.like(100)) #January

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

levelplot(r_stack_diff,col.regions=matlab.like(100),

          ylab=NULL,xlab=NULL,

          par.settings = list(axis.text = list(font = 2, cex = 2),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),

) #January

dev.off()

###

#scales=list(log="e",x=list(cex=.3),y=list(cex=.3)) #change font size

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

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

LC_subset <- c("LC1","LC5","LC6","LC7","LC9","LC11")  

LC_names <- c("LC1_forest", "LC5_shrub", "LC6_grass", "LC7_crop", "LC9_urban","LC11_barren")

LC_s <- subset(s_raster,LC_subset)

names(LC_s) <- LC_names

plot(LC_s)

avl<-c(0,10,1,10,20,2,20,30,3,30,40,4,40,50,5,50,60,6,60,70,7,70,80,8,80,90,9,90,100,10)#Note that category 1 does not include 0!!

stat_list <- extract_diff_by_landcover(r_stack_diff,s_raster,LC_subset,LC_names,avl)

#r_subset_name <- "elev_s"

#r_names <- c("elev_s")

#stat_list_elev <- extract_diff_by_landcover(r_stack_diff,s_raster,LC_subset,LC_names,avl)

#write_out_raster_fun(s_raster,out_suffix=out_prefix,out_dir=out_dir,NA_flag_val=-9999,file_format=".rst")

#show correlation with LST by day over the year, ok writeout s_raster of coveriate??

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

#reorder the list

order_list<- c(1,4,2,5,3,6)

## Now create plots

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

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

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

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

par(mfrow=layout_m)    

#funciton plot

for (i in 1:length(stat_list$avg)){

  #i=i+1

  index <- order_list[i]

  zones_stat <- as.data.frame(stat_list$avg[[index]])

  zones_stat$zones <- 0:10

  plot(zones_stat$zones,zones_stat[,1],type="b",ylim=c(-4.5,6),

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

  #mtext("difference between mod3 and mod6 (degree C)",line=3,side=2,cex=1.2,font=2) #Add ylab with distance 3 from box

  #mtext("land cover percent classes",side=1,cex=1.2,line=3,font=2)

  lines(zones_stat$zones,zones_stat[,2],col="red",lty="dashed",pch=2) #shrub

  points(zones_stat$zones,zones_stat[,2],col="red",lty="dashed",pch=2) #shrub

  lines(zones_stat$zones,zones_stat[,3],col="green",lty="dotted",pch=3) #grass

  points(zones_stat$zones,zones_stat[,3],col="green",lty="dotted",pch=3) #grass

  lines(zones_stat$zones,zones_stat[,4],col="blue",lty="dashed",pch=4) #crop

  points(zones_stat$zones,zones_stat[,4],col="blue",lty="dashed",pch=4) #crop

  lines(zones_stat$zones,zones_stat[,5],col="darkgreen",lty="dashed",pch=5)

  points(zones_stat$zones,zones_stat[,5],col="darkgreen",lty="dashed",pch=5)

  lines(zones_stat$zones,zones_stat[,6],col="purple",lty="dashed",pch=6)

  points(zones_stat$zones,zones_stat[,6],col="purple",lty="dashed",pch=6)

  breaks_lab<-zones_stat$zones

  #make it slanted...

  tick_lab<-c("0","1-10","","20-30","","40-50","","60-70","","80-90","90-100") #Not enough space for  

  #tick_lab<-c("0","10-20","30-40","60-70","80-90","90-100")

  axis(side=1,las=1,tick=TRUE,

       at=breaks_lab,labels=tick_lab, cex.axis=1.2,font=2) #reduce number of labels to Jan and June

  #text(tick_lab, par(\u201cusr\u201d)[3], labels = tick_lab, srt = 45, adj = c(1.1,1.1), xpd = TRUE, cex=.9)

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

  box()

  legend("topleft",legend=names(zones_stat)[-7], 

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

        lty=1,pch=1:7)

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

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

}

#### Figure 11: Spatial lag profiles  

#date_selected <- "20100901"

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

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

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

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

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

p1 <- list_moranI_plots[[1]]

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

#Figure 12: Monthly climatology, Daily deviation and bias

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

#par(mfrow=layout_m)    

#plot(temp_day_s,col=  temp.colors(25),nr=layout_m,nc=layout_m,

#     cex=1.5) #use nr to choose number of columns in layout!!

p1 <- levelplot(temp_day_s,col.region=temp.colors(25),layer=1)

p2 <- levelplot(temp_day_s,col.region=temp.colors(25),layer=2)

p3 <- levelplot(temp_day_s,col.region=temp.colors(25),layer=3)

grid.arrange(p1,p2,p3,ncol=3)

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

####### Figure 13 : Disussion section

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

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

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

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

#Use data from accuraccy with no monthly hold out

##Now plot Figure 13a

plot(TMax_avgm,type="b",ylim=c(-10,50),col="black",xlab="month",ylab="tmax (degree C)")

lines(1:12,LST_avgm$LST,type="b",col="black",lty="dashed",pch=2)

lines(1:12,LST_avgm_min$LST,type="b",col="black",lty="dashed",pch=3)

lines(1:12,LST_avgm_max$LST,type="b",col="black",lty="dashed",pch=4)

text(1:12,LST_avgm$LST,month.abb,cex=0.8,pos=2)

legend("topleft",legend=c("TMax","LST","LST min","LST max"), cex=1, col=c("black","black","black","black"),

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

##Now plot Figure 13b

plot(TMax_avgm,type="b",ylim=c(-7,50),col="black",xlab="month",ylab="tmax (degree C)")

lines(1:12,LST_avgm$LST,type="b",col="black",lty="dashed",pch=2)

lines(LST_avgm_forest,type="b",col="black",lty="dotted",pch=3)

lines(LST_avgm_grass,type="b",col="black",lty="dotdash",pch=4)

legend("topleft",legend=c("TMax","LST","LST_forest","LST_grass"), cex=0.8, col=c("black","black","black","black","black"),

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

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

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

#get the  list of all prediction for a specific method

#list_lf <-extract_list_from_list_obj(load_obj(list_raster_obj_files$gam_CAI)$method_mod_obj,"dailyTmax") #getting objet

list_lf_gam_CAI <-extract_list_from_list_obj(load_obj(list_raster_obj_files[["gam_CAI"]])$method_mod_obj,"dailyTmax") #getting objet

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

nb_lag <- 10

list_filters<-lapply(1:nb_lag,FUN=autocor_filter_fun,f_type="queen") #generate lag 10 filters

list_param_stat_moran_CAI <- list(filter=list_filters[[10]],lf_list=list_lf_gam_CAI)

#tt <- stat_moran_std_raster_fun(1,list_param=list_param_stat_moran)

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

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

#list_param_stat_moran <- list(filter=list_filters[[10]],lf_list=list_lf_gam_fss)

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

#list_param_stat_moran <- list(filter=list_filters[[10]],lf_list=list_lf_gam_CAI)

#tt_CAI <- 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_CAI,file=file.path(out_dir,"moran_std_tt_CAI.RData"))

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

tt_fss2 <- load_obj("moran_std_tt_fss2.RData")

tx<- do.call(rbind,tt_fss2)

dates<-unique(raster_obj$tb_diagnostic_v$date)

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

dates_proc <- as.data.frame(dates_proc)

dates_proc$index <- 1:365

tx2 <- merge(tx,dates_proc,by="index")

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

names(tx2) <- c("index","moranI","std","pred_mod","date","month")

t<-melt(tx2,

          #measure=mod_var, 

          id=c("pred_mod","month"),

          na.rm=F)

#t$value<-as.numeric(t$value) #problem with char!!!

tb_mod_m_avg <-cast(t,pred_mod+month~variable,mean) #monthly mean for every model

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

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

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

# end of pasted

x<-subset(tb_mod_m_avg,pred_mod=="mod7")

##Now plot Figure 13c

plot(x$month,x$moranI,ylim=c(0.72,0.92),ylab="Lag 10 Moran's I autocorrelation",

     xlab="Month",

     type="b",col="black",lty="solid",pch=1)

plot(x$month,x$std,type="b",col="black",lty="dashed",pch=2,axes=F,

     xlab="",ylab="")

  #axis(4,xlab="",ylab="elevation(m)")  

axis(4,cex=1.2)

legend("topleft",legend=c("Moran's I","Stand Dev."), cex=0.8, col=c("black","black"),

       lty=1:2,pch=1:2,bty="n")

title("Spatial variation in FSS surfaces")

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

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

x_cor_tmp1 <-layerStats(r_tmp,stat="pearson",na.rm=T)

##Now plot Figure 13d

     ylab="Pearson corelation",xlab="month",main="Correlation between LST-elevation and LST-Forest")

lines(1:12,x_cor_tmp2[[1]][1:12,13],type="b",lty="dashed",pch=2)

legend("topleft",legend=c("Correlation LST-Elevation","Correlation LST-Forest"), cex=1.1, col=c("black","black"),

#closing file for figure 13