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

    trans_data<-extract(r_stack, transect)

    if (disp==FALSE){

    }

    #Plot layer values for specific transect

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

      } 

    }

    title(title_plot[i])

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

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

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"

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

s_raster <- brick(infile_covariates)

names(s_raster)<-covar_names

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

                          rep("gwr",7))    

#Check input covariates and model formula:

#strsplit(list_formulas,"~")

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

#Testing siginificance between models

res_pix<-960

col_mfrow<-1

row_mfrow<-1

    width=col_mfrow*res_pix,height=row_mfrow*res_pix)

par(mfrow=c(1,1))

plot(ghcn_dat_WGS84,add=T)

title_text <-c("Elevation (m) and meteorological"," stations in Oregon")

legend("topleft",legend=title_text,cex=2.1,bty="n")

par(mar = c(0,0,0,0)) # remove margin

#opar <- par(fig=c(0.9,0.95,0.8, 0.85), new=TRUE)

#opar <- par(fig=c(0.85,0.95,0.8, 0.9), new=TRUE)

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

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

lst_md<-raster(ref_rast_name)

projection(lst_md)<-projection(s_raster)

lst_mm_09<-subset(s_raster,"mm_09")

y_var_name <-"dailyTmax"

index<-244 #index corresponding to Sept 1

date_selected <- "20109101"

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

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

names_layers<-methods_names

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

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

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

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

  print(p)

  dev.off()

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

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

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

#r_stack<-rast_pred

#title_plot2

#rast_pred2

debug(plot_transect_m2)

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

dev.off()

###################### END OF SCRIPT #######################