#Toolkit version number
r.tool.name = "Model Diagnostics"
r.diag.toolkit.version = "5.4"
tkd.version = "tkdfileversion10"
vsd.version="VSDfileversion10"
svn.version.diag = "$ProjectRevision: SVN 13255 (2026-06-03) $"
svn.prj.name.diag = "$ProjectName: /cerc_cerc/Products/Utilities/ModelEvaluationToolkit/Tags/ModelEvaluationToolkit_5_4_1 $"


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

#Load common functions
argv = commandArgs(trailingOnly = FALSE)
install_dir = dirname(substring(argv[grep("--file=", argv)], 8))
source(paste(install_dir, "CommonFunctions.r", sep=.Platform$file.sep))
#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 Diagnostics tool.")
  quit("no")
} else if (exists("r.common.functions.version")){
  if (r.common.functions.version != r.diag.toolkit.version){
    stop("The file CommonFunctions.r is incompatible with this version of the Model Diagnostics tool.")
    quit("no")
  }
}

#Set up input and output
diag.file = commandArgs(trailingOnly=TRUE)[[1]]
# Determine if settings file is TKD or VSD
diag.file.type=tolower(substr(diag.file, nchar(diag.file)-3, nchar(diag.file)))
# Check settings file type is valid and get output.file.path.root 
diag.output.root = substr(diag.file,0,nchar(diag.file)-4)
if (diag.file.type ==".vsd"){
   diag.file.type="VSD"
} else if (diag.file.type==".tkd"){
   diag.file.type="TKD"
} else {
   stop("The settings file is not a recognised file type. Settings files for the Model Diagnostics tool must either be a TKD file or a VSD file.")
}

setwd(dirname(diag.file))
log.file = file(paste(diag.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.diag.toolkit.version,": ",r.tool.name," Tool \n",sep=""))
if (length(grep("$ProjectRevision: ",svn.version.diag,fixed=TRUE))>0) {
  svn.version.diag.numberonly=sub("$ProjectRevision: ","",svn.version.diag,fixed=TRUE)
  svn.version.diag.numberonly=sub(" $","",svn.version.diag.numberonly,fixed=TRUE)
  cat(paste("Internal version ",svn.version.diag.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="")) 

#Check length of output directory to ensure it is not too long (exceeds 200 characters) which could cause problems later
if(nchar(diag.output.root)>200){
  if (diag.file.type=="TKD"){
    stop(paste("The TKD file path has too many characters: ", nchar(diag.output.root),". This may cause issues with saving outputs later on. \nPlease reduce the length of the TKD file path to be less than 200 characters.", sep=""))
    quit("no")
  } else if (diag.file.type=="VSD"){
    stop(paste("The VSD file path has too many characters: ", nchar(diag.output.root),". This may cause issues with saving outputs later on. \nPlease reduce the length of the VSD file path to be less than 200 characters.", sep=""))
    quit("no")
  }
}

#Read in settings file
if (!file.exists(diag.file)) {
  if (diag.file.type=="TKD"){
  stop(paste("The .tkd file specified does not exist: ",diag.file, sep=""))
  quit("no")   
  } else if (diag.file.type=="VSD"){
  stop(paste("The .vsd file specified does not exist: ",diag.file, sep=""))
  quit("no")  
  }
}

# If settings file is a VSD file, remove commented out lines
# Read in data from diag file 
if(diag.file.type=="VSD"){
  VSD.data=readLines(diag.file)
  VSD.data=VSD.data[-which(startsWith(VSD.data,"#"))]  
  tkd.data=VSD.data
} else {
  tkd.data = readLines(diag.file)
}

#getting the data from tkd
tkd.data = sub("PATH = \"","",as.character(tkd.data),fixed=TRUE)
tkd.data = sub("\"","",as.character(tkd.data),fixed=TRUE)

#Check that the TKD/VSD file is for the right version
if(diag.file.type=="TKD") {  
  tkd.version.in = as.character(tkd.data[1])
  if (tkd.version.in != tkd.version){
    stop(paste("The TKD file version (",tkd.version.in,") does not match the expected version(",tkd.version,")",sep=""))
    quit("no")
  }
} else if (diag.file.type =="VSD"){
  vsd.version.in = as.character(tkd.data[1])
  if (vsd.version.in != vsd.version){
    stop(paste("The VSD file version (",vsd.version.in,") does not match the expected version(",vsd.version,")",sep=""))
    quit("no")
  }
}

#LOAD R WORKSPACE
nProject = which(tkd.data=="Project")
data.input.workspace = tkd.data[nProject+2]
cat("Loading workspace...\n")
if (data.input.workspace=="") {
  stop("No workspace selected.")
  quit("no")
}
data.loaded=tryCatch(load(data.input.workspace),error=function(e) return(NULL))
if (is.null(data.loaded)) {
  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 Diagnostics tool.")
  quit("no")
} else if (exists("r.data.toolkit.version")) {
  if (r.data.toolkit.version != r.diag.toolkit.version) {
    stop("The selected workspace is incompatible with this version of the Model Diagnostics tool.")
    quit("no")
  }
}

#Maintain compatibility with workspaces generated before the change to user-formatted station names in output (CR0117)
if (!("station.output.name" %in% names(all.data))) {
  all.data$station.output.name=all.data$station
}
if (!("station.output.name" %in% names(stations))) {
  stations$station.output.name=stations$name
}

#Finding line numbers for the data
nProject = which(tkd.data=="Project")
nRSettings = which(tkd.data=="RSettings")
nStation.Pollutant.Model = which(tkd.data=="Station, pollutant and model")
nData.to.Plot = which(tkd.data=="Data to Plot")
nPlot.Options = which(tkd.data=="Plot options")
nOutput = which(tkd.data=="Output")

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

#OTHER DATA
AllStations=tkd.data[nStation.Pollutant.Model+2]=="1"
SelectedStations=tkd.data[nStation.Pollutant.Model+2]=="0"
selected.stations=tkd.data[nStation.Pollutant.Model+3]
if(SelectedStations & selected.stations==""){
  stop("No station name entered.")
  quit("no")
}

AllPols=tkd.data[nStation.Pollutant.Model+4]=="1"
SelectedPols=tkd.data[nStation.Pollutant.Model+4]=="0"
selected.pollutants=tkd.data[nStation.Pollutant.Model+5]
if(SelectedPols & selected.pollutants ==""){
  stop("No pollutant name entered.")
  quit("no")
}

AllModels=tkd.data[nStation.Pollutant.Model+6]=="1"
SelectedModels=tkd.data[nStation.Pollutant.Model+6]=="0"
selected.models=tkd.data[nStation.Pollutant.Model+7]
if (SelectedModels & selected.models=="") {
  stop("No modelled data label entered.")
  quit("no")
}


bDateSubset = tkd.data[nStation.Pollutant.Model+8]=="Y"
start.date.string = tkd.data[nStation.Pollutant.Model+9]
end.date.string = tkd.data[nStation.Pollutant.Model+10]

#Convert date strings to date values
if (bDateSubset) {
  start.date.req = as.POSIXct(start.date.string,format="%Y,%m,%d,%H",tz="GMT")
  end.date.req = as.POSIXct(end.date.string,format="%Y,%m,%d,%H",tz="GMT")
  if (end.date.req<start.date.req) {
    stop("The end date required must be later than or the same as the start date required.")
    quit("no")    
  }
}

#AVERAGING TIME SETTINGS
nData.to.Plot = which(tkd.data=="Data to Plot")
bNoAvg = tkd.data[nData.to.Plot+2]=="0"
bConcAvg = tkd.data[nData.to.Plot+2]=="1"
bForecastIndex = tkd.data[nData.to.Plot+2]=="2"
ref.values = tkd.data[nData.to.Plot+3]
b8HourRollingMean = tkd.data[nData.to.Plot+4]=="8hrrollmean"
b8HourMean = tkd.data[nData.to.Plot+4]=="8hrmean"
bDailyMean = tkd.data[nData.to.Plot+4]=="Dailymean"
bDailyMax = tkd.data[nData.to.Plot+4]=="Dailymax"
bDailyMax8HourRollingMean = tkd.data[nData.to.Plot+4]=="Dailymax_8hrrollingmean"
data.threshold = as.numeric(tkd.data[nData.to.Plot+5])
index.scales.file = tkd.data[nData.to.Plot+6]
bPlotAlert = tkd.data[nData.to.Plot+7]=="Y"
alert.thresholds.file = tkd.data[nData.to.Plot+8]
bPlotRefValues = ref.values!=""
if (bPlotRefValues) {
  ref.values = strsplit(ref.values,",")[[1]]
  if (any(is.na(as.numeric(ref.values)))) {
    stop("Invalid reference value(s) provided.")
    quit("no")
  } else if (length(ref.values) > 10) {
    stop("Too many reference values entered.")
    quit("no")
  }
}

if (bForecastIndex) {
  if (index.scales.file=="") {
    stop("No index scale file selected.")
    quit("no")
  }
  
}

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

#DIAGNOSTIC TOOLS
bTimeVariation =tkd.data[nPlot.Options+2]=="Y"
bTimePlot = as.logical(tkd.data[nPlot.Options+3]=="Y")
bTimePlotMarkers = as.logical(tkd.data[nPlot.Options+4]=="Y")
bTimePlotLimits = as.logical(tkd.data[nPlot.Options+5]=="Y")
timeplot.ymin = as.numeric(tkd.data[nPlot.Options+6])
timeplot.ymax = as.numeric(tkd.data[nPlot.Options+7])
bScatterPlot = tkd.data[nPlot.Options+8]=="Y" & tkd.data[nPlot.Options+9]=="1"
bFreqScatterPlot = tkd.data[nPlot.Options+8]=="Y" & tkd.data[nPlot.Options+9]=="0"
scatter.filters = c("No filtering","Weekday","Month","Season")
scatter.filter = scatter.filters[as.numeric(tkd.data[nPlot.Options+10])+1]
scatter.selected = (scatter.filter == scatter.filters)
bNoFilter = scatter.selected[1]
bFilterWeekday = scatter.selected[2]
bFilterMonth = scatter.selected[3]
bFilterSeason = scatter.selected[4]
bScatterLog = as.logical(tkd.data[nPlot.Options+11]=="Y")
bScatterBestLine=as.logical(tkd.data[nPlot.Options+12]=="Y")
bScatterLimits = as.logical(tkd.data[nPlot.Options+13]=="Y")
scatter.xymin = as.numeric(tkd.data[nPlot.Options+14])
scatter.xymax = as.numeric(tkd.data[nPlot.Options+15])
bPolarPlots = as.logical(tkd.data[nPlot.Options+16]=="Y")
bPolar.filter.by.stability=tkd.data[nPlot.Options+17]=="Stability"
bPolar.filter.by.weekday=tkd.data[nPlot.Options+17]=="Weekday"
bPolar.filter.by.month=tkd.data[nPlot.Options+17]=="Month"
bPolar.filter.by.season=tkd.data[nPlot.Options+17]=="Season"
bPolar.conc.plot = as.logical(tkd.data[nPlot.Options+18]=="Y")
bPolar.radial.axis.windspeed=tkd.data[nPlot.Options+19]=="windspeed"
bPolar.radial.axis.hour=tkd.data[nPlot.Options+19]=="hour"
bPolar.radial.axis.month=tkd.data[nPlot.Options+19]=="month"
Polar.conc.plot.min.bin=as.numeric(tkd.data[nPlot.Options+20])
bPollution.rose.plot=as.logical(tkd.data[nPlot.Options+21]=="Y")
bPolar.wind.freq.plot=as.logical(tkd.data[nPlot.Options+22]=="Y")
Polar.wind.plot.min.bin=as.numeric(tkd.data[nPlot.Options+23])
if (!bTimeVariation & !bScatterPlot & !bFreqScatterPlot & !bTimePlot & !bPolar.conc.plot & !bPollution.rose.plot & !bPolar.wind.freq.plot) {
  stop("No diagnostic graph output selected.")
  quit("no")
}
if (bPolarPlots) {
  if (!exists("met.data")){
    stop("You have chosen to output Polar Plots however the workspace you have specified does not contain met data.\n")
    quit("no")
  }
  if (!bPolar.conc.plot & !bPollution.rose.plot & !bPolar.wind.freq.plot) {
    stop("Polar plots option selected but no polar plot graphs selected.")
    quit("no")
  }
  if (!bNoAvg){
    stop(paste("Polar plots can only plot raw concentration data, not calculated statistics or forecast indices.\n"))
    quit("no")
  }
}
# Set met data parameter
if (!exists("MetDataFormat")){
  stop("Cannot find parameter 'MetDataFormat' in R workspace. Please check the workspace.\n")
  quit("no")
}
PolarPlotsMetFormat=MetDataFormat
# Check the filtering of the polar plots is compatibility with the meteorological  data type 
if (bPolarPlots & PolarPlotsMetFormat %in% c("wrf","ADMSmetnc","met") & bPolar.filter.by.stability){
  stop("You have chosen to filter the Polar Plots by atmospheric stability however the meteorological data format does not contain sufficient information for calculating atmospheric stability.\nPlease change the filtering option for the Polar Plots.")
  quit("no")
}
# Check time plot limits are valid, if they have been requested
if (bTimePlot & bTimePlotLimits) {
  if (is.na(timeplot.ymin)) {
    stop("The minimum time plot y-axis limit is invalid. Ensure the value is numeric.\n") # Interface already forces it to be numeric, this is to protect against users modifying TKD files directly.
    quit("no")
  }
  if (is.na(timeplot.ymax)) {
    stop("The maximum time plot y-axis limit is invalid. Ensure the value is numeric.\n") # Interface already forces it to be numeric, this is to protect against users modifying TKD files directly.
    quit("no")
  }
  if (timeplot.ymin >= timeplot.ymax) {
    stop("The maximum time plot y-axis limit must be greater than the minimum y-axis limit.\n") # There is no check in the interface that forces min to be less than max, checking here.
    quit("no")
  }
}
# Check scatter plot limits are valid, if they have been requested
if ((bScatterPlot | bFreqScatterPlot) & bScatterLimits) {
  if (is.na(scatter.xymin)) {
    stop("The minimum scatter plot x and y axes limits are invalid. Ensure the value is numeric.\n") # Interface already forces it to be numeric, this is to protect against users modifying TKD files directly.
    quit("no")
  }
  if (is.na(scatter.xymax)) {
    stop("The maximum scatter plot x and y axes limits are invalid. Ensure the value is numeric.\n") # Interface already forces it to be numeric, this is to protect against users modifying TKD files directly.
    quit("no")
  }
  if (scatter.xymin >= scatter.xymax) {
    stop("The maximum scatter plot x and y axes limits must be greater than the minimum x and y axes limits.\n") # There is no check in the interface that forces min to be less than max, checking here.
    quit("no")
  }
}

output.types = c("JPG","PNG","PDF")
output.type = (output.types[as.numeric(tkd.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(tkd.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]

#R CODE:

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

#Create directories based on output plots
if (bTimeVariation){
  time.var.dir=(paste(dirname(diag.output.root),.Platform$file.sep,"Time_var_plots",sep=""))
  if (!dir.exists(time.var.dir)){
    dir.create(time.var.dir) 
  }
}
if (bTimePlot){
  time.plot.dir=paste(dirname(diag.output.root),.Platform$file.sep,"Time_plots",sep="")
  if (!dir.exists(time.plot.dir)){
    dir.create(time.plot.dir) 
  }
} 
if(bScatterPlot | bFreqScatterPlot){
  scatter.plot.dir=paste(dirname(diag.output.root),.Platform$file.sep,"Scatter_plots",sep="")
  if (!dir.exists(scatter.plot.dir)){
    dir.create(scatter.plot.dir) 
  }
}
if (bPolarPlots & bPolar.conc.plot){
  polar.conc.plot.dir=paste(dirname(diag.output.root),.Platform$file.sep,"Polar_conc_plots",sep="")
  if (!dir.exists(polar.conc.plot.dir)){
    dir.create(polar.conc.plot.dir) 
  }
}
if (bPolarPlots & bPollution.rose.plot){
  pollution.rose.plot.dir=paste(dirname(diag.output.root),.Platform$file.sep,"Pollution_rose_plots",sep="")
  if (!dir.exists(pollution.rose.plot.dir)){
    dir.create(pollution.rose.plot.dir) 
  }
}
if (bPolarPlots & bPolar.wind.freq.plot){
  polar.wind.freq.dir=paste(dirname(diag.output.root),.Platform$file.sep,"Polar_ws_freq_plots",sep="")
  if (!dir.exists(polar.wind.freq.dir)){
    dir.create(polar.wind.freq.dir) 
  }
}
diag.file.name=basename(diag.output.root)
  

#get list of selected stations from the tki file and compare with stations in the data
unavailable.stations=c()
available.stations=unique(all.data$station) #list of stations present in the data 
if(SelectedStations) {
  selected.stations.vector=gsub(", ",",",selected.stations)
  selected.stations.vector=unlist(strsplit(selected.stations.vector,","))#create a vector of selected stations from tki file
  selected.stations=selected.stations.vector
  #get list of selected stations that are not available
  unavailable.stations=selected.stations[which(!(tolower(selected.stations) %in% tolower(available.stations)))]#Find specified stations that are not present in the data 
  #if there are any unavailable stations then issue an error message and stop later
  if (!(length(unavailable.stations)==0)) { #write error message to log if there are unavailable stations, do not stop yet, check mod and pols first
    unavailable.stations.str=paste(unavailable.stations,collapse=",")
    available.stations.str=paste(available.stations,collapse=",")
    cat(paste("Error: No data for selected stations: ", unavailable.stations.str,
              ". Data are available for the following stations: ",available.stations.str, ".\n", sep=""))
  } else {
    for (i in 1:length(selected.stations)){
      selected.stations[i]=available.stations[which(tolower(available.stations) %in% tolower(selected.stations[i]))]
    }
  }
}

#get list of selected models from tki file and compare with models in the data 
unavailable.models=c()
modelled.data.names = names(all.data)[which(!(names(all.data) %in% c("date","obs", "station", "pollutant","station.output.name","type","network","auto.download")))] #find modelled data names from data
if (SelectedModels) {
  selected.models.vector=gsub(", ",",",selected.models)
  selected.models.vector=unlist(strsplit(selected.models.vector,","))#create a vector of selected models from tki file
  selected.models=selected.models.vector
  #get list of selected models that are not available
  unavailable.models=selected.models[which(!(tolower(selected.models) %in% tolower(modelled.data.names)))]#Find specified models that are not present in the data 
  if (!(length(unavailable.models)==0)) { #write error message to log if there are unavailable models, do not stop yet, check pols
    unavailable.models.str=paste(unavailable.models,collapse=",")
    available.models.str=paste(modelled.data.names,collapse=",")
    cat(paste("Error: No data for selected models: ", unavailable.models.str,
              ". Data are available for the following models: ",available.models.str, ".\n", sep=""))
  } else {
    for (i in 1:length(selected.models)){
      selected.models[i]=modelled.data.names[which(tolower(modelled.data.names) %in% tolower(selected.models[i]))]
    }
  }
}

#get list of selected pollutants
unavailable.pollutants=c()
available.pollutants=pollutants$pollutant #list pols availabel in data
if(SelectedPols) {
  selected.pollutants.vector=gsub(", ",",",selected.pollutants)
  selected.pollutants.vector=unlist(strsplit(selected.pollutants.vector,","))#create a vector of selected pols from tki file
  selected.pollutants=selected.pollutants.vector
  #get list of selected pollutants that are not available 
  unavailable.pollutants=selected.pollutants[which(!(tolower(selected.pollutants) %in% tolower(available.pollutants)))]#Find specified pols that are not present in the data 
  if (!(length(unavailable.pollutants)==0)) {  #write error message to log if there are unavailable pollutants
    unavailable.pollutants.str=paste(unavailable.pollutants,collapse=",")
    available.pollutants.str=paste(available.pollutants,collapse=",")
    cat(paste("Error: No data for selected pollutants: ", unavailable.pollutants.str,
              ". Data are available for the following pollutants: ",available.pollutants.str, ".\n", sep=""))
  } else {
    for (i in 1:length(selected.pollutants)){
      selected.pollutants[i]=available.pollutants[which(tolower(available.pollutants) %in% tolower(selected.pollutants[i]))]
    }
  }
}

#stop if there are any unavailable pols, mods or stations. Error messages will have already been written in messages above
if ((!(length(unavailable.pollutants)==0))|(!(length(unavailable.models)==0))|(!(length(unavailable.stations)==0))){
  stop()
  quit("no")
}

#Subset by date if required
if (bDateSubset) {
  cat(paste("Selecting data from",format(start.date.req,format="%d/%m/%y %H:%M",tz="GMT"),"to",format(end.date.req,format="%d/%m/%y %H:%M",tz="GMT"),"inclusive\n"))
  all.data.for.plotting=subset(all.data,(date >= start.date.req) & (date <= end.date.req))
  if (length(all.data.for.plotting$date)==0) {
    stop("No data for selected date period")
  }
} else {
  all.data.for.plotting=all.data
}

#if all pollutants,stations or models are selected, set the data as all available
if (AllStations){
  selected.stations=available.stations
}
if (AllModels){
  selected.models=modelled.data.names
}
if (AllPols){
  selected.pollutants=available.pollutants
}
#Reduce full dataset to only those pollutants and stations selected
all.data.for.plotting = subset(all.data.for.plotting,
                               station %in% selected.stations &
                                 pollutant %in% selected.pollutants)
if (length(all.data.for.plotting$date)==0) {
  stop("No data for selected stations and pollutants")
}

# If producing polar plots, check that the met data date range overlaps 
if (bPolarPlots){
  #subset met data to data.for.plotting range
  names(met.data)=tolower(names(met.data))
  if (PolarPlotsMetFormat %in% c("wrf", "ADMSmetnc")) {
    #reduce dataset to only those selected stations
    met.data=subset(met.data,station %in% selected.stations)
  }
  #reduce dataset to only overlapping hours
  met.data.overlap=subset(met.data, date>=range(all.data.for.plotting$date)[1] & date<=range(all.data.for.plotting$date)[2])
  #If there isn't enough met data for polar plots, disable this option
  if (length(met.data.overlap$date)==0){
    cat("\nWarning: Meteorological data is not available for any of the analysis period, so polar plots will not be produced.\n")
    cat(paste("         Analysis range: ",paste(range(all.data.for.plotting$date), collapse=": "),"\n"))
    cat(paste("         Met data range: ",paste(range(met.data$date), collapse=": "),"\n\n"))
    bPolarPlots = FALSE
  } else {
    met.data = met.data.overlap
    if (range(met.data$date)[1]!=range(all.data.for.plotting$date)[1] |
        range(met.data$date)[2]!=range(all.data.for.plotting$date)[2]){
      cat("\nWarning: Polar plots have been selected but meteorological data is only available for part of the analysis period,\n")
      cat("         therefore the analysis period will be truncated to match the meteorological data period for all plot types.\n")
      cat(paste("         Analysis period: ",paste(range(all.data.for.plotting$date), collapse=": "),"\n"))
      cat(paste("         Met data period: ",paste(range(met.data.overlap$date), collapse=": "),"\n"))
      cat(paste("         New analysis period: ",paste(range(met.data$date), collapse=": "),"\n\n"))
    }  
  }
}
if (!bTimeVariation & !bScatterPlot & !bFreqScatterPlot & !bTimePlot & !bPolarPlots) {
  stop("No diagnostic graph output to plot.")
  quit("no")
}

#Load the required packages
argv = commandArgs(trailingOnly = FALSE)
install_dir = dirname(substring(argv[grep("--file=", argv)], 8))
packages.required = c("openair","hexbin","nlme","reshape2","tibble","stringr", "ggplot2")
install.load.packages(packages.required)
cat("\n")

# Openair versions 3.1 or higher is required to run the Toolkit due to plotting changes
if (packageVersion("openair") < "3.1.0") {
  stop("Openair v3.1.0 or higher is required to run the Toolkit. \nTo update packages go to 'R Support > Update R packages' in the Toolkit interface.\n")
  quit("no")
}

if (bForecastIndex){ #Add the index scales data to the pollutants data frame if needed
  pollutants = process.index.scales(index.scales.file, pollutants)
}

#set value of original plot parameters
bFreqScatterPlotAll=bFreqScatterPlot
bScatterPlotAll=bScatterPlot
bTimePlotAll=bTimePlot
bTimeVariationAll=bTimeVariation
bPolarPlotsAll=bPolarPlots

num.selected.models=length(selected.models)
num.selected.pollutants=length(selected.pollutants)
num.selected.stations=length(selected.stations)
for(m in 1:num.selected.models) {
  single.model=selected.models[m]
  data.for.plotting.model=subset(all.data.for.plotting, select=c("date","obs","station","station.output.name","pollutant","type",single.model))
  cat(paste("Processing data for model ",m," of ",num.selected.models,": ",single.model,"...\n", sep=""))
  
  if (nrow(data.for.plotting.model)==0 ){
    cat(paste("  Warning: No data to graph for model", single.model, "- skipping this model\n"))
    next
  } else {
    for(j in 1:num.selected.pollutants) {
      
      single.pollutant=selected.pollutants[j]
      cat(paste("  Processing data for pollutant ",j," of ",num.selected.pollutants,": ",single.pollutant,"...\n", sep=""))
      single.pollutant.lower=tolower(single.pollutant)
      single.pollutant.upper=toupper(single.pollutant)
      data.for.plotting.pollutant=subset(data.for.plotting.model, pollutant %in% single.pollutant.lower)
      
      if (nrow(data.for.plotting.pollutant)==0) {
        cat(paste("    Warning: No data to graph for pollutant", single.pollutant,"for model",single.model, "- skipping this pollutant for this model\n"))
        next
      } else  {
        for(s in 1:num.selected.stations) {
          #reset values of b output plot parameters to original values  
          bFreqScatterPlot=bFreqScatterPlotAll
          bScatterPlot=bScatterPlotAll
          bTimePlot=bTimePlotAll
          bTimeVariation=bTimeVariationAll
          bPolarPlots=bPolarPlotsAll
          single.station=selected.stations[s]
          cat(paste("    Processing data for station ",s," of ",num.selected.stations,": ",single.station,"...\n", sep=""))
          single.station.upper=toupper(single.station)
          single.station.lower=tolower(single.station)
          data.for.plotting=subset(data.for.plotting.pollutant, station %in% single.station.lower)
          # data.for.plotting.complete retains NAs
          data.for.plotting.complete = data.for.plotting
          data.for.plotting.complete[is.na(data.for.plotting.complete[["obs"]]),single.model] = NA
          data.for.plotting.complete[is.na(data.for.plotting.complete[[single.model]]),"obs"] = NA
          # remove NAs from beginning and end of data
          valid.indices = which(!is.na(data.for.plotting.complete$obs))
          data.for.plotting.complete = tryCatch(data.for.plotting.complete[valid.indices[1]:valid.indices[length(valid.indices)],], 
                                                 error=function(e) return(NA)) # if no valid data, return NA, and error will be caught in data.for.plotting check below
          data.for.plotting=data.for.plotting[which(!is.na(data.for.plotting[[single.model]])),]
          if (nrow(data.for.plotting)==0) {
            cat(paste("      Warning: No modelled data for pollutant", single.pollutant,"for model",single.model,"for station",single.station, "- skipping this pollutant for this model for this station\n"))
            next
          }
          data.for.plotting=data.for.plotting[which(!is.na(data.for.plotting[["obs"]])),]
          if (nrow(data.for.plotting)==0) {
            cat(paste("      Warning: No observed data for pollutant", single.pollutant,"for model",single.model,"for station",single.station, "- skipping this pollutant for this model for this station\n"))
            next
          } else {
            if (bPolarPlots) {
              # We want all the graphs to cover the same period, so merge in the met data here before doing anything else
              # Get station met data
              if (PolarPlotsMetFormat %in% c("wrf", "ADMSmetnc")){
                station.met.data=subset(met.data, station %in% single.station)
              } else {
                station.met.data=met.data
              }
              station.met.data = subset(station.met.data,!is.na(ws) & !is.na(wd))
              if (length(station.met.data$date)==0){
                cat(paste("      Warning: No valid met data for station",single.station, "- skipping this pollutant for this model for this station\n"))
                next
              } else {
                data.for.plotting = merge(data.for.plotting,station.met.data,by="date")
              }
            }
            data.start.date.for.graphs = format(range(data.for.plotting$date)[1],"%d/%m/%y %H:%M")
            data.end.date.for.graphs = format(range(data.for.plotting$date)[2], "%d/%m/%y %H:%M")
            
            #Change column headings for openair functions
            names(data.for.plotting)[names(data.for.plotting)==single.model] = "mod"
            names(data.for.plotting.complete)[names(data.for.plotting.complete)==single.model] = "mod"
            pol.units = pollutants$output.units[pollutants$pollutant==single.pollutant.lower]
            single.station.output.name = unique(data.for.plotting$station.output.name)
            #Apply settings to data for graphs
            if (bFreqScatterPlot | bScatterPlot) {
              #If a filtering option has been selected, add new column containing filter
              if (bNoFilter) {
                filter.type = "default"
                filter.name=", No Filtering"
              } else if (bFilterWeekday) {
                filter.type="weekday"
                filter.name=", Filtered by Weekday"
              } else if (bFilterMonth) {
                filter.type="month"
                filter.name=", Filtered by Month"
              } else if (bFilterSeason) {
                filter.type="season"
                filter.name=", Filtered by Season"
              }
            }
            
            if (bNoAvg) {
              avg.name = ", Raw Data"
              #Set data for plotting for each graph
              time.variation.data = data.for.plotting
              scatter.plot.data = data.for.plotting
              #remove NAs
              scatter.plot.data = subset(scatter.plot.data, !(is.na(obs) | is.na(mod)) )
              time.plot.data = data.for.plotting.complete
              polar.plot.data = data.for.plotting
              
            } else if (bConcAvg) {
              #Check that the averaging time is suitable for applying further averaging, i.e. 1 hour mean
              pollutant.data = subset(pollutants,pollutant==single.pollutant.lower)
              if (!(pollutant.data$statistic=="mean" & pollutant.data$avg.time %in% c("hour","1 hour"))) {
                stop("Input data must be hourly in order to plot calculated statistics.")
                quit("no")
              }
              #Apply averaging time to data
              if (b8HourRollingMean | bDailyMax8HourRollingMean) {
                avgTime = 8
                avg.data.for.graph.obs = rollingMean(data.for.plotting, pollutant="obs", width=avgTime, new.name = "rolling.obs",align="right",data.thresh=data.threshold)
                avg.data.for.graph.obs = subset(avg.data.for.graph.obs,select=c("date","station","rolling.obs"))
                names(avg.data.for.graph.obs)[names(avg.data.for.graph.obs)=="rolling.obs"]="obs"
                avg.data.for.graph.mod = rollingMean(data.for.plotting, pollutant="mod", width=avgTime, new.name = "rolling.mod",align="right",data.thresh=data.threshold)
                avg.data.for.graph.mod = subset(avg.data.for.graph.mod,select=c("date","station","rolling.mod"))
                names(avg.data.for.graph.mod)[names(avg.data.for.graph.mod)=="rolling.mod"]="mod"      
                avg.data.for.graph = merge(avg.data.for.graph.obs,avg.data.for.graph.mod)
                
                # Apply 8-hour rolling mean to complete dataset
                complete.obs = rollingMean(data.for.plotting.complete, pollutant="obs", width=avgTime, new.name = "rolling.obs", align="right", data.thresh=data.threshold)
                complete.obs = subset(complete.obs, select=c("date","station","rolling.obs"))
                names(complete.obs)[names(complete.obs)=="rolling.obs"] = "obs"
                complete.mod = rollingMean(data.for.plotting.complete, pollutant="mod", width=avgTime, new.name = "rolling.mod", align="right", data.thresh=data.threshold)
                complete.mod = subset(complete.mod, select=c("date","station","rolling.mod"))
                names(complete.mod)[names(complete.mod)=="rolling.mod"] = "mod"
                avg.data.for.graph.complete = merge(complete.obs, complete.mod)
                
                if (bDailyMax8HourRollingMean) {
                  avg.name = ", Daily Maximum of 8-Hour Rolling Mean"
                  avg.data.for.graph = timeAverage(avg.data.for.graph,avg.time="day",statistic="max",data.thresh=data.threshold)
                  avg.data.for.graph.complete = timeAverage(avg.data.for.graph.complete, avg.time="day", statistic="max", data.thresh=data.threshold)
                } else {
                  avg.name = ", 8-Hour Rolling Mean"
                }
              } else {
                if (b8HourMean) {
                  avgTime = "8 hour"
                  avgStatistic = "mean"
                  avg.name = ", 8-Hour Mean"
                } else if (bDailyMean) {
                  avgTime = "day"
                  avgStatistic = "mean"
                  avg.name = ", Daily Mean"
                } else if (bDailyMax) {
                  avgTime = "day"
                  avgStatistic = "max"
                  avg.name = ", Daily Maximum"
                }
                avg.data.for.graph = timeAverage(data.for.plotting,avg.time=avgTime,statistic=avgStatistic,data.thresh=data.threshold)
                # data retaining NA values
                avg.data.for.graph.complete = timeAverage(data.for.plotting.complete,avg.time=avgTime,statistic=avgStatistic,data.thresh=data.threshold)
              }
              
              #remove NAs
              data.for.plotting=subset(data.for.plotting, !(is.na(obs) | is.na(mod)) )
              avg.data.for.graph=subset(avg.data.for.graph, !(is.na(obs) | is.na(mod)) )
              # remove NAs from beginning and end of data
              valid.indices = which(!is.na(avg.data.for.graph.complete$obs))
              # if no valid data, return original data, and error will be caught in time.plot.data check below
              avg.data.for.graph.complete = tryCatch(avg.data.for.graph.complete[valid.indices[1]:valid.indices[length(valid.indices)],], 
                                                     error=function(e) return(avg.data.for.graph.complete)) 
              #Set data for plotting for each graph
              if (bTimeVariation) time.variation.data = data.for.plotting
              if (bScatterPlot | bFreqScatterPlot) scatter.plot.data = avg.data.for.graph
              if (bTimePlot) time.plot.data = avg.data.for.graph.complete
              if (bPolarPlots) polar.plot.data = data.for.plotting
              
            } else if (bForecastIndex) {
              bScatterLog = FALSE
              avg.name = ""
              index.lower.bounds = find.index.lower.bounds(pollutants)
              pollutant.to.index =  which(pollutants$pollutant==single.pollutant.lower)
              names(data.for.plotting)[names(data.for.plotting)=="station"] = "site"
              names(data.for.plotting.complete)[names(data.for.plotting.complete)=="station"] = "site"
              data.for.index.graph = find.data.for.indexing(data.for.plotting, pollutants, pollutant.to.index, index.lower.bounds, data.threshold)
              data.for.index.graph = data.for.index.graph[!is.na(data.for.index.graph$index.obs)& !is.na(data.for.index.graph$index.mod),]
              data.for.index.graph.complete = find.data.for.indexing(data.for.plotting.complete, pollutants, pollutant.to.index, index.lower.bounds, data.threshold)
              # remove NAs from beginning and end of data
              valid.indices = which(!is.na(data.for.index.graph.complete$index.obs))
              # if no valid data, return original data, and error will be caught in time.plot.data check below
              data.for.index.graph.complete = tryCatch(data.for.index.graph.complete[valid.indices[1]:valid.indices[length(valid.indices)],], 
                                                      error=function(e) return(data.for.index.graph.complete))
              #If scatter plot is selected, calculate frequencies of each x,y combination
              
              if (bScatterPlot | bFreqScatterPlot) {
                #If a filtering option has been selected, add new column containing filter
                data.for.scatter.only = cutData(data.for.index.graph, type=filter.type)
                names(data.for.scatter.only)[names(data.for.scatter.only)==filter.type]="filter"
                data.for.scatter.only = as.data.frame(table(data.for.scatter.only$index.obs, data.for.scatter.only$index.mod,data.for.scatter.only$filter))
                names(data.for.scatter.only) = c("index.obs","index.mod","filter","count")
                data.for.scatter.only$index.obs=as.numeric(as.character(data.for.scatter.only$index.obs))
                data.for.scatter.only$index.mod=as.numeric(as.character(data.for.scatter.only$index.mod))
                
                scatter.plot.data = subset(data.for.scatter.only,count>0)
                filter.type = "filter"
              }
              
              #Set data for plotting for graphs
              if (bTimeVariation) time.variation.data = data.for.plotting   
              if (bTimePlot) time.plot.data = data.for.index.graph.complete
              if (bPolarPlots) polar.plot.data = data.for.plotting
            }
            
            
            if (bTimeVariation) {
              #Need at least one data point for a time variation plot
              if (nrow(time.variation.data)<1){
                bTimeVariation=FALSE
                cat(paste("      Warning: Insufficient data for time variation plot - skipping this plot for this pollutant, model and station.\n",sep=""))
              }
            }
            if (bScatterPlot | bFreqScatterPlot){
              #Need at least two data points for a scatter plot
              if (nrow(scatter.plot.data)<1){
                bScatterPlot=FALSE
                bFreqScatterPlot=FALSE
                cat(paste("      Warning: Insufficient data for scatter plot - skipping this plot for this pollutant, model and station.\n",sep=""))
              }
            }
            if(bTimePlot) {
              #Need at least 2 data points for a time plot
              if (nrow(na.omit(time.plot.data))<2) {
                bTimePlot=FALSE
                cat(paste("      Warning: Insufficient data for time plot - skipping this plot for this pollutant, model and station.\n",sep=""))
              }
            }
            if(bPolarPlots) {
              if (nrow(polar.plot.data)==0){
                bPolarPlots=FALSE
                cat(paste("      Warning: Insufficient data for polar plots - skipping these plots for this pollutant, model and station.\n",sep=""))
              }
            }
            
            #remove zeros if plotting to log_10
            if (bScatterLog) {
              if (any(scatter.plot.data[,c("obs","mod")] == 0)) {
                scatter.plot.data = subset(scatter.plot.data, (obs!=0 & mod!=0))
              }
            }
            
            #Diagnostic graphs
            
            #Set graphical parameters
            #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
              max.title.char.length = 90
              max.title.char.length.scatter = 80
            } else if (bPaperMedium) {
              paper.height = round(170*sqrt(2))
              paper.width = 170   
              max.title.char.length = 90
              max.title.char.length.scatter = 65
            } else if (bPaperSmall) { # set square plots to width*width, mimics A6 but a bit smaller
              paper.height = round(125*sqrt(2))
              paper.width = 125   
              max.title.char.length = 65
              max.title.char.length.scatter = 50
            }
            if (!bForecastIndex) {
              if (bPlotRefValues) {
                ref.selected.x = list(v=as.numeric(ref.values),lty=3)
                ref.selected.y = list(h=as.numeric(ref.values),lty=3)
              } else {
                ref.selected.x = NULL
                ref.selected.y = NULL
              }
              graph.obs = "obs"
              graph.mod = "mod"
              if (bNoAvg) {
                pollutant.max = max(data.for.plotting$obs,data.for.plotting$mod,na.rm=TRUE)
                pollutant.min = min(data.for.plotting$obs,data.for.plotting$mod,na.rm=TRUE)
              } else if (bConcAvg) {
                pollutant.max = max(avg.data.for.graph$obs,avg.data.for.graph$mod,na.rm=TRUE)
                pollutant.min = min(avg.data.for.graph$obs,avg.data.for.graph$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(" (",pol.units,")",sep="")
              graph.z=NA
              scatter.cex=NULL
              scatter.mod.line=TRUE
              
              if (bFreqScatterPlot) {
                # check the number of counts in bins in case too few or need to change colour scale
                if (bScatterLog) {
                  scatter.dummy.data = scatter.plot.data
                  scatter.dummy.data$obs = log10(scatter.dummy.data$obs)
                  scatter.dummy.data$mod = log10(scatter.dummy.data$mod)            
                } else {
                  scatter.dummy.data = scatter.plot.data
                }
                # determine range of data and no. hexbins to modify colour key if necessary
                # determine the maximum count (in bins) i.e. max label of colour scale
                if (!bNoFilter) {
                  if (bFilterWeekday) {
                    max.count = rep(NA, 7)
                    unique.matches = c("Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday")
                  } else if (bFilterMonth) {
                    max.count = rep(NA, 12)
                    unique.matches = c("January","February","March","April","May","June","July","August","September","October","November","December")
                  } else if (bFilterSeason) {
                    max.count = rep(NA, 4)
                    unique.matches = c("December","January","February","March","April","May","June","July","August","September","October","November")
                  }
                  for (i in 1:length(max.count)) {
                    if (bFilterWeekday) {
                      filter.subset = subset(scatter.dummy.data,weekdays(scatter.dummy.data$date)==unique.matches[i])
                    } else if (bFilterMonth) {
                      filter.subset = subset(scatter.dummy.data,months(scatter.dummy.data$date)==unique.matches[i])
                    } else if (bFilterSeason) {
                      filter.subset = subset(scatter.dummy.data, subset = months(scatter.dummy.data$date) %in% unique.matches[c((((i-1)*3)+1):(((i-1)*3)+3))] )
                    }
                    max.count[i] = max(hexbin(filter.subset$obs, filter.subset$mod, xbnds=range(filter.subset[,c("obs","mod")]), ybnds=range(filter.subset[,c("obs","mod")]))@count)
                    rm(filter.subset)
                  }
                  hexbin.count = max(max.count, na.rm=TRUE)
                } else { #bNoFilter
                  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)
                # determine if linear or non-linear colour scale should be used
                # hexbin.count greater than 17 will use a log10 colour scale, otherwise, apply no transformation
                if (hexbin.count <= 17) {
                  scatter.trans = "identity"
                } else {
                  scatter.trans = "log10"
                }
                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.method="hexbin"   
                if (bScatterLog) {
                  scatter.xbnds=log10(range(scatter.plot.data[,c("obs","mod")]))
                } else {
                  scatter.xbnds=range(scatter.plot.data[,c("obs","mod")])            
                }
                scatter.titletext="Frequency Scatter"
              } else if (bScatterPlot) {
                scatter.cols="black"
                scatter.pch=16 # filled dots
                scatter.method="scatter"  
                scatter.xbnds=NULL
                scatter.trans=FALSE 
                scatter.titletext="Scatter"
              }
            } else if (bForecastIndex) {
              pol.display = paste(" (Forecast Indices)",sep="")
              if (bPlotAlert) {
                alert.thresholds = read.csv(alert.thresholds.file,header=TRUE,sep=",",colClasses="character")
                names(alert.thresholds)=tolower(names(alert.thresholds))
                expected.alert.col.names = c("index","name")
                for (i in 1:length(expected.alert.col.names)) {
                  if (length(which(names(alert.thresholds)==expected.alert.col.names[i]))==0) {
                    stop(paste("Missing column in alert thresholds file: ",expected.alert.col.names[i]))
                    quit("no")
                  }
                } 
                if (any(is.na(as.numeric(alert.thresholds$index)))) {
                  stop("Invalid reference value(s) provided.")
                  quit("no")
                }
                ref.selected.x = list(v=as.numeric(alert.thresholds$index),lty=3)
                ref.selected.y = list(h=as.numeric(alert.thresholds$index),lty=3)
              } else {
                ref.selected.x = NULL
                ref.selected.y = NULL
              }
              graph.obs = "index.obs"
              graph.mod="index.mod"
              max.val.x = 1.05*max(as.numeric(data.for.index.graph$index.obs),na.rm=TRUE)
              max.val.y = 1.05*max(as.numeric(data.for.index.graph$index.mod),na.rm=TRUE)
              max.val=max(max.val.x,max.val.y)
              min.val=0
              if (bScatterPlot | bFreqScatterPlot){
                if (any(scatter.plot.data[,"count"] > 1)) {
                  graph.z="count"   
                  scatter.pch=19
                } else {
                  #Don't have counts so output conventional scatter plot
                  graph.z=NA
                  scatter.pch=16
                }
                scatter.cols="increment"
                if (bNoFilter) {
                  scatter.cex=3
                  filter.type="default"
                } else {
                  scatter.cex=1
                }
                scatter.mod.line=FALSE
                scatter.method="scatter"
                scatter.xbnds=NULL
                scatter.trans=FALSE 
                scatter.titletext="Frequency Scatter"
              }
            }
            
            #Produce time variation graph
            if (bTimeVariation) {
              cat("      Producing time variation graph...\n")
              #Remove data points where either modelled or observed data are missing
              time.variation.data = subset(time.variation.data,!is.na(mod) & !is.na(obs))
              if (bJPG) {
                jpeg(paste(time.var.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimeVariation.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
              } else if (bPNG) {
                png(paste(time.var.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimeVariation.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
              } else if (bPDF) {
                pdf(paste(time.var.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimeVariation.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
              }
              #Construct titles (max 2 rows)
              titletext1a = paste("Time Variation Analysis: ",diag.name,", ",single.model,sep="")
              titletext1b = paste("Station ",single.station.output.name,", Pollutant ",single.pollutant.upper," (",pol.units,")",sep="")
              titletext2 =  paste("Date range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep="")
              subtext=""
              if (nchar(titletext1a)+nchar(titletext1b)+nchar(titletext2)+4<max.title.char.length) {
                #Fits on 1 row
                titletext = paste(titletext1a,titletext1b,titletext2,sep=", ")
              } else {
                if (nchar(titletext1a)+nchar(titletext1b)+2<max.title.char.length) {
                  #Fits on 2 rows
                  titletext1 = paste(titletext1a,titletext1b,sep=", ")
                  titletext = paste(titletext1,titletext2,sep="\n")
                } else {
                  #Needs 3 rows, so date range must go in subtitle
                  titletext = paste(titletext1a,titletext1b,sep="\n")
                  subtext = titletext2
                }
              }
              tv.object = timeVariation(
                time.variation.data,
                pollutant=c("mod","obs"),
                name.pol=c("Modelled","Observed"),
                ylab=paste(single.pollutant.upper," (",pol.units,")",sep=""),
                title=titletext,
                caption=subtext,
                plot=FALSE)
              plot(plot(tv.object) + patchwork::plot_annotation(theme = ggplot2::theme(plot.title = ggplot2::element_text(hjust=0.5),
                                                                                       plot.caption = ggplot2::element_text(hjust=0.5))))
              graphics.off()
            }
            
            #Produce scatter plot
            if (bFreqScatterPlot | bScatterPlot) {
              cat("      Producing scatter plot...\n")
            }
            if (!bForecastIndex & bFreqScatterPlot) {
              # check frequency scatter plot has .gt. 3 count in each bin
              if (hexbin.count <= 3) {
                cat(paste("        Warning: Too few data for frequency scatter plot output for pollutant - choose conventional scatter plot instead\n"))
                bFreqScatterPlot = FALSE
              }
            }
            if (bFreqScatterPlot | bScatterPlot) {
              
              if (bJPG) {
                jpeg(paste(scatter.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_ScatterPlot.jpg",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
              } else if (bPNG) {
                png(paste(scatter.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_ScatterPlot.png",sep=""),width=paper.width,height=paper.width,units="mm",res=300)
              }  else if (bPDF) {
                pdf(paste(scatter.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_ScatterPlot.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title) 
              }
              if(!is.null(ref.selected.y)){
                if (any(as.numeric(ref.selected.y$h) > max.val)) {
                  cat("        Warning - one or more reference values lie outside of the data limits and will not be plotted on the scatter plot.\n")
                }
              }
              #Construct titles (max 2 rows)
              titletext1a = paste(scatter.titletext," Plot: ",diag.name,", ",single.model,sep="")
              titletext1b = single.station.output.name 
              titletext2 = paste(substring(avg.name,3)," ",single.pollutant.upper,pol.display,filter.name,sep="")
              if (nchar(titletext1a)+nchar(titletext1b)+2<max.title.char.length.scatter) {
                #Can fit station name on first row
                titletext = paste(titletext1a,", ",titletext1b," \n",titletext2,sep="")
              } else {
                #Put station name on second row
                titletext = paste(titletext1a,"\n",titletext1b,", ",titletext2,sep="")
              }
              if (length(scatter.plot.data[[graph.obs]]) > 0){
                plot_object = suppressMessages(scatterPlot(
                  scatter.plot.data,
                  x=graph.obs,y=graph.mod,z=graph.z,
                  xlab="Observed",ylab="Modelled",
                  log.x=bScatterLog, log.y=bScatterLog,
                  cols=scatter.cols,pch=scatter.pch,cex=scatter.cex,
                  type=filter.type,method=scatter.method,
                  mod.line=scatter.mod.line,ref.x=ref.selected.x,ref.y=ref.selected.y,
                  main=titletext,
                  sub=paste("Date range: ", data.start.date.for.graphs, " to ", data.end.date.for.graphs, sep = ""),
                  xbnds=scatter.xbnds, ybnds=scatter.xbnds,
                  trans=scatter.trans,
                  linear=bScatterBestLine,
                  plot=FALSE           # return ggplot object, don't auto-plot
                )$plot +                       # extract the ggplot object and add to it
                  coord_fixed(ratio = 1))
                
                # Use the custom axes limits, if defined
                if (bScatterLimits) {
                  min.val = scatter.xymin
                  max.val = scatter.xymax
                  
                  # Compute data min and max for comparison
                  data.vals = c(scatter.plot.data[[graph.obs]], scatter.plot.data[[graph.mod]])
                  data.min = min(data.vals, na.rm = TRUE)
                  data.max = max(data.vals, na.rm = TRUE)
                  
                  if (min.val > data.min) {
                    cat(paste0("        Warning: The minimum axes limit is greater than the minimum value in the data. Some data will not be plotted\n          Min axes limit: ", min.val, "\n          Min data value: ", data.min, "\n"))
                  }
                  if (max.val < data.max) {
                    cat(paste0("        Warning: The maximum axes limit is less than the maximum value in the data. Some data will not be plotted\n          Max axes limit: ", max.val, "\n          Max data value: ", data.max, "\n"))
                  }
                  
                }
                
                # Set the axes based on the value of bScatterLog
                if (bScatterLog) {
                  plot_object = suppressMessages(plot_object + 
                                                   scale_x_log10(limits = c(min.val, max.val)) + 
                                                   scale_y_log10(limits = c(min.val, max.val)))
                } else {
                  plot_object = suppressMessages(plot_object + 
                                                   scale_x_continuous(limits = c(min.val, max.val)) + 
                                                   scale_y_continuous(limits = c(min.val, max.val)))
                }
                
                suppressMessages(plot(plot_object))              # render the plot
              }else{
                cat("        Insufficient data to create scatter plot - skipping this station. \n")
              }
              graphics.off()
            }
            
            #Produce time plot
            if (bTimePlot) {   
              cat("      Producing time plot...\n")
              if (bJPG) {
                jpeg(paste(time.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimePlot.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
              } else if (bPNG) {
                png(paste(time.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimePlot.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
              } else if (bPDF) {   
                pdf(paste(time.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimePlot.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
              }
              if(!is.null(ref.selected.y)){
                if (any(as.numeric(ref.selected.y$h) > max.val)) {
                  cat("        Warning - one or more reference values lie outside of the data limits and will not be plotted on the time plot.\n")
                }
              }
              if (bPaperSmall) {
                titletext = paste(paste("Time Plot: ",diag.name,", ",single.model,", ",single.station.output.name,sep=""),paste(substring(avg.name,3)," ",single.pollutant.upper,pol.display,sep=""), sep="\n")
              } else {
                titletext = paste("Time Plot: ",diag.name,", ",single.model,", Station ",single.station.output.name,", Pollutant ",single.pollutant.upper,pol.display,avg.name,sep="")
                if (nchar(titletext)>65) {
                  titletext = paste(paste("Time Plot: ",diag.name,", ",single.model,", ",single.station.output.name,sep=""),paste(substring(avg.name,3)," ",single.pollutant.upper,pol.display,sep=""), sep="\n")
                }
              }
              if (bTimePlotMarkers) {
                pch.value = c(16,16)
              } else {
                pch.value = c(NA,NA)
              }
              if (bTimePlotLimits) { # Use custom y-axis limits if defined.
                ylim.value = c(timeplot.ymin, timeplot.ymax)
                
                # Extract values from data to plot
                y.vals = as.numeric(unlist(time.plot.data[, c(graph.mod, graph.obs)]))
                
                data.min = min(y.vals, na.rm = TRUE)
                data.max = max(y.vals, na.rm = TRUE)
                
                # Check lower limit
                if (data.min < ylim.value[1]) {
                  cat(paste0("        Warning: The minimum y-axis limit is greater than the minimum value in the data. Some data will not be plotted\n          Min axis limit: ", ylim.value[1], "\n          Min data value: ", data.min, "\n"))
                }
                
                # Check upper limit
                if (data.max > ylim.value[2]) {
                  cat(paste0("        Warning: The maximum y-axis limit is greater than the minimum value in the data. Some data will not be plotted\n          Max axis limit: ", ylim.value[2], "\n          Max data value: ", data.max, "\n"))
                }
              } else {
                ylim.value = c(0,max.val)
              }
              plot_object = timePlot(
                time.plot.data,
                pollutant=c(graph.mod,graph.obs),
                name.pol=c("Modelled","Observed"),
                group=TRUE,ylab=paste(single.pollutant.upper,pol.display,sep=""),
                cols=c("red","blue"),ref.y=ref.selected.y,lty=c(1,1),ylim=ylim.value,
                pch=pch.value,
                main=titletext,
                sub=paste("Date Range: ",data.start.date.for.graphs, " to ",data.end.date.for.graphs,sep=""),
                plot=FALSE
              )$plot
              plot(plot_object)
              graphics.off()
            }
            
            #Produce polar plots
            if (bPolarPlots) {
              names(polar.plot.data)[names(polar.plot.data)=="mod"]="Modelled"
              names(polar.plot.data)[names(polar.plot.data)=="obs"]="Observed"
              if (bPolar.filter.by.stability) {
                filter.type="h_lmo"
              } else if (bPolar.filter.by.season) {
                filter.type="season"
              } else if (bPolar.filter.by.weekday) {
                filter.type="weekday"
              } else if (bPolar.filter.by.month) {
                filter.type="month"
              } else {
                filter.type="default"
              }  
              if (bPolar.conc.plot) {
                cat("      Producing concentration polar plot...\n")
                if (bJPG) {
                  jpeg(paste(polar.conc.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PolarConcPlot.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
                } else if (bPNG) {
                  png(paste(polar.conc.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PolarConcPlot.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
                } else if (bPDF) {   
                  pdf(paste(polar.conc.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PolarConcPlot.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
                }
                #Construct titles (max 2 rows)
                titletext1a = paste("Polar concentration plot: ",diag.name,", ",single.model,sep="")
                titletext1b = paste("Station ",single.station.output.name,", Pollutant ",single.pollutant.upper," (",pol.units,")",sep="")
                titletext2 = paste("Date range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep="")
                subtext.date=""
                if (nchar(titletext1a)+nchar(titletext1b)+nchar(titletext2)+4<max.title.char.length) {
                  #Fits on 1 row
                  titletext = paste(titletext1a,titletext1b,titletext2,sep=", ")
                } else {
                  if (nchar(titletext1a)+nchar(titletext1b)+2<max.title.char.length) {
                    #Fits on 2 rows
                    titletext1 = paste(titletext1a,titletext1b,sep=", ")
                    titletext = paste(titletext1,titletext2,sep="\n")
                  } else {
                    #Needs 3 rows, so date range must go in subtitle
                    titletext = paste(titletext1a,titletext1b,sep="\n")
                    subtext.date = titletext2
                  }
                }
                if (bPolar.radial.axis.windspeed) {
                  radial.axis="ws"
                  axis.text="wind speed (m/s)"
                } else if (bPolar.radial.axis.hour) {
                  polar.plot.data$hour=as.numeric(format(polar.plot.data$date,format="%H"))
                  radial.axis="hour"
                  axis.text="hour"
                } else if (bPolar.radial.axis.month) {
                  polar.plot.data$month=as.numeric(format(polar.plot.data$date,format="%m"))
                  radial.axis="month"
                  axis.text="month"
                }
                subtext = paste("Minimum number of data points per bin = ",Polar.conc.plot.min.bin,"; Radial axis = ",axis.text,sep="")
                if (nchar(subtext.date)>0) {
                  subtext = paste(subtext,subtext.date,sep="\n")
                }
                polar.plot = tryCatch(polarPlot(
                  polar.plot.data,
                  x=radial.axis,
                  pollutant=c("Modelled","Observed"),
                  type=filter.type,
                  main=titletext,
                  key.header="",key.footer="",
                  min.bin=Polar.conc.plot.min.bin,
                  units=axis.text,
                  sub=subtext
                ),error=function(e) return(NA))
                if(class(polar.plot)!="openair") {
                  cat(paste("        Error plotting polar concentration plot, likely to be caused by insufficient data. Try reducing the minimum allowed number of data points per bin.\n"))
                }
                graphics.off()
              }
              if (bPollution.rose.plot) {
                cat("      Producing pollution rose plot...\n")
                if (bJPG) {
                  jpeg(paste(pollution.rose.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PollutionRose.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
                } else if (bPNG) {
                  png(paste(pollution.rose.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PollutionRose.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
                } else if (bPDF) {   
                  pdf(paste(pollution.rose.plot.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PollutionRose.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
                }
                #Construct titles (max 2 rows)
                titletext1 = paste("Pollution rose plot: ",diag.name,", ",single.model,sep="")
                titletext2 = paste("Station ",single.station.output.name,", Pollutant ",single.pollutant.upper," (",pol.units,")",sep="")
                subtext.date=paste("Date range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep="")
                if (nchar(titletext1)+nchar(titletext2)+2<max.title.char.length) {
                  #Fits on 1 row
                  titletext = paste(titletext1,titletext2,sep=", ")
                } else {
                  #Fits on 2 rows
                  titletext = paste(titletext1,titletext2,sep="\n")
                }
                polar.plot.data.melt=melt(polar.plot.data,measure.vars=c("Modelled","Observed"),variable.name="Quantity")
                if (filter.type=="default") {polrose.filter.type="Quantity"} else {polrose.filter.type=c(filter.type,"Quantity")}
                pollution.rose = tryCatch(pollutionRose(
                  polar.plot.data.melt,
                  pollutant="value",
                  type=polrose.filter.type,
                  main=titletext,
                  key.footer="", key.header="",
                  annotate=FALSE,
                  sub=subtext.date
                ),error=function(e) return(NA))
                if(class(pollution.rose)!="openair") {
                  cat(paste("        Error plotting pollution rose plot, likely to be caused by insufficient data.\n"))
                }
                graphics.off()
              }  
              if (bPolar.wind.freq.plot) {
                cat("      Producing polar wind speed frequency plot...\n")
                if (bJPG) {
                  jpeg(paste(polar.wind.freq.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PolarWindFreq.jpg",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
                } else if (bPNG) {
                  png(paste(polar.wind.freq.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PolarWindFreq.png",sep=""),width=paper.height,height=paper.width,units="mm",res=300)
                } else if (bPDF) {   
                  pdf(paste(polar.wind.freq.dir,"/",diag.file.name,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_PolarWindFreq.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
                }
                titletext = paste(
                  "Polar wind speed frequency plot: ",diag.name,", ",
                  data.start.date.for.graphs," to ",data.end.date.for.graphs,sep="")
                if (nchar(titletext) > max.title.char.length) {
                  titletext = paste(
                    paste("Polar wind speed frequency plot:",diag.name),
                    paste("Date range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep=""),
                    sep="\n")
                }
                polar.freq = tryCatch(polarFreq(
                  polar.plot.data,
                  type=filter.type,
                  main=titletext,
                  border.col="white",
                  min.bin=Polar.wind.plot.min.bin,
                  key.footer="", key.header="",
                  sub=paste("Minimum number of data points per bin = ",Polar.wind.plot.min.bin,sep="")
                ),error=function(e) return(NA))
                if(class(polar.freq)!="openair") {
                  cat(paste("        Error plotting polar wind speed frequency plot, likely to be caused by insufficient data. Try reducing the minimum allowed number of data points per bin.\n"))
                }
                graphics.off()
                }
            }
          } #end of test for data to plot for this station, for this pollutant, for this model
        } #end of loop over stations 
      } #end of test for data to plot for this pollutant, for this model
    } #end of loop over pollutants
  } #end of test for data to plot for this model
} #end of loop over models

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

