#Toolkit version number
r.tool.name = "Model Evaluation"
r.eval.toolkit.version = "4.0"
tkm.version = "tkmfileversion1"
svn.version.eval = "$ProjectRevision: SVN 1942 (2016-01-29) $"
svn.prj.name.eval = "$ProjectName: /cerc/Products/Utilities/ModelEvaluationToolkit/Tags/ModelEvaluationToolkit_4_0_0 $"

#Suppress warnings that we don't want to worry the user with
options(warn=-1)

RHC.n = 26 # Robust highest concentration: number of data points that represent the highest values (tail) of the concentration distribution
# RHC = X(n) + (\bar{X} - X(n)) * ln( (3n-1)/2 ) 
#  where X are the (n-1) largest concentrations and X(n) is the n^th largest concentration

#Set up input and output
settings.eval.file = commandArgs(trailingOnly=TRUE)[[1]]
eval.output.root = sub(".tkm","",settings.eval.file,ignore.case=TRUE)
log.file = file(paste(eval.output.root,".log",sep=""),open="wt")
sink(file=log.file,type="output",split=TRUE)
sink(file=log.file,type="message")
cat(paste("Model Evaluation Toolkit version ",r.eval.toolkit.version,": ",r.tool.name," Tool \n",sep=""))
if (length(grep("$ProjectRevision: ",svn.version.eval,fixed=TRUE))>0) {
  svn.version.eval.numberonly=sub("$ProjectRevision: ","",svn.version.eval,fixed=TRUE)
  svn.version.eval.numberonly=sub(" $","",svn.version.eval.numberonly,fixed=TRUE)
  cat(paste("Internal version ",svn.version.eval.numberonly," \n",sep=""))
} else {
  cat("CERC INTERNAL USE ONLY\n")
}
cat(paste("Running in",R.version$version.string,"\n"))
cat(paste(date(),"\n\n",sep=""))

#Read in settings file
if (!file.exists(settings.eval.file)) {
   stop(paste("The .tkm file specified does not exist:",settings.eval.file))
   quit("no")   
}

#getting the data from tkm
tkm.data = readLines(settings.eval.file)
tkm.data = sub("PATH = \"","",as.character(tkm.data),fixed=TRUE)
tkm.data = sub("\"","",as.character(tkm.data),fixed=TRUE)

#Check that the tkm file is for the right version
tkm.version.in = as.character(tkm.data[1])
if (tkm.version.in != tkm.version){
  stop(paste("The TKM file version (",tkm.version.in,") does not match the expected version(",tkm.version,")",sep=""))
  quit("no")
}

#LOAD R WORKSPACE
nProject = which(tkm.data == "Project")
data.input.workspace = tkm.data[nProject+2]
cat("Loading workspace...\n")
if (data.input.workspace=="") {
   stop("No workspace selected.")
   quit("no")
}
if (grepl(",",data.input.workspace)) {
   stop("Workspace file path must not contain commas")
   quit("no")
}
data.loaded=tryCatch(load(data.input.workspace),error=function(e) return("error"))
if (data.loaded=="error"){
  stop("The workspace could not be imported.")
  quit("no")
}

#check R script compatibility
if (!exists("r.data.toolkit.version")){
  stop("The selected workspace is incompatible with this version of the Model Evaluation tool.")
quit("no")
} else if (exists("r.data.toolkit.version")){
  if (r.data.toolkit.version != r.eval.toolkit.version){
    stop("The selected workspace is incompatible with this version of the Model Evaluation tool.")
    quit("no")
  }
}

#Finding line numbers for the data
nProject = which(tkm.data == "Project")
nRSettings = which(tkm.data=="RSettings")
nModStatPoll = which(tkm.data == "Models, Stations and Pollutants")
nGroupFilt = which(tkm.data == "Grouping and Filtering")
nConcEval = which(tkm.data == "Concentration Evaluation")
nForeIndEval = which(tkm.data == "Forecast index Evaluation")
nOutput = which(tkm.data == "Output")

#PROJECT
project.name = tkm.data[nProject+3]
if (project.name==""){
  stop("No project name entered")
  quit("no")
}
if (grepl(",",project.name)) {
  stop("Project name must not contain commas")
  quit("no")
}

#R SETTINGS
mirror=tkm.data[nRSettings+2]
update.packages = tkm.data[nRSettings+3]=="Y"
mirror.data = tryCatch(url(mirror),error=function(e) return(NA))
if (is.na(mirror.data)) {
  stop("Invalid CRAN mirror web address. Please refer to the help menu for a complete list of available CRAN mirrors.")
  quit("no")
}

#OTHER DATA
bAllModels = tkm.data[nModStatPoll+2]=="0"
bSingleModel = tkm.data[nModStatPoll+2]=="1"
single.model = tkm.data[nModStatPoll+3]
if (bSingleModel & single.model=="") {
   stop("No model name entered.")
   quit("no")
}
bAllStations = tkm.data[nModStatPoll+4]=="0"
bSingleStationType = tkm.data[nModStatPoll+4]=="1"
single.station.type = tkm.data[nModStatPoll+5]
if (bSingleStationType & single.station.type=="") {
   stop("No station type entered.")
   quit("no")
}
bAllPols = tkm.data[nModStatPoll+6]=="0"
bSinglePollutant = tkm.data[nModStatPoll+6]=="1"
single.pollutant = tkm.data[nModStatPoll+7]
if (bSinglePollutant & single.pollutant=="") {
   stop("No pollutant name entered.")
   quit("no")
}

#GROUPING and FILTERING
Grouping.and.filtering = tkm.data[nGroupFilt+2]


if (Grouping.and.filtering == "0"){
  #group by station, filter by pollutant and model
  bGrpStnName = TRUE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = TRUE
  bFilterStnType = FALSE
}else if (Grouping.and.filtering == "1"){
  #group by station, filter by pollutant, station type and model
  bGrpStnName = TRUE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = TRUE
  bFilterStnType = TRUE
}else if (Grouping.and.filtering == "2"){
  #group by station, filter by model and station type
  bGrpStnName = TRUE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = FALSE
  bFilterStnType = TRUE
}else if (Grouping.and.filtering == "3"){
  #group by station name, filter by model
  bGrpStnName = TRUE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = FALSE
  bFilterStnType = FALSE
}else if (Grouping.and.filtering == "4"){
  #group by station type, filter by model and pollutant
  bGrpStnName = FALSE
  bGrpStnType = TRUE 
  bGrpPollutant = FALSE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = TRUE
  bFilterStnType = FALSE
}else if (Grouping.and.filtering == "5"){
  #group by station type, filter by model 
  bGrpStnName = FALSE
  bGrpStnType = TRUE 
  bGrpPollutant = FALSE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = FALSE
  bFilterStnType = FALSE
}else if (Grouping.and.filtering == "6"){
  #group by model, filter by pollutant and station type
  bGrpStnName = FALSE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = TRUE
  bFilterModel = FALSE
  bFilterPollutant = TRUE
  bFilterStnType = TRUE
}else if (Grouping.and.filtering == "7"){
  #group by model, filter by pollutant
  bGrpStnName = FALSE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = TRUE
  bFilterModel = FALSE
  bFilterPollutant = TRUE
  bFilterStnType = FALSE
}else if (Grouping.and.filtering == "8"){
  #group by model, filter by station type
  bGrpStnName = FALSE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = TRUE
  bFilterModel = FALSE
  bFilterPollutant = FALSE
  bFilterStnType = TRUE
}else if (Grouping.and.filtering == "9"){
  #group by model, no filtering
  bGrpStnName = FALSE
  bGrpStnType = FALSE 
  bGrpPollutant = FALSE
  bGrpModel = TRUE
  bFilterModel = FALSE
  bFilterPollutant = FALSE
  bFilterStnType = FALSE
}else if (Grouping.and.filtering == "10"){
  #group by pollutant, filter by model and station type
  bGrpStnName = FALSE
  bGrpStnType = FALSE 
  bGrpPollutant = TRUE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = FALSE
  bFilterStnType = TRUE
}else if (Grouping.and.filtering == "11"){
  #Group by pollutant, filter by model
  bGrpStnName = FALSE
  bGrpStnType = FALSE 
  bGrpPollutant = TRUE
  bGrpModel = FALSE
  bFilterModel = TRUE
  bFilterPollutant = FALSE
  bFilterStnType = FALSE
}

#CONCENTRATION EVALUATION
bAssessConcs = tkm.data[nConcEval+1]=="Y"
output.averaging.times.file = tkm.data[nConcEval+2]
data.threshold.conc = as.numeric(tkm.data[nConcEval+3])
bOpenairModStats = tkm.data[nConcEval+4]=="Y"
bBoxWhisker = tkm.data[nConcEval+5]=="Y"
bBoxOutliers = tkm.data[nConcEval+6]=="Y"
bBoxLog = tkm.data[nConcEval+7]=="Y"
bScatterPlot = tkm.data[nConcEval+8]=="Y" & tkm.data[nConcEval+10]=="1"
bFreqScatterPlot = tkm.data[nConcEval+8]=="Y" & tkm.data[nConcEval+10]=="0"
bScatterLog = tkm.data[nConcEval+9]=="Y"
bQQPlot = tkm.data[nConcEval+11]=="Y"
bQQLog = tkm.data[nConcEval+12]=="Y"
bTimePlot = tkm.data[nConcEval+13]=="Y"
bTargetPlot = tkm.data[nConcEval+14]=="Y"
bTargetPlot.DELTA.3.3 = tkm.data[nConcEval+15]=="Y"
TargetPlot.DELTA.3.3.uncertainty.file = tkm.data[nConcEval+16]
if (bAssessConcs) {
   #Check output averaging times file has been selected and that it contains no commas
   if (output.averaging.times.file=="") {
      stop("No output averaging times file selected.")
      quit("no")
   }
   if (grepl(",",output.averaging.times.file)) {
      stop("Output averaging times file path must not contain commas")
      quit("no")
   }

   # check some output has been selected
   if (!bTargetPlot & !bTargetPlot.DELTA.3.3 & !bBoxWhisker & !bFreqScatterPlot & !bScatterPlot & !bQQPlot & !bTimePlot) {
      stop("No concentration evaluation graph option selected.")
      quit("no")
   }
   # check compatibility of output options with group/filter options
   # note that (bGrpModel | bFilterModel) = TRUE earlier so don't need to repeat bFilterModel check
   if (bTargetPlot & ( (bGrpStnName & !bFilterPollutant) | (bGrpStnType & !bFilterPollutant) | (bGrpModel & !bFilterPollutant) )) {
      cat("   Warning: Grouping and/or filtering options not suitable for target plot (DELTA 1.2) - this plot will not be output\n")
      bTargetPlot = FALSE  
   }
   if (bTargetPlot.DELTA.3.3 & ( (bGrpStnName & !bFilterPollutant) | (bGrpStnType & !bFilterPollutant)  | (bGrpModel & !bFilterPollutant) )) {
      cat("   Warning: Grouping and/or filtering options not suitable for target plot (DELTA 3.3) - this plot will not be output\n")
      bTargetPlot.DELTA.3.3 = FALSE  
   }
   if (bBoxWhisker & ( (bGrpStnName & !bFilterPollutant) | (bGrpStnType & !bFilterPollutant) | (bGrpPollutant) | (bGrpModel & (!bFilterPollutant | bFilterStnType)) )) {
      cat("   Warning: Grouping and/or filtering options not suitable for box and whisker plot - this plot will not be output\n")
      bBoxWhisker = FALSE
   }
   if (bFreqScatterPlot & ( (bGrpStnName & !bFilterPollutant) | (bGrpStnType) | (bGrpPollutant) | (bGrpModel) )) {
      cat("   Warning: Grouping and/or filtering options not suitable for frequency scatter plot - this plot will not be output\n")
      bFreqScatterPlot = FALSE
   }
   if (bScatterPlot & ( (bGrpStnName & !bFilterPollutant) | (bGrpStnType) | (bGrpPollutant) | (bGrpModel) )) {
      cat("   Warning: Grouping and/or filtering options not suitable for scatter plot - this plot will not be output\n")
      bScatterPlot = FALSE
   }
   if (bQQPlot & ( (bGrpStnName & !bFilterPollutant) | (bGrpStnType) | (bGrpPollutant) | (bGrpModel) )) {
      cat("   Warning: Grouping and/or filtering options not suitable for quantile-quantile plot - this plot will not be output\n")
      bQQPlot = FALSE
   }
   if (bTimePlot & !(bGrpStnName & bFilterPollutant & bFilterModel & !bFilterStnType)) {
     cat("   Warning: Grouping and/or filtering options not suitable for time plot averaged over all stations - this plot will not be output\n")
     bTimePlot = FALSE
   }
   # check uncertainty file has been entered for DELTA 3.3 target plot
   if (bTargetPlot.DELTA.3.3) {
     if (TargetPlot.DELTA.3.3.uncertainty.file=="") {
      stop("No uncertainties file selected.")
      quit("no")
     }
     if (grepl(",",TargetPlot.DELTA.3.3.uncertainty.file)) {
      stop("Uncertainties file path must not contain commas")
      quit("no")
     }
   }
}

#FORECAST INDEX EVALUATION
bAssessForecasts = tkm.data[nForeIndEval+1]=="Y"
index.scale.file = tkm.data[nForeIndEval+2]
thresholds.file = tkm.data[nForeIndEval+3]
data.threshold.index = as.numeric(tkm.data[nForeIndEval+4])

if (bAssessForecasts) {
  if (index.scale.file=="") {
    stop("No index scale file selected.")
    quit("no")
  }
  
  if (thresholds.file=="") {
      stop("No alert threshold file selected.")
      quit("no")
   }
}

bForecastAccuracy = tkm.data[nForeIndEval+5]=="Y"
bOddsRatioSkillScore = tkm.data[nForeIndEval+6]=="Y"
bForecastPerformanceMetrics = tkm.data[nForeIndEval+7]=="Y"
if (bAssessForecasts & !bForecastAccuracy & !bOddsRatioSkillScore & !bForecastPerformanceMetrics) {
   stop("No forecast assessment graph output selected.")
   quit("no")
}
# regarding compatibility of output options with group/filter options
# note that (bGrpModel | bFilterModel) = TRUE earlier so don't need any other checks

#OUTPUT
output.types = c("JPG","PNG","PDF")
output.type = output.types[as.numeric(tkm.data[nOutput+2])+1]
output.selected = (output.type == output.types)
bJPG = output.selected[1]
bPNG = output.selected[2]
bPDF = output.selected[3]
output.paper.sizes = c("Large (A4)","Medium","Small")
output.paper.size = output.paper.sizes[as.numeric(tkm.data[nOutput+3])+1]
output.paper.selected = (output.paper.size == output.paper.sizes)
bPaperLarge = output.paper.selected[1]
bPaperMedium = output.paper.selected[2]
bPaperSmall = output.paper.selected[3]
boutput.eval.stats.as.csv = tkm.data[nOutput+4]=="Y"

if (!bAssessForecasts & !bAssessConcs) {
   stop("No output selected.")
   quit("no")
}

#re-define overwritten settings after data read-in
boutput.stats.as.csv = boutput.eval.stats.as.csv

#Load common functions
source_local("CommonFunctions.r")
#Check CommonFunctions.r is compatible with this script
if (!exists("r.common.functions.version")){
  stop("The file CommonFunctions.r is incompatible with this version of the Model Evaluation tool.")
  quit("no")
} else if (exists("r.common.functions.version")){
  if (r.common.functions.version != r.eval.toolkit.version){
    stop("The file CommonFunctions.r is incompatible with this version of the Model Evaluation tool.")
    quit("no")
  }
}

#FUNCTIONS
CRMSE = function(x) {sqrt(mean(((x$mod - mean(x$mod)) - (x$obs - mean(x$obs)))^2))}
StDev = function (x) {sqrt(mean((x - mean(x))^2))}
MBias = function(x) {mean(x$mod-x$obs)}
NMSE = function(x) {mean((x$mod - x$obs)^2) / (mean(x$obs)*mean(x$mod))}
cor.data = function(x) {cor(x$obs,x$mod)}
Fac2 = function(x) { y = subset(x, x$obs>0)
   Fac2 = nrow( subset(y, abs(y$mod/y$obs)>=0.5 & abs(y$mod/y$obs)<=2.0 )) / nrow(y) }
Fb = function(x) { (mean(x$mod) - mean(x$obs))/(0.5*(mean(x$obs)+mean(x$mod))) }
Fs = function(x) { (StDev(x$mod) - StDev(x$obs)) /(0.5*(StDev(x$obs)+StDev(x$mod))) }
RHC = function(x) { y = sort(x, decreasing=TRUE)
   RHC = y[RHC.n] + ( (mean(y[1:(RHC.n-1)])-y[RHC.n])*log((3*RHC.n-1)/2) )}
boxplot.stats.wrap = function(x) {boxplot.stats(x)$stats}
length.unique = function(x) {length(unique(x)[,1])}

#Functions for drawing a table
make.table = function(nr, nc) { 
  plot(c(0, nc*6 + 1), c(0, -(nr + 1)), 
  type="n", xlab="", ylab="", axes=FALSE)
}
get.r = function(i, nr) i-floor((i-1)/nr)*nr 
get.c = function(i, nr) 1+floor((i-1)/nr) 
draw.title = function(title, nc) { 
  text((nc*6 + 1)/2, 0, title) 
}
draw.sample.cell = function(id, key, i, nr) { 
  r = get.r(i, nr) 
  c = get.c(i, nr)
  text(6*(c - 1) + 2, -r, id, col="violetred",cex=0.7) 
  text(6*c-1, -r, key,cex=0.7) 
  rect(6*(c - 1) + .5, -(r - .5), 6*c + .5, -(r + .5), border="grey") 
}

#R CODE:

#Set embedded title for PDF files
pdf.title = paste("Output from DS-LOCAL Model Evaluation Toolkit Version",r.eval.toolkit.version)

#Having written the settings and loaded the data, check whether only a single 
# pollutant and / or model exists
modelled.data.names = names(all.data)[names(all.data)!="date" & names(all.data)!="obs" & names(all.data)!="station" & names(all.data)!="pollutant" & names(all.data)!="type"]
if (bAllModels & length(modelled.data.names)==1) {
   bAllModels = FALSE
   bSingleModel = TRUE
   single.model = as.character(modelled.data.names)
}
if (bAllPols & length(unique(pollutants$pollutant))==1) {
   bAllPols = FALSE
   bSinglePollutant = TRUE
   single.pollutant = as.character(unique(pollutants$pollutant))
}

#MODEL
#Check model name exists and reduce data set if required, prior to any further work
if (bSingleModel) {
   if (!any(single.model == modelled.data.names)) {
      models.available = paste(modelled.data.names,collapse=", ")
      stop(paste("No data for selected model ",single.model,". Data are available for the following models: ",models.available,"\n",sep=""))
      quit("no")
   } else {
      all.data = subset(all.data,select=c("date","obs","station","pollutant","type",single.model))
      modelled.data.names = single.model
   }
}

#POLLUTANTS
#Reduce to a subset of pollutants if required
if (bSinglePollutant) {
   pollutants.new = subset(pollutants,pollutant==tolower(single.pollutant))
   if (length(pollutants.new$pollutant)==0) {
      pollutants.available = toupper(paste(unique(pollutants$pollutant),collapse=", "))
      stop(paste("No data for selected pollutant ",toupper(single.pollutant),". Data are available for the following pollutants: ",pollutants.available,"\n",sep=""))
      quit("no")
   }   
   pollutants = pollutants.new
}
# output averaging time and statistic
# check whether data is to be averaged over different time period or statistic than the modelled data was output:
if (bAssessConcs) {
    # process output.averaging.times.file and append to pollutants data frame
    pollutants = process.output.averaging(output.averaging.times.file, pollutants)
    # store title strings for output averaging times and statistic
    pollutants$output.avg.time.and.statistic = vector(mode="character", length = length(pollutants$pollutant)) # preallocate
    for (j in 1:length(pollutants$pollutant)) {
      if (pollutants$output.avg.time.hours[j]==24) {
        if (pollutants$output.statistic[j]=="max") {
          pollutants$output.avg.time.and.statistic[j] = "daily maximum"
        } else {
          pollutants$output.avg.time.and.statistic[j] = paste("daily", pollutants$output.statistic[j], sep=" ")
        }
      } else if (pollutants$output.avg.time.hours[j]==1) {
        if (pollutants$output.statistic[j]=="max") {
          pollutants$output.avg.time.and.statistic[j] = "hourly maximum"
        } else {
          pollutants$output.avg.time.and.statistic[j] = paste("hourly", pollutants$output.statistic[j], sep=" ")  
        }
      } else {
        if (pollutants$output.statistic[j]=="max") {
          pollutants$output.avg.time.and.statistic[j] = paste(paste(pollutants$output.avg.time.hours[j],"-hour", sep=""), "maximum", sep=" ")
        } else {
          pollutants$output.avg.time.and.statistic[j] = paste(paste(pollutants$output.avg.time.hours[j],"-hour", sep=""),pollutants$output.statistic[j],sep=" ")
        }
      }
      if (pollutants$output.daily.max[j]=="yes") {
        pollutants$output.avg.time.and.statistic[j]=paste("daily maximum of",pollutants$output.avg.time.and.statistic[j])
      }
    }
}
#Read and process the uncertainty values file
if (bAssessConcs & bTargetPlot.DELTA.3.3) {
  Target.plot.uncertainties =  read.csv(file=TargetPlot.DELTA.3.3.uncertainty.file)
  names(Target.plot.uncertainties) = tolower(names(Target.plot.uncertainties))
  # make sure each column exists
  if (is.null(Target.plot.uncertainties$k) | is.null(Target.plot.uncertainties$ur) | is.null(Target.plot.uncertainties$lv)| is.null(Target.plot.uncertainties$alpha) | 
      is.null(Target.plot.uncertainties$target.units) | is.null(Target.plot.uncertainties$target.avg.time.hours) | is.null(Target.plot.uncertainties$target.statistic) ) {
    stop("You have chosen to create a target plot (DELTA 3.3) but you have not included all of the required information in the uncertainties file.")
    quit("no")
  }
  Target.plot.uncertainties$pollutant = tolower(as.character(Target.plot.uncertainties$pollutant))
  Target.plot.uncertainties$target.units = tolower(as.character(Target.plot.uncertainties$target.units))
  Target.plot.uncertainties$target.statistic = tolower(as.character(Target.plot.uncertainties$target.statistic)) 
 
 # Initialise plotting target plot to true
  pollutants$bTargetPlotDELTA3.3 = TRUE
  # check each pollutant in the pollutant data frame (which may be a single pollutant) is included in the user defined averaging file
  for (j in 1:length(pollutants$pollutant)) {   
    if (length(which(Target.plot.uncertainties$pollutant==pollutants$pollutant[j]))==0 ) {
      cat(paste("Pollutant",toupper(pollutants$pollutant[j]),"is not in the uncertainties file, no target plot (DELTA 3.3) will be produced for this pollutant.\n"))
      pollutants$bTargetPlotDELTA3.3[j] = FALSE
    }
    if (length(which(Target.plot.uncertainties$pollutant==pollutants$pollutant[j])) > 1 ) {
      stop(paste("Pollutant",toupper(pollutants$pollutant[j]),"is listed more than once in the uncertainties file."))
      quit("no")
     }
   }
   # now merge with pollutants file
   pollutants = merge(pollutants,Target.plot.uncertainties,by="pollutant",all.x=TRUE)
   # check units
  for (j in 1:length(pollutants$pollutant)) { 
    if (pollutants$bTargetPlotDELTA3.3[j]) {
      if (pollutants$output.units[j] != pollutants$target.units[j]) {
        stop(paste("Pollutant",toupper(pollutants$pollutant[j])," does not have the correct units in the output averaging times file for a target plot (DELTA 3.3)."))
        quit("no")
      }
      if (pollutants$output.avg.time.hours[j] != pollutants$target.avg.time.hours[j]) {
        stop(paste("Pollutant",toupper(pollutants$pollutant[j])," does not have the correct averaging time in the output averaging times file for a target plot (DELTA 3.3)."))
        quit("no")
      }
      if (pollutants$output.statistic[j] != pollutants$target.statistic[j]) {
        stop(paste("Pollutant",toupper(pollutants$pollutant[j])," does not have the correct statistic in the output averaging times file for a target plot (DELTA 3.3)."))
        quit("no")
      }
    }
  }
  # Check there is at least one pollutant to plot the target plot (DELTA 3.3) for
  if (all(!pollutants$bTargetPlotDELTA3.3)) {
    stop("No valid pollutants for creating target plot (DELTA 3.3)")
    quite("no")
  }
}
#Read and process the index scales file
if (bAssessForecasts) {
   
   #Read in index scales file and add to pollutants
   pollutants = process.index.scales(index.scale.file, pollutants)   
   
   #Find lower bounds of index
   index.lower.bounds = find.index.lower.bounds(pollutants)
   
}

#STATIONS
if (bSingleStationType) {
   stations.new = subset(stations,type==tolower(single.station.type))
   if (length(stations.new$name)==0) {
      station.types.available = toupper(paste(unique(stations$type),collapse=", "))
      stop(paste("No data for selected station type ",toupper(single.station.type),". Data are available for the following station types: ",station.types.available,"\n",sep=""))
      quit("no")
   }
   stations = stations.new   
}
# rename stations column "name" to "station" for subsequent merging 
names(stations)[names(stations)=="name"]="station"
stations.station = subset(stations,select=c("station","type"))
station.types.available = sort(as.character(unique(stations$type)))

#GROUPING AND FILTERING BY
#Amend plotting settings if only one value, for titles and layout using filtering options
if (bSingleModel) {
   bFilterModel = FALSE
}
if (bSinglePollutant) {
   bFilterPollutant = FALSE
}
# predefine some strings and variables 
if (bGrpStnName) {
   groupby.col.label = "station"
   groupby.col.openair = "site"
   groupby.key.string = "Station" 
   groupby.title.string = "ALL STATIONS"  
} else if (bGrpStnType) {
   groupby.col.label = "type"
   groupby.key.string = "Station type"
   groupby.col.openair = "type"
   groupby.title.string = "ALL STATION TYPES" 
} else if (bGrpPollutant) { 
   groupby.col.label = "pollutant"
   groupby.key.string = "Pollutant"
   groupby.col.openair = "pollutant"
   groupby.title.string = "ALL POLLUTANTS" 
} else if (bGrpModel) {
   groupby.col.label = "model"
   groupby.key.string = "Model"
   groupby.col.openair = "model"
   groupby.title.string = "ALL MODELS" 
}
if (bGrpModel | !bFilterModel) {
   if (!bFilterStnType & !bFilterPollutant) {
      # length of the plotting loop
      filterby.length = 1
      filterby.list = "all data"
      filterby.col.label = groupby.col.label  
   } else if (bFilterStnType & !bFilterPollutant) {
      # length of the plotting loop
      filterby.length = length(station.types.available)
      filterby.list = as.data.frame(station.types.available)
      names(filterby.list) = "type"
      filterby.col.label = "type"    
   } else if (bFilterPollutant & !bFilterStnType) { 
      # length of the plotting loop
      filterby.length = length(pollutants$pollutant)
      filterby.list = as.data.frame(as.character(pollutants$pollutant))
      names(filterby.list) = "pollutant"
      filterby.col.label = "pollutant"   
   } else if (bFilterPollutant & bFilterStnType) {
      # length of the plotting loop
      filterby.length = length(station.types.available) * length(pollutants$pollutant)
      filterby.list = as.data.frame(merge(station.types.available,as.character(pollutants$pollutant)))
      names(filterby.list) = c("type","pollutant")
      filterby.col.label = c("type","pollutant")
   }
} else { #bFilterModel
   if (!bFilterStnType & !bFilterPollutant) {
      # length of the plotting loop
      filterby.length = length(modelled.data.names)
      filterby.list = as.data.frame(as.character(modelled.data.names))
      names(filterby.list) = "model"
      filterby.col.label = "model"  
   } else if (bFilterStnType & !bFilterPollutant) {
      # length of the plotting loop
      filterby.length = length(modelled.data.names) * length(station.types.available)
      filterby.list = as.data.frame(merge(station.types.available,as.character(modelled.data.names)))
      names(filterby.list) = c("type","model")
      filterby.col.label = c("type","model")    
   } else if (bFilterPollutant & !bFilterStnType) { 
      # length of the plotting loop
      filterby.length = length(modelled.data.names) * length(pollutants$pollutant)
      filterby.list = as.data.frame(merge(as.character(pollutants$pollutant),as.character(modelled.data.names))) 
      names(filterby.list) = c("pollutant","model")
      filterby.col.label = c("pollutant","model")   
   } else if (bFilterPollutant & bFilterStnType) {
      # length of the plotting loop
      filterby.length = length(modelled.data.names) * length(station.types.available) * length(pollutants$pollutant)
      filterby.list = as.data.frame(merge(station.types.available,as.character(pollutants$pollutant)))
      names(filterby.list) = c("type","pollutant")
      filterby.list = merge(filterby.list,as.character(modelled.data.names))
      names(filterby.list) = c("type","pollutant","model")      
      filterby.col.label = c("type","pollutant","model")
   }
}

#READ AND PROCESS THRESHOLD LEVELS FILE (separator here is ",")
#header = "index, name"
if (bAssessForecasts) {
   cat("Importing alert threshold details...\n")
   thresholds = read.csv(thresholds.file,header=TRUE,sep=",")
} 

#Load the required packages
packages.required = c("openair","lattice","nlme","hexbin")
update.install.load.packages(packages.required,mirror,update.packages)

#Loop over model and pollutants
if (bSingleStationType) {
   all.data = subset(all.data,type==tolower(single.station.type))
}
station.types = as.data.frame(table(stations$type))
names(station.types)=c("type","num.stations")

#Remove single rows of data to stop rolling average errors
cat("Removing single rows of data for rolling averages to prevent errors...\n")
all.data$keep = TRUE
total.num.stns.removed = 0
for (j in 1:length(pollutants$pollutant)) {
  pol = pollutants$pollutant[j]
  num.stns.removed = 0
  bCheckRollingForConcs = FALSE
  bCheckRollingForFcsts = FALSE
  if (bAssessConcs) {
    if (pollutants$output.statistic[j]=="rolling mean" & 
        !(pollutants$mod.statistic[j]==pollutants$output.statistic[j] & pollutants$mod.avg.time.hours[j]==pollutants$output.avg.time.hours[j])) {
      bCheckRollingForConcs = TRUE
    }
  }
  if (bAssessForecasts) {
    if (pollutants$index.statistic[j]=="rolling mean" & 
        !(pollutants$mod.statistic[j]==pollutants$index.statistic[j] & pollutants$mod.avg.time.hours[j]==pollutants$index.avg.time.hours[j])) {
      bCheckRollingForFcsts = TRUE
    }
  }
  if (bCheckRollingForConcs | bCheckRollingForFcsts) {
    for (i in 1:length(stations$station)) {
      stn = stations$station[i]
      num.records = nrow(subset(all.data,pollutant==pol & station==stn))
      if (num.records==1) {
        all.data$keep[all.data$pollutant==pol & all.data$station==stn] = FALSE
        cat(paste("  Station ",stn," removed from dataset for ",toupper(pol),"\n",sep=""))
        num.stns.removed = num.stns.removed+1
      }
    }
    cat(paste("  ",num.stns.removed," station(s) removed from dataset for ",toupper(pol),"\n",sep=""))
    total.num.stns.removed = total.num.stns.removed + num.stns.removed
  }
}
all.data = subset(all.data,keep==TRUE)
all.data = all.data[-which(names(all.data)=="keep")]
cat(paste("  In total, ",total.num.stns.removed," station(s) removed from dataset \n",sep=""))

# Loop through each model and pollutant to initially compare obs and mod data
#all comparisons between obs and modelled data are by date-time / station / pollutant / model combination
bDoneFirstModelConc = FALSE
bDoneFirstModelIndex = FALSE
for (m in 1:length(modelled.data.names)) {
   cat(paste("Processing model ",modelled.data.names[m],"...\n",sep=""))
   #Get all the data for this model
   all.data.onemodel = subset(all.data,select=c("date","station","type","pollutant","obs",modelled.data.names[m]))
   names(all.data.onemodel)[names(all.data.onemodel)==modelled.data.names[m]] = "mod"
   all.data.onemodel$model = modelled.data.names[m]
   
   bDoneFirstPolIndex = FALSE
   bDoneFirstPolConc = FALSE
   for (j in 1:length(pollutants$pollutant)) {
      bNullIndex = FALSE
      bNullConc = FALSE
      #Get all the data for this pollutant
      all.data.onepol = subset(all.data.onemodel,pollutant==pollutants$pollutant[j])
      if (length(all.data.onepol$station)==0) next   
      cat(paste("   Processing pollutant ",toupper(pollutants$pollutant[j]),"...\n",sep=""))
      all.data.onepol = subset(all.data.onepol,select=c("date","station","type","obs","mod","model"))
      #rename 'station' to 'site' so that openair recognises it
      names(all.data.onepol)[names(all.data.onepol)=="station"] = "site"
      if (bAssessForecasts) {
         cat(paste("      Performing forecast assessment...\n",sep=""))
         
         #Obtain indices for data
         all.data.for.indexing = find.data.for.indexing(all.data.onepol,pollutants,j,index.lower.bounds,data.threshold.index)

         #EVALUATE THE MODEL INDEX FORECAST PERFORMANCE
         all.data.for.indexing$index.absdiff = abs(all.data.for.indexing$index.mod - all.data.for.indexing$index.obs)
   
         #add the station type and pollutant back into the data frame for completeness
         names(all.data.for.indexing)[names(all.data.for.indexing)=="site"] = "station"
         all.data.for.indexing = merge(all.data.for.indexing, stations[ ,c("station", "type")], by="station")
         names(all.data.for.indexing)[names(all.data.for.indexing)=="station"] = "site"
         all.data.for.indexing$pollutant = pollutants$pollutant[j]
         all.data.for.indexing$model = modelled.data.names[m]
   
         #Save the indexed data including missing data
         fcst.data.for.output.onepol = all.data.for.indexing
   
         #Remove missing data
         all.data.for.indexing = subset(all.data.for.indexing,!is.na(obs.output.units) & !is.na(mod.output.units))

         if (length(all.data.for.indexing$date)==0) {

            cat(paste("      Warning: No valid forecast index data for ",modelled.data.names[m]," ",toupper(pollutants$pollutant[j])," with ", data.threshold.index, "% data capture threshold\n",sep=""))
            bNullIndex = TRUE

         } else {
            #Determine the statistics about accuracy of index prediction for each station
            fcst.accuracy.onepol = data.frame(station = unique(all.data.for.indexing$site), pollutant = pollutants$pollutant[j], model=modelled.data.names[m])
            #Need station type for filtering or grouping later
            if (bFilterStnType | bGrpStnType) {
               fcst.accuracy.onepol = as.data.frame(merge(fcst.accuracy.onepol,stations[ ,c("station","type")], by="station"))
            }
            scores = table(subset(all.data.for.indexing,select=c("site","index.absdiff")),useNA="no")
            fcst.accuracy.onepol$num.valid.indices = rowSums(scores)
            if (any(dimnames(scores)[[2]]=="0")) {
               num.same = scores[,name="0"]
            } else {
               num.same = 0
            }
            if (any(dimnames(scores)[[2]]=="1")) {
               num.oneapart = scores[,name="1"]
            } else {
               num.oneapart = 0
            }
            fcst.accuracy.onepol$indices.same = num.same
            fcst.accuracy.onepol$indices.oneapart = num.oneapart
      
            #Find the alert level for both obs and mod for each station
            all.data.for.indexing$alert.level.obs=findInterval(all.data.for.indexing$index.obs,thresholds$index)
            all.data.for.indexing$alert.level.mod=findInterval(all.data.for.indexing$index.mod,thresholds$index)
      
            bDoneFirstAlert=FALSE
            for (k in 1:length(thresholds$name)) {
      
               #Work out whether an event was forecast AND observed, forecast but not observed etc
               all.data.for.indexing$mod.obs = 0
               all.data.for.indexing$mod.obs[all.data.for.indexing$alert.level.mod == k & all.data.for.indexing$alert.level.obs == k] = 1
               all.data.for.indexing$mod.noobs = 0
               all.data.for.indexing$mod.noobs[all.data.for.indexing$alert.level.mod == k & all.data.for.indexing$alert.level.obs != k] = 1
               all.data.for.indexing$nomod.obs = 0
               all.data.for.indexing$nomod.obs[all.data.for.indexing$alert.level.mod != k & all.data.for.indexing$alert.level.obs == k] = 1
               all.data.for.indexing$nomod.noobs = 0
               all.data.for.indexing$nomod.noobs[all.data.for.indexing$alert.level.mod != k & all.data.for.indexing$alert.level.obs != k] = 1
      
               #Work out the components of the odds ratio skill score
               fcst.onealert = data.frame(alert.name=thresholds$name[k],pollutant=pollutants$pollutant[j],station=unique(all.data.for.indexing$site),model=modelled.data.names[m])
               if (any(all.data.for.indexing$alert.level.obs==k)) {
                  fcst.onealert$num.obs.alerts = table(subset(all.data.for.indexing,select=c("site","alert.level.obs")))[,name=as.character(k)]
               } else {
                  fcst.onealert$num.obs.alerts = 0
               }
               if (any(all.data.for.indexing$mod.obs==1)) {
                  fcst.onealert$a = table(subset(all.data.for.indexing,select=c("site","mod.obs")))[,name="1"]
               } else {
                  fcst.onealert$a = 0
               }
               if (any(all.data.for.indexing$mod.noobs==1)) {
                  fcst.onealert$b = table(subset(all.data.for.indexing,select=c("site","mod.noobs")))[,name="1"]
               } else {
                  fcst.onealert$b = 0
               }
               if (any(all.data.for.indexing$nomod.obs==1)) {
                  fcst.onealert$c = table(subset(all.data.for.indexing,select=c("site","nomod.obs")))[,name="1"]
               } else {
                  fcst.onealert$c = 0
               }
               if (any(all.data.for.indexing$nomod.noobs==1)) {
                  fcst.onealert$d = table(subset(all.data.for.indexing,select=c("site","nomod.noobs")))[,name="1"]
               } else {
                  fcst.onealert$d = 0
               }
      
               if (bDoneFirstAlert) {
                  fcst.onepol = rbind(fcst.onepol,fcst.onealert)
               } else {
                  fcst.onepol = fcst.onealert
                  bDoneFirstAlert = TRUE
               }

            } #End of loop over thresholds
      
            #Need station type for filtering or grouping later
            if (bFilterStnType | bGrpStnType) {
               fcst.onepol = as.data.frame(merge(fcst.onepol,stations[ ,c("station","type")], by="station"))
            }

         } #End of check for null data   

      } #End of forecast assessment code
                  
      if (bAssessConcs) {
         cat(paste("      Performing concentration assessment...\n",sep=""))
   
         # ASSESS CONCS AVERAGING
         # we want the output averaging time and statistic for the data in the concentration evaluation,
         # for the statistics and the plots.
         if ( (pollutants$mod.avg.time.hours[j]==pollutants$output.avg.time.hours[j]) & identical(pollutants$mod.statistic[j],pollutants$output.statistic[j])) {
             # no change to data
             all.data.onepol.averaged = all.data.onepol
         } else {
             all.data.onepol.averaged = average.data(all.data.onepol, pollutants, j, stations.station, data.threshold.conc) # function defined in DataInput.rgg
             all.data.onepol.averaged$model = modelled.data.names[m]
         } 

         # keep averaged data for plots.  NB includes all data inc NAs.
         all.data.onepol.averaged$pollutant = pollutants$pollutant[j]
         all.data.onepol.averaged$output.avg.time.and.statistic = pollutants$output.avg.time.and.statistic[j]

         #remove all lines where one or other or both of obs and mod is missing
         all.data.onepol = subset(all.data.onepol.averaged,!is.na(obs) & !is.na(mod))

         if (length(all.data.onepol$site)==0) {

            cat(paste("      Warning: No valid concentration data for ",modelled.data.names[m]," ",toupper(pollutants$pollutant[j])," with ", data.threshold.conc, "% data capture threshold\n",sep=""))
            bNullConc = TRUE

         } else {

            #store the data inc. NAs for plotting - revert back to columns for model labels
            names(all.data.onepol.averaged)[names(all.data.onepol.averaged)=="mod"] = modelled.data.names[m]
            all.data.onepol.averaged = subset(all.data.onepol.averaged, select=names(all.data.onepol.averaged)[names(all.data.onepol.averaged)!="model"])
            # if not first model, append other model data on to data frame, 
            # for each pollutant; can then rbind at end 
            if (bDoneFirstModelConc) {
               all.data.onepol.averaged = merge(all.data.onepol.averaged, all.data.averaged[ all.data.averaged$pollutant==pollutants$pollutant[j], ], all.y=TRUE)   
            }

            #Calculate statistics
            # define the columns needed by the grouping and filtering option chosen
            # pollutant and model are singular (subsetted) and type relates to site so want just one column:
            if (bGrpStnName) {
               conc.col.string = groupby.col.openair # "site"
            } else {
               if (bGrpStnType | bFilterStnType) {
                  conc.col.string = "type"         
               } else {
                  conc.col.string = "pollutant"         
               }
            }      
      
            # allocate the identifiers in the first columns
            if (bGrpStnName) {
               conc.onepol.grouped = data.frame(station=sort(unique(all.data.onepol$site)), pollutant=pollutants$pollutant[j], output.avg.time.and.statistic=pollutants$output.avg.time.and.statistic[j], model=modelled.data.names[m])         
               conc.onepol.grouped$station = as.character(conc.onepol.grouped$station) # 'station' as character to allow row for "all" later
               if (bFilterStnType) {
                  conc.onepol.grouped = as.data.frame(merge(stations[ ,c("station","type")], conc.onepol.grouped, by="station"))
               } 
            } else {
               if (bGrpStnType | bFilterStnType) {
                  conc.onepol.grouped = data.frame(type=sort(unique(all.data.onepol$type)), pollutant=pollutants$pollutant[j], output.avg.time.and.statistic=pollutants$output.avg.time.and.statistic[j], model=modelled.data.names[m])         
                  conc.onepol.grouped$type = as.character(conc.onepol.grouped$type) # 'type' as character to allow row for "all" later
               } else {
                  conc.onepol.grouped = data.frame(pollutant=pollutants$pollutant[j], output.avg.time.and.statistic=pollutants$output.avg.time.and.statistic[j], model=modelled.data.names[m])         
               }
            } 
            conc.onepol.grouped$model = as.character(conc.onepol.grouped$model) # 'model' as character to allow row for "all" later
            # then apply validation / BOOT-like statistics over the group/filter combinations
            conc.onepol.grouped$num.valid.values = as.integer(table(all.data.onepol[ ,c(conc.col.string)])) #nb table sorts alphabetically
            conc.onepol.grouped$obs.mean = as.numeric(t(as.data.frame(lapply(split(all.data.onepol$obs, all.data.onepol[,conc.col.string]), FUN=mean))))
            conc.onepol.grouped$mod.mean = as.numeric(t(as.data.frame(lapply(split(all.data.onepol$mod, all.data.onepol[,conc.col.string]), FUN=mean))))
            conc.onepol.grouped$SDO = as.numeric(t(as.data.frame(lapply(split(all.data.onepol$obs, all.data.onepol[,conc.col.string]), FUN=StDev))))
            conc.onepol.grouped$SDM = as.numeric(t(as.data.frame(lapply(split(all.data.onepol$mod, all.data.onepol[,conc.col.string]), FUN=StDev))))
            conc.onepol.grouped$MB= gapply(all.data.onepol,FUN=MBias,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$NMSE = gapply(all.data.onepol,FUN=NMSE,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$R = gapply(all.data.onepol,FUN=cor.data,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$Fac2 = gapply(all.data.onepol,FUN=Fac2,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$Fb = gapply(all.data.onepol,FUN=Fb,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$Fs = gapply(all.data.onepol,FUN=Fs,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$obs.max = gapply(all.data.onepol,which="obs",FUN=max,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$mod.max = gapply(all.data.onepol,which="mod",FUN=max,groups=all.data.onepol[,conc.col.string])
            conc.onepol.grouped$obs.RHC = tapply(all.data.onepol$obs, all.data.onepol[,conc.col.string], RHC)
            conc.onepol.grouped$mod.RHC = tapply(all.data.onepol$mod, all.data.onepol[,conc.col.string], RHC)  
            #Calculate box and whisker plot statistics
            #Observed statistics
            obs.data.box.stats = t(as.data.frame(lapply(split(all.data.onepol[,"obs"],all.data.onepol[,conc.col.string]),FUN=boxplot.stats.wrap)))
            conc.onepol.grouped$obs.lower.whisker=obs.data.box.stats[,1]
            conc.onepol.grouped$obs.lower.quartile=obs.data.box.stats[,2]
            conc.onepol.grouped$obs.median=obs.data.box.stats[,3]
            conc.onepol.grouped$obs.upper.quartile=obs.data.box.stats[,4]
            conc.onepol.grouped$obs.upper.whisker=obs.data.box.stats[,5]
            #Modelled statistics
            mod.data.box.stats = t(as.data.frame(lapply(split(all.data.onepol[,"mod"],all.data.onepol[,conc.col.string]),FUN=boxplot.stats.wrap)))
            conc.onepol.grouped$mod.lower.whisker=mod.data.box.stats[,1]
            conc.onepol.grouped$mod.lower.quartile=mod.data.box.stats[,2]
            conc.onepol.grouped$mod.median=mod.data.box.stats[,3]
            conc.onepol.grouped$mod.upper.quartile=mod.data.box.stats[,4]
            conc.onepol.grouped$mod.upper.whisker=mod.data.box.stats[,5]
            # target plot statistics
            conc.onepol.grouped$CRMSE = gapply(all.data.onepol,FUN=CRMSE,groups=all.data.onepol[,conc.col.string])
            if (bTargetPlot.DELTA.3.3) {
              if (pollutants$bTargetPlotDELTA3.3[j]) {
                conc.onepol.grouped$RMSu = pollutants$k[j]*pollutants$ur[j]*sqrt((1-pollutants$alpha[j])*(conc.onepol.grouped$obs.mean^2+conc.onepol.grouped$SDO^2)+pollutants$alpha[j]*pollutants$lv[j]^2)
                conc.onepol.grouped$NMSD = (conc.onepol.grouped$SDM - conc.onepol.grouped$SDO)/conc.onepol.grouped$SDO
              } else {
                conc.onepol.grouped$RMSu = NA
                conc.onepol.grouped$NMSD = NA
              }
            }
                 
            # and create a row for "all" data for each pollutant - i.e. over the whole model
            if (!all(c(groupby.col.label,filterby.col.label) == "pollutant" | c(groupby.col.label,filterby.col.label) == "model")) { # ignore if only grouping by pollutant or model with same filters
               end.row = length(conc.onepol.grouped[,1])      
               end.row.data = conc.onepol.grouped[1,] # preallocate
               end.row.data[1, ] = rep(NA,length(end.row.data[1,])) # preallocate with NAs
               conc.onepol.grouped = rbind(conc.onepol.grouped, end.row.data) # preallocate with NAs
               if (any(names(conc.onepol.grouped)=="station")) {
                  conc.onepol.grouped[end.row+1, "station"] = "all"   # can't use conc.col.string as use "site" in all.data.onepol not "station"
               } 
               if (any(names(conc.onepol.grouped)=="type")) {
                  conc.onepol.grouped[end.row+1, "type"] = "all"          
               }     
               conc.onepol.grouped[end.row+1, c("pollutant","output.avg.time.and.statistic","model")] = c(pollutants$pollutant[j],pollutants$output.avg.time.and.statistic[j],modelled.data.names[m])
               conc.onepol.grouped[end.row+1,"num.valid.values"] = length(all.data.onepol$obs)
               conc.onepol.grouped[end.row+1,"obs.mean"] = mean(all.data.onepol$obs)
               conc.onepol.grouped[end.row+1,"mod.mean"] = mean(all.data.onepol$mod)
               conc.onepol.grouped[end.row+1,"SDO"] = StDev(all.data.onepol$obs)
               conc.onepol.grouped[end.row+1,"SDM"] = StDev(all.data.onepol$mod)
               conc.onepol.grouped[end.row+1,"MB"] = MBias(all.data.onepol)
               conc.onepol.grouped[end.row+1,"NMSE"] = NMSE(all.data.onepol)
               conc.onepol.grouped[end.row+1,"R"] = cor.data(all.data.onepol)
               conc.onepol.grouped[end.row+1,"Fac2"] = Fac2(all.data.onepol)
               conc.onepol.grouped[end.row+1,"Fb"] = Fb(all.data.onepol)
               conc.onepol.grouped[end.row+1,"Fs"] = Fs(all.data.onepol)
               conc.onepol.grouped[end.row+1,"obs.max"] = max(all.data.onepol$obs)
               conc.onepol.grouped[end.row+1,"mod.max"] = max(all.data.onepol$mod)
               conc.onepol.grouped[end.row+1,"obs.RHC"] = RHC(all.data.onepol$obs)
               conc.onepol.grouped[end.row+1,"mod.RHC"] = RHC(all.data.onepol$mod)
               obs.data.box.stats = boxplot.stats.wrap(all.data.onepol$obs)
               conc.onepol.grouped[end.row+1,"obs.lower.whisker"]=obs.data.box.stats[1]
               conc.onepol.grouped[end.row+1,"obs.lower.quartile"]=obs.data.box.stats[2]
               conc.onepol.grouped[end.row+1,"obs.median"]=obs.data.box.stats[3]
               conc.onepol.grouped[end.row+1,"obs.upper.quartile"]=obs.data.box.stats[4]
               conc.onepol.grouped[end.row+1,"obs.upper.whisker"]=obs.data.box.stats[5]
               mod.data.box.stats = boxplot.stats.wrap(all.data.onepol$mod)
               conc.onepol.grouped[end.row+1,"mod.lower.whisker"]=mod.data.box.stats[1]
               conc.onepol.grouped[end.row+1,"mod.lower.quartile"]=mod.data.box.stats[2]
               conc.onepol.grouped[end.row+1,"mod.median"]=mod.data.box.stats[3]
               conc.onepol.grouped[end.row+1,"mod.upper.quartile"]=mod.data.box.stats[4]
               conc.onepol.grouped[end.row+1,"mod.upper.whisker"]=mod.data.box.stats[5]
               conc.onepol.grouped[end.row+1,"CRMSE"] = CRMSE(all.data.onepol)
               if (bTargetPlot.DELTA.3.3) {
                 if (pollutants$bTargetPlotDELTA3.3[j]) {
                   conc.onepol.grouped[end.row+1,"RMSu"] = pollutants$k[j]*pollutants$ur[j]*sqrt((1-pollutants$alpha[j])*(mean(all.data.onepol$obs)^2+StDev(all.data.onepol$obs)^2)+pollutants$alpha[j]*pollutants$lv[j]^2)
                   conc.onepol.grouped[end.row+1,"NMSD"] = (StDev(all.data.onepol$mod) - StDev(all.data.onepol$obs))/StDev(all.data.onepol$obs)
                 } else {
                   conc.onepol.grouped[end.row+1,"RMSu"] = NA
                   conc.onepol.grouped[end.row+1,"NMSD"] = NA
                 }
               }
            }
                  
            # and add target plot stats (applies to every row)
            conc.onepol.grouped$CRMSE.sign.DELTA.1.2[conc.onepol.grouped$SDO > conc.onepol.grouped$SDM]= conc.onepol.grouped$CRMSE[conc.onepol.grouped$SDO > conc.onepol.grouped$SDM]
            conc.onepol.grouped$CRMSE.sign.DELTA.1.2[conc.onepol.grouped$SDO < conc.onepol.grouped$SDM]= -conc.onepol.grouped$CRMSE[conc.onepol.grouped$SDO < conc.onepol.grouped$SDM]
            conc.onepol.grouped$CRMSE.sign.DELTA.1.2[conc.onepol.grouped$SDO == conc.onepol.grouped$SDM]= 0
            conc.onepol.grouped$T.DELTA.1.2 = sqrt((conc.onepol.grouped$MB/conc.onepol.grouped$SDO)^2+(conc.onepol.grouped$CRMSE/conc.onepol.grouped$SDO)^2)
            if (bTargetPlot.DELTA.3.3) {
              if (pollutants$bTargetPlotDELTA3.3[j]) {
                conc.onepol.grouped$CRMSE.sign.DELTA.3.3 = sign(abs(conc.onepol.grouped$NMSD) - sqrt(2*(1-conc.onepol.grouped$R)))*conc.onepol.grouped$CRMSE
                conc.onepol.grouped$T.DELTA.3.3 = sqrt((conc.onepol.grouped$MB/(2*conc.onepol.grouped$RMSu))^2+(conc.onepol.grouped$CRMSE/(2*conc.onepol.grouped$RMSu))^2)
              } else {
                conc.onepol.grouped$CRMSE.sign.DELTA.3.3 = NA
                conc.onepol.grouped$T.DELTA.3.3 = NA
              }
            }
         }  #End of check for null concentration data
         
      } #End of concentration assessment code
      
      if (bAssessConcs & !bNullConc) {
         #Store summary information for each pollutant   
         summary.onepol.conc = station.types[order(station.types$type), ]
         summary.onepol.conc$pollutant = pollutants$pollutant[j]
         summary.onepol.conc$model = modelled.data.names[m]
         temp = as.data.frame(gapply(all.data.onepol,which="site",FUN=length.unique,groups=all.data.onepol$type))
         names(temp)="num.valid.stations"
         temp$type = row.names(temp)
         summary.onepol.conc = merge(summary.onepol.conc,temp,by="type",all.x=TRUE)
         if (bDoneFirstPolConc) {
            summary.by.pol.and.type.one.model.conc = rbind(summary.by.pol.and.type.one.model.conc,summary.onepol.conc)
            all.data.averaged.one.model = rbind(all.data.averaged.one.model, all.data.onepol.averaged)
            rm(all.data.onepol.averaged, all.data.onepol)
            conc.grouped.one.model = rbind(conc.grouped.one.model,conc.onepol.grouped)
         } else {
            summary.by.pol.and.type.one.model.conc = summary.onepol.conc
            # keep averaged data for plots
            # NB "site" not "station"
            all.data.averaged.one.model = all.data.onepol.averaged 
            rm(all.data.onepol.averaged, all.data.onepol)
            # and concentration summary information by station and for all
            conc.grouped.one.model = conc.onepol.grouped
            bDoneFirstPolConc = TRUE
         }
      }
      if (bAssessForecasts & !bNullIndex) {
         #Store summary information for each pollutant   
         summary.onepol.index = station.types[order(station.types$type), ]
         summary.onepol.index$pollutant = pollutants$pollutant[j]
         summary.onepol.index$model = modelled.data.names[m]
         temp.data =  merge(fcst.onepol[!is.na(fcst.onepol$num.obs.alerts),c("station","num.obs.alerts")], stations)
         temp = as.data.frame(gapply(temp.data,which="station",FUN=length.unique,groups=temp.data$type))
         names(temp)="num.valid.stations"
         temp$type = row.names(temp)
         summary.onepol.index = merge(summary.onepol.index,temp,by="type",all.x=TRUE)
         if (bDoneFirstPolIndex) {
            summary.by.pol.and.type.one.model.index = rbind(summary.by.pol.and.type.one.model.index,summary.onepol.index)
            fcst.allstns.one.model = rbind(fcst.allstns.one.model,fcst.onepol)
            fcst.data.for.output.one.model = rbind(fcst.data.for.output.one.model,fcst.data.for.output.onepol)
            fcst.accuracy.allstns.one.model = rbind(fcst.accuracy.allstns.one.model,fcst.accuracy.onepol)
         } else {
            summary.by.pol.and.type.one.model.index = summary.onepol.index
            fcst.allstns.one.model = fcst.onepol
            fcst.data.for.output.one.model = fcst.data.for.output.onepol
            fcst.accuracy.allstns.one.model = fcst.accuracy.onepol
            bDoneFirstPolIndex = TRUE
         }
      }                  
         
   } #End of loop over pollutants

   if (bAssessConcs) {
      if (bDoneFirstPolConc) {
         if (bDoneFirstModelConc) {
            summary.by.pol.and.type.conc = rbind(summary.by.pol.and.type.conc,summary.by.pol.and.type.one.model.conc)
            all.data.averaged = all.data.averaged.one.model
            rm(all.data.averaged.one.model, all.data.onemodel) 
            conc.grouped = rbind(conc.grouped,conc.grouped.one.model)
         } else {
            summary.by.pol.and.type.conc = summary.by.pol.and.type.one.model.conc
            # keep averaged data for plots
            # NB "site" not "station"
            all.data.averaged = all.data.averaged.one.model
            rm(all.data.averaged.one.model, all.data.onemodel) 
            # and concentration summary information by station and for all
            conc.grouped = conc.grouped.one.model
            bDoneFirstModelConc = TRUE
         }
      }
   }
   if (bAssessForecasts) {
      if (bDoneFirstPolIndex) {
         if (bDoneFirstModelIndex) {
            summary.by.pol.and.type.index = rbind(summary.by.pol.and.type.index,summary.by.pol.and.type.one.model.index)
            fcst.allstns = rbind(fcst.allstns,fcst.allstns.one.model)
            fcst.data.for.output = rbind(fcst.data.for.output,fcst.data.for.output.one.model)
            fcst.accuracy.allstns = rbind(fcst.accuracy.allstns,fcst.accuracy.allstns.one.model)
         } else {
            summary.by.pol.and.type.index = summary.by.pol.and.type.one.model.index
            fcst.allstns = fcst.allstns.one.model
            fcst.data.for.output = fcst.data.for.output.one.model
            fcst.accuracy.allstns = fcst.accuracy.allstns.one.model
            bDoneFirstModelIndex = TRUE
         }
      }
   }
           
} #end of loop over models

if (bAssessConcs & !bDoneFirstModelConc) {
   cat(paste("Warning: No valid concentration data with ", data.threshold.conc, "% data capture threshold - no concentration evaluation plots will be output\n",sep=""))
   bAssessConcs = FALSE
}

if (bAssessForecasts & !bDoneFirstModelIndex) {
   cat(paste("Warning: No valid forecast index data with ", data.threshold.index, "% data capture threshold - no forecast index evaluation plots will be output\n",sep=""))
   bAssessForecasts = FALSE
}

if (!bAssessConcs & !bAssessForecasts) {
   stop("No valid data to graph - stopping")
   quit("no")
}

# once collated for all pollutants, make names(all.data.averaged)==names(all.data)
if (bAssessConcs) {
   names(all.data.averaged)[names(all.data.averaged)=="site"] = "station"
}

#Get summary data by filter variable for plots
#create summary data frame
if (bAssessConcs) {
   summary.by.filter.conc = as.data.frame(filterby.list)
   names(summary.by.filter.conc) = filterby.col.label
}
if (bAssessForecasts) {
   summary.by.filter.index = as.data.frame(filterby.list)
   names(summary.by.filter.index) = filterby.col.label
}
#When grouping by model with no filters or single pollutant, single model chosen
# need to amend the filterby.col.label column in order to 'group' over all data
if (filterby.col.label == groupby.col.label) {
   if (bAssessConcs) {summary.by.pol.and.type.conc$station = rep("all data", length(summary.by.pol.and.type.conc[,1]))}
   if (bAssessForecasts) {summary.by.pol.and.type.index$station = rep("all data", length(summary.by.pol.and.type.index[,1]))}
   filterby.col.label = "station"
}
# picking a single pollutant removes the filtering by pollutant, which is needed for the box plot
if (bSinglePollutant) {
   if (bAssessConcs) {summary.by.filter.conc$pollutant = tolower(single.pollutant)}
   if (bAssessForecasts) {summary.by.filter.index$pollutant = tolower(single.pollutant)}
}
if (bGrpStnName) {
   # create a list of the variables to be applied to group in gapply
   if (filterby.list == list("all data")) {
      if (bAssessConcs) {filterby.data.as.list.conc = list(rep("all data",length(summary.by.pol.and.type.conc[,1])))}
      if (bAssessForecasts) {filterby.data.as.list.index = list(rep("all data",length(summary.by.pol.and.type.index[,1])))}
   } else {
      if (length(filterby.list) == 1) {
         if (bAssessConcs) {filterby.data.as.list.conc = eval(parse(text=paste("list(summary.by.pol.and.type.conc$",filterby.col.label[1],")",sep="")))}
         if (bAssessForecasts) {filterby.data.as.list.index = eval(parse(text=paste("list(summary.by.pol.and.type.index$",filterby.col.label[1],")",sep="")))}
      } else if (length(filterby.list)==2) {
         if (bAssessConcs) {filterby.data.as.list.conc = eval(parse(text=paste("list(summary.by.pol.and.type.conc$",filterby.col.label[1],",summary.by.pol.and.type.conc$",filterby.col.label[2],")",sep="")))}
         if (bAssessForecasts) {filterby.data.as.list.index = eval(parse(text=paste("list(summary.by.pol.and.type.index$",filterby.col.label[1],",summary.by.pol.and.type.index$",filterby.col.label[2],")",sep="")))}
      } else if (length(filterby.list)==3) {
         if (bAssessConcs) {filterby.data.as.list.conc = eval(parse(text=paste("list(summary.by.pol.and.type.conc$",filterby.col.label[1],",summary.by.pol.and.type.conc$",filterby.col.label[2],",summary.by.pol.and.type.conc$",filterby.col.label[3],")",sep="")))}
         if (bAssessForecasts) {filterby.data.as.list.index = eval(parse(text=paste("list(summary.by.pol.and.type.index$",filterby.col.label[1],",summary.by.pol.and.type.index$",filterby.col.label[2],",summary.by.pol.and.type.index$",filterby.col.label[3],")",sep="")))}
      }
   }
   #apply gapply to summary.by.pol.and.type depending on filters
   if (!bFilterStnType) {
      #collate all station types together so SUM
      if (bAssessConcs) {temp.b.conc = aggregate(summary.by.pol.and.type.conc$num.stations, by=filterby.data.as.list.conc,FUN=sum,na.rm=TRUE)} 
      if (bAssessForecasts) {temp.b.index = aggregate(summary.by.pol.and.type.index$num.stations, by=filterby.data.as.list.index,FUN=sum,na.rm=TRUE)} 
      if (bAssessConcs) {temp.a.conc = aggregate(summary.by.pol.and.type.conc$num.valid.stations, by=filterby.data.as.list.conc,FUN=sum,na.rm=TRUE)} 
      if (bAssessForecasts) {temp.a.index = aggregate(summary.by.pol.and.type.index$num.valid.stations, by=filterby.data.as.list.index,FUN=sum,na.rm=TRUE)} 
   } else {
      #find maximums in each station type or model, over pollutant/model groups
      if (bAssessConcs) {temp.b.conc = aggregate(summary.by.pol.and.type.conc$num.stations, by=filterby.data.as.list.conc,FUN=max,na.rm=TRUE)}
      if (bAssessForecasts) {temp.b.index = aggregate(summary.by.pol.and.type.index$num.stations, by=filterby.data.as.list.index,FUN=max,na.rm=TRUE)}
      if (bAssessConcs) {temp.a.conc = aggregate(summary.by.pol.and.type.conc$num.valid.stations, by=filterby.data.as.list.conc,FUN=max,na.rm=TRUE)}
      if (bAssessForecasts) {temp.a.index = aggregate(summary.by.pol.and.type.index$num.valid.stations, by=filterby.data.as.list.index,FUN=max,na.rm=TRUE)}
   }
   if (bAssessConcs) {
      names(temp.b.conc) = c(filterby.col.label,"num.stations")
      names(temp.a.conc) = c(filterby.col.label,"num.valid.stations")
      summary.by.filter.conc = merge(summary.by.filter.conc, temp.a.conc, by=filterby.col.label)
      summary.by.filter.conc = merge(summary.by.filter.conc, temp.b.conc, by=filterby.col.label)
      # re-order to filterby.list (alphabetical) order following merge
      for (j in 1:length(filterby.list)) {
         summary.by.filter.conc = summary.by.filter.conc[order(summary.by.filter.conc[, filterby.col.label[j]]), ]
      }
   }
   if (bAssessForecasts) {
      names(temp.b.index) = c(filterby.col.label,"num.stations")
      names(temp.a.index) = c(filterby.col.label,"num.valid.stations")
      summary.by.filter.index = merge(summary.by.filter.index, temp.a.index, by=filterby.col.label)
      summary.by.filter.index = merge(summary.by.filter.index, temp.b.index, by=filterby.col.label)
      # re-order to filterby.list (alphabetical) order following merge
      for (j in 1:length(filterby.list)) {
         summary.by.filter.index = summary.by.filter.index[order(summary.by.filter.index[, filterby.col.label[j]]), ]
      }
   }
} else if (bGrpStnType) {
   if (bAssessConcs) {
      #need to remove NA lines before finding length of unique 'type'
      summary.by.pol.and.type.conc = summary.by.pol.and.type.conc[ !is.na(summary.by.pol.and.type.conc$num.valid.stations), ]
      summary.by.filter.conc[ ,c("num.valid.stations")] = as.data.frame(gapply(summary.by.pol.and.type.conc, which="type",FUN=length.unique, groups=summary.by.pol.and.type.conc[ ,filterby.col.label]))
      summary.by.filter.conc[ ,c("num.stations")] = as.data.frame(rep(length(station.types.available), filterby.length))
   }
   if (bAssessForecasts) {
      #need to remove NA lines before finding length of unique 'type'
      summary.by.pol.and.type.index = summary.by.pol.and.type.index[ !is.na(summary.by.pol.and.type.index$num.valid.stations), ]
      summary.by.filter.index[ ,c("num.valid.stations")] = as.data.frame(gapply(summary.by.pol.and.type.index, which="type",FUN=length.unique, groups=summary.by.pol.and.type.index[ ,filterby.col.label]))
      summary.by.filter.index[ ,c("num.stations")] = as.data.frame(rep(length(station.types.available), filterby.length))
   }
} else if (bGrpPollutant) {
   if (bAssessConcs) {
      #need to remove NA lines before finding length of unique 'pollutant'
      summary.by.pol.and.type.conc = summary.by.pol.and.type.conc[ !is.na(summary.by.pol.and.type.conc$num.valid.stations), ]
      summary.by.filter.conc[ ,c("num.valid.stations")] = as.data.frame(gapply(summary.by.pol.and.type.conc, which="pollutant",FUN=length.unique, groups=summary.by.pol.and.type.conc[ ,filterby.col.label]))
      summary.by.filter.conc[ ,c("num.stations")] = as.data.frame(rep(length(unique(pollutants$pollutant)), filterby.length))
   }
   if (bAssessForecasts) {
      #need to remove NA lines before finding length of unique 'pollutant'
      summary.by.pol.and.type.index = summary.by.pol.and.type.index[ !is.na(summary.by.pol.and.type.index$num.valid.stations), ]
      summary.by.filter.index[ ,c("num.valid.stations")] = as.data.frame(gapply(summary.by.pol.and.type.index, which="pollutant",FUN=length.unique, groups=summary.by.pol.and.type.index[ ,filterby.col.label]))
      summary.by.filter.index[ ,c("num.stations")] = as.data.frame(rep(length(unique(pollutants$pollutant)), filterby.length))
   }
} else if (bGrpModel) {
   if (bAssessConcs) {
      #need to remove NA lines before finding length of unique 'model'
      summary.by.pol.and.type.conc = summary.by.pol.and.type.conc[ !is.na(summary.by.pol.and.type.conc$num.valid.stations), ]
      summary.by.filter.conc[ ,c("num.valid.stations")] = as.data.frame(gapply(summary.by.pol.and.type.conc, which="model",FUN=length.unique, groups=summary.by.pol.and.type.conc[ ,filterby.col.label]))
      summary.by.filter.conc[ ,c("num.stations")] = as.data.frame(rep(length(names(all.data)[names(all.data)!="date" & names(all.data)!="obs" & names(all.data)!="station" & names(all.data)!="type" & names(all.data)!="pollutant"]), filterby.length))
   }
   if (bAssessForecasts) {
      #need to remove NA lines before finding length of unique 'model'
      summary.by.pol.and.type.index = summary.by.pol.and.type.index[ !is.na(summary.by.pol.and.type.index$num.valid.stations), ]
      summary.by.filter.index[ ,c("num.valid.stations")] = as.data.frame(gapply(summary.by.pol.and.type.index, which="model",FUN=length.unique, groups=summary.by.pol.and.type.index[ ,filterby.col.label]))
      summary.by.filter.index[ ,c("num.stations")] = as.data.frame(rep(length(names(all.data)[names(all.data)!="date" & names(all.data)!="obs" & names(all.data)!="station" & names(all.data)!="type" & names(all.data)!="pollutant"]), filterby.length))
   }
}
if (bAssessConcs) {rm(temp.a.conc, temp.b.conc)}
if (bAssessForecasts) {rm(temp.a.index, temp.b.index)}
if (filterby.col.label == "station") {
   filterby.col.label = groupby.col.label
}

#Set up graphics page size
if (bPDF & (bPaperMedium | bPaperSmall)) {
   bPaperLarge = TRUE
   bPaperMedium = FALSE
   bPaperSmall = FALSE
   cat("   Warning: PDF output will be to A4\n")   
}
# note "height" is longest length, for square output use "width" twice
if (bPaperLarge) {  
   paper.height = 297
   paper.width = 210
} else if (bPaperMedium) {
   paper.height = round(170*sqrt(2))
   paper.width = 170   
} else if (bPaperSmall) { # set square plots to width*width, mimics A6 but a bit smaller
   paper.height = round(125*sqrt(2))
   paper.width = 125   
}

if (bAssessForecasts) {
   #Summarise the data over the "groupby" variable, with appropriate filtering
   # Index stats table
   temp.accuracy.allstns = fcst.accuracy.allstns
   # Use ave() function to allow more than one column in the grouping
   temp.accuracy.allstns$sum.num.valid.indices = ave(temp.accuracy.allstns$num.valid.indices, temp.accuracy.allstns[ ,c(groupby.col.label, filterby.col.label)], FUN = sum)
   #do percentage calculations based on number of indices / total number valid indices
   temp.accuracy.allstns$perc.indices.same = 100*ave(temp.accuracy.allstns$indices.same, temp.accuracy.allstns[ ,c(groupby.col.label, filterby.col.label)], FUN=sum)/temp.accuracy.allstns$sum.num.valid.indices
   temp.accuracy.allstns$perc.indices.oneapart = 100*ave(temp.accuracy.allstns$indices.oneapart, temp.accuracy.allstns[ ,c(groupby.col.label, filterby.col.label)], FUN=sum)/temp.accuracy.allstns$sum.num.valid.indices
   # create data frame
   fcst.accuracy = temp.accuracy.allstns[!duplicated(temp.accuracy.allstns[ ,c(groupby.col.label, filterby.col.label)]), c(groupby.col.label, filterby.col.label,"sum.num.valid.indices","perc.indices.same","perc.indices.oneapart")] 
   #remove duplicate columns when no filtering occurs
   fcst.accuracy = fcst.accuracy[ ,!colnames(fcst.accuracy)==paste(groupby.col.label,".1",sep="")]
   rm(temp.accuracy) 
   #rename columns of data frame to match those used up to v2.2 of the toolkit
   names(fcst.accuracy)[names(fcst.accuracy)=="sum.num.valid.indices"] = "num.valid.indices"   
 
   #Alert stats table
   temp.fcst = fcst.allstns
   temp.fcst$sum.num.obs.alerts = ave(temp.fcst$num.obs.alerts, temp.fcst[ ,c(groupby.col.label, filterby.col.label,"alert.name")], FUN=sum) 
   temp.fcst$sum.a = ave(temp.fcst$a, temp.fcst[ ,c(groupby.col.label, filterby.col.label,"alert.name")], FUN=sum)
   temp.fcst$sum.b = ave(temp.fcst$b, temp.fcst[ ,c(groupby.col.label, filterby.col.label,"alert.name")], FUN=sum)
   temp.fcst$sum.c = ave(temp.fcst$c, temp.fcst[ ,c(groupby.col.label, filterby.col.label,"alert.name")], FUN=sum)
   temp.fcst$sum.d = ave(temp.fcst$d, temp.fcst[ ,c(groupby.col.label, filterby.col.label,"alert.name")], FUN=sum)
   fcst = temp.fcst[!duplicated(temp.fcst[ ,c(groupby.col.label, filterby.col.label,"alert.name")]), c(groupby.col.label, filterby.col.label,"alert.name","sum.num.obs.alerts","sum.a","sum.b","sum.c","sum.d")] 
   #remove duplicate columns when no filtering occurs
   fcst = fcst[ ,!colnames(fcst)==paste(groupby.col.label,".1",sep="")]
   rm(temp.fcst)
   #rename columns of data frame to match those used up to v2.2 of the toolkit
   names(fcst)[names(fcst)=="sum.num.obs.alerts"] = "num.obs.alerts"
   names(fcst)[names(fcst)=="sum.a"] = "a" 
   names(fcst)[names(fcst)=="sum.b"] = "b" 
   names(fcst)[names(fcst)=="sum.c"] = "c"
   names(fcst)[names(fcst)=="sum.d"] = "d"
    
   #Calculate forecast assessment metrics
   cat("Calculating forecast assessment metrics...\n")
   #Forecast performance metrics:
   #Frequency bias (FB) [NA implies no events observed]
   fcst$FB=NA
   fcst$FB = (fcst$a + fcst$b)/(fcst$a + fcst$c)
   fcst$FB[(fcst$a + fcst$c) == 0] = NA
   #Probability of event detection (POD) [NA implies no events observed]
   fcst$POD=NA
   fcst$POD = fcst$a /(fcst$a + fcst$c)
   fcst$POD[(fcst$a + fcst$c) == 0] = NA
   #Probability of a correct forecast (PCF) [NA implies no data - shouldn't happen]
   fcst$PCF=NA
   fcst$PCF = (fcst$a + fcst$d)/(fcst$a + fcst$b + fcst$c + fcst$d)
   fcst$PCF[(fcst$a + fcst$b + fcst$c + fcst$d) == 0] = NA
   #False alarm ratio (FAR) [NA implies no events predicted]
   fcst$FAR=NA
   fcst$FAR = fcst$b /(fcst$a + fcst$b)
   fcst$FAR[(fcst$a + fcst$b) == 0] = NA
   #Probability of false detection (POFD) [NA implies no non-events]
   fcst$POFD=NA
   fcst$POFD = fcst$b /(fcst$b + fcst$d)
   fcst$POFD[(fcst$b + fcst$d) == 0] = NA
   #Overall skill
   fcst$total.skill = 0.25*(fcst$PCF+fcst$POD+(1-fcst$FAR)+(1-fcst$POFD))
   #Odds ratio skill score (ORSS)
   fcst$OR = NA
   fcst$OR = (fcst$a*fcst$d)/(fcst$b*fcst$c)
   #Meaningless if a+c=0 (no obs alerts) or a+b=0 (no mod alerts)
   fcst$OR[fcst$a+fcst$c==0 | fcst$a+fcst$b==0] = NA
   fcst$ORSS = NA
   fcst$ORSS = (fcst$OR-1)/(fcst$OR+1)
   #ORSS score is perfect (1) if b*c=0
   fcst$ORSS[fcst$b*fcst$c==0]=1
   #Meaningless if a+c=0 (no obs alerts) or a+b=0 (no mod alerts)
   fcst$ORSS[fcst$a+fcst$c==0 | fcst$a+fcst$b==0] = NA

   #Output data to csv
   if (boutput.stats.as.csv) {
      cat("Exporting forecast index data, stats and alert stats to CSV...\n")
      write.csv(fcst,file=paste(eval.output.root,"_forecast_alert_stats.csv",sep=""),row.names = FALSE)
      #Rename 'site' column as 'station' to be consistent with all other output
      names(fcst.data.for.output)[names(fcst.data.for.output)=="site"]="station"
      write.csv(fcst.data.for.output,file=paste(eval.output.root,"_forecast_index_data.csv",sep=""),row.names = FALSE)
      write.csv(fcst.accuracy,file=paste(eval.output.root,"_forecast_index_stats.csv",sep=""),row.names = FALSE)
   }

   if (bPDF) {
      cat("Sending forecast assessment charts to pdf...\n")
   } else if (bJPG) {
      cat("Sending forecast assessment charts to jpg...\n")
   } else if (bPNG) {
      cat("Sending forecast assessment charts to png...\n")
   }
   max.x=max(fcst$num.obs.alerts)

   if (bGrpStnName) {
      if (bAllStations) {
          text.stntype = "ALL STATIONS"
      } else if (bSingleStationType) {
          text.stntype = toupper(single.station.type)
      }
   } else if (bGrpStnType) {
      text.stntype = "ALL STATION TYPES"
   } else if (bGrpPollutant) {
      if (bAllPols) {
         text.stntype = "ALL POLLUTANTS"
      } else if (bSinglePollutant) {
         text.stntype = toupper(single.pollutant)
      }
   }

   #Forecast index accuracy: stacked bar chart with number of valid forecasts marked on
   if (bForecastAccuracy) {      
      if (bPDF) {
         pdf(paste(eval.output.root,"_forecast_accuracy.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
      }
      for (j in 1:filterby.length) {
         # data
         if (!bFilterPollutant & !bFilterStnType & !bFilterModel) {
            fcst.onepol.allstations = fcst.accuracy
         } else {
            # multiple filters
            subset.string = ""
            for (s in 1:(bFilterPollutant + bFilterStnType + bFilterModel)) {
               subset.string = paste(subset.string,filterby.col.label[s],"=='",filterby.list[j,s],"' & ",sep="")
            }
            subset.string = substr(subset.string, 1, nchar(subset.string)-3)  #remove final " & " from subset string
            fcst.onepol.allstations = subset(fcst.accuracy,eval(parse(text=subset.string)))             
         }
         # strings for plot titles
         if (bFilterStnType) {
            fcst.string.stn.j = paste(toupper(filterby.list[j, 'type'])) 
         } else if (bSingleStationType) {
            fcst.string.stn.j = toupper(single.station.type)
         } else if (bGrpStnType) {
            fcst.string.stn.j = "ALL STATION TYPES"
         } else {
            fcst.string.stn.j = "ALL STATIONS"            
         }
         if (bFilterPollutant) {
            fcst.string.pol.j = paste(toupper(filterby.list[j, 'pollutant'])) 
         } else if (bSinglePollutant) {
            fcst.string.pol.j = toupper(single.pollutant)         
         } else {
            fcst.string.pol.j = "ALL POLLUTANTS"            
         }
         if (bFilterModel) {
            filterby.string.j = paste(filterby.list[j, 'model'],fcst.string.stn.j,fcst.string.pol.j)
         } else if (bSingleModel) {
            filterby.string.j = paste(single.model,fcst.string.stn.j,fcst.string.pol.j)
         } else {
            filterby.string.j = paste(fcst.string.stn.j,fcst.string.pol.j)         
         }
         if (bFilterPollutant | bSinglePollutant) {
            if (bFilterStnType | bSingleStationType) {
               filterby.filename.j = paste(tolower(fcst.string.stn.j), "_",tolower(fcst.string.pol.j), sep="")
            } else if (!bFilterStnType) {
               filterby.filename.j = tolower(fcst.string.pol.j)
            }
         } else if (!bFilterPollutant) {
            if (bFilterStnType | bSingleStationType) {
               filterby.filename.j = tolower(fcst.string.stn.j)
            } else if (!bFilterStnType) {
               filterby.filename.j = "all"            
            }
         }
         if (bFilterModel) {
            filterby.filename.j = paste(filterby.list[j, 'model'],"_",filterby.filename.j,sep="")
         } else if (bSingleModel) {
            filterby.filename.j = paste(single.model,"_",filterby.filename.j,sep="")
         }

         if (length(fcst.onepol.allstations[ ,1])==0) {
            cat(paste("   Warning: no data to graph for",tolower(filterby.string.j),"- skipping this selection\n"))
            next
         }
         max.num.indices = max(fcst.onepol.allstations$num.valid.indices)
         sorted.data = fcst.onepol.allstations[order(fcst.onepol.allstations$num.valid.indices),]
         sorted.data$sort.order = c(1:length(sorted.data[ ,1]))
         sorted.data.matrix = matrix(c(sorted.data$perc.indices.same,sorted.data$perc.indices.oneapart),
            nrow=length(sorted.data$sort.order),ncol=2)
         if (bJPG) {
            jpeg(paste(eval.output.root,"_",filterby.filename.j,"_forecast_accuracy.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
         } else if (bPNG) {
            png(paste(eval.output.root,"_",filterby.filename.j,"_forecast_accuracy.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
         }
         #Set up the three drawing panels
         layout(matrix(c(1,2,3),nrow=3,ncol=1,byrow=TRUE),heights=c(8,3,2))
         #Put the graph in the top panel
         par(mar=c(4,4,5,5),las=1,col.axis="violetred")
         barplot(t(sorted.data.matrix),border=NA,col=c("lightgreen","lightgrey"),
            xlab=paste(groupby.key.string,"number (ordered by the number of forecast indices valid for comparison)"),
            ylab="Percentage of forecast indices valid for comparison",
            ylim=c(0,120),
            names.arg=sorted.data$sort.order,
            main=list(paste("Forecast Index Accuracy:",filterby.string.j),cex=1.5),
            legend.text=c("modelled = observed",expression("modelled = observed" %+-% "1")),
            axes=FALSE,space=0.5,
            args.legend=c(x="top",bty="n",horiz=TRUE)
         )
         box(which="plot",col="grey")
         par(col.axis="black")
         axis(2, at=c(0,50,100))
         abline(a=100,b=0,lty=2,col="grey")
         if (bPaperMedium | bPaperSmall) { 
            mtext("Percentage of forecast indices valid for comparison where the modelled index is either equal to or 1 band apart from the observed index",
               3,line=0.5,font=3,cex=0.7)
         } else {
            mtext("Percentage of forecast indices valid for comparison where the modelled index is either equal to or 1 band apart from the observed index",
               3,line=0.5,font=3)
         }
         par(new=TRUE)
         #Add the number of valid indices for each station
         plot(sorted.data$num.valid.indices,axes=FALSE,
            xlab=NA,ylab=NA,
            ylim=c(0,1.2*max.num.indices),yaxs="i",
            xlim=c(0.75,length(sorted.data$num.valid.indices)+0.25),
            col="blue")
         par(col.axis="blue")
         axis(4,col="blue",col.text="blue",at=floor(seq(from=0,to=max.num.indices, by = max.num.indices/5)))
         mtext("Number of forecast indices valid for comparison",4,line=3,col="blue",las=0,cex=0.7)
         #Put the legend in the middle panel
         par(mar=c(1,0,1,1))
         nr=10
         nc=floor((length(sorted.data$sort.order)-1)/nr)+1
         make.table(nr,nc)
         draw.title("Key to numbers",nc)
         for (n in 1:length(sorted.data$sort.order)) {
            draw.sample.cell(sorted.data$sort.order[n], toupper(eval(parse(text=paste("sorted.data$",groupby.col.label,"[n]",sep="")))), n, nr)
         }
         #Put the summary table in the bottom panel
         par(mar=c(1,0,3,0))
         make.table(1,2)
         draw.title("Summary",2)
         draw.sample.cell("Date range",paste(data.start.date.for.graphs," to ",data.end.date.for.graphs),1,1)
         draw.sample.cell(paste("Valid",paste(tolower(groupby.key.string),"s",sep="")),paste(summary.by.filter.index$num.valid.stations[j],"out of",summary.by.filter.index$num.stations[j]),2,1)
         mtext(airsheds.text,1,line=0,col="black",las=0,cex=0.8)
         if (bJPG | bPNG) {
            graphics.off()
         }
      }
      if (bPDF) {
         graphics.off()
      }
   }      
   #Odds ratio skill score (ORSS)
   if (bOddsRatioSkillScore) {
      if (bPDF) {
         pdf(paste(eval.output.root,"_forecast_ORSS.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
      }
      for (j in 1:filterby.length) {
         # data
         if (!bFilterPollutant & !bFilterStnType & !bFilterModel) {
            fcst.onepol = fcst
         } else {
            # multiple filters
            subset.string = ""
            for (s in 1:(bFilterPollutant + bFilterStnType + bFilterModel)) {
               subset.string = paste(subset.string,filterby.col.label[s],"=='",filterby.list[j,s],"' & ",sep="")
            }
            subset.string = substr(subset.string, 1, nchar(subset.string)-3)  #remove final " & " from subset string
            fcst.onepol = subset(fcst,eval(parse(text=subset.string)))             
         }
         # strings for plot titles
         if (bFilterStnType) {
            fcst.string.stn.j = paste(toupper(filterby.list[j, 'type'])) 
         } else if (bSingleStationType) {
            fcst.string.stn.j = toupper(single.station.type)
         } else if (bGrpStnType) {
            fcst.string.stn.j = "ALL STATION TYPES"
         } else {
            fcst.string.stn.j = "ALL STATIONS"            
         }
         if (bFilterPollutant) {
            fcst.string.pol.j = paste(toupper(filterby.list[j, 'pollutant'])) 
         } else if (bSinglePollutant) {
            fcst.string.pol.j = toupper(single.pollutant)         
         } else {
            fcst.string.pol.j = "ALL POLLUTANTS"            
         }
         if (bFilterModel) {
            filterby.string.j = paste(filterby.list[j, 'model'],fcst.string.stn.j)
         } else if (bSingleModel) {
            filterby.string.j = paste(single.model,fcst.string.stn.j)
         } else {
            filterby.string.j = fcst.string.stn.j        
         }
         if (bFilterPollutant | bSinglePollutant) {
            if (bFilterStnType | bSingleStationType) {
               filterby.filename.j = paste(tolower(fcst.string.stn.j), "_",tolower(fcst.string.pol.j), sep="")
            } else if (!bFilterStnType) {
               filterby.filename.j = tolower(fcst.string.pol.j)
            }
         } else if (!bFilterPollutant) {
            if (bFilterStnType | bSingleStationType) {
               filterby.filename.j = tolower(fcst.string.stn.j)
            } else if (!bFilterStnType) {
               filterby.filename.j = "all"            
            }
         }
         if (bFilterModel) {
            filterby.filename.j = paste(filterby.list[j, 'model'],"_",filterby.filename.j,sep="")
         } else if (bSingleModel) {
            filterby.filename.j = paste(single.model,"_",filterby.filename.j,sep="")
         }

         if (length(fcst.onepol[ ,1])==0) {
            cat(paste("   Warning: no data to graph for",tolower(fcst.string.stn.j),tolower(fcst.string.pol.j),"- skipping this selection\n"))
            next
         }
         for (k in 1:length(thresholds$name)) {
            fcst.onepol.onealert = subset(fcst.onepol,alert.name==thresholds$name[k])
            sorted.data = fcst.onepol.onealert[order(fcst.onepol.onealert$num.obs.alerts),]
            sorted.data$sort.order = c(1:length(sorted.data[ ,1]))
            max.alerts = max(max(sorted.data$num.obs.alerts),max(sorted.data$a+sorted.data$b))
            if (bJPG) {
               jpeg(paste(eval.output.root,"_",filterby.filename.j,"_",thresholds$name[k],"_forecast_ORSS.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
            } else if (bPNG) {
               png(paste(eval.output.root,"_",filterby.filename.j,"_",thresholds$name[k],"_forecast_ORSS.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
            }
            #Set up the three drawing panels
            layout(matrix(c(1,2,3),nrow=3,ncol=1,byrow=TRUE),heights=c(8,3,2))
            #Put the graph in the top panel
            par(mar=c(4,4,5,5),las=1,col.axis="violetred",xpd=FALSE)
            barplot(sorted.data$ORSS,border=NA,
               xlab=paste(groupby.key.string,"number (ordered by the number of observed alert threshold exceedences)"),
               ylab="Odds Ratio Skill Score (ORSS)",
               ylim=c(-1,1.3),
               names.arg=sorted.data$sort.order,
               main=list(paste("Odds Ratio Skill Score (ORSS):",filterby.string.j,toupper(thresholds$name[k]),fcst.string.pol.j),cex=1.5),
               axes=FALSE,space=0.5,
            )
            par(col.axis="black")
            axis(2, at=c(-1,0,1), labels=c("-1 (Bad)","0 (Random)","1 (Good)"),las=0)
            abline(a=1,b=0,col="green")
            abline(a=0,b=0,col="grey")
            abline(a=-1,b=0,col="red")
            mtext("No bars are shown where ORSS is invalid (the number of observed or modelled alert threshold exceedences is zero)",
               3,line=0.5,font=3,cex=0.7)
            #Add the number of observed alert threshold exceedences for each station
            if (max.alerts==0) {
              axis.limits = c(0,1.15)
              axis.ticks = c(0,1)
            } else {
              axis.limits = c(0,1.15*max.alerts)
              axis.ticks=floor(seq(from=0,to=max.alerts, by = max.alerts/5))
            }  
            par(new=TRUE,xpd=TRUE)
            plot(sorted.data$num.obs.alerts,axes=FALSE,
               xlab=NA,ylab=NA,
               xlim=c(0.75,length(sorted.data$sort.order)+0.25),
               ylim=axis.limits,yaxs="i",
               col="blue")
            par(col.axis="blue")
            axis(4,col="blue",col.text="blue",at=axis.ticks)
            mtext("Number of alert threshold exceedences",4,line=3,col="blue",las=0, cex=0.7)
            #Add the number of modelled alert threshold exceedences for each station
            par(new=TRUE,xpd=TRUE)
            plot(sorted.data$a+sorted.data$b,axes=FALSE,
               xlab=NA,ylab=NA,
               xlim=c(0.75,length(sorted.data$sort.order)+0.25),
               ylim=axis.limits,yaxs="i",
               col="blue",pch=4)
            legend(x="topright",c("observed","modelled"),pch=c(1,4),bty="n",text.col="blue",col="blue",horiz=TRUE)
            #Put the legend in the middle panel
            par(mar=c(1,0,1,1))
            nr=10
            nc=floor((length(sorted.data$sort.order)-1)/nr)+1
            make.table(nr,nc)
            draw.title("Key to numbers",nc)
            for (n in 1:length(sorted.data$sort.order)) {
               draw.sample.cell(sorted.data$sort.order[n], toupper(eval(parse(text=paste("sorted.data$",groupby.col.label,"[n]",sep="")))), n, nr)
            }
            #Put the summary table in the bottom panel
            par(mar=c(1,0,3,0))
            make.table(1,2)
            draw.title("Summary",2)
            draw.sample.cell("Date range",paste(data.start.date.for.graphs," to ",data.end.date.for.graphs),1,1)
            draw.sample.cell(paste("Valid",paste(tolower(groupby.key.string),"s",sep="")),paste(summary.by.filter.index$num.valid.stations[j],"out of",summary.by.filter.index$num.stations[j]),2,1)
            mtext(airsheds.text,1,line=0,col="black",las=0,cex=0.8)
            if (bJPG | bPNG) {
               graphics.off()
            }            
         }
      }
      if (bPDF) {
         graphics.off()
      }
   }

   #Function to draw one of the performance metrics graphs
   # - plots scatter plot of y as a function of x (black circles)
   # - title.text is the title for the graph
   # - a green line will be plotted at y=good
   # - a red line will be plotted at y=bad
   performance.graph = function(x,y,title.text,good,bad) {
            par(mar=c(2,3,2,0),las=1,xpd=TRUE,col.axis="black")
            plot(x,y,
               xlab=NA,
               ylab=NA,
               ylim=c(0,1),yaxs="i",
               main=title.text,
               axes=FALSE,type="n"
            )
            box(which="plot",col="grey")
            axis(2, at=c(0,0.5,1))
            par(xpd=FALSE)
            abline(a=0.5,b=0,lty=3,col="grey")
            abline(v=c(1:length(x)),col="grey")
            par(col.axis="violetred")
            axis(1,at=floor(seq(from=1,to=length(sorted.data$sort.order),by=(length(sorted.data$sort.order)-1)/10)))
            abline(a=good,b=0,col="green")
            abline(a=bad,b=0,col="red")
            par(new=TRUE,xpd=TRUE)
            plot(x,y,
               xlab=NA,
               ylab=NA,
               ylim=c(0,1),yaxs="i",
               axes=FALSE
            )
   }

   #Forecast performance metrics
   if (bForecastPerformanceMetrics) {
      #This document has one page per pollutant per threshold
      if (bPDF) {
         pdf(paste(eval.output.root,"_forecast_performance.pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
      }
      for (j in 1:filterby.length) {
         # data
         if (!bFilterPollutant & !bFilterStnType & !bFilterModel) {
            fcst.onepol = fcst
         } else {
            # multiple filters
            subset.string = ""
            for (s in 1:(bFilterPollutant + bFilterStnType + bFilterModel)) {
               subset.string = paste(subset.string,filterby.col.label[s],"=='",filterby.list[j,s],"' & ",sep="")
            }
            subset.string = substr(subset.string, 1, nchar(subset.string)-3)  #remove final " & " from subset string
            fcst.onepol = subset(fcst,eval(parse(text=subset.string)))             
         }
         # strings for plot titles
         if (bFilterStnType) {
            fcst.string.stn.j = paste(toupper(filterby.list[j, 'type'])) 
         } else if (bSingleStationType) {
            fcst.string.stn.j = toupper(single.station.type)
         } else if (bGrpStnType) {
            fcst.string.stn.j = "ALL STATION TYPES"
         } else {
            fcst.string.stn.j = "ALL STATIONS"            
         }
         if (bFilterPollutant) {
            fcst.string.pol.j = paste(toupper(filterby.list[j, 'pollutant'])) 
         } else if (bSinglePollutant) {
            fcst.string.pol.j = toupper(single.pollutant)         
         } else {
            fcst.string.pol.j = "ALL POLLUTANTS"            
         }
         if (bFilterModel) {
            filterby.string.j = paste(filterby.list[j, 'model'],fcst.string.stn.j)
         } else if (bSingleModel) {
            filterby.string.j = paste(single.model,fcst.string.stn.j)
         } else {
            filterby.string.j = fcst.string.stn.j        
         }
         if (bFilterPollutant | bSinglePollutant) {
            if (bFilterStnType | bSingleStationType) {
               filterby.filename.j = paste(tolower(fcst.string.stn.j), "_",tolower(fcst.string.pol.j), sep="")
            } else if (!bFilterStnType) {
               filterby.filename.j = tolower(fcst.string.pol.j)
            }
         } else if (!bFilterPollutant) {
            if (bFilterStnType | bSingleStationType) {
               filterby.filename.j = tolower(fcst.string.stn.j)
            } else if (!bFilterStnType) {
               filterby.filename.j = "all"            
            }
         }
         if (bFilterModel) {
            filterby.filename.j = paste(filterby.list[j, 'model'],"_",filterby.filename.j,sep="")
         } else if (bSingleModel) {
            filterby.filename.j = paste(single.model,"_",filterby.filename.j,sep="")
         }
          
         if (length(fcst.onepol[ ,1])==0) {
            cat(paste("   Warning: no data to graph for",tolower(fcst.string.stn.j), tolower(fcst.string.pol.j),"- skipping this selection\n"))
            next
         }
         for (k in 1:length(thresholds$name)) {
            fcst.onepol.onealert = subset(fcst.onepol,alert.name==thresholds$name[k])
            sorted.data = fcst.onepol.onealert[order(fcst.onepol.onealert$num.obs.alerts),]
            sorted.data$sort.order = c(1:length(sorted.data[ ,1]))
            max.alerts = max(max(sorted.data$num.obs.alerts),max(sorted.data$a+sorted.data$b))
            if (bJPG) {
               jpeg(paste(eval.output.root,"_",filterby.filename.j,"_",thresholds$name[k],"_forecast_performance.jpg",sep=""),width=paper.width,height=paper.height,units="mm",res=300)
            } else if (bPNG) {
               png(paste(eval.output.root,"_",filterby.filename.j,"_",thresholds$name[k],"_forecast_performance.png",sep=""),width=paper.width,height=paper.height,units="mm",res=300)
            }
            #Set the outer margins
            par(oma=c(0,0,3,0))
            #Set up the 7 equally-sized drawing panels 
            #(5 metrics graphs, 1 graph of number of alert threshold exceedences, 1 table (key))
            layout(matrix(c(1,2,3,4,5,6,7),nrow=7,ncol=1,byrow=TRUE),heights=c(1,1,1,1,1,1.5,1))
            #Draw first graph
            performance.graph(sorted.data$sort.order,sorted.data$PCF,"Probability of a correct forecast (PCF)",1,0)
            mtext(airsheds.text,3,line=4,col="black",las=0,cex=0.5)
            #Add title for the whole page
            titletext = paste("Forecast Performance Metrics:",filterby.string.j,toupper(thresholds$name[k]),fcst.string.pol.j)
            if ((nchar(titletext)>90 & bPaperLarge) | (nchar(titletext)>75 & bPaperMedium) | (nchar(titletext)>65 & bPaperSmall)) {
               titletext = paste("Forecast Performance Metrics:",paste(fcst.string.stn.j,toupper(thresholds$name[k]),fcst.string.pol.j),sep="\n")
            }
            mtext(side=3,outer=TRUE,font=2,titletext)
            #Draw the other 3 metrics graphs
            performance.graph(sorted.data$sort.order,sorted.data$POD,"Probability of detection (POD)",1,0)
            performance.graph(sorted.data$sort.order,sorted.data$FAR,"False alarm ratio (FAR)",0,1)
            performance.graph(sorted.data$sort.order,sorted.data$POFD,"Probability of false detection (POFD)",0,1)
            #Now draw the graph of the number of observed and modelled alert threshold exceedences
            if (max.alerts==0) {
              axis.limits = c(0,1.3)
              axis.ticks = c(0,1)
            } else {
              axis.limits = c(0,1.3*max.alerts)
              axis.ticks=floor(seq(from=0,to=max.alerts, by = max.alerts/5))
            }  
            par(mar=c(1,3,2,0))
            plot(sorted.data$num.obs.alerts,axes=FALSE,
               xlab=NA,ylab=NA,
               xlim=c(0.75,length(sorted.data$sort.order)+0.25),
               ylim=axis.limits,yaxs="i",
               col="blue",type="n")
            box(which="plot",col="grey")
            par(xpd=FALSE)
            abline(v=c(1:length(sorted.data$sort.order)),col="grey")
            par(xpd=TRUE)
            par(new=TRUE)
            plot(sorted.data$num.obs.alerts,axes=FALSE,
               xlab=NA,ylab=NA,
               xlim=c(0.75,length(sorted.data$sort.order)+0.25),
               ylim=axis.limits,yaxs="i",
               col="blue")
            par(col.axis="blue")
            axis(2,col="blue",col.text="blue",at=axis.ticks)
            par(col.axis="violetred")
            axis(1,at=floor(seq(from=1,to=length(sorted.data$sort.order),by=(length(sorted.data$sort.order)-1)/10)))
            mtext(paste(groupby.key.string, "number (ordered by the number of observed alert threshold exceedences)"),1,line=2,cex=0.7)
            par(new=TRUE)
            plot(sorted.data$a+sorted.data$b,axes=FALSE,
               xlab=NA,ylab=NA,
               xlim=c(0.75,length(sorted.data$sort.order)+0.25),
               ylim=axis.limits,yaxs="i",
               main="Number of alert threshold exceedences",
               col="blue",pch=4)
            par(xpd=FALSE)
            abline(a=max.alerts,b=0,col="blue",lty=3)
            legend(x="topright",c("observed","modelled"),pch=c(1,4),bty="n",text.col="blue",col="blue",horiz=TRUE)
            #Draw the key
            par(mar=c(0,0,3,0))
            nr=10
            nc=floor((length(sorted.data$sort.order)-1)/nr)+1
            make.table(nr,nc)
            draw.title("Key to numbers",nc)
            for (n in 1:length(sorted.data$sort.order)) {
               draw.sample.cell(sorted.data$sort.order[n], toupper(eval(parse(text=paste("sorted.data$",groupby.col.label,"[n]",sep="")))), n, nr)
            }
            #Draw the summary table
            par(mar=c(0,0,3,0))
            make.table(1,2)
            draw.title("Summary",2)
            draw.sample.cell("Date range",paste(data.start.date.for.graphs," to ",data.end.date.for.graphs),1,1)    
            draw.sample.cell(paste("Valid",paste(tolower(groupby.key.string),"s",sep="")),paste(summary.by.filter.index$num.valid.stations[j],"out of",summary.by.filter.index$num.stations[j]),2,1)
            if (bJPG | bPNG) {
               graphics.off()
            }
        }
      }
      if (bPDF) {
         graphics.off()
      }
   }
} #End of forecast assessment graphics

if (bAssessConcs) {
  
   #Output data to csv
   if (boutput.stats.as.csv) {
     
     if (bOpenairModStats) {
       #Calculate COE and IOA using openair modStats function
       cat("Calculating additional statistics using openair modStats function...\n")
       if (bGrpStnName) {
         filter.text = c("station","pollutant")
       } else if (bGrpStnType) {
         filter.text = c("type","pollutant")
       } else if (bGrpPollutant) {
         if (bFilterStnType) {
           filter.text = c("type","pollutant")
         } else {
           filter.text = "pollutant"
         }
       } else if (bGrpModel) {
         if (bFilterStnType) {
           filter.text = c("type","pollutant")
         } else {
           filter.text = "pollutant"
         }
       }
       for (i in 1:length(modelled.data.names)) {
         extra.stats.onemodel = modStats(all.data.averaged,mod=modelled.data.names[i],obs="obs",type=filter.text)  
         if (bGrpStnName) {
           #Add extra "all" station
           extra.stats.all.stations = modStats(all.data.averaged,mod=modelled.data.names[i],obs="obs",type="pollutant")
           extra.stats.all.stations$station="all"
           extra.stats.onemodel=rbind(extra.stats.onemodel,extra.stats.all.stations)
         }
         #Rename columns to avoid conflicts with existing stats
         names(extra.stats.onemodel)[names(extra.stats.onemodel) != filter.text]=paste("modstats.",names(extra.stats.onemodel)[names(extra.stats.onemodel) != filter.text],sep="")
         #Add the model name
         extra.stats.onemodel$model = modelled.data.names[i]
         if (i==1) {
           extra.stats=extra.stats.onemodel
         } else {
           extra.stats=rbind(extra.stats,extra.stats.onemodel)
         }
       }
       merge.by.columns = c(filter.text,"model")
       #Merge the extra stats with the existing stats
       conc.grouped=merge(conc.grouped,extra.stats,by=merge.by.columns)
     }
       
     cat("Exporting concentration stats to CSV...\n")
     write.csv(conc.grouped,file=paste(eval.output.root,"_conc_stats.csv",sep=""),row.names = FALSE)
   }
  
   #Target plot
   if (bTargetPlot) {

      #Set the plot title text
      if (bAllPols) text.pol = "ALL POLLUTANTS" else text.pol = toupper(paste(unique(conc.grouped$output.avg.time.and.statistic), single.pollutant))
      if (bAllStations) {
         text.stntype = "all stations" 
      } else if (bSingleStationType) {
         text.stntype = single.station.type
      }
      if (bAllModels) {
         text.filters = paste(toupper(text.stntype),", ",text.pol,sep="")
      } else {
         text.filters = paste(single.model,", ",toupper(text.stntype),", ",text.pol,sep="")
      }
      title.text=paste("Target plot: ",text.filters,", from ", data.start.date.for.graphs," to ",data.end.date.for.graphs, sep="")
      if ((nchar(title.text)>90 & bPaperLarge) | (nchar(title.text)>75 & bPaperMedium) | (nchar(title.text)>65 & bPaperSmall)) {
         title.text = paste(paste("Target plot: ",text.filters,sep=""),paste("from ", data.start.date.for.graphs," to ",data.end.date.for.graphs),sep="\n")
      }  
      
      #setup data frame for target plot with a vector of appropriate strings for the titles  
      target.plot.data = conc.grouped[ ,c(groupby.col.label,filterby.col.label,"SDO","MB","CRMSE.sign.DELTA.1.2")]
      target.plot.data$polstring = paste(conc.grouped$output.avg.time.and.statistic, toupper(conc.grouped$pollutant))
      # remove 'all' stations row used for statistics
      target.plot.data = target.plot.data[!apply(target.plot.data[,c(filterby.col.label,groupby.col.label)], MARGIN=1, FUN=function(x) any(x=="all")), ]

      #Set the formula for plotting
      # filtering by station type acheived by colouring graph symbols
      # filtering by model and pollutant leads to panels or tiles 
      # If only one model and/or only one pollutant, title indicates this so remove filtering      
      if (bGrpPollutant) {
         if (bSingleModel & !bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO)")   
         if (bSingleModel & bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | type")             
         if (!bSingleModel & !bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | model") 
         if (!bSingleModel & bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | model +type")
      } else if (bGrpModel | !bFilterModel) {
         if (bSinglePollutant & !bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO)")
         if (bSinglePollutant & bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | type")          
         if (!bSinglePollutant & !bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | polstring")  
         if (!bSinglePollutant & bFilterStnType) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | polstring + type")   
      } else {
         if (bSingleModel & !bSinglePollutant) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | polstring")
         if (bSinglePollutant & !bSingleModel) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | model")
         if (!bSinglePollutant & !bSingleModel) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO) | model + polstring")
         if (bSinglePollutant & bSingleModel) xyformula = as.formula("(MB/SDO) ~ (CRMSE.sign.DELTA.1.2/SDO)")
      }
      
      if (bGrpPollutant) {
         groupby.for.target = "polstring" # use string including averaging time in the key
      } else {
         groupby.for.target = groupby.col.label
      }
      
      #plot key
      max.chars = max(nchar(as.character(target.plot.data[,groupby.for.target])))
      target.key.columns = floor((paper.width-30)/(max.chars+16))
      if (target.key.columns > 8) {target.key.columns = 8} # don't allow too many columns
      
      #set a rainbow palette except yellow (hard to see) and with grey and black added
      palette(c(rainbow(6)[c(1,3:6)],"black","gray"))
      # determine colour scheme, with and without colour by station type
      if (!bFilterStnType | bGrpPollutant | bGrpModel) {
         target.par.settings = list(superpose.symbol=list(pch=1:25, col = palette()),strip.background=list(col="white"))
      } else { #bFilterStnType, use colour scheme on station types
         target.plot.data$type.colours = rep(NA, length(target.plot.data[ ,1]))
         target.plot.data$type.pch = rep(NA, length(target.plot.data[ ,1]))
         for (k in 1:length(station.types.available)) {
            # loop over the 7 distinct colours in palette() and give 1:25 symbols in turn
            if (k%%length(palette())==0) {
               target.plot.data$type.colours[target.plot.data$type==station.types.available[k]] = palette()[length(palette())] 
            } else {
               target.plot.data$type.colours[target.plot.data$type==station.types.available[k]] = palette()[k%%length(palette())] 
            }
            target.plot.data$type.pch[target.plot.data$type==station.types.available[k]] = c(1:sum(target.plot.data$type==station.types.available[k]))%%25        
         }
         target.plot.data$type.pch[target.plot.data$type.pch==0] = 25
         target.par.settings = list(superpose.symbol=list(col=target.plot.data$type.colours, pch=target.plot.data$type.pch), strip.background=list(col="white"))
      }
      
      if (bPDF) {    
         cat("Sending concentration assessment target plots to pdf...\n")
         pdf(paste(eval.output.root,"_conc_targetplot.pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
      } else if (bJPG) {
         cat("Sending concentration assessment target plots to jpg...\n")
         jpeg(paste(eval.output.root,"_conc_targetplot.jpg",sep=""),width=paper.width,height=paper.height,units="mm",res=300)
      } else if (bPNG) {
         cat("Sending concentration assessment target plots to png...\n")
         png(paste(eval.output.root,"_conc_targetplot.png",sep=""),width=paper.width,height=paper.height,units="mm",res=300)
      }
     
      #Target plots
      plot(xyplot(
         xyformula,
         target.plot.data,
         groups=eval(parse(text=groupby.for.target)),auto.key=list(columns=target.key.columns,space="bottom",cex=0.5),
         main=list(title.text,cex=0.8),
         sub=list(airsheds.text,cex=0.6),
         xlim=c(-2.1,2.1),ylim=c(-2.1,2.1),
         par.settings=target.par.settings,
         xlab=expression(CRMSE/sigma[obs]),
         ylab=expression(BIAS/sigma[obs]),
         scales=list(at=c(-2,-1,0,1,2)),aspect="iso",
         panel=function(x,y,...){
            panel.grid(h=-1,v=-1)
            panel.abline(h=0)
            panel.abline(v=0)
            panel.curve(sqrt(1-x^2),from=-1,to=1,n=1001)
            panel.curve(-sqrt(1-x^2),from=-1,to=1,n=1001)
            panel.text(x=0.2,y=1.1,"T=1",cex=0.7)
            panel.curve(sqrt(0.8^2-x^2),from=-0.8,to=0.8,col="forestgreen",n=1001)
            panel.curve(-sqrt(0.8^2-x^2),from=-0.8,to=0.8,col="forestgreen",n=1001)
            panel.curve(sqrt(0.65^2-x^2),from=-0.65,to=0.65,col="green",n=1001)
            panel.curve(-sqrt(0.65^2-x^2),from=-0.65,to=0.65,col="green",n=1001)
            panel.curve(sqrt(0.3^2-x^2),from=-0.3,to=0.3,lty=2,n=1001)
            panel.curve(-sqrt(0.3^2-x^2),from=-0.3,to=0.3,lty=2,n=1001)
            panel.text(x=-1,y=-1.8,expression(sigma[mod]>sigma[obs]),cex=0.8)
            panel.text(x=1,y=-1.8,expression(sigma[obs]>sigma[mod]),cex=0.8)
            panel.xyplot(x,y,...)
         }))
      graphics.off() #relevant to all output types as target plots are always output onto one page
      # reset colour palette
      palette("default")
   }

   #Target plot DELTA 3.3
   if (bTargetPlot.DELTA.3.3) {

      #Set the plot title text
      if (bAllPols) text.pol = "ALL POLLUTANTS" else text.pol = toupper(paste(unique(conc.grouped$output.avg.time.and.statistic), single.pollutant))
      if (bAllStations) {
         text.stntype = "all stations" 
      } else if (bSingleStationType) {
         text.stntype = single.station.type
      }
      if (bAllModels) {
         text.filters = paste(toupper(text.stntype),", ",text.pol,sep="")
      } else {
         text.filters = paste(single.model,", ",toupper(text.stntype),", ",text.pol,sep="")
      }
      title.text=paste("Target plot (DELTA 3.3): ",text.filters,", from ", data.start.date.for.graphs," to ",data.end.date.for.graphs, sep="")
      if ((nchar(title.text) >90 & bPaperLarge) | (nchar(title.text) >75 & bPaperMedium) | (nchar(title.text) >65 & bPaperSmall)) {
         title.text = paste(paste("Target plot (DELTA 3.3): ",text.filters, sep=""),paste("from ", data.start.date.for.graphs," to ",data.end.date.for.graphs),sep="\n")
      }  
      
      #setup data frame for target plot with a vector of appropriate strings for the titles  
      target.plot.data = conc.grouped[ ,c(groupby.col.label,filterby.col.label,"RMSu","MB","CRMSE.sign.DELTA.3.3")]
      target.plot.data$polstring = paste(conc.grouped$output.avg.time.and.statistic, toupper(conc.grouped$pollutant))
      # remove 'all' stations row used for statistics
      target.plot.data = target.plot.data[!apply(target.plot.data[,c(filterby.col.label,groupby.col.label)], MARGIN=1, FUN=function(x) any(x=="all")), ]
      # remove any where RMSu hasn't been calculated
      target.plot.data = subset(target.plot.data,!is.na(RMSu))

      #Set the formula for plotting
      # filtering by station type acheived by colouring graph symbols
      # filtering by model and pollutant leads to panels or tiles 
      # If only one model and/or only one pollutant, title indicates this so remove filtering      
      if (bGrpPollutant) {
         if (bSingleModel & !bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu))")      
         if (bSingleModel & bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | type")  
         if (!bSingleModel & !bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | model")   
         if (!bSingleModel & bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | model + type") 
      } else if (bGrpModel | !bFilterModel) {
         if (bSinglePollutant & !bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu))")      
         if (bSinglePollutant & bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | type")      
         if (!bSinglePollutant & !bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | polstring")    
         if (!bSinglePollutant & bFilterStnType) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | polstring + type")   
      } else {
         if (bSingleModel & !bSinglePollutant) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | polstring")
         if (bSinglePollutant & !bSingleModel) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | model")
         if (!bSinglePollutant & !bSingleModel) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu)) | model + polstring")
         if (bSinglePollutant & bSingleModel) xyformula = as.formula("(MB/(2*RMSu)) ~ (CRMSE.sign.DELTA.3.3/(2*RMSu))")
      }

      if (bGrpPollutant) {
         groupby.for.target = "polstring" # use string including averaging time in the key
      } else {
         groupby.for.target = groupby.col.label
      }
      
      #plot key
      max.chars = max(nchar(as.character(target.plot.data[,groupby.for.target])))
      target.key.columns = floor((paper.width-30)/(max.chars+16))
      if (target.key.columns > 8) {target.key.columns = 8} # don't allow too many columns
      
      #set a rainbow palette except yellow (hard to see) and with grey and black added
      palette(c(rainbow(6)[c(1,3:6)],"black","gray"))
      # determine colour scheme, with and without colour by station type
      if (!bFilterStnType | bGrpPollutant | bGrpModel) {
         target.par.settings = list(superpose.symbol=list(pch=1:25, col = palette()),strip.background=list(col="white"))
      } else { #bFilterStnType, use colour scheme on station types
         target.plot.data$type.colours = rep(NA, length(target.plot.data[ ,1]))
         target.plot.data$type.pch = rep(NA, length(target.plot.data[ ,1]))
         for (k in 1:length(station.types.available)) {
            # loop over the 7 distinct colours in palette() and give 1:25 symbols in turn
            if (k%%length(palette())==0) {
               target.plot.data$type.colours[target.plot.data$type==station.types.available[k]] = palette()[length(palette())] 
            } else {
               target.plot.data$type.colours[target.plot.data$type==station.types.available[k]] = palette()[k%%length(palette())] 
            }
            target.plot.data$type.pch[target.plot.data$type==station.types.available[k]] = c(1:sum(target.plot.data$type==station.types.available[k]))%%25        
         }
         target.plot.data$type.pch[target.plot.data$type.pch==0] = 25
         target.par.settings = list(superpose.symbol=list(col=target.plot.data$type.colours, pch=target.plot.data$type.pch), strip.background=list(col="white"))
      }
      
      if (bPDF) {    
         cat("Sending concentration assessment target plots to pdf...\n")
         pdf(paste(eval.output.root,"_conc_targetplotD33.pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
      } else if (bJPG) {
         cat("Sending concentration assessment target plots to jpg...\n")
         jpeg(paste(eval.output.root,"_conc_targetplotD33.jpg",sep=""),width=paper.width,height=paper.height,units="mm",res=300)
      } else if (bPNG) {
         cat("Sending concentration assessment target plots to png...\n")
         png(paste(eval.output.root,"_conc_targetplotD33.png",sep=""),width=paper.width,height=paper.height,units="mm",res=300)
      }
     
      #Target plots
      plot(xyplot(
         xyformula,
         target.plot.data,
         groups=eval(parse(text=groupby.for.target)),auto.key=list(columns=target.key.columns,space="bottom",cex=0.5),
         main=list(title.text,cex=0.8),
         sub=list(airsheds.text,cex=0.6),
         xlim=c(-2.1,2.1),ylim=c(-2.1,2.1),
         par.settings=target.par.settings,
         xlab=expression(CRMSE/(2*RMSu)),
         ylab=expression(BIAS/(2*RMSu)),
         scales=list(at=c(-2,-1,0,1,2)),aspect="iso",
         panel=function(x,y,...){
            panel.grid(h=-1,v=-1)
            panel.abline(h=0)
            panel.abline(v=0)
            panel.text(x=0.2,y=1.1,"T=1",cex=0.7)
            panel.curve(x^1,from=-2,to=-sqrt(0.5),col="grey50",n=1001)
            panel.curve(x^1,from=sqrt(0.5),to=2,col="grey50",n=1001)
            panel.curve(-x,from=-2,to=-sqrt(0.5),col="grey50",n=1001)
            panel.curve(-x,from=sqrt(0.5),to=2,col="grey50",n=1001)
            panel.curve(sqrt(1^2-x^2),from=-1,to=1,col="black",n=1001)
            panel.curve(-sqrt(1^2-x^2),from=-1,to=1,col="black",n=1001)
            panel.curve(sqrt(0.5^2-x^2),from=-0.5,to=0.5,col="black",lty=2,n=1001)
            panel.curve(-sqrt(0.5^2-x^2),from=-0.5,to=0.5,col="black",lty=2,n=1001)
            panel.text(x=-0.05,y=1.8,expression(Bias > 0),cex=0.8,col="grey50")
            panel.text(x=-0.05,y=-1.8,expression(Bias < 0),cex=0.8,col="grey50")
            panel.text(x=-1.8,y=0.1,expression(R),cex=0.8,col="grey50")
            panel.text(x=1.8,y=0.1,expression(SD),cex=0.8,col="grey50")
            panel.xyplot(x,y,...)
         }))
      graphics.off() #relevant to all output types as target plots are always output onto one page
      # reset colour palette
      palette("default")
   }

   #Box and whisker plot
   if (bBoxWhisker) {
      if (bPDF) {
         cat("Sending concentration assessment box and whisker plots to pdf...\n")
         pdf(paste(eval.output.root,"_conc_boxplot.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
      } else if (bJPG) {
         cat("Sending concentration assessment box and whisker plots to jpg...\n")
      } else if (bPNG) {
         cat("Sending concentration assessment box and whisker plots to png...\n")
      }
      #Loop over pollutants
      for (j in 1:length(pollutants$pollutant)) {
         #to make conditioning work, need "mod" as column and rows corresponding to models
         bDoneFirstModel = FALSE
         for (m in (1:length(modelled.data.names))) {
            temp.data.onemodel = subset(all.data.averaged, pollutant==pollutants$pollutant[j], select=c("date","obs","station","type","pollutant","output.avg.time.and.statistic",modelled.data.names[m]))      
            names(temp.data.onemodel)[names(temp.data.onemodel)==modelled.data.names[m]] = "mod"
            if (length(temp.data.onemodel$mod)==0) {
              next
            }
            temp.data.onemodel$model = modelled.data.names[m]
            if (bDoneFirstModel) {
               all.data.by.pol = rbind(all.data.by.pol,temp.data.onemodel)
            } else {
               all.data.by.pol = temp.data.onemodel
               bDoneFirstModel = TRUE
            }         
            rm(temp.data.onemodel)
         }
         rm(temp)
         #Check there is some valid data and all.data.by.pol exists
         if (!bDoneFirstModel) {
            cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
            next
         }

         all.data.by.pol = subset(all.data.by.pol, !is.na(obs) & !is.na(mod))
         
         #Check that there is data to plot and return an error if not
         if (length(all.data.by.pol$pollutant)==0) {
            cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
            next
         }
         
         #can't plot zeros on log scale
         if (bBoxLog & any(all.data.by.pol[,c("obs","mod")] == 0)) {
            all.data.by.pol = subset(all.data.by.pol, (obs!=0 & mod!=0))
         }   
      
         #Determine if outputs are to be filtered by station type
         if (!bFilterStnType) {
            if (bSingleStationType) {
              station.type = single.station.type
            } else {
              station.type = "all stations"
            }
         } else if (bFilterStnType) {
            station.type = sort(unique(all.data.by.pol$type))
         }
         #Find maximum y values for axis based on outliers
         if (bBoxOutliers) {
            box.outline = TRUE
            pollutant.max = max(all.data.by.pol$obs,all.data.by.pol$mod,na.rm=TRUE)
            pollutant.min = min(all.data.by.pol$obs,all.data.by.pol$mod,na.rm=TRUE)
         } else {
            box.outline = FALSE
            #Find the maximum concentration within range of data to be shown on graph
            conc.pol = subset(conc.grouped, pollutant==pollutants$pollutant[j])
            #remove 'all' stations statistics from box plot data
            conc.pol = conc.pol[eval(parse(text=paste("!conc.pol$",groupby.col.label," == 'all'",sep=""))), ] 
            pollutant.max = max(conc.pol$obs.upper.whisker, conc.pol$mod.upper.whisker)
            pollutant.min = min(conc.pol$obs.lower.whisker, conc.pol$mod.lower.whisker)
         }
         #Set up graph labels
         pollutant.units = (subset(pollutants, pollutant==pollutants$pollutant[j]))$output.units
         if (bGrpModel) {
            x.axis.label = paste("Stations",sep="")         
         } else {
            x.axis.label = paste(groupby.key.string,"s",sep="")
         }
         y.axis.label = paste(toupper(pollutants$pollutant[j])," (",pollutant.units,") ")

         #Loop over station types
         for (i in 1:length(station.type)) {
         
            #Set the plot title text and subset of data to use
            if (bSingleStationType) {
               all.data.by.pol.and.type=all.data.by.pol
               text.stntype.title = single.station.type
               summary.stats=subset(summary.by.filter.conc,pollutant==pollutants$pollutant[j])
            } else if (!bFilterStnType) {
               all.data.by.pol.and.type=all.data.by.pol
               text.stntype.title = station.type
               summary.stats=subset(summary.by.filter.conc,pollutant==pollutants$pollutant[j])
            } else if (bFilterStnType) {
               all.data.by.pol.and.type=subset(all.data.by.pol, type==station.type[i])
               summary.stats=subset(summary.by.filter.conc,pollutant==pollutants$pollutant[j] & type==station.type[i])
               text.stntype.title = station.type[i]
            }
            text.stntype.filename = gsub(" ", "", text.stntype.title)  # remove any white space for the filename
            models.available.by.pol.and.type = unique(all.data.by.pol.and.type$model)
            group.var.available.by.pol.and.type = unique(all.data.by.pol.and.type[,groupby.col.label])
            
            number.stations = length(group.var.available.by.pol.and.type)
            #Set the size of outliers
            if (number.stations>=100) {
               outlier.size = 0.25
            } else {
               outlier.size = 1.25 - (number.stations/100)
            }
            
            #set log/linear options
            if (bBoxLog) {
               #y-axis range and labels
               log.y.max = 10^ceiling(log10(1.05*pollutant.max))
               log.y.min = 10^floor(log10(0.95*pollutant.min))
               log.y.labels = as.vector(log.y.max)
               n=2
               while (log.y.labels[n-1] > min(1,log.y.min)) {
                  log.y.labels[n] = log.y.labels[n-1]/10
                  n = n+1
               }  
            }
            
            #If JPG or PNG is selected, set up the page
            if (bJPG) {
               jpeg(paste(eval.output.root,"_",pollutants$pollutant[j],"_",text.stntype.filename,"_conc_boxplot.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
            } else if (bPNG) {
               png(paste(eval.output.root,"_",pollutants$pollutant[j],"_",text.stntype.filename,"_conc_boxplot.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
            }   
            
            if (bGrpModel) {
               boxplotobsformula = as.formula(paste("obs~pollutant"))
               boxplotmodformula = as.formula(paste("mod~pollutant"))
               # predefine formatting for placement of bars forced to one central point
               at.sequence.centered = seq(from=1,to=1,length=1)         
            } else {
               # predefine formatting for placement of bars because of number of models
               at.sequence.centered = seq(from=1,to=(length(group.var.available.by.pol.and.type)),length=(length(group.var.available.by.pol.and.type)))
               boxplotobsformula = as.formula(paste("obs~",groupby.col.label))
               boxplotmodformula = as.formula(paste("mod~",groupby.col.label))
            }
            number.models = length(modelled.data.names)
            if (number.models > 2) {  # place bars at +- 0.2 from each groupby.col.label
               delta.at.models = 0.4/number.models
            } else {  # is small number ,place bars at +- 0.1
               delta.at.models = 0.2/number.models
            }
            #predefine at-sequence depending on if a model-station combination has valid data
            at.entry = at.sequence.centered-(delta.at.models*number.models/2)
            #Aesthetics
            box.model.fill.colours = rep(c("coral","coral2","coral3","coral4","indianred1","indianred2","indianred3","indianred4"),ceiling(number.models/8))
            box.model.outline.colours = rep(c("darkorange3","indianred3","brown3","brown4","pink","coral3","red","red4"),ceiling(number.models/8))
            
            title.text=paste("Box and Whisker Plot: ",toupper(project.name),", ",toupper(text.stntype.title),", ",toupper(pollutants$output.avg.time.and.statistic[j])," ",toupper(pollutants$pollutant[j]),sep="")  
            if ((nchar(title.text)>90 & bPaperLarge) | (nchar(title.text)>75 & bPaperMedium) | (nchar(title.text)>65 & bPaperSmall)) {
               title.text = paste(paste("Box and Whisker plot: ",toupper(project.name),sep=""),paste(toupper(text.stntype.title),", ",toupper(pollutants$output.avg.time.and.statistic[j])," ",toupper(pollutants$pollutant[j]),sep=""),sep="\n")
            } 
            
            if (bPaperLarge) {
               box.cex.axis = 0.7
            } else {
               box.cex.axis = 0.9
            }
            
            #If large number of stations, resort to numbering
            if ((number.stations >= 100 & bPaperLarge) | (number.stations >= 65 & bPaperMedium) | (number.stations >= 45 & bPaperSmall)) {
               #Set up 3 drawing panels
               layout(matrix(c(1,2,3),nrow=3,ncol=1,byrow=TRUE),heights=c(8,3,2))
               # boxplot
               par(mar=c(4,4,5,4),xpd=FALSE)
               #Draw first boxplot to set plot area
               if (bBoxLog) {
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,log="y",col="white",border="white",xlab=x.axis.label,ylab=y.axis.label,main=title.text,boxwex=0.5,outline=box.outline,
                     outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis, cex.lab=1.0,axes=FALSE,ylim=c(log.y.min,log.y.max),las=2)
               } else {
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,col="white",border="white",xlab=x.axis.label,ylab=y.axis.label,main=title.text,boxwex=0.5,outline=box.outline,
                     outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis, cex.lab=1.0,axes=FALSE,ylim=c(0, 1.2*pollutant.max),las=2)               
               }
               axis(side=1,at=unique(5*floor(c(1:number.stations)/5)+1),labels=as.character(unique(5*floor(c(1:number.stations)/5)+1)),las=0,tick=TRUE) # make numered x-axis rather than station names
               #Draw horizontal grid lines
               grid(nx=NA,ny=NULL,lty=3,col="grey")
               #Draw vertical grid lines
               abline(v=c(1:length(group.var.available.by.pol.and.type)),lty=3,col="grey")
               par(new=TRUE)
               if (bBoxLog) {
                  #Plot observed data boxplot
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,log="y",col="lightblue",border="blue",outline=box.outline,boxwex=0.5,ylim=c(log.y.min,log.y.max),outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                     at=at.entry,las=2,outwex=0.5,xaxt="n",yaxt="n")
                  #Plot modelled data boxplot
                  for (m in 1:number.models) {
                     #Reduce size of 'at' if some stations/entries are invalid for this model
                     at.entry.mod = at.entry[match(unique(all.data.by.pol.and.type[all.data.by.pol.and.type[,"model"]==models.available.by.pol.and.type[m], groupby.col.label]), unique(all.data.by.pol.and.type[ , groupby.col.label]))]
                     #plot
                     boxplot(boxplotmodformula,data=all.data.by.pol.and.type,subset=all.data.by.pol.and.type$model==models.available.by.pol.and.type[m],log="y",col=box.model.fill.colours[m],boxwex=0.5,add=TRUE,border=box.model.outline.colours[m],outline=box.outline,ylim=c(log.y.min,log.y.max),outpch=4,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                        at=at.entry.mod + (m*delta.at.models),outwex=0.5,las=2,xaxt="n",yaxt="n",width=rep(1,length(at.entry.mod)))
                  }
                  #add y-axis and labels
                  axis(side=2, at=log.y.labels, lab=as.character(log.y.labels), las=2, cex.axis=box.cex.axis) 
               } else {
                  #Plot observed data boxplot
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,col="lightblue",border="blue",outline=box.outline,boxwex=0.5,ylim=c(0,1.2*pollutant.max),outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                     at=at.entry,las=2,outwex=0.5,xaxt="n")
                  #Plot modelled data boxplot
                  for (m in 1:number.models) {
                     #Reduce size of 'at' if some stations/entries are invalid for this model
                     at.entry.mod = at.entry[match(unique(all.data.by.pol.and.type[all.data.by.pol.and.type[,"model"]==models.available.by.pol.and.type[m], groupby.col.label]), unique(all.data.by.pol.and.type[ , groupby.col.label]))]
                     #plot
                     boxplot(boxplotmodformula,data=all.data.by.pol.and.type,subset=all.data.by.pol.and.type$model==models.available.by.pol.and.type[m],col=box.model.fill.colours[m],boxwex=0.5,add=TRUE,border=box.model.outline.colours[m],outline=box.outline,ylim=c(0,1.2*pollutant.max),outpch=4,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                     at=at.entry.mod + (m*delta.at.models) ,outwex=0.5,las=2,xaxt="n",width=rep(1,length(at.entry.mod)))
                  }
               }
               legend(x="topleft",c("observed",models.available.by.pol.and.type),pch=22,col=c("blue",box.model.outline.colours[1:number.models]),pt.bg=c("lightblue",box.model.fill.colours[1:number.models]),bg="white")       
   
               #Key to stations   
               #Draw the key
               par(mar=c(0,0,3,0),xpd=FALSE)
               nr=10
               nc=floor((number.stations-1)/nr)+1
               make.table(nr,nc)
               draw.title("Key to numbers",nc)
               for (n in 1:number.stations) {
                  draw.sample.cell(n, toupper(group.var.available.by.pol.and.type[n]), n, nr)
               }
               
               #Draw the summary table
               par(mar=c(0,0,3,0))
               make.table(1,2)
               draw.title("Summary",2)
               draw.sample.cell("Date range",paste(data.start.date.for.graphs," to ",data.end.date.for.graphs),1,1)        
               draw.sample.cell(paste("Valid",paste(tolower(groupby.key.string),"s",sep="")),paste(max(summary.stats$num.valid.stations),"out of",max(summary.stats$num.stations)),2,1)
               mtext(airsheds.text,1,line=0,col="black",las=0,cex=0.8)

            } else {             
               #Set up 2 drawing panels
               if (bPaperLarge) {
                  layout(matrix(c(1,2),nrow=2,ncol=1,byrow=TRUE),heights=c(8,2))
               } else {  # medium or small, increase Summary proportion of page for aesthetics
                  layout(matrix(c(1,2),nrow=2,ncol=1,byrow=TRUE),heights=c(3,1))
               }
               par(mar=c(4,4,5,4),xpd=FALSE)
               #Draw first boxplot to set plot area
               if (bBoxLog) {
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,log="y",yaxt="n",col="white",border="white",xlab=x.axis.label,ylab=y.axis.label,main=title.text,boxwex=0.5,outline=box.outline,
                     outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis, cex.lab=1.0,ylim=c(log.y.min,log.y.max),las=2)
               } else {
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,col="white",border="white",xlab=x.axis.label,ylab=y.axis.label,main=title.text,boxwex=0.5,outline=box.outline,
                     outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis, cex.lab=1.0,ylim=c(0, 1.2*pollutant.max),las=2)               
               }
               #Draw horizontal grid lines
               grid(nx=NA,ny=NULL,lty=3,col="grey")
               #Draw vertical grid lines
               abline(v=c(1:length(group.var.available.by.pol.and.type)),lty=3,col="grey")
               par(new=TRUE)
               if (bBoxLog) {              
                  #Plot observed data boxplot
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,log="y",yaxt="n",col="lightblue",border="blue",outline=box.outline,boxwex=0.5,ylim=c(log.y.min,log.y.max),outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                     at=at.entry,las=2,outwex=0.5,xaxt="n")
                  #Plot modelled data boxplot
                  for (m in 1:number.models) {
                     #Reduce size of 'at' if some stations/entries are invalid for this model
                     at.entry.mod = at.entry[match(unique(all.data.by.pol.and.type[all.data.by.pol.and.type[,"model"]==models.available.by.pol.and.type[m], groupby.col.label]), unique(all.data.by.pol.and.type[ , groupby.col.label]))]
                     #plot
                     boxplot(boxplotmodformula,data=all.data.by.pol.and.type,subset=all.data.by.pol.and.type$model==models.available.by.pol.and.type[m],log="y",yaxt="n",col=box.model.fill.colours[m],boxwex=0.5,add=TRUE,border=box.model.outline.colours[m],outline=box.outline,ylim=c(log.y.min,log.y.max),outpch=4,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                        at=at.entry.mod + (m*delta.at.models),outwex=0.5,las=2,xaxt="n",width=rep(1,length(at.entry.mod)))
                  }
                  #add y-axis and labels
                  axis(side=2, at=log.y.labels, lab=as.character(log.y.labels), las=2, cex.axis=box.cex.axis) 
               } else {
                  #Plot observed data boxplot
                  boxplot(boxplotobsformula,data=all.data.by.pol.and.type,col="lightblue",border="blue",outline=box.outline,boxwex=0.5,ylim=c(0,1.2*pollutant.max),outpch=21,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                     at=at.entry,las=2,outwex=0.5,xaxt="n")
                  #Plot modelled data boxplot
                  for (m in 1:number.models) {
                     #Reduce size of 'at' if some stations/entries are invalid for this model
                     at.entry.mod = at.entry[match(unique(all.data.by.pol.and.type[all.data.by.pol.and.type[,"model"]==models.available.by.pol.and.type[m], groupby.col.label]), unique(all.data.by.pol.and.type[ , groupby.col.label]))]
                     #plot
                     boxplot(boxplotmodformula,data=all.data.by.pol.and.type,subset=all.data.by.pol.and.type$model==models.available.by.pol.and.type[m],col=box.model.fill.colours[m],boxwex=0.5,add=TRUE,border=box.model.outline.colours[m],outline=box.outline,ylim=c(0,1.2*pollutant.max),outpch=4,outcex=outlier.size,lty=1,cex.axis=box.cex.axis,
                        at=at.entry.mod + (m*delta.at.models),outwex=0.5,las=2,xaxt="n",width=rep(1,length(at.entry.mod)))
                  }
               }
               legend(x="topleft",c("observed",models.available.by.pol.and.type),pch=22,col=c("blue",box.model.outline.colours[1:number.models]),pt.bg=c("lightblue",box.model.fill.colours[1:number.models]),bg="white")       
   
               #If there is only one station on the graph, make this the x-axis label
               if (number.stations==1) {
                  if (bGrpModel) {
                     axis(side=1,at=at.sequence.centered - (delta.at.models*number.models/2),labels=c("all stations"),cex.axis=box.cex.axis,las=0)                  
                  } else {
                     axis(side=1,at=at.sequence.centered - (delta.at.models*number.models/2),labels=tolower(sort(group.var.available.by.pol.and.type)),cex.axis=box.cex.axis,las=2)
                  }
               }
               
               #Draw the summary table
               par(mar=c(1,0,3,0),xpd=TRUE)
               make.table(1,2)
               draw.title("Summary",2)
               draw.sample.cell("Date range",paste(data.start.date.for.graphs," to ",data.end.date.for.graphs),1,1)        
               draw.sample.cell(paste("Valid",paste(tolower(groupby.key.string),"s",sep="")),paste(max(summary.stats$num.valid.stations),"out of",max(summary.stats$num.stations)),2,1)
               mtext(airsheds.text,1,line=0,col="black",las=0,cex=0.8)

            }
            if (bJPG | bPNG) {
               graphics.off()
            }
            rm(all.data.by.pol.and.type)
         } #end loop over station types
         rm(all.data.by.pol)                        
      } #end loop over pollutants
      if (bPDF) {
         graphics.off()
      }
   }   

   #Scatter Plot, noting only one of FreqScatterPlot or bScatterPlot may be true 
   if (bFreqScatterPlot | bScatterPlot) {
      # predefine strings associated with conventional or frequency scatter plot
      if (bFreqScatterPlot) {
         scatter.string = "frequency"
         scatter.title.string = "Frequency Scatter Plot: "
      } else if (bScatterPlot) {
         scatter.string = "conventional"
         scatter.title.string = "Scatter Plot: "
      }
     
      if (bPDF) {
         cat(paste("Sending concentration assessment",scatter.string,"scatter plots to pdf...\n"))
         pdf(paste(eval.output.root,"_conc_",scatter.string,"_scatter_plot.pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
      } else if (bJPG) {
         cat(paste("Sending concentration assessment",scatter.string,"scatter plots to jpg...\n"))
      } else if (bPNG) {
         cat(paste("Sending concentration assessment",scatter.string,"scatter plots to png...\n"))
      }
      
      #Loop over pollutants
      for (j in 1:length(pollutants$pollutant)) {
         #to make conditioning work, need "mod" as column and rows corresponding to models
         bDoneFirstModel = FALSE
         for (m in (1:length(modelled.data.names))) {
            temp.data.onemodel = subset(all.data.averaged, pollutant==pollutants$pollutant[j], select=c("date","obs","station","type","pollutant","output.avg.time.and.statistic",modelled.data.names[m]))      
            names(temp.data.onemodel)[names(temp.data.onemodel)==modelled.data.names[m]] = "mod"
            if (length(temp.data.onemodel$mod)==0) {
              next
            }
            temp.data.onemodel$model = modelled.data.names[m]
            if (bDoneFirstModel) {
               all.data.by.pol = rbind(all.data.by.pol,temp.data.onemodel)
            } else {
               all.data.by.pol = temp.data.onemodel
               bDoneFirstModel = TRUE
            }         
            rm(temp.data.onemodel)
         }
         rm(temp)  
         #Check that there is valid data to plot and all.data.by.pol exists
         if (!bDoneFirstModel) {
            cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
            next
         }
          
         all.data.by.pol = subset(all.data.by.pol, !is.na(obs) & !is.na(mod))

         #Check that there is data to plot and return an error if not
         if (length(all.data.by.pol$pollutant)==0) {
            cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
            next
         }   
         
         #can't plot zeros on log scale
         if (bScatterLog & any(all.data.by.pol[,c("obs","mod")] == 0)) {
            all.data.by.pol = subset(all.data.by.pol, (obs!=0 & mod!=0))
         } 
        
         #setup pollutants
         pollutant.units = (subset(pollutants, pollutant==pollutants$pollutant[j]))$output.units         
         
         #Set graphical parameters
         graph.obs = "obs"
         graph.mod = "mod"
         pollutant.max = max(all.data.by.pol$obs,all.data.by.pol$mod,na.rm=TRUE)
         pollutant.min = min(all.data.by.pol$obs,all.data.by.pol$mod,na.rm=TRUE)
         if (bScatterLog) {
            max.val = 10^ceiling(log10(1.05*pollutant.max))
            min.val = 10^floor(log10(0.95*pollutant.min))
         } else {   
            # define linear axis extent and labels
            max.val = 1.05*pollutant.max
            min.val = 0
         }
         pol.display = paste(" (",pollutant.units,")",sep="")
         graph.z=NA
         if (bFreqScatterPlot) {
            # determine range of data and no. hexbins to modify colour key if necessary
            if (bScatterLog) {
               # if log-log scale, need to determine hexbin extents in log-log transford
               scatter.dummy.data = all.data.by.pol
               scatter.dummy.data$obs = log10(scatter.dummy.data$obs)
               scatter.dummy.data$mod = log10(scatter.dummy.data$mod)
            } else {
               scatter.dummy.data = all.data.by.pol
            }
            # determine the maximum count (in bins) i.e. max label of colour scale
            if (bFilterStnType) {
               station.types.available.by.pol = unique(scatter.dummy.data$type)
               max.count = rep(NA, length(station.types.available.by.pol))
               for (i in 1:length(station.types.available.by.pol)) {
                  type.subset = subset(scatter.dummy.data,type==station.types.available.by.pol[i])[ ,c("obs","mod")]
                  max.count[i] = max(hexbin(type.subset$obs, type.subset$mod, xbnds=range(type.subset[,c("obs","mod")]), ybnds=range(type.subset[,c("obs","mod")]))@count)
                  rm(type.subset)
               }
               hexbin.count = max(max.count, na.rm=TRUE)
            } else {
               hexbin.count = max(hexbin(scatter.dummy.data$obs, scatter.dummy.data$mod, xbnds=range(scatter.dummy.data[,c("obs","mod")]), ybnds=range(scatter.dummy.data[,c("obs","mod")]))@count)
            }
            rm(scatter.dummy.data)  
            # check frequency scatter plot has .gt. 3 count in each bin to avoid plotting errors
            if (hexbin.count <= 3) {
               cat(paste("   Warning: Too few data for frequency scatter plot output for pollutant", pollutants$pollutant[j], "- choose conventional scatter plot instead\n"))
               next
            } 
            # .gt. 32 count in each bin for non-linear key of 16 colours
            # .le. 32 counts, use non-linear key and reduce number of colours
            # .le. 16 counts, use linear key and reduce number of colours 
            if (hexbin.count <= 17) {
               scatter.trans = function(x) {x}
               scatter.inv = function(x) {x}
               scatter.colorcut = seq(0, 1, length = hexbin.count) # linear 1:1 scale
            } else if (hexbin.count <= 986) {  #floor(log(hexbin.count)/log(1.5)).le.17, first hexbin.count to guarantee 17 unique labels on 16 colours = 986
               scatter.trans = function(x) {log(x)}
               scatter.inv = function(x) {exp(x)} 
               scatter.colorcut = seq(0, 1, length = floor(log(hexbin.count)/log(1.5))) # log scale, ensure second level rounds up to integer 2       
            } else {
               scatter.trans = function(x) {log(x)}
               scatter.inv = function(x) {exp(x)}        
               scatter.colorcut = seq(0, 1, length = 17) # default  
            }
            # other graphical settings
            if (bScatterLog) {
               scatter.xbnds=range(log10(all.data.by.pol[,c("obs","mod")]))
            } else {
               scatter.xbnds=range(all.data.by.pol[,c("obs","mod")])            
            }
            scatter.cols=c(hsv(seq(0.675,0.400,by=-0.275/4), 0.6, 0.75), hsv(seq(0.350,0.2,by=-0.2/3),0.35,0.9), hsv(seq(0.2,0.1,by=-0.1/3),0.5,1.0), hsv(seq(0.1,0.0,by=-0.1/4), 0.75, 0.9), hsv(seq(1.0, 0.975,by=-0.025/1),0.75, 0.6)) # colour palette with reduced saturation, to mimic "hue", 17 colours to match default levels used by hexbin 
            scatter.pch=1
            scatter.cex=NULL
            scatter.mod.line=TRUE
            scatter.method="hexbin"     
         } else if (bScatterPlot) {
            scatter.trans = NULL
            scatter.inv = NULL      
            scatter.colorcut = NULL
            scatter.xbnds=NULL
            scatter.cols="black" 
            scatter.pch=16
            scatter.cex=NULL
            scatter.mod.line=TRUE
            scatter.method="scatter"          
         }
         
         # R doesn't like the column header "type" as it is also an intrinsic function
         scatter.plot.data = all.data.by.pol
         names(scatter.plot.data)[names(scatter.plot.data)=="type"] = "station.type"   
                 
         #Determine if outputs are to be filtered by station type
         # setup options related to filtering or not
         if (bAllStations & bFilterStnType) {
            # use scatterPlot's intrinsic filter function to display each station type on one page, one page per pollutant:
            filter.type=c("model","station.type")
            # titles
            text.stntype.title = "ALL STATIONS"
            text.stntype.filename = "allstationsbytype"
         } else {
            # scatterPlot(type = ...)          
            filter.type="model" 
            # titles
            if (bSingleStationType) {
              #Check that there is data to plot and return an error if not
              if (length(scatter.plot.data$pollutant)==0) {
                cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
                next
              }   
              text.stntype.title = single.station.type
              text.stntype.filename = single.station.type
            } else if (bAllStations & !bFilterStnType) {
              text.stntype.title = "ALL STATIONS"
              text.stntype.filename = "allstations"
            }        
         }
         if (bSingleModel) {
            text.filters = paste(single.model,", ",toupper(text.stntype.title),", ",toupper(pollutants$output.avg.time.and.statistic[j])," ",toupper(pollutants$pollutant[j]),pol.display,sep="")
         } else {
            text.filters = paste(toupper(text.stntype.title),", ",toupper(pollutants$output.avg.time.and.statistic[j])," ",toupper(pollutants$pollutant[j]),pol.display,sep="")
         }
         # setup output page
         if (bJPG) {
            jpeg(paste(eval.output.root,"_",pollutants$pollutant[j],"_",text.stntype.filename,"_conc_",scatter.string,"_scatter_plot.jpg",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
         } else if (bPNG) {
            png(paste(eval.output.root,"_",pollutants$pollutant[j],"_",text.stntype.filename,"_conc_",scatter.string,"_scatter_plot.png",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
         } 
         #title
         titletext=paste(scatter.title.string,toupper(project.name),", ",text.filters,sep="")  
         if (nchar(titletext)>65) {  # split onto two rows
            titletext=paste(paste(scatter.title.string,toupper(project.name), sep=""), text.filters, sep="\n")  
         }         
         # Draw frequency scatter plot
         scatterPlot(
            scatter.plot.data,
            x=graph.obs,y=graph.mod,z=graph.z,
            xlim=c(min.val,max.val),ylim=c(min.val,max.val),
            log.x=bScatterLog, log.y=bScatterLog,
            xlab="Observed",ylab="Modelled",
            cols=scatter.cols,pch=scatter.pch,cex=scatter.cex,
            type=filter.type, 
            method=scatter.method,
            mod.line=scatter.mod.line,
            main=titletext,
            sub=paste("Date range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep=""),
            aspect="iso",
            xbnds=scatter.xbnds, ybnds=scatter.xbnds,
            trans=scatter.trans, inv=scatter.inv,
            colorcut=scatter.colorcut
         )
                
         if (bJPG | bPNG) {
            graphics.off()
         }           
         rm(all.data.by.pol, scatter.plot.data)
	  } # next pollutant
	     
      if (bPDF) {
         graphics.off()
      }
       
   }  # endif bFreqScatterPlot | bScatterPlot

   #Quantile-Quantile
   if (bQQPlot) {  
      if (bPDF) {
         cat("Sending concentration assessment quantile-quantile plots to pdf...\n")
         pdf(paste(eval.output.root,"_conc_qq_plot.pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
      } else if (bJPG) {
         cat("Sending concentration assessment quantile-quantile plots to jpg...\n")
      } else if (bPNG) {
         cat("Sending concentration assessment quantile-quantile plots to png...\n")
      }
      
      #Loop over pollutants
      for (j in 1:length(pollutants$pollutant)) {
         #to make conditioning work, need "mod" as column and rows corresponding to models
         bDoneFirstModel = FALSE
         for (m in (1:length(modelled.data.names))) {
            temp.data.onemodel = subset(all.data.averaged, pollutant==pollutants$pollutant[j], select=c("date","obs","station","type","pollutant","output.avg.time.and.statistic",modelled.data.names[m]))      
            names(temp.data.onemodel)[names(temp.data.onemodel)==modelled.data.names[m]] = "mod"
            if (length(temp.data.onemodel$mod)==0) {
              next
            }
            temp.data.onemodel$model = modelled.data.names[m]
            if (bDoneFirstModel) {
               all.data.by.pol = rbind(all.data.by.pol,temp.data.onemodel)
            } else {
               all.data.by.pol = temp.data.onemodel
               bDoneFirstModel = TRUE
            }         
            rm(temp.data.onemodel)
         }
         rm(temp)  
         #Check that there is valid data to plot and all.data.by.pol exists
         if (!bDoneFirstModel) {
            cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
            next
         } 
           
         all.data.by.pol = subset(all.data.by.pol, !is.na(obs) & !is.na(mod))
         
         #Check that there is data to plot and return an error if not
         if (length(all.data.by.pol$pollutant)==0) {
            cat(paste("   Warning: no data to graph for pollutant", pollutants$pollutant[j], "- skipping this pollutant\n"))
            next
         }   
  
         #can't plot zeros on log scale
         if (bQQLog & any(all.data.by.pol[,c("obs","mod")] == 0)) {
            all.data.by.pol = subset(all.data.by.pol, (obs!=0 & mod!=0))
         }
         
         #setup pollutants
         pollutant.units = (subset(pollutants, pollutant==pollutants$pollutant[j]))$output.units         
         
         #Set graphical parameters
         graph.obs = "obs"
         graph.mod = "mod"
         pollutant.max = max(all.data.by.pol$obs,all.data.by.pol$mod,na.rm=TRUE)
         pollutant.min = min(all.data.by.pol$obs,all.data.by.pol$mod,na.rm=TRUE)
         if (bQQLog) {
            max.val = 10^ceiling(log10(1.05*pollutant.max))
            min.val = 10^floor(log10(0.95*pollutant.min))
         } else {   
            # define linear axis extent and labels
            max.val = 1.05*pollutant.max
            min.val = 0 
         }   
         pol.display = paste(" (",pollutant.units,")",sep="")
         graph.z=NA
         QQ.cols="black"
         QQ.pch=16
         scatter.cex=NULL
         scatter.mod.line=TRUE
         QQ.method="scatter"
         filter.type="default"
         
         # sort QQ data (obs, mod) separately
         # has to be filtered by model, so needs ordering within 'model' group
         # if filtered by type, each subset needs ordering within the 'type' of station
         bDoneFirstModel = FALSE
         for (m in 1:length(modelled.data.names)) {
            data.bymodel = subset(all.data.by.pol, model==unique(all.data.by.pol$model)[m]) 
            if (bAllStations & bFilterStnType) {
               for (k in 1:length(unique(all.data.by.pol$type))) {
                  data.bytype = subset(data.bymodel, type==unique(all.data.by.pol$type)[k], select=c("obs","mod","type","model"))
                  # sort QQ data (obs, mod) separately within each 'type'
                  data.bytype$obs = sort(data.bytype$obs, decreasing = TRUE)
                  data.bytype$mod = sort(data.bytype$mod, decreasing = TRUE)              
                  if (k == 1) {
                     qq.plot.data.bymodel = data.bytype
                  } else {
                     qq.plot.data.bymodel = rbind(qq.plot.data.bymodel, data.bytype)
                  }              
               }
            } else {
               # sort QQ data (obs, mod) separately within each 'type'
               qq.plot.data.bymodel = data.bymodel[ , c("obs", "mod", "model")]
               qq.plot.data.bymodel$obs = sort(data.bymodel$obs, decreasing = TRUE)
               qq.plot.data.bymodel$mod = sort(data.bymodel$mod, decreasing = TRUE)              
            }                
            if (!bDoneFirstModel) {
               qq.plot.data = qq.plot.data.bymodel
               bDoneFirstModel = TRUE
            } else {
               qq.plot.data = rbind(qq.plot.data, qq.plot.data.bymodel)            
            } 
            rm(qq.plot.data.bymodel, data.bymodel, data.bytype)           
         }         
         # setup plot options:         
         if (bAllStations & bFilterStnType) {
            # use scatterPlot's intrinsic filter function to display each station type on one page, one page per pollutant.
            # R doesn't like the column header "type" as it is also an intrinsic function
            names(qq.plot.data)[names(qq.plot.data)=="type"] = "station.type"
            filter.type=c("model","station.type")
            # titles
            text.stntype.title = "ALL STATIONS"
            text.stntype.filename = "allstationsbytype"
         } else {
            filter.type="model" 
            # titles
            if (bSingleStationType) {
              text.stntype.title = single.station.type
              text.stntype.filename = single.station.type
            } else if (bAllStations & !bFilterStnType) {
              text.stntype.title = "ALL STATIONS"
              text.stntype.filename = "allstations"
            } 
         }
         if (bSingleModel) {
            text.filters = paste(single.model,", ",toupper(text.stntype.title),", ",toupper(pollutants$output.avg.time.and.statistic[j])," ",toupper(pollutants$pollutant[j]),pol.display,sep="")
         } else {
            text.filters = paste(toupper(text.stntype.title),", ",toupper(pollutants$output.avg.time.and.statistic[j])," ",toupper(pollutants$pollutant[j]),pol.display,sep="")
         }
         # setup output page
         if (bJPG) {
            jpeg(paste(eval.output.root,"_",pollutants$pollutant[j],"_",text.stntype.filename,"_conc_qq_plot.jpg",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
         } else if (bPNG) {
            png(paste(eval.output.root,"_",pollutants$pollutant[j],"_",text.stntype.filename,"_conc_qq_plot.png",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
         } 
         #title
         titletext=paste("Quantile-Quantile Plot: ",toupper(project.name),", ",text.filters,sep="")  
         if (nchar(titletext)>65) {  # split onto two rows
            titletext=paste(paste("Quantile-Quantile Plot: ",toupper(project.name), sep=""), text.filters, sep="\n")  
         }
         
         # Draw quantile-quantile plot
         scatterPlot(
            qq.plot.data,
            x=graph.obs,y=graph.mod,z=graph.z,
            xlim=c(min.val,max.val),ylim=c(min.val,max.val),
            log.x=bQQLog, log.y=bQQLog,
            xlab="Observed",ylab="Modelled",
            cols=QQ.cols,pch=QQ.pch,cex=scatter.cex,
            type=filter.type, 
            method=QQ.method,
            mod.line=scatter.mod.line,
            main=titletext,
            sub=paste("Date range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep=""),
            aspect="iso"
         )
                
         if (bJPG | bPNG) {
            graphics.off()
         }           
         
         rm(all.data.by.pol, qq.plot.data)   
	   } # next pollutant
	     
	   if (bPDF) {
         graphics.off()
      }
       
   }  # endif bQQPlot
   
   if (bTimePlot) {
     # note "height" is longest length, for square output use "width" twice
     if (bPaperLarge) {  
       paper.height = 297
       paper.width = 210
     } else if (bPaperMedium) {
       paper.height = round(170*sqrt(2))
       paper.width = 170   
     } else if (bPaperSmall) { # set square plots to width*width, mimics A6 but a bit smaller
       paper.height = round(125*sqrt(2))
       paper.width = 125   
     }
     if (bPDF) {
       cat("Sending concentration assessment time plots to pdf...\n")
       pdf(paste(eval.output.root,"_conc_timeplot.pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
     } else if (bJPG) {
       cat("Sending concentration assessment time plots to jpg...\n")
       jpeg(paste(eval.output.root,"_conc_timeplot.jpg",sep=""),width=300,height=150,units="mm",res=300)
     } else if (bPNG) {
       cat("Sending concentration assessment time plots to png...\n")
       png(paste(eval.output.root,"_conc_timeplot.png",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
     }
     
     timeplot.data = aggregate(subset(all.data.averaged,select=c("obs",modelled.data.names)),
                               list(pollutant=all.data.averaged$pollutant,date=all.data.averaged$date),
                               FUN=mean,na.rm=TRUE)
     all.data.averaged.na.omit = na.omit(all.data.averaged)
     stations.per.datapoint = aggregate(subset(all.data.averaged.na.omit,select=c("obs",modelled.data.names)),
                                        list(pollutant=all.data.averaged.na.omit$pollutant,date=all.data.averaged.na.omit$date),
                                        FUN=length)
     names(stations.per.datapoint)[names(stations.per.datapoint)=="obs"]="num.obs.stations"
     for (i in 1:length(modelled.data.names)) {
       names(stations.per.datapoint)[names(stations.per.datapoint)==modelled.data.names[i]] = paste("num.",modelled.data.names[i],".stations",sep="")
     }         
     timeplot.data.plot = merge(timeplot.data,stations.per.datapoint,by=c("date","pollutant"))
     timeplot.data.csv = timeplot.data.plot
     titletext=paste("Time Plot: ",toupper(project.name),"Average concentration over all stations")  
     for (i in 1:length(pollutants$pollutant)) {
       if (pollutants$pollutant[i]=="nox") {next}
       if (is.null(pollutants$i4[i]) | is.null(pollutants$i7[i])) {next} #No index information available
       moderate.threshold=pollutants$i4[i]
       high.threshold=pollutants$i7[i]
       conc.values=subset(timeplot.data.plot,select=c("obs",modelled.data.names),pollutant==pollutants$pollutant[i])
       max.conc=range(conc.values)[2]
       if (moderate.threshold < 1.5*max.conc) {timeplot.data.plot$moderate.threshold[timeplot.data.plot$pollutant==pollutants$pollutant[i]]=moderate.threshold}
       if (high.threshold < 1.5*max.conc) {timeplot.data.plot$high.threshold[timeplot.data.plot$pollutant==pollutants$pollutant[i]]=high.threshold}       
       timeplot.data.csv$moderate.threshold[timeplot.data.csv$pollutant==pollutants$pollutant[i]]=moderate.threshold
       timeplot.data.csv$high.threshold[timeplot.data.csv$pollutant==pollutants$pollutant[i]]=high.threshold
     }
     if (all(is.null(pollutants$i4)) & all(is.null(pollutants$i7))) {
       pollutants.list = c(modelled.data.names,"obs","num.obs.stations")
       cols.list = c(rainbow(length(modelled.data.names),start=2/6,v=0.7),"black","slategrey")
       names.list = c(modelled.data.names,"Observed","Number of monitoring stations available")
       lty.list = c(rep(1,times=length(modelled.data.names)),1,1)
       lwd.list = c(rep(2,times=length(modelled.data.names)),2,1)
     } else {
       pollutants.list = c(modelled.data.names,"obs","num.obs.stations","moderate.threshold","high.threshold")
       cols.list = c(rainbow(length(modelled.data.names),start=2/6,v=0.7),"black","slategrey","darkorange","red")
       names.list = c(modelled.data.names,"Observed","Number of monitoring stations available","Moderate","High")
       lty.list = c(rep(1,times=length(modelled.data.names)),1,1,2,2)
       lwd.list = c(rep(2,times=length(modelled.data.names)),2,1,1.5,1.5)
     }
     if (length(pollutants$pollutant==4)) {
       plot.layout=c(2,2)
     } else {
       plot.layout=NULL
     }
     timePlot(
       timeplot.data.plot,
       pollutant=pollutants.list,
       group=TRUE,
       type="pollutant",
       key=TRUE,
       y.relation="free",
       cols=cols.list,
       name.pol=names.list,
       lty=lty.list,
       lwd=lwd.list,
       date.format="%d/%m",
       layout=plot.layout,
       key.columns=2,
       main=titletext)
     graphics.off()
     #Output data to csv
     if (boutput.stats.as.csv) {
       cat("Exporting aggregated time series to CSV...\n")
       write.csv(timeplot.data.csv,file=paste(eval.output.root,"_conc_agg.csv",sep=""),row.names = FALSE)
     }
     
   }

} #End of concentration assessment graphics

#The openair rollingMean() function seems to generate a spurious 'Rplots.pdf' file, so remove it here
if (file.exists("Rplots.pdf")) {bDum = file.remove("Rplots.pdf")}

cat("\nFinished!\n")
