#MODIFIED ON: 11/15/2013            

plot_transect_m2<-function(list_trans,r_stack,title_plot,disp=FALSE,m_layers){

  #This function creates plot of transects for stack of raster images.

  #Arguments:

  #list_trans: list of files containing the transects lines in shapefile format

  #r_stack: raster stack containing the information to extect

  #title_plot: plot title

  #disp: display and save from X11 if TRUE or plot to png file if FALSE

  #m_layers: index for layerers containing alternate units to be drawned on a differnt scale

  #RETURN:

  #list containing transect information

  nb<-length(list_trans)

  t_col<-rainbow(nb)

  t_col<-c("red","green","black")

  lty_list<-c("dashed","solid","dotted")

  list_trans_data<-vector("list",nb)

  #For scale 1

  for (i in 1:nb){

    trans_file<-list_trans[[i]][1]

    filename<-sub(".shp","",trans_file)             #Removing the extension from file.

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

    trans_data<-extract(r_stack, transect)

    if (disp==FALSE){

      png(file=paste(list_trans[[i]]),".png",sep="")

    }

    #Plot layer values for specific transect

    for (k in 1:ncol(trans_data[[1]])){

      y<-trans_data[[1]][,k]

      x<-1:length(y)

      m<-match(k,m_layers)

      if (k==1 & is.na(m)){

        plot(x,y,type="l",xlab="transect distance from coastal origin (km)", ylab=" maximum temperature (degree C)",

             ,cex=1.2,col=t_col[k])

        #axis(2)

      }

      if (k==1 & !is.na(m)){

        plot(x,y,type="l",col=t_col[k],lty="dotted",axes=F) #plotting fusion profile

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

        axis(4,cex=1.2)

      }

      if (k!=1 & is.na(m)){

        #par(new=TRUE)              # new plot without erasing old

        lines(x,y,type="l",xlab="",ylab="",col=t_col[k],axes=F) #plotting fusion profile

        #axis(2,xlab="",ylab="tmax (in degree C)")

      }

      if (k!=1 & !is.na(m)){

        par(new=TRUE)              # key: ask for new plot without erasing old

        plot(x,y,type="l",col=t_col[k],xlab="",ylab="",lty="dotted",axes=F) #plotting fusion profile

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

        axis(4,cex=1.2)

      } 

    }

    title(title_plot[i])

    legend("topleft",legend=layerNames(r_stack)[1:2], 

           cex=1.2, col=t_col,lty=1,bty="n")

    legend("topright",legend=layerNames(r_stack)[3], 

           cex=1.2, col=t_col[3],lty="dotted",bty="n")

    if (disp==TRUE){

      savePlot(file=paste(list_trans[[i]][2],".png",sep=""),type="png")

    }

    if (disp==FALSE){

      dev.off()

    }

    list_trans_data[[i]]<-trans_data

  }

  names(list_trans_data)<-names(list_trans)

  return(list_trans_data)

}

ref_rast_name<- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day244_rescaled.rst"                     #This is the shape file of outline of the study area. #local raster name defining resolution, exent, local projection--. set on the fly??

infile_reg_outline <- "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/OR83M_state_outline.shp"  #input region outline defined by polygon: Oregon

ref_rast_name <-"/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day244_rescaled.rst"  #local raster name defining resolution, exent: oregon

ref_rast_d001 <-"/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/mean_day001_rescaled.rst"

transect_list <-c("/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/transect_OR_1.shp",

                  "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/transect_OR_2.shp",

                  "/data/project/layers/commons/data_workflow/inputs/region_outlines_ref_files/transect_OR_3.shp")

met_stations_obj <- load_obj(file.path(in_dir1,met_obj_file_1))

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

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

### CREATE FIGURE MEAN DAILY AND MEAN MONTHLY: AAG 2013  ####

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("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 January 1")

plot(lst_mm_01)

plot(interp_area,add=TRUE)

title("Mean for month of January")

dev.off()

y_var_name <-"dailyTmax"

index<-244 #index corresponding to Sept 1

lf1 <- raster_prediction_obj_1$method_mod_obj[[index]][[y_var_name]] #select relevant raster images for the given dates

lf4 <- raster_prediction_obj_4$method_mod_obj[[index]][[y_var_name]]

lf7 <- raster_prediction_obj_7$method_mod_obj[[index]][[y_var_name]]

date_selected <- "20109101"

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

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

names_layers<-methods_names

lf <- (list(lf1,lf4[1:7],lf7[1:7]))

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

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

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

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

  #lf2 <- raster_prediction_obj_2$method_mod_obj[[index]][[y_var_name]]

  #lf2 #contains the models for gam

  #s.range <- c(min(minValue(pred_temp_s)), max(maxValue(pred_temp_s)))

  #s.range <- s.range+c(5,-5)

  #col.breaks <- pretty(s.range, n=200)

  #lab.breaks <- pretty(s.range, n=100)

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

  #max_val<-s.range[2]

  #min_val <-s.range[1]

  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="Interpolated Surfaces Model Comparison baseline 1", 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),

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

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

  print(p)

  dev.off()

}

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

#Figure 7: Spatial transects for one day (maps)

#######Figure 7a: Map of transects

nb_transect<-3

list_transect2<-vector("list",nb_transect)

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

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

                                             # the largest range.

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

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

layers_names<-layerNames(rast_pred2)<-c("mod1","mod6","elev")

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

transect_list

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

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

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

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

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

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

names(list_transect2)<-c("transect_OR1","transect_OR2","transect_OR3")

#X11(width=9,height=9)

#png(paste("fig3_elevation_transect1_path_CAI_fusion_",date_selected,out_prefix,".png", sep=""))

#plot(elev)

#k<-1  #transect to plot

#trans_file<-list_transect2[[k]][[1]]

#filename<-sub(".shp","",trans_file)             #Removing the extension from file.

#transect<-readOGR(".", filename)                 #reading shapefile 

#plot(transect,add=TRUE)

#title("Transect Oregon")

#dev.off()

layerNames(rast_pred2)<-layers_names

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

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

#r_stack<-rast_pred

m_layers_sc<-c(3)

#title_plot2

#rast_pred2

debug(plot_transect_m2)

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

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

par(mfrow=c(1,2))

dev.off()