#Toolkit version number
r.tool.name = "Model Diagnostics"
r.diag.toolkit.version = "4.0"
tkd.version = "tkdfileversion1"
svn.version.diag = "$ProjectRevision: SVN 1942 (2016-01-29) $"
svn.prj.name.diag = "$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)

#Set up input and output
diag.file = commandArgs(trailingOnly=TRUE)[[1]]
diag.output.root = sub(".tkd","",diag.file,ignore.case=TRUE)
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="")) 

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

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

#Check that the tkd file is for the right version
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")
}

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

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

#R SETTINGS
mirror=tkd.data[nRSettings+2]
update.packages = tkd.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
single.station = tkd.data[nStation.Pollutant.Model+2]
if (single.station=="") {
   stop("No station name entered.")
   quit("no")
}
single.station.upper = toupper(single.station)
single.station.lower = tolower(single.station)

single.pollutant = tkd.data[nStation.Pollutant.Model+3]
if (single.pollutant=="") {
   stop("No pollutant name entered.")
   quit("no")
}
single.pollutant.upper = toupper(single.pollutant)
single.pollutant.lower = tolower(single.pollutant)

single.model = tkd.data[nStation.Pollutant.Model+4]
if (single.model=="") {
   stop("No modelled data label entered.")
   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]
if (ref.values!="") {
  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")
bScatterPlot = tkd.data[nPlot.Options+4]=="Y" & tkd.data[nPlot.Options+5]=="1"
bFreqScatterPlot = tkd.data[nPlot.Options+4]=="Y" & tkd.data[nPlot.Options+5]=="0"
scatter.filters = c("No filtering","Weekday","Month","Season")
scatter.filter = scatter.filters[as.numeric(tkd.data[nPlot.Options+6])+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+7]=="Y")

if (!bTimeVariation & !bScatterPlot & !bFreqScatterPlot & !bTimePlot) {
   stop("No diagnostic graph output selected.")
   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]



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

#R CODE:

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

#Get the data for the selected station, pollutant and model
if (all(names(all.data)!=single.model)) {
   modelled.data.labels = paste(names(all.data)[names(all.data)!="date" & names(all.data)!="obs" & names(all.data)!="station" & names(all.data)!="pollutant" & names(all.data)!="type"], collapse=", ")
   stop(paste("No data for selected modelled data ",single.model,".  This workspace contains modelled data labelled: ",modelled.data.labels,"\n",sep=""))
   quit("no")
}
data.for.plotting = subset(all.data,pollutant==single.pollutant.lower & station==single.station.lower, select=c("date","obs","station","pollutant","type",single.model))
#Check there are some valid (not NA) modelled data values
if (length(which(!is.na(data.for.plotting[,single.model])))==0) {
   data.station = subset(all.data,station==single.station.lower)
   if (length(which(!is.na(data.station[,single.model]))) > 0) {
      pols = toupper(paste(unique(data.station[!is.na(data.station[,single.model]),"pollutant"]),collapse=" "))
      stop(paste("No data for selected station ",single.station.upper," for selected pollutant ",single.pollutant.upper," in the ",single.model," data. This station has data for these pollutants: ",pols,"\n",sep=""))
   } else {
      stop(paste("No data for selected station",single.station.upper,"in the",single.model,"data."))
   }
   quit("no")
}
#Change column headings for openair functions
names(data.for.plotting)[names(data.for.plotting)==single.model] = "mod"
pol.units = pollutants$output.units[pollutants$pollutant==single.pollutant.lower]   

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

#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

} 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$mod.statistic=="mean" & pollutant.data$mod.avg.time.hours=="1")) {
      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)
      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)
      } 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)
   }
   
   #Set data for plotting for each graph
   time.variation.data = data.for.plotting
   scatter.plot.data = avg.data.for.graph
   #remove NAs
   scatter.plot.data = subset(scatter.plot.data, !(is.na(obs) | is.na(mod)) )
   time.plot.data = avg.data.for.graph

} else if (bForecastIndex) {

   bScatterLog = FALSE
   
   avg.name = ""
   #Call functions to create data with indices applied
   pollutants = process.index.scales(index.scales.file, pollutants)
   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"
   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),]
   #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
   time.variation.data = data.for.plotting   
   time.plot.data = data.for.index.graph
}
#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
} 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 (!bForecastIndex) {
   ref.selected.x = list(v=ref.values,lty=3)
   ref.selected.y = list(h=ref.values,lty=3)
   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
      # >32 count in each bin for non-linear key of 16 colours
      # <=32 counts, use non-linear key and reduce number of colours
      # <=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  
      }
      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.inv=FALSE
      scatter.colorcut=NULL 
      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 (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.inv=FALSE
   scatter.colorcut=NULL  
   scatter.titletext="Frequency Scatter"
}

#Produce time variation graph
if (bTimeVariation) {
   cat("Producing time variation graph...\n")
   if (bJPG) {
      jpeg(paste(diag.output.root,"_",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(diag.output.root,"_",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(diag.output.root,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimeVariation.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
   }
   titletext = paste("Time Variation Analysis: ",diag.name,", ",single.model,", Station ",single.station.upper,", Pollutant ",single.pollutant.upper," (",pol.units,"), ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep="")
   if ((nchar(titletext) >90 & bPaperLarge) | (nchar(titletext) >75 & bPaperMedium) | (nchar(titletext) >65 & bPaperSmall)) {
      titletext = paste(paste("Time variation analysis: ",diag.name,", ",single.model,", Station ",single.station.upper,", Pollutant ",single.pollutant.upper," (",pol.units,")",sep=""),paste(data.start.date.for.graphs," to ",data.end.date.for.graphs,sep=""),sep="\n")
   }
   timeVariation(
      time.variation.data,
      pollutant=c("mod","obs"),
      name.pol=c("Modelled","Observed"),
      ylab=paste(single.pollutant.upper," (",pol.units,")",sep=""),
      main=titletext)
   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(diag.output.root,"_",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(diag.output.root,"_",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(diag.output.root,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_ScatterPlot.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title) 
   }
   if (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")
      }
   }
   if (bPaperSmall) {
      titletext = paste(paste(scatter.titletext," Plot: ",diag.name,", ",single.model,", ",single.station.upper,sep=""),paste(substring(avg.name,3)," ",single.pollutant.upper,pol.display,filter.name,sep=""), sep="\n")
   } else {
      titletext = paste(scatter.titletext," Plot: ",diag.name,", ",single.model,", Station ",single.station.upper,", Pollutant ",single.pollutant.upper,pol.display,avg.name,filter.name,sep="")
      if (nchar(titletext)>65) {
         titletext = paste(paste(scatter.titletext," Plot: ",diag.name,", ",single.model,", ",single.station.upper,sep=""),paste(substring(avg.name,3)," ",single.pollutant.upper,pol.display,filter.name,sep=""), sep="\n")
      }
   }
   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),
      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=""),
      aspect="iso", 
      xbnds=scatter.xbnds, ybnds=scatter.xbnds,
      trans=scatter.trans, inv=scatter.inv, 
      colorcut=scatter.colorcut
   )
   graphics.off()
}
#Produce time plot
if (bTimePlot) {   
   cat("Producing time plot...\n")
   if (bJPG) {
      jpeg(paste(diag.output.root,"_",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(diag.output.root,"_",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(diag.output.root,"_",single.model,"_",single.pollutant.upper,"_",single.station.upper,"_TimePlot.pdf",sep=""),paper="a4r",width=0,height=0,title=pdf.title)
   }
   if (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.upper,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.upper,", Pollutant ",single.pollutant.upper,pol.display,avg.name,sep="")
     if (nchar(titletext)>65) {
       titletext = paste(paste("Time Plot: ",diag.name,", ",single.model,", ",single.station.upper,sep=""),paste(substring(avg.name,3)," ",single.pollutant.upper,pol.display,sep=""), sep="\n")
     }
   }
   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=c(0,max.val),
      main=titletext,
      sub=paste("Date Range: ",data.start.date.for.graphs," to ",data.end.date.for.graphs,sep=""))   
   graphics.off()
}

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