# R BATCH --no-save < /project/inf1w/gilliam/ENSEMBLE_CASE/aq_stats_presentations/boxplot.cities.R 
 require(ncdf)
 require(maps)
 require(mapdata)
 require(date)
 require(fields)
 source("/project/inf39w/gilliam/AMET_v12/R/MET_amet.stats-lib.R")

####################################################################################
# Process time series of daily model and observed ozone from ensembles
dir_to_look<-"/project/inf39w/gilliam/ENSEMBLE_CASE/aq_stats_presentations/houston"
str_look<-"*.csv"

files<-system(paste("ls -l ",dir_to_look,"/",str_look,sep=""),intern=T)
 
 ens_id<-c("ENS_Output","EM_Ctl","EM_P1","EM_P2","EM_N1","EM_N2","NMM_Ctl","NMM_P1",
           "NMM_P2","NMM_N1","NMM_N2","ETA_Ctl1","ETA_Ctl2","ETA_P1",
           "ETA_P2","ETA_N1","ETA_N2")
 ens_name_amet<-paste("CMAQ_WRF_Ensemble_",ens_id,sep="")
 
 nens<-length(ens_id)
 ndays<-4
 date<-array(NA,c(ndays))
 obs<-array(NA,c(ndays))
 mod<-array(NA,c(nens,ndays))
 ens_id<-array(NA,c(nens))
 
 
 nf<-length(files)
 
 for(f in 1:nf) {
  files[f]
  parts<-unlist(strsplit(files[f], " "))
  file<-parts[length(parts)] 
  writeLines(paste("WRF file:",file))
  
  a<-read.csv(file[1],fill=T, header=F)
  nl<-dim(a)[1]
  
  for(l in 1:nl){
   for(ens in 1:nens){
     if(a[l,1] == ens_name_amet[ens]){
       writeLines(paste("Project record found, will grab data",ens_name_amet[ens]))
       obs[1:ndays]<- as.numeric(as.character(a[(l+2):(l+1+ndays),2]))
       mod[ens,1:ndays]<-as.numeric(as.character(a[(l+2):(l+1+ndays),3]))
     }
   
   } 
  }
 
 }
pdf(outname_pdf, width=8, height=8)                                             # Set output device with options

{
   par(mai=c(1,1,0.5,0.5))                                                         # Set plot margins
   boxplot_info <- boxplot(mod[,1],plot=F)
   boxplot(mod[,1],range=0, border="transparent", col=colors, ylim=c(y.axis.min, y.axis.max), xlab="", ylab=y_label, names=labels, cex.axis=1, cex.lab=1, cex=1, las=1, boxwex=.4) # Create boxplot
   x.axis.max <- length(run_names)
}
####################################################################################


####################################################################################
# Process spatial daily ozone CSV files from AMET for various Ensemble analyses
dir<-"/project/inf39w/gilliam/ENSEMBLE_CASE/aq_stats_presentations/spatial/BASE_RUN2"
 
ens_id<-c("ens_mean","em_ctl","em_p1","em_p2","em_n1","em_n2","nmm_ctl","nmm_p1","nmm_p2","nmm_n1","nmm_n2","eta_ctl1","eta_ctl2","eta_p1","eta_n1","eta_n2")
ens_id<-c("em_ctl","em_p1","em_p2","em_n1","em_n2","nmm_ctl","nmm_p1","nmm_p2","nmm_n1","nmm_n2","eta_ctl1","eta_p1","eta_p2","eta_n1","eta_n2")
# Note this is the one used for paper with missing eta_p2 member
ens_id<-c("em_ctl","em_p1","em_p2","em_n1","em_n2","nmm_ctl","nmm_p1","nmm_p2","nmm_n1","nmm_n2","eta_ctl1","eta_p1","eta_n1","eta_n2")
files<-paste(dir,"/spatial_",ens_id,".csv",sep="")

nf<-length(files)

####################################################################### 
for(f in 1:nf){
 writeLines(paste("Processing file:",files[f])) 
 a<-read.csv(files[f],fill=T, header=T,skip=4)
 ############### 
 if(f ==1 ) {
   statx<-a[,2]
   latx<-a[,3]
   lonx<-a[,4]
   obdx<-a[,5]
   obsx<- a[,9]
   modx<- a[,10]
   
   stat_id<-sort(unique(statx))
   dates<-sort(unique(obdx))

   ns<-length(stat_id)
   nd<-length(dates)
   obs<-array(NA,c(ns,nd))
   mod<-array(NA,c(ns,nd,nf))   
   lats<-array(NA,c(ns))   
   lons<-array(NA,c(ns))   
 }
 ############### 
 for(s in 1:ns){
   for(d in 1:nd){
     ind<-which(stat_id[s] == a[,2] & dates[d] == a[,5])
     # get data if a hit is found 
     if(length(ind)> 0){
       obs[s,d]<-a[ind,9]
       mod[s,d,f]<-a[ind,10]
       lats[s]<-a[ind,3]  
       lons[s]<-a[ind,4]   
     }
   }
 }
 ############### 
 
} 
obs<-ifelse(obs <= 0, NA, obs)
mod<-ifelse(mod <=0, NA, mod)
####################################################################### 

# GFS or Baseline run file read in. This is done seprately since it's not part of the ensemble

fileBaseline<-paste(dir,"/spatial_gfs.csv",sep="")
modBaseline<-array(NA,c(ns,nd))   
 writeLines(paste("Processing file:",files[f])) 
 a<-read.csv(fileBaseline,fill=T, header=T,skip=4)
 ############### 
 for(s in 1:ns){
   for(d in 1:nd){
     ind<-which(stat_id[s] == a[,2] & dates[d] == a[,5])
     # get data if a hit is found 
     if(length(ind)> 0){
       modBaseline[s,d]<-a[ind,10]
     }
   }
 }
############### 
####################################################################### 
##################  ENS Scatterplots of all Sites   ##################
# 
####################################################################### 
# 1. obs vs. model baseline all days
rng<-c(0,round(max(obs[,],modBaseline[,],na.rm=T)))
m1<-round(mbias(modBaseline,obs,na.rm=T),2)
m2<-round(rmserror(modBaseline,obs,na.rm=T),2)

plot(obs[,],modBaseline[,],ylim=rng,xlim=rng, ylab="Base Model (ppb)", xlab="Observations (ppb)")
text(3,rng[2],paste("RMSE:",m2),pos=4)
text(3,rng[2]-5,paste("ME:",m1),pos=4)

# 2. obs vs. ENS median, range, best
med3rd<-c(1,2,3,4,5,11,12,13,14,15)
ens_median<-array(NA,c(ns,nd))
ens_median3rd<-array(NA,c(ns,nd))
ens_range<-array(NA,c(ns,nd))
ens_best<-array(NA,c(ns,nd))
ens_besti<-array(NA,c(ns,nd))
for(s in 1:ns){
 for(d in 1:nd){
   ens_median[s,d]<-median(mod[s,d,],na.rm=T)
   ens_median3rd[s,d]<-median(mod[s,d,med3rd],na.rm=T)
   ens_range[s,d]<-abs(diff(range(mod[s,d,],na.rm=T)))
   
   diffmo<-mod[s,d,]-obs[s,d]
   bestvalabs<-min(abs(diffmo),na.rm=T)
   check<-which(bestvalabs == diffmo)
   bestval<-bestvalabs
   if(length(check) > 0){
   }
   else {
    check<-which(bestvalabs*-1 == diffmo)
    bestval<-bestvalabs*-1
   }
   ens_best[s,d]<-mod[s,d,check[1]]
   ens_besti[s,d]<-check[1]
 }
}

################################################
# Obs vs. Best/Closest ENS member
m1<-round(mbias(ens_best,obs,na.rm=T),2)
m2<-round(rmserror(ens_best,obs,na.rm=T),2)
rng<-c(0,round(max(obs[,],ens_best[,],na.rm=T)))
levs<-seq(1,15,by=1)
levcols<-(tim.colors(length(levs)-1))
levcols<-gray(levs/15)
pcols<-levcols[cut(ens_besti,br=levs,labels=FALSE,include.lowest=T,right=T)]


plot(obs[,],ens_best[,],ylim=rng,xlim=rng, ylab="Closest ENS(ppb)", xlab="Observations (ppb)")
text(3,rng[2],paste("RMSE:",m2),pos=4)
text(3,rng[2]-5,paste("ME:",m1),pos=4)

# Histogram of best Ensemble member match for 4 day daily Max 8hr ozone
hist(ens_besti)
################################################

################################################
# Obs vs. median ENS
m1<-round(mbias(ens_median,obs,na.rm=T),2)
m2<-round(rmserror(ens_median,obs,na.rm=T),2)
rng<-c(0,round(max(obs[,],ens_median[,],na.rm=T)))

plot(obs[,],ens_median[,],ylim=rng,xlim=rng, ylab="Median ENS(ppb)", xlab="Observations (ppb)")
text(3,rng[2],paste("RMSE:",m2),pos=4)
text(3,rng[2]-5,paste("ME:",m1),pos=4)
################################################

################################################
# Obs vs. median of top best 2/3 of the ENS
m1<-round(mbias(ens_median3rd,obs,na.rm=T),2)
m2<-round(rmserror(ens_median3rd,obs,na.rm=T),2)
rng<-c(0,round(max(obs[,],ens_median3rd[,],na.rm=T)))

plot(obs[,],ens_median3rd[,],ylim=rng,xlim=rng, ylab="Median ENS(ppb)", xlab="Observations (ppb)")
text(3,rng[2],paste("RMSE:",m2),pos=4)
text(3,rng[2]-5,paste("ME:",m1),pos=4)
################################################

####################################################################### 
##################  Ranked O3 Histogram   ##################
# 
####################################################################### 
# Rank Hist for All Days combined.
rnk_hista<-array(NA,c(ns*nd))
count<-1
for(s in 1:ns){
 for(d in 1:nd){
  vec<-sort(mod[s,d,])
  if(is.na(obs[s,d])) { count<-count+1; next }
  if(obs[s,d] < vec[1]) { rnk_hista[count]<-1 }
  if(obs[s,d] > vec[nf]){ rnk_hista[count]<-nf+1 }
  for(ens in 1:(nf-1)) {  
   if(obs[s,d] > vec[ens] & obs[s,d] < vec[ens+1]) {rnk_hista[count]<-ens+1 }
  }
  count<-count+1
 }
}
figure<-paste(dir,"/ENS_RankedHist_O3_Alldays.png",sep="") 
writeLines(paste("Plotting:",figure))
png(file=figure, width=800, height=800)
par(cex.axis=1.5)
par(cex.lab=1.5)
hist(rnk_hista, main ="Ranked Histogram of Ozone",freq=F,xlab="Rank")
dev.off()

# Rank Hist for each day 
rnk_histd<-array(NA,c(ns,nd))
for(s in 1:ns){
 for(d in 1:nd){
  vec<-sort(mod[s,d,])
  if(is.na(obs[s,d])) { next }
  if(obs[s,d] < vec[1]) { rnk_histd[s,d]<-1 }
  if(obs[s,d] > vec[nf]){ rnk_histd[s,d]<-nf+1 }
  for(ens in 1:(nf-1)) {  
   if(obs[s,d] > vec[ens] & obs[s,d] < vec[ens+1]) {rnk_histd[s,d]<-ens+1 }
  }
 }
}
for(d in 1:nd){
 figure<-paste(dir,"/ENS_RankedHist_O3_",dates[d],".png",sep="") 
 writeLines(paste("Plotting:",figure))
 png(file=figure, width=800, height=800)
 par(cex.axis=1.5)
 par(cex.lab=1.5)
 hist(rnk_histd[,d], main ="Ranked Histogram of Ozone",freq=F,xlab="Rank")
 dev.off()
}
####################################################################### 
####################################################################### 

####################################################################### 
##################   Box Plots for Defined Cities #####################
####################################################################### 
# Box plots for defined Cities
ens_small_id<-c("ctl","p1","p2","n1","n2")
cities<-c("LosAngeles","Sacramento","Phoenix","Denver","Dallas","Houston","Chicago","Columbus","Tampa","Charlotte","DC-Balt","NewYork","Cleveland")
cityabrv<-c("LA","SAC","PHX","DEN","DFW","HOU","CHI","COL","TMP","CHT","DCB","NYC","CLV")
nc<-length(cities)
files.city<-paste(dir,"/sites_stats_",cityabrv,".csv",sep="")

c_ens_rmse<-array(NA,c(nc,nd,nf))
c_base_rmse<-array(NA,c(nc,nd))

box.range<-range(obs,mod,modBaseline,na.rm=T)
ind<-1:ns
for(c in 1:nc){
 a<-read.csv(files.city[c],fill=T, header=T,skip=9)
 ### Find the index of sites in defined city
 cities.to.get<-a[,2]
 nctg<-length(cities.to.get)
 ind<-array(NA,nctg)
 ens_memb_range<-array(NA,c(nctg,nd,nf+2))
 for(cc in 1:nctg){
   indtg<-which(cities.to.get[cc] == stat_id)
   if(length(indtg > 0)){
     ind[cc]<-indtg
   }  
 }
 ######
 for(d in 1:nd){
 
   obs_box<-obs[ind,d]
   baseline_box<-modBaseline[ind,d]
   mod_box<-matrix(mod[ind,d,])
   obsx<-array(0.9,length(obs_box))
   basex<-array(1.1,length(baseline_box))
   # Compute number and percentage of obs hits that fall into model range
   range.mod<-range(mod_box)
   obs.check<-ifelse( obs_box > range.mod[1] & obs_box < range.mod[2],1,0)
   obs.check2<-sum(ifelse(is.na(obs.check),1,0),na.rm=T)
   num.hits<-sum(obs.check, na.rm=T)
   perc.hits<-round(100*num.hits/(length(obs_box)-obs.check2))

   med_obs<-median(obs_box,na.rm=T)
   med_modbase<-median(baseline_box,na.rm=T)

  #################################
  #    PLOT SECTION
  #################################
  # Obs, Base vs. Each Ensemble
  figure<-paste(dir,"/ENSJRG_",cityabrv[c],"_",dates[d],".png",sep="") 
  writeLines(paste("Plotting:",figure))
  png(file=figure, width=500, height=500)

  ylims<-range(obs_box,baseline_box,mod_box,na.rm=T)
  plot(17,1,col="white",xlim=c(1,17),ylim=ylims,ylab="",xlab="",axes=FALSE)
  par(new=T)
  plot(array(1,c(nctg)),obs_box,col="black",xlim=c(1,17),ylim=ylims,ylab="",xlab="",axes=FALSE)
  par(new=T)
  plot(array(2,c(nctg)),baseline_box,col="brown",xlim=c(1,17),ylim=ylims,ylab="",xlab="",axes=FALSE)
  c_base_rmse[c,d]<-round(rmserror(baseline_box,obs_box,na.rm=T),1)
  text(2,ylims[2]+2,c_base_rmse[c,d],pos=1,cex=0.85)
  text(1,ylims[1]+1,"obs",pos=1,cex=0.8)
  text(2,ylims[1]+1,"base",pos=1,cex=0.8)
  box()
  text(5,ylims[1]+3,"EM",pos=1,cex=1)
  for(e in 1:5){
   c_ens_rmse[c,d,e]<-round(rmserror(mod[ind,d,e],obs_box,na.rm=T),1)
   par(new=T)
   plot(array(e+2,c(nctg)),mod[ind,d,e],col="red",xlim=c(1,17),ylim=ylims,ylab="",xlab="",axes=FALSE)
   text(e+2,ylims[2]+2,c_ens_rmse[c,d,e],pos=1,cex=0.85)
   text(e+2,ylims[1]+1,ens_small_id[e],pos=1,cex=0.8)
  }
  text(10,ylims[1]+3,"NMM",pos=1,cex=1)
  for(e in 1:5){
   c_ens_rmse[c,d,e+5]<-round(rmserror(mod[ind,d,e+5],obs_box,na.rm=T),1)
   par(new=T)
   plot(array(e+7,c(nctg)),mod[ind,d,e+5],col="green",xlim=c(1,17),ylim=ylims,ylab="",xlab="",axes=FALSE)
   text(e+7,ylims[2]+2,c_ens_rmse[c,d,e+5],pos=1,cex=0.85)
   text(e+7,ylims[1]+1,ens_small_id[e],pos=1,cex=0.8)
  }
  text(15,ylims[1]+3,"ETA",pos=1,cex=1)
  for(e in 1:4){
   c_ens_rmse[c,d,e+10]<-round(rmserror(mod[ind,d,e+10],obs_box,na.rm=T),1)
   par(new=T)
   plot(array(e+12,c(nctg)),mod[ind,d,e+10],col="blue",xlim=c(1,17),ylim=ylims,ylab="",xlab="",axes=FALSE)
   text(e+12,ylims[2]+2,c_ens_rmse[c,d,e+10],pos=1,cex=0.85)
   text(e+12,ylims[1]+1,ens_small_id[e],pos=1,cex=0.8)
  }
  axis(2,col="black",tick=T)
  title(paste(cities[c],"Area Ozone (ppb) on",dates[d]))
  dev.off()
  #################################
  #################################
}}


  #################################
  #################################
  #################################
  # Plot Barchart NOTE: Copy and paste this in plot section above
  figure<-paste(dir,"/BOX.ALL_",cityabrv[c],"_",dates[d],".png",sep="") 
  writeLines(paste("Plotting:",figure))
  png(file=figure, width=200, height=400)
   par(mai=c(.1,.35,0.1,0.1))
   boxplot(mod_box,col="gray", ylim=box.range, range=0,xlab=cities[c], ylab="8-hr MaxOzone(pbb)", border="black", pars = list(boxwex = 0.8, staplewex = 0.5, outwex = 0.5),varwidth=T)
   points(0.9,med_obs,pch=19,col="red",cex=2)
   points(1.1,med_modbase,pch=19,col="blue",cex=2)
   points(obsx,obs_box,pch=20,col="red")
   points(obsx,obs_box,pch=21,col="black")
   points(0.9,med_obs,pch=18,col="red")
   points(basex,baseline_box,pch=20,col="blue")
   points(basex,baseline_box,pch=21,col="black")
   hit.str<-paste("Obs Hits: ",num.hits," of ",length(obs_box)-obs.check2," or ~",perc.hits,"%",sep="")
   text(.45,box.range[2]-10,hit.str,pos=4)
   text(.45,box.range[2]-5,paste("Date:",dates[d]),pos=4)
   text(.45,box.range[2],paste("City:",cities[c]),pos=4)
  dev.off()
  #################################
  #################################

####################################################################### 
##################   Ranked Histograms for Defined Cities #############
####################################################################### 
# Box plots for defined Cities
ens_small_id<-c("ctl","p1","p2","n1","n2")
cities<-c("LosAngeles","Sacramento","Phoenix","Denver","Dallas","Houston","Chicago","Columbus","Tampa","Charlotte","DC-Balt","NewYork","Cleveland")
cityabrv<-c("LA","SAC","PHX","DEN","DFW","HOU","CHI","COL","TMP","CHT","DCB","NYC","CLV")
nc<-length(cities)
files.city<-paste(dir,"/sites_stats_",cityabrv,".csv",sep="")

c_ens_rmse<-array(NA,c(nc,nd,nf))
c_base_rmse<-array(NA,c(nc,nd))

box.range<-range(obs,mod,modBaseline,na.rm=T)
ind<-1:ns
for(c in 1:nc){
 a<-read.csv(files.city[c],fill=T, header=T,skip=9)
 ### Find the index of sites in defined city
 cities.to.get<-a[,2]
 nctg<-length(cities.to.get)
 ind<-array(NA,nctg)
 ens_memb_range<-array(NA,c(nctg,nd,nf+2))
 for(cc in 1:nctg){
   indtg<-which(cities.to.get[cc] == stat_id)
   if(length(indtg > 0)){
     ind[cc]<-indtg
   }  
 }
 #################################
 #  RANKED HIST PLOT SECTION
 #################################
 ni<-length(ind)
 rnk_hista<-array(NA,c(ni*(nd-1)))
 count<-1
 for(i in 1:ni){
  for(d in 2:nd){
   vec<-sort(mod[ind[i],d,])
   if(is.na(obs[ind[i],d])) { count<-count+1; next }
   if(obs[ind[i],d] < vec[1]) { rnk_hista[count]<-1 }
   if(obs[ind[i],d] > vec[nf]){ rnk_hista[count]<-nf+1 }
   for(ens in 1:(nf-1)) {  
    if(obs[ind[i],d] > vec[ens] & obs[ind[i],d] < vec[ens+1]) {rnk_hista[count]<-ens+1 }
   }
   count<-count+1
   }
  }
 writeLines(paste("Number of samples for",cities[c],"is",length(rnk_hista)))
 figure<-paste(dir,"/RankedHist_O3_Alldays_",cityabrv[c],".png",sep="") 
 writeLines(paste("Plotting:",figure))
 png(file=figure, width=800, height=800)
 par(cex.axis=1.5)
 par(cex.lab=1.5)
 hist(rnk_hista, main =paste("Ranked Histogram of Ozone:",cities[c]),freq=F,breaks=nf+1,xlab="Rank")
 dev.off()
 #################################
 #################################
}
# END RANKED HIST
####################################################################### 


####################################################################### 
##################   Spatial Hit Percentage By Day ####################
####################################################################### 
range.lat<-range(lats,na.rm=T)
range.lon<-range(lons,na.rm=T)

hit.ozone<-array(NA,c(ns,nd))
min.diff<-array(NA,c(ns,nd))
range.ozone<-array(NA,c(ns,nd))

for(d in 1:nd){
 for(s in 1:ns){
   
   range.mod<-range(mod[s,d,],na.rm=T)
   hit.ozone[s,d]<-NA

   if( is.na(obs[s,d]) ) { next;}
   hit.ozone[s,d]<-0
   if( obs[s,d] > range.mod[1] & obs[s,d] < range.mod[2]) {  hit.ozone[s,d]<-1 }
   min.diff[s,d]<-min(abs(mod[s,d,] - obs[s,d]), na.rm=T)
   range.ozone[s,d]<-range.mod[2]-range.mod[1]
   writeLines(paste("Site",s,ns,"Range Ozone",range.mod[1],range.mod[2],obs[s,d],hit.ozone[s,d]))
   
 }

 figure<-paste(dir,"/spatial.hit.",dates[d],".pdf",sep="")
 pdf(file=figure, width = 11, height = 8.5)
 writeLines(paste("Plotting figure",figure))
 m<-map('usa',plot=FALSE)
 map("state", xlim=c( min(lons)-5,max(lons)+5),ylim=c( min(lats)-5,max(lats)+5),resolution=0, boundary = TRUE, lty = 1)
 box()
 points(lons,lats,pch=as.character(hit.ozone[,d]),cex= 1+range.ozone[,d]/max(range.ozone[,d],na.rm=T))
 dev.off() 
}
#### Generate 4 day cummulative hits and csv output data.frame
hit.total<-array(NA,c(ns))
mean.min.diff<-array(NA,c(ns))
mean.range<-array(NA,c(ns))
 for(s in 1:ns){
   hit.total[s]<-sum(hit.ozone[s,],na.rm=T)
   mean.min.diff[s]<-mean(min.diff[s,],na.rm=T)
   mean.range[s]<-mean(range.ozone[s,],na.rm=T)
 }
google.data<- data.frame(stat_id,lats,lons,hit.total,mean.min.diff,mean.range)
####################################################################### 

####################################################################### 
##################   Spatial Error (MAE) and SDev of Error ############
####################################################################### 
range.lat<-range(lats,na.rm=T)
range.lon<-range(lons,na.rm=T)
range.lon[1]<-range.lon[1]-5
range.lat[1]<-range.lat[1]-2
range.lat[2]<-range.lat[2]+2

mea.ozone<-array(NA,c(ns))
bias.ozone<-array(NA,c(ns))
std.ozone<-array(NA,c(ns))
range.mea.ozone<-array(NA,c(ns))

mea.ozone.d<-array(NA,c(ns,nd))
bias.ozone.d<-array(NA,c(ns,nd))
std.ozone.d<-array(NA,c(ns,nd))
range.mea.ozone.d<-array(NA,c(ns,nd))

for(s in 1:ns){
   all.diff<-matrix(mod[s,,]-obs[s,])
   range.diff<-range(abs(all.diff),na.rm=T) 
   mea.ozone[s]<-mean(abs(all.diff),na.rm=T)
   bias.ozone[s]<-mean(all.diff,na.rm=T)
   std.ozone[s]<-sd(abs(all.diff),na.rm=T)   
   range.mea.ozone[s]<-  range.diff[2]-range.diff[1]
   for(d in 1:nd){
    all.diff<-matrix(mod[s,d,]-obs[s,d])
    range.diff<-range(abs(all.diff),na.rm=T) 
    mea.ozone.d[s,d]<-mean(abs(all.diff),na.rm=T)
    bias.ozone.d[s,d]<-mean(all.diff,na.rm=T)
    std.ozone.d[s,d]<-sd(abs(all.diff),na.rm=T)   
    range.mea.ozone.d[s,d]<-  range.diff[2]-range.diff[1]
   }
}

#####################################################################################################
# Mean MEA
field1<-mea.ozone
field2<-std.ozone
title<-paste("Mean Ensemble MAE (pbb)")
levs<-c(seq(0,20,by=2),30,40)
#####
levcols<-(tim.colors(length(levs)-1))

symb_size<-0.5 +  (field2/mean(field2,na.rm=T))
pdf(file=paste("spatial.mean.ens.mae.pdf",sep=""), width = 11, height = 8.5)

 m<-map('usa',plot=FALSE)
 map("state", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col=1, lty = 1,lwd=1)
 map("county", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col="gray", lty = 1, lwd=0.25, add=T)
 box()

 pcols<-levcols[cut(field1,br=levs,labels=FALSE,include.lowest=T,right=T)]
 points(lons,lats,pch=20, cex=symb_size, col=pcols)
legend(range.lon[1],range.lat[1],levs,col=levcols,pch=20,xjust=0,yjust=0, pt.cex=0.75, cex=0.75)
 title(main=title,cex.main = 0.90, line=0.25)
dev.off() 
#####################################################################################################
#####################################################################################################
# StDev of Ensemble MEA
field1<-std.ozone
field2<-1
title<-paste("StDev of Ensemble MAE (pbb)")
levs<-c(seq(0,20,by=2),30,40)
#####
levcols<-(tim.colors(length(levs)-1))

symb_size<-0.5 +  (field2/mean(field2,na.rm=T))
pdf(file=paste("spatial.sdev.of.ens.mae.pdf",sep=""), width = 11, height = 8.5)

 m<-map('usa',plot=FALSE)
 map("state", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col=1, lty = 1,lwd=1)
 map("county", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col="gray", lty = 1, lwd=0.25, add=T)
 box()

 pcols<-levcols[cut(field1,br=levs,labels=FALSE,include.lowest=T,right=T)]
 points(lons,lats,pch=20, cex=symb_size, col=pcols)
legend(range.lon[1],range.lat[1],levs,col=levcols,pch=20,xjust=0,yjust=0, pt.cex=0.75, cex=0.75)
 title(main=title,cex.main = 0.90, line=0.25)
dev.off() 
#####################################################################################################
#####################################################################################################
# Range of Ensemble MEA
field1<-range.mea.ozone
field2<-1
title<-paste("Range of Ensemble MAE (pbb)")
levs<-c(seq(0,50,by=5))
#####
levcols<-(tim.colors(length(levs)-1))

symb_size<-0.5 +  (field2/mean(field2,na.rm=T))
pdf(file=paste("spatial.range.of.ens.mae.pdf",sep=""), width = 11, height = 8.5)

 m<-map('usa',plot=FALSE)
 map("state", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col=1, lty = 1,lwd=1)
 map("county", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col="gray", lty = 1, lwd=0.25, add=T)
 box()

 pcols<-levcols[cut(field1,br=levs,labels=FALSE,include.lowest=T,right=T)]
 points(lons,lats,pch=20, cex=symb_size, col=pcols)
legend(range.lon[1],range.lat[1],levs,col=levcols,pch=20,xjust=0,yjust=0, pt.cex=0.75, cex=0.75)
 title(main=title,cex.main = 0.90, line=0.25)
dev.off() 
#####################################################################################################
#####################################################################################################
#####################################################################################################
 for(d in 1:nd){
#####################################################################################################
# Mean MEA
field1<-mea.ozone.d[,d]
field2<-std.ozone.d[,d]
title<-paste("Mean Ensemble MAE (pbb)",dates[d])
levs<-c(seq(0,20,by=2),30,40)
#####
levcols<-(tim.colors(length(levs)-1))

symb_size<-0.5 +  (field2/mean(field2,na.rm=T))
pdf(file=paste("spatial.",dates[d],".mean.ens.mae.pdf",sep=""), width = 11, height = 8.5)

 m<-map('usa',plot=FALSE)
 map("state", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col=1, lty = 1,lwd=1)
 map("county", xlim=range.lon,ylim=range.lat,resolution=0, boundary = TRUE, col="gray", lty = 1, lwd=0.25, add=T)
 box()

 pcols<-levcols[cut(field1,br=levs,labels=FALSE,include.lowest=T,right=T)]
 points(lons,lats,pch=20, cex=symb_size, col=pcols)
legend(range.lon[1],range.lat[1],levs,col=levcols,pch=20,xjust=0,yjust=0, pt.cex=0.75, cex=0.75)
 title(main=title,cex.main = 0.90, line=0.25)
dev.off() 
#####################################################################################################
}
#####################################################################################################
#####################################################################################################

####################################################################### 
#######  AMET-BASED Ranked Histogram of Meteorology Timeseries  #######
#######################################################################
site<-"KDCA"
plotdir<-"/project/inf39w/gilliam/ENSEMBLE_CASE/aq_stats_presentations/spatial/BASE_RUN2"
varid<-c("T","Q","WS","WD")
varlab<-c("2-m Temp","2-m WVMR","10-m WS","10-m WD")
dir<-"/project/inf39w/gilliam/AMET_v12/output/" 
ens_id<-c("ens_em_ctl","ens_em_p1","ens_em_p2","ens_em_n1","ens_em_n2","ens_nmm_ctl","ens_nmm_p1","ens_nmm_p2","ens_nmm_n1","ens_nmm_n2","ens_eta_ctl1","ens_eta_p1","ens_eta_p2","ens_eta_n1","ens_eta_n2")
files<-paste(dir,ens_id,"/timeseries/",ens_id,".",site,".txt",sep="")

nf<-length(files)

####################################################################### 
for(f in 1:nf){
 writeLines(paste("Processing file:",files[f])) 
 a<-read.csv(files[f],fill=T, header=T,skip=1)
 ############### 
 if(f ==1 ) {
   nt<-dim(a)[1]
   ens<-array(NA,c(nt,nf,4))
   obs<-array(NA,c(nt,4))
   dates<-a[,1]
   obs[,1]<-a[,3]
   obs[,2]<-a[,6]
   obs[,3]<-a[,8]
   obs[,4]<-a[,10]
 }
 ############### 
 ens[,f,1]<-a[,2]
 ens[,f,2]<-a[,5]
 ens[,f,3]<-a[,7]
 ens[,f,4]<-a[,9]
 ############### 
} 
obs<-ifelse(obs <= 0, NA, obs)
ens<-ifelse(ens <=0, NA, ens)

 rnk_hista<-array(NA,c(nt,4))
 for(v in 1:4){
  count<-1
  for(tt in 1:nt){
   vec<-sort(ens[tt,,v])
   if(is.na(obs[tt,v])) { count<-count+1; next }
   if(obs[tt,v] < vec[1]) { rnk_hista[tt,v]<-1 }
   if(obs[tt,v] > vec[nf]){ rnk_hista[tt,v]<-nf+1 }
   for(e in 1:(nf-1)) {  
    if(obs[tt,v] > vec[e] & obs[tt,v] < vec[e+1]) {rnk_hista[tt,v]<-e+1 }
   }
  }
 }
 for(v in 1:4){
  figure<-paste(plotdir,"/RnkHist_",varid[v],".",site,".png",sep="") 
  writeLines(paste("Plotting:",figure))
  png(file=figure, width=800, height=800)
  par(cex.axis=1.5)
  par(cex.lab=1.5)
  hist(rnk_hista[,v], main =paste("Ranked Histogram of",varid[v],"at",site),freq=T,xlab="Rank")
  dev.off()
 }

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



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