############################  Script for manuscript analyses,tables and figures: MULTI TIME SCALE  ##############################

#This script uses the worklfow code applied to the Oregon case study. Single and Multi timescale procedures are compared.

#Single and Multitime scale (CAI and FSS)procedures are tested with different covariate combinations

#using GAM,GWR, Kriging. Results conurrently reside on NCEAS ATLAS server.

#Accuracy methods are added in the the function scripts to evaluate results.

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

#MODIFIED ON: 12/25/2013            

#Version: 3

function_analyses_paper1 <-"contribution_of_covariates_paper_interpolation_functions_10222013.R" #first interp paper

function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_12252013.R"

source(file.path(script_path,function_analyses_paper1)) #source all functions used in this script 1.

source(file.path(script_path,function_analyses_paper2)) #source all functions used in this script 2.

#Singlte time scale -direct methods: gam, kriging, gwr (no monthly holdout but 30%daily hold out)

#Multi time scale - CAI: gam, kriging, gwr

#Multi time scale - FSS: gam, kriging, gwr

### Multi time scale: Results from monthly hold out 0-70%

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

#Testing siginificance between models                    

## Transects image location in OR  

layout_m<-c(1,1)

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

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

vy_text  <-c(473000,297045.9,112611.5)

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

  #plot(transect2,add=TRUE)

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

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

#reorder the list

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

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

  index <- order_list[i]

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

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

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

#The list of files for 365 days in 2010...

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

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

#The list of files for 365 days in 2010...

#y_var_name <- dailyTmax

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

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

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

index <- 289

temp_day_s <- stack(lf_clim_gam_CAI[[index]]$mod7,lf_delta_gam_CAI[[index]]$mod7,lf_daily_gam_CAI[[index]]$mod7)

names(temp_day_s) <- c("Monhtly Climatology","Daily Devation","Daily Prediction")

layout_m <- c(1,3)

png(paste("Figure12_paper_climatology_daily_deviation_daily_prediction_model7_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)

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

dev.off()

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

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

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

#Get Tmax and LST averages

## open device for plottging

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

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

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

par(mfrow=layout_m)    

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

LST_avgm_min <- aggregate(LST~month,data=data_month_all,min)

LST_avgm_max <- aggregate(LST~month,data=data_month_all,max)

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

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

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

              lty=c("solid","dashed","dashed","dashed"),pch=1:4,bty="n")

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)

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

#lines(LST_avgm_urban,type="b",col="black",lty="longdash",pch=4)

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

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

nb_lag <- 10

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

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

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

dates<-unique(raster_obj$tb_diagnostic_v$date)

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

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

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

elev <- subset(s_raster,"elev_s")

LC1 <- subset(s_raster,"LC1")

#### Now elev

r_tmp1 <-stack(LST_s,elev)

names(r_tmp1) <- 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_tmp1 <-layerStats(r_tmp,stat="pearson",na.rm=T)

r_tmp2 <-stack(LST_s,LC1)

names(r_tmp2) <- 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","LC1")

x_cor_tmp2 <-layerStats(r_tmp2,stat="pearson",na.rm=T)

##Now plot Figure 13d

plot(1:12,x_cor_tmp1[[1]][1:12,13],ylim=c(-0.8,0.5),

     type="b",lty="solid",pch=1,

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

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

#closing file for figure 13

dev.off()

#DATE MODIFIED: 12/25/2013            

#Version: 2

#plot_mean_nobs_r_stack_by_month