#Toolkit version number
r.tool.name = "Data Input"
r.data.toolkit.version = "4.0"
tki.version = "tkifileversion1"
svn.version.data = "$ProjectRevision: SVN 1942 (2016-01-29) $"
svn.prj.name.data = "$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)
#Flag to enable extra output for development testing; should always be false in checked in versions
btesting = FALSE

#Set up input and output
settings.data.file = commandArgs(trailingOnly=TRUE)[[1]]
output.file.path.root = sub(".tki","",settings.data.file,ignore.case=TRUE)
log.file = file(paste(output.file.path.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.data.toolkit.version,": ",r.tool.name," Tool \n",sep=""))
if (length(grep("$ProjectRevision: ",svn.version.data,fixed=TRUE))>0) {
  svn.version.data.numberonly=sub("$ProjectRevision: ","",svn.version.data,fixed=TRUE)
  svn.version.data.numberonly=sub(" $","",svn.version.data.numberonly,fixed=TRUE)
  cat(paste("Internal version ",svn.version.data.numberonly," \n",sep=""))
} else {
  cat("CERC INTERNAL USE ONLY\n")
}
cat(paste("Running in ",R.version$version.string," \n",sep=""))
cat(paste(date(),"\n\n",sep=""))

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

tki.data = readLines(settings.data.file)
tki.data = sub("PATH = \"","",as.character(tki.data),fixed=TRUE)
tki.data = sub("\"","",as.character(tki.data),fixed=TRUE)
#Check that the settings file is in the right format
tki.version.in = as.character(tki.data[1])
if (tki.version.in != tki.version) {
  stop(paste("tki file version (",tki.version.in,") does not match the expected version (",tki.version,")", sep=""))
  quit("no")   
}   

#SETTINGS
nSetup = which(tki.data=="Settings")
bNewData=tki.data[nSetup+2]=="0"
bAddData=tki.data[nSetup+2]=="1"
existing.workspace.file =tki.data[nSetup+3]
pollutants.file = tki.data[nSetup+4] 
stations.file = tki.data[nSetup+5]
boutput.stats.as.csv = tki.data[nSetup+6]=="Y"
if (bAddData & existing.workspace.file=="") {
  stop("No workspace selected to add new model data set to.")
  quit("no")
}
if (bNewData) {
   if (pollutants.file=="") {
      stop("No pollutant definitions file selected.")
      quit("no")
   }
   if (stations.file=="") {
     stop("No station definitions file selected.")
     quit("no")
   }
}

#R SETTINGS
nRSettings = which(tki.data=="RSettings")
mirror=tki.data[nRSettings+2]
update.packages = tki.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")
}


#MODELLED DATA
netCDF.formats = c("AIRSHEDS","CAMS Ensemble","CMAQ","CAMS Ensemble CERC archive","Prevair CERC archive")
nMod = which(tki.data=="Concentration data")
import.model.name = tki.data[nMod+2]
pollutants.mod.file = tki.data[nMod+3]
mod.format=tki.data[nMod+4]
bfile.mod=tki.data[nMod+5]=="file"
bdir.mod=tki.data[nMod+5]=="dir"
file.path.mod=tki.data[nMod+6]
dir.path.mod=tki.data[nMod+7]
data.separator.mod=tki.data[nMod+11]
missing.data.mod=tki.data[nMod+12]
#NETCDF option
bmod.ncdf = any(mod.format==netCDF.formats)
mod.ncdf.format = mod.format
bmod.ncdf.select.file = bfile.mod
mod.ncdf.file = file.path.mod
bmod.ncdf.select.dir = bdir.mod
mod.ncdf.dir = dir.path.mod
start.date.archive=tki.data[nMod+8]
end.date.archive=tki.data[nMod+9]
forecast.day.archive=tki.data[nMod+10]
mod.ncdf.missing.data = missing.data.mod
#ADMS option
bmod.adms.pst = mod.format=="ADMS PST"
bmod.pst.select.file =  bfile.mod
mod.pst.file = file.path.mod
bmod.pst.select.dir = bdir.mod
mod.pst.dir = dir.path.mod
mod.pst.missing.data = "-999"
#CSV option
bmod.csv = mod.format=="CSV"
bmod.csv.select.file =  bfile.mod
mod.csv.file = file.path.mod
bmod.csv.select.dir = bdir.mod
mod.csv.dir = dir.path.mod
mod.csv.separator = data.separator.mod
mod.csv.missing.data = missing.data.mod
#Check all required information has been entered
if (import.model.name=="") {
   stop("No model name given.")
   quit("no") 
}
if (grepl("[^[:alnum:]]",import.model.name) & !(grepl("[.]",import.model.name) | grepl("[-]",import.model.name) | grepl("[_]",import.model.name) | grepl("[()]",import.model.name) )) {
   stop("Modelled data label must be alpha-numeric.")
   quit("no") 
}
if (pollutants.mod.file=="") {
  stop("No modelled data pollutants file selected.")
  quit("no")
}
if (bfile.mod & file.path.mod=="") {
  stop("No modelled data file selected.")
  quit("no")
}
if (bdir.mod & dir.path.mod=="") {
  stop("No modelled data directory selected.")
  quit("no")
}
if (bmod.ncdf) {
   if (bmod.ncdf.select.dir & mod.ncdf.format == "DELTA") {
      stop("The option to select a whole directory of netCDF files is not permitted for DELTA format files, which contain no date information.")
      quit("no")
   }
}

#OBSERVED DATA
nObs=which(tki.data=="Observed data")
benter.obs.data = tki.data[nObs+1]=="Y"
pollutants.obs.file = tki.data[nObs+2]
obs.format=tki.data[nObs+3]
bfile.obs=tki.data[nObs+4]=="file"
bdir.obs=tki.data[nObs+4]=="dir"
file.path.obs=tki.data[nObs+5]
dir.path.obs=tki.data[nObs+6]
data.separator.obs=tki.data[nObs+7]
missing.data.obs=tki.data[nObs+8]
#CSV option
bobs.csv = obs.format=="CSV"
bobs.csv.select.file = bfile.obs
obs.csv.file = file.path.obs
bobs.csv.select.dir =  bdir.obs
obs.csv.dir = dir.path.obs
obs.csv.separator = data.separator.obs
obs.csv.missing.data = missing.data.obs
#KCL option
bobs.kcl = obs.format=="KCL"
#AURN option
bobs.aurn = obs.format=="AURN"
#Error handling:
if (bNewData) {
   if (!benter.obs.data) {
      stop("Observed data must be entered when a new workspace is being created.")
      quit("no")
   }
   if (pollutants.obs.file=="") {
      stop("No observed data pollutants file selected.")
      quit("no")
   }
   if ((bmod.ncdf & mod.ncdf.format=="DELTA") & !bobs.csv) {
      stop("When importing DELTA-format modelled data, the observed data must also be of DELTA-format.")
      quit("no")
   }
   if (bobs.csv) {
      if (bobs.csv.select.file & obs.csv.file=="") {
         stop("No observed CSV data file selected.")
         quit("no")
      }
      if (bobs.csv.select.dir & obs.csv.dir=="") {
         stop("No observed CSV data directory selected.")
         quit("no")
      }
   }
} else if (bAddData) {
   if (benter.obs.data) {
      stop("No observed data can be entered when modelled data is being added to an existing workspace.")
      quit("no")
   }
}

#FUNCTIONS:
# only need to define if new workspace
if (bNewData) {
   #Explicitly define the modulus function, fmod
   fmod = function (k,m)
   {
       j = floor(k/m)
       a = k-(m*j)
       return(a) 
   }
   
   #Modification of GEOmap function to handle greenwich meridian
   cerc.lambert.cc.xy = function (phi, lam, PROJ.DATA)
   {
       lam = fmod(lam, 360)
       phi0 = PROJ.DATA$LAT0
       lam0 = PROJ.DATA$LON0
       phi1 = PROJ.DATA$LAT1
       phi2 = PROJ.DATA$LAT2
       FE = PROJ.DATA$FE
       FN = PROJ.DATA$FN
       phi = phi * pi/180
       lam = lam * pi/180
       phi1 = phi1 * pi/180
       phi2 = phi2 * pi/180
       phi0 = phi0 * pi/180
       lam0 = lam0 * pi/180
       R = MAPconstants()$A.MAPK
       if (phi1==phi2) {
          #1 standard parallel
          n = sin (phi1)
       } else {
          #2 standard parallels
          n = log(cos(phi1)/cos(phi2))/log(tan(pi/4 + phi2/2)/tan(pi/4 +
              phi1/2))
       }
       F = cos(phi1) * (tan(pi/4 + phi1/2))^n/n
       rho0 = R * F/(tan(pi/4 + phi0/2))^n
       rho = R * F/((tan(pi/4 + phi/2))^n)
       lam.diff = lam - lam0
       for (i in 1:length(lam)) {
          if (lam.diff[i] > pi) {lam.diff[i]  = lam.diff[i] - 2*pi}
          if (lam.diff[i] < -1*pi) {lam.diff[i]  = lam.diff[i] + 2*pi}
       }
       theta = n * lam.diff
       x = rho * sin(theta) + FE
       y = rho0 - rho * cos(theta) + FN
       return(list(x = x, y = y))
   }
   
   #Conversion of data between different units
   convert.units = function (in.data, in.units, out.units, conv.ugm3.ppb)
   {
       #First convert data to ug/m3
       if (in.units == "ppb") {
          data.ugm3 = in.data / conv.ugm3.ppb
       } else if (in.units == "ppm") {
          data.ugm3 = 1000 * in.data / conv.ugm3.ppb
       } else if (in.units == "mol/mol") {
          data.ugm3 = 1000000000 * in.data / conv.ugm3.ppb
       } else if (in.units == "ng/m3") {
          data.ugm3 = in.data / 1000
       } else if (in.units == "ug/m3") {
          data.ugm3 = in.data
       } else if (in.units == "mg/m3") {
          data.ugm3 = 1000 * in.data
       } else if (in.units == "g/m3") {
          data.ugm3 = 1000000 * in.data
       } else if (in.units == "kg/m3") {
          data.ugm3 = 1000000000 * in.data
       } else {
          stop("Unrecognised units: must be either ppb, ppm, mol/mol, ng/m3, ug/m3, mg/m3, g/m3 or kg/m3")
          quit("no")
       }           
       #Now convert to output units
       if (out.units == "ppb") {
          out.data = data.ugm3 * conv.ugm3.ppb
       } else if (out.units == "ppm") {
          out.data = data.ugm3 * conv.ugm3.ppb / 1000
       } else if (out.units == "mol/mol") {
          out.data = data.ugm3 * conv.ugm3.ppb / 1000000000
       } else if (out.units == "ng/m3") {
          out.data = 1000 * data.ugm3
       } else if (out.units == "ug/m3") {
          out.data = data.ugm3
       } else if (out.units == "mg/m3") {
          out.data = data.ugm3 / 1000
       } else if (out.units == "g/m3") {
          out.data = data.ugm3 / 1000000
       } else if (out.units == "kg/m3") {
          out.data = data.ugm3 / 1000000000
       } else {
          stop("Unrecognised units: must be either ppb, ppm, mol/mol, ng/m3, ug/m3, mg/m3, g/m3 or kg/m3")
          quit("no")
       }
       return(out.data)
   }
   
   #Conversion of time from HHMMSS format to HH, MM, SS
   convert.time = function (HHMMSS,mode)
   {
       HH = floor(HHMMSS/10000)
       MM = floor((HHMMSS-HH*10000)/100)
       SS = HHMMSS - HH*10000 - MM*100
       if (mode=="string") {
          out.time = paste(HH,MM,SS)
       } else if (mode=="seconds") {
          out.time = SS + MM*60 + HH*3600
       } else {
          stop("Unidentified mode in convert.time: must be either string or seconds")
          quit("no")
       }
       return(out.time)
   }
   
   #Strip out station name from filename
   strip.name = function (filename)
   {
       return(sub("[.][^.]*$", "", basename(filename)))
   }
   
   #Open file of csv data and add file name as station if no station column
   open.csv.and.add.station = function (filename,csv.separator)
   {
      file.data = read.csv(filename,header=TRUE,sep=csv.separator,check.names=FALSE)
      #Check for duplicate or invalid column headers
      if (!all(names(file.data) == make.names(names(file.data),unique=TRUE))) {
         stop(paste("Invalid or duplicate column headers in CSV file ", filename))
         quit("no")
      }
      if (length(which(tolower(names(file.data))=="station"))==0) {
         #There is no 'station' column, so set the station name to be the filename
         file.data$station = strip.name(filename)
      }
      return(file.data)
   }
   
   #Check whether two ranges intersect; returns FALSE if they don't.
   branges.intersect = function(r1start, r1end, r2start, r2end)
   {
       start = max(r1start,r2start) # max of starts 
       end = min(r1end,r2end) # min of ends 
       return (start <= end) 
   }
   
   #Load in a source file from the local directory
   source_local = function(fname){
     argv = commandArgs(trailingOnly = FALSE)
     base_dir = dirname(substring(argv[grep("--file=", argv)], 8))
     source(paste(base_dir, fname, sep=.Platform$file.sep))
   }
   
   #Update, install and load all required R packages
   update.install.load.packages = function(packages.required,mirror,update.packages) {
     
     cat("Setting the CRAN mirror...\n")
     r = getOption("repos")
     r["CRAN"] = mirror
     options(repos=r)
     if (update.packages) {
       cat("Updating installed R packages...\n")
       update.packages(ask=FALSE)
     }
     cat("Loading required R packages...\n")
     is.installed = c(packages.required[1], require(packages.required[1],character.only=TRUE,quietly=TRUE))
     for (i in 2:length(packages.required)) {
       is.installed = rbind(is.installed, c(packages.required[i], require(packages.required[i],character.only=TRUE,quietly=TRUE)))
     }
     #Add column names to the array of packages that have tried to load
     colnames(is.installed) = c("package", "installed?")
     
     #Check if any packages did not load 
     if (any(is.installed[,"installed?"]=="FALSE")) {
       is.not.installed = subset(is.installed, is.installed[,"installed?"]=="FALSE")
       packages.to.install = is.not.installed[1,"package"]
       if (nrow(is.not.installed) < 2) {
         cat(paste("The following package is missing and will be installed: ", packages.to.install,"\n", sep = "\n"))
       } else if (nrow(is.not.installed)>2) {
         for (p in 2:(nrow(is.not.installed)-1)) {
           packages.to.install = paste(packages.to.install,is.not.installed[p,"package"],sep=", ")
         }
         packages.to.install = paste(packages.to.install, is.not.installed[nrow(is.not.installed),"package"],sep=" and ")
         cat(paste("The following packages are missing and will be installed: ", packages.to.install,"\n", sep ="\n"))
       } else {
         packages.to.install = paste(packages.to.install, is.not.installed[nrow(is.not.installed),"package"],sep=" and ")
         cat(paste("The following packages are missing and will be installed: ", packages.to.install,"\n", sep ="\n"))
       }
       cat("Installing required R packages...\n")
       install.packages(is.not.installed[,1])
       cat("Loading newly installed R packages...\n")
       for (i in 1:length(is.not.installed[,1])) {require(is.not.installed[i,1],character.only=TRUE,quietly=TRUE)}
     }
     
   }
   
} # end if (bNewData)

#R CODE:

#If adding model data to existing workspace, load workspace and settings 
if (bAddData) {
   #Firstly store user-defined settings that are to be used in this run through the tool
   r.data.toolkit.version.ADD = r.data.toolkit.version
   output.file.path.root.ADD = output.file.path.root
   bNewData.ADD = bNewData
   bAddData.ADD = bAddData
   mirror.ADD = mirror
   update.packages.ADD = update.packages
   import.model.name.ADD = import.model.name
   pollutants.mod.file.ADD = pollutants.mod.file
   existing.workspace.file.ADD = existing.workspace.file
   bmod.ncdf.ADD = bmod.ncdf 
   mod.ncdf.format.ADD = mod.ncdf.format
   bmod.ncdf.select.file.ADD = bmod.ncdf.select.file
   mod.ncdf.file.ADD = mod.ncdf.file
   bmod.ncdf.select.dir.ADD = bmod.ncdf.select.dir
   mod.ncdf.dir.ADD = mod.ncdf.dir
   mod.ncdf.missing.data.ADD = mod.ncdf.missing.data
   bmod.adms.pst.ADD = bmod.adms.pst
   bmod.pst.select.file.ADD = bmod.pst.select.file
   mod.pst.file.ADD = mod.pst.file
   bmod.pst.select.dir.ADD = bmod.pst.select.dir
   mod.pst.dir.ADD = mod.pst.dir
   mod.pst.missing.data.ADD = mod.pst.missing.data
   bmod.csv.ADD = bmod.csv
   bmod.csv.select.file.ADD = bmod.csv.select.file
   mod.csv.file.ADD = mod.csv.file
   bmod.csv.select.dir.ADD = bmod.csv.select.dir  
   mod.csv.dir.ADD=mod.csv.dir
   mod.csv.separator.ADD=mod.csv.separator
   mod.csv.missing.data.ADD=mod.csv.missing.data
   boutput.stats.as.csv.ADD=boutput.stats.as.csv
   
   #load workspace
   cat("Loading workspace...\n")
   data.loaded=tryCatch(load(existing.workspace.file.ADD),error=function(e) return("error"))
   if (data.loaded=="error") {
      stop("The workspace could not be imported.")
      quit("no")
   }
   
   #check consistency with version of Data Input tool
   if (r.data.toolkit.version.ADD != r.data.toolkit.version) {
      stop(paste("Toolkit version (",r.data.toolkit.version.ADD,") does not match version used to create workspace file (",r.data.toolkit.version,")", sep=""))
      quit("no")
   }

   #check model name is not duplicated
   existing.modelled.data.names = names(all.data)[names(all.data)!="date" & names(all.data)!="obs" & names(all.data)!="station" & names(all.data)!="pollutant" & names(all.data)!="type"]
   if (any(import.model.name.ADD == existing.modelled.data.names)) {
      model.names = paste(existing.modelled.data.names,collapse=", ")
      stop(paste("Modelled data label ",import.model.name.ADD," already exists in workspace. Modelled data labels in existing workspace are: ",model.names, sep=""))
      quit("no")
   }

   # re-set settings from GUI   
   #Firstly store user-defined settings that are to be used in this run through the tool
   r.data.toolkit.version = r.data.toolkit.version.ADD
   output.file.path.root = output.file.path.root.ADD
   bNewData = bNewData.ADD
   bAddData = bAddData.ADD
   mirror = mirror.ADD
   update.packages = update.packages.ADD
   import.model.name = import.model.name.ADD
   pollutants.mod.file = pollutants.mod.file.ADD
   existing.workspace.file = existing.workspace.file.ADD
   bmod.ncdf = bmod.ncdf.ADD 
   mod.ncdf.format = mod.ncdf.format.ADD
   bmod.ncdf.select.file = bmod.ncdf.select.file.ADD
   mod.ncdf.file = mod.ncdf.file.ADD
   bmod.ncdf.select.dir = bmod.ncdf.select.dir.ADD
   mod.ncdf.dir = mod.ncdf.dir.ADD
   mod.ncdf.missing.data = mod.ncdf.missing.data.ADD
   bmod.adms.pst = bmod.adms.pst.ADD
   bmod.pst.select.file = bmod.pst.select.file.ADD
   mod.pst.file = mod.pst.file.ADD
   bmod.pst.select.dir = bmod.pst.select.dir.ADD
   mod.pst.dir = mod.pst.dir.ADD
   mod.pst.missing.data = mod.pst.missing.data.ADD
   bmod.csv = bmod.csv.ADD
   bmod.csv.select.file = bmod.csv.select.file.ADD
   mod.csv.file = mod.csv.file.ADD
   bmod.csv.select.dir = bmod.csv.select.dir.ADD  
   mod.csv.dir=mod.csv.dir.ADD
   mod.csv.separator=mod.csv.separator.ADD
   mod.csv.missing.data=mod.csv.missing.data.ADD
   boutput.stats.as.csv=boutput.stats.as.csv.ADD
   match.vars = grep(".ADD",ls(), ignore.case=FALSE,value=TRUE)
   rm(list=c(match.vars,"match.vars"))
}

allowed.units = c("ppb", "ppm", "mol/mol", "ng/m3", "ug/m3", "mg/m3", "g/m3", "kg/m3")
allowed.statistics = c("max","mean","rolling mean")
if (bNewData) {

   #Set embedded title for PDF files
   pdf.title = paste("Output from DS-LOCAL Model Evaluation Toolkit Version",r.data.toolkit.version)
   
   #READ AND PROCESS POLLUTANT DEFINITION FILE (separator here is ",")
   #header = "pollutant, output.units, conv.ugm3.ppb, min.allowed, max.allowed"
   cat("Importing pollutant details...\n")
   pollutants=read.csv(pollutants.file,header=TRUE,sep=",")
   names(pollutants) = tolower(names(pollutants))
   pol.expected.col.names = c("pollutant","output.units","conv.ugm3.ppb","min.allowed","max.allowed")
   for (j in 1:length(pol.expected.col.names)) {
      if (length(which(names(pollutants)==pol.expected.col.names[j]))==0) {
         stop(paste("Missing column in pollutant definitions file: ",pol.expected.col.names[j]))
         quit("no")
      }
   }
   for (j in 1:length(names(pollutants))) {
      if (length(which(pol.expected.col.names==names(pollutants)[j]))==0) {
         stop(paste("Pollutant definitions file contains unexpected column: ", names(pollutants)[j]))
         quit("no")
      }
   }
   #Check that fields that should be numeric are numeric
   if (!is.numeric(pollutants$conv.ugm3.ppb)) {
      stop("The pollutant definitions file contains a non-numeric entry in the 'conv.ugm3.ppb' column")
      quit("no")
   }
   if (!is.numeric(pollutants$min.allowed)) {
     stop("The pollutant definitions file contains a non-numeric entry in the 'min.allowed' column")
     quit("no")
   }
   if (!is.numeric(pollutants$max.allowed)) {
     stop("The pollutant definitions file contains a non-numeric entry in the 'max.allowed' column")
     quit("no")
   }
   #Convert pollutant names from factors to strings
   pollutants$pollutant=tolower(as.character(pollutants$pollutant))
   #Error handling:
   #Make sure that the units entered are allowed
   pollutants$output.units = tolower(as.character(pollutants$output.units))
   for (j in 1:length(pollutants$pollutant)) {
     if (length(which(allowed.units==pollutants$output.units[j]))==0) {
       stop(paste("Concentration units for",toupper(pollutants$pollutant[j]),"must be either ppb, ppm, ng/m3, ug/m3, mg/m3 or g/m3"))
       quit("no")
     } 
   }
   
   #READ AND PROCESS OBS DATA POLLUTANT DEFINITION FILE (separator here is ",")
   #header = "pollutant, obs.alias, obs.units, obs.avg.time.hours, obs.statistic
   cat("Importing obs data pollutant details...\n")
   pollutants.obs=read.csv(pollutants.obs.file,header=TRUE,sep=",")
   names(pollutants.obs) = tolower(names(pollutants.obs))
   pol.obs.expected.col.names = c("pollutant","obs.alias","obs.units","obs.avg.time.hours","obs.statistic")
   for (j in 1:length(pol.obs.expected.col.names)) {
     if (length(which(names(pollutants.obs)==pol.obs.expected.col.names[j]))==0) {
       stop(paste("Missing column in obs data pollutant definitions file: ",pol.obs.expected.col.names[j]))
       quit("no")
     }
   }
   for (j in 1:length(names(pollutants.obs))) {
     if (length(which(pol.obs.expected.col.names==names(pollutants.obs)[j]))==0) {
       stop(paste("Obs data pollutant definitions file contains unexpected column: ", names(pollutants.obs)[j]))
       quit("no")
     }
   }
   #Check that fields that should be numeric are numeric
   if (!is.numeric(pollutants.obs$obs.avg.time.hours)) {
     stop("The obs data pollutant definitions file contains a non-numeric entry in the 'obs.avg.time.hours' column")
     quit("no")
   }
   #Convert pollutant names from factors to strings
   pollutants.obs$pollutant=tolower(as.character(pollutants.obs$pollutant))
   pollutants.obs$obs.alias=tolower(as.character(pollutants.obs$obs.alias))
   #Error handling:
   #Make sure that the units entered are allowed
   pollutants.obs$obs.units = tolower(as.character(pollutants.obs$obs.units))
   for (j in 1:length(pollutants.obs$pollutant)) {
     if (length(which(allowed.units==pollutants.obs$obs.units[j]))==0) {
       stop(paste("Obs pollutant data concentration units for",toupper(pollutants.obs$pollutant[j]),"must be either ppb, ppm, ng/m3, ug/m3, mg/m3 or g/m3"))
       quit("no")
     } 
   }

   #Make sure the statistics entered are allowed
   pollutants.obs$obs.statistic=tolower(as.character(pollutants.obs$obs.statistic))
   for (j in 1:length(pollutants.obs$pollutant)) {
     if (length(which(allowed.statistics==pollutants.obs$obs.statistic[j]))==0) {
       stop(paste("Obs data averaging statistic for",toupper(pollutants.obs$pollutant[j]),"must be either max, mean or rolling mean"))
       quit("no")
     } 
   }
   #Combine the main pollutants file with the obs pollutant file
   pollutants = merge(pollutants,pollutants.obs,by="pollutant")
   if (length(pollutants$pollutant)==0) {
     stop("The obs pollutant data file does not contain any of the same pollutants as the main pollutants data file")
     quit("no")      
   }
   
   #READ AND PROCESS THE STATION DETAILS FILE (separator here is ",")
   #header = "name, type"
   cat("Importing station details...\n")
   stations = read.csv(stations.file,header=TRUE,sep=",")
   names(stations) = tolower(names(stations))
   if (bmod.ncdf & any(mod.ncdf.format==c("CMAQ","AIRSHEDS","CAMS Ensemble","CAMS Ensemble CERC archive","Prevair CERC archive"))) {
      stn.expected.col.names = c("name","type","latitude","longitude")
   } else {
      stn.expected.col.names = c("name","type")
   }
   for (i in 1:length(stn.expected.col.names)) {
      if (length(which(names(stations)==stn.expected.col.names[i]))==0) {
         stop(paste("Missing column in station information file: ",stn.expected.col.names[i]))
         quit("no")
      }
   }
   stations$name = tolower(as.character(stations$name))
   stations$type = tolower(as.character(stations$type))

} else {
   #Removed the old modelled data from pollutants data
   pollutants = subset(pollutants,select=c("pollutant", "output.units", "conv.ugm3.ppb","obs.alias","obs.units","obs.avg.time.hours","obs.statistic"))
} # end if bNewData

#Load the required packages
packages.required = c("openair","latticeExtra","plyr")
if (bmod.ncdf) packages.required = c(packages.required,"ncdf","GEOmap","fields","maps","akima")
update.install.load.packages(packages.required,mirror,update.packages)

#Get all modelled data
cat("Importing modelled data...\n")

#READ AND PROCESS MOD DATA POLLUTANT DEFINITION FILE (separator here is ",")
#header = "pollutant, mod.alias, mod.units, mod.avg.time.hours, mod.statistic"
cat("Importing modelled data pollutant details...\n")
pollutants.mod=read.csv(pollutants.mod.file,header=TRUE,sep=",")
names(pollutants.mod) = tolower(names(pollutants.mod))
pol.mod.expected.col.names = c("pollutant","mod.alias","mod.units","mod.avg.time.hours","mod.statistic")
for (j in 1:length(pol.mod.expected.col.names)) {
    if (length(which(names(pollutants.mod)==pol.mod.expected.col.names[j]))==0) {
       stop(paste("Missing column in mod data pollutant definitions file: ",pol.mod.expected.col.names[j]))
       quit("no")
    }
}
for (j in 1:length(names(pollutants.mod))) {
    if (length(which(pol.mod.expected.col.names==names(pollutants.mod)[j]))==0) {
       stop(paste("Mod data pollutant definitions file contains unexpected column: ", names(pollutants.mod)[j]))
       quit("no")
    }
}
#Check that fields that should be numeric are numeric
if (!is.numeric(pollutants.mod$mod.avg.time.hours)) {
   stop("The mod data pollutant definitions file contains a non-numeric entry in the 'mod.avg.time.hours' column")
   quit("no")
}
#Convert pollutant names from factors to strings
pollutants.mod$pollutant=tolower(as.character(pollutants.mod$pollutant))
pollutants.mod$mod.alias=tolower(as.character(pollutants.mod$mod.alias))
#Error handling:
#Make sure that the units entered are allowed
pollutants.mod$mod.units = tolower(as.character(pollutants.mod$mod.units))
for (j in 1:length(pollutants.mod$pollutant)) {
    if (length(which(allowed.units==pollutants.mod$mod.units[j]))==0) {
       stop(paste("Mod data concentration units for",toupper(pollutants.mod$pollutant[j]),"must be either ppb, ppm, ng/m3, ug/m3, mg/m3 or g/m3"))
       quit("no")
    } 
}
#Make sure the statistics entered are allowed
pollutants.mod$mod.statistic=tolower(as.character(pollutants.mod$mod.statistic))
for (j in 1:length(pollutants.mod$pollutant)) {
    if (length(which(allowed.statistics==pollutants.mod$mod.statistic[j]))==0 ) {
       stop(paste("Mod data averaging statistics for",toupper(pollutants.mod$pollutant[j]),"must be either max, mean or rolling mean"))
       quit("no")
    } 
}
#Combine the main pollutants data with the mod pollutant data
pollutants = merge(pollutants,pollutants.mod,by="pollutant")
if (length(pollutants$pollutant)==0) {
   stop("The mod pollutant data file does not contain any of the required pollutants")
   quit("no")      
}
#Make sure that the observed averaging time and statistic is either same as the modelled data or 1 hour mean
not.same = subset(pollutants, obs.statistic!=mod.statistic | obs.avg.time.hours!=mod.avg.time.hours)
not.same.not.1hrmean = subset(not.same, obs.avg.time.hours!=1 | obs.statistic!="mean")
if (length(not.same.not.1hrmean$pollutant) !=0 ) {
   stop("The observed data must either be a 1 hour mean or the same averaging time and statistic as the modelled data.")
   quit("no")
}

airsheds.text = "" #initialise
if (bmod.ncdf) {
   mod.missing.data = mod.ncdf.missing.data
   if (bmod.ncdf.select.dir) {
     if (any(mod.ncdf.format==c("CAMS Ensemble CERC archive","Prevair CERC archive"))) {
       start.date.archive.dateformat = as.POSIXct(start.date.archive,format="%Y%m%d",tz="GMT")
       end.date.archive.dateformat = as.POSIXct(end.date.archive,format="%Y%m%d",tz="GMT")
       archive.dirs.dates = seq(from = start.date.archive.dateformat, to = end.date.archive.dateformat,by="day")  
       archive.dirs.dates.text = format(archive.dirs.dates,format="%Y%m%d")
       if (forecast.day.archive == "0") {
         first.timestep = 1
         last.timestep = 24
       } else if (forecast.day.archive == "1") {
         first.timestep = 25
         last.timestep = 48
       } else if (forecast.day.archive == "2") {
         first.timestep = 49
         last.timestep = 72
       } else if (forecast.day.archive == "3") {
         first.timestep = 73
         last.timestep = 96
       }
       archive.timesteps = formatC(seq(from=first.timestep,to=last.timestep),width=2,flag="0")
       archive.dirs = paste(mod.ncdf.dir,.Platform$file.sep,archive.dirs.dates.text,sep="")
       if (mod.ncdf.format=="Prevair CERC archive") {
         #Prevair has one file per directory
         mod.ncdf.files = vector(mode="character",length=length(archive.dirs))
       } else if (mod.ncdf.format=="CAMS Ensemble CERC archive") {
         #CAMS Ens has one file per forecast hour
         mod.ncdf.files = vector(mode="character",length=length(archive.dirs)*length(archive.timesteps)*2)
       }
       for (t in 1:length(archive.dirs)) {
         if (mod.ncdf.format=="Prevair CERC archive") {
           mod.ncdf.files = paste(archive.dirs,.Platform$file.sep,"out.",archive.dirs.dates.text,".nc",sep="")
         } else if (mod.ncdf.format=="CAMS Ensemble CERC archive") {
           iStart1 = (t-1)*length(archive.timesteps)+1
           iEnd1 = iStart1+length(archive.timesteps)-1
           mod.ncdf.files[iStart1:iEnd1] = paste(archive.dirs[t],.Platform$file.sep,"HRES_ENS_",archive.dirs.dates.text[t],"00+0",archive.timesteps,".nc",sep="")
           iStart2 = length(archive.dirs)*length(archive.timesteps)+(t-1)*length(archive.timesteps)+1
           iEnd2 = iStart2+length(archive.timesteps)-1
           mod.ncdf.files[iStart2:iEnd2] = paste(archive.dirs[t],.Platform$file.sep,"HRES_ENS_",archive.dirs.dates.text[t],"00+0",archive.timesteps,"_addspecies.nc",sep="")
         }
       }
       if (mod.ncdf.format=="CAMS Ensemble CERC archive") {first.timestep=0}
     } else {
       first.timestep=0
       mod.ncdf.files = list.files(path=mod.ncdf.dir,pattern=".nc",full.names=TRUE)
     }
   } else if (bmod.ncdf.select.file) {
     first.timestep=0
     mod.ncdf.files = as.vector(mod.ncdf.file)
   }
   if (mod.ncdf.format == "CAMS Ensemble CERC archive") {
     bDoneFirstFile.a = FALSE
     bDoneFirstFile.b = FALSE
   } else {
     bDoneFirstFile = FALSE
   }
   for (f in 1:length(mod.ncdf.files)) {
      modelled.nc = tryCatch(open.ncdf(mod.ncdf.files[f]),error=function(e) return(NA))
      if (is.na(modelled.nc)) next
      cat(paste("   Importing",basename(mod.ncdf.files[f]),"...\n"))
      if (mod.ncdf.format == "CMAQ") {
         #Check that the data are gridded and the projection is lambert conformal conic
         missing.atts = ""
         FTYPE.att = att.get.ncdf(modelled.nc,0,"FTYPE")
         if (FTYPE.att$hasatt) {FTYPE = FTYPE.att$value} else {missing.atts = paste(missing.atts,"FTYPE")}
         GDTYP.att = att.get.ncdf(modelled.nc,0,"GDTYP")
         if (GDTYP.att$hasatt) {GDTYP = GDTYP.att$value} else {missing.atts = paste(missing.atts,"GDTYP")}
         #Stop if any expected attributes are missing
         if (nchar(missing.atts)!=0) {
            stop(paste("CMAQ data is not in the expected format; the following attributes are missing:",missing.atts))
            quit("no")
         }
         if (FTYPE != 1) {
            stop("CMAQ data must be gridded")
            quit("no")
         }
         if (GDTYP != 2) {
            stop("CMAQ data grid projection must be Lambert Conformal Conic")
            quit("no")
         }
         
         #Get all the required attributes
         missing.atts = ""
         P_ALP.att = att.get.ncdf(modelled.nc,0,"P_ALP")
         if (P_ALP.att$hasatt) {P_ALP = P_ALP.att$value} else {missing.atts = paste(missing.atts,"P_ALP")}
         P_BET.att = att.get.ncdf(modelled.nc,0,"P_BET")
         if (P_BET.att$hasatt) {P_BET = P_BET.att$value} else {missing.atts = paste(missing.atts,"P_BET")}
         XCENT.att = att.get.ncdf(modelled.nc,0,"XCENT")
         if (XCENT.att$hasatt) {XCENT = XCENT.att$value} else {missing.atts = paste(missing.atts,"XCENT")}
         YCENT.att = att.get.ncdf(modelled.nc,0,"YCENT")
         if (YCENT.att$hasatt) {YCENT = YCENT.att$value} else {missing.atts = paste(missing.atts,"YCENT")}
         XORIG.att = att.get.ncdf(modelled.nc,0,"XORIG")
         if (XORIG.att$hasatt) {XORIG = XORIG.att$value} else {missing.atts = paste(missing.atts,"XORIG")}
         YORIG.att = att.get.ncdf(modelled.nc,0,"YORIG")
         if (YORIG.att$hasatt) {YORIG = YORIG.att$value} else {missing.atts = paste(missing.atts,"YORIG")}
         XCELL.att = att.get.ncdf(modelled.nc,0,"XCELL")
         if (XCELL.att$hasatt) {XCELL = XCELL.att$value} else {missing.atts = paste(missing.atts,"XCELL")}
         YCELL.att = att.get.ncdf(modelled.nc,0,"YCELL")
         if (YCELL.att$hasatt) {YCELL = YCELL.att$value} else {missing.atts = paste(missing.atts,"YCELL")}
         SDATE.att = att.get.ncdf(modelled.nc,0,"SDATE")
         if (SDATE.att$hasatt) {SDATE = SDATE.att$value} else {missing.atts = paste(missing.atts,"SDATE")}
         STIME.att = att.get.ncdf(modelled.nc,0,"STIME")
         if (STIME.att$hasatt) {STIME = STIME.att$value} else {missing.atts = paste(missing.atts,"STIME")}
         TSTEP.att = att.get.ncdf(modelled.nc,0,"TSTEP")
         if (TSTEP.att$hasatt) {TSTEP = TSTEP.att$value} else {missing.atts = paste(missing.atts,"TSTEP")}
         VARLIST.att = att.get.ncdf(modelled.nc,0,"VAR-LIST")
         if (VARLIST.att$hasatt) {VARLIST = VARLIST.att$value} else {missing.atts = paste(missing.atts,"VAR-LIST")}
         #Stop if any expected attributes are missing
         if (nchar(missing.atts)!=0) {
            stop(paste("CMAQ data is not in the expected format; the following attributes are missing:",missing.atts))
            quit("no")
         }
         
         #Get all the required dimensions
         missing.dims = ""
         dim.names = names(modelled.nc$dim)
         if (length(which(dim.names =="TSTEP"))==0) {missing.dims = paste(missing.dims,"TSTEP")} else {num.timesteps = modelled.nc$dim$TSTEP$len}
         if (length(which(dim.names =="COL"))==0) {missing.dims = paste(missing.dims,"COL")} else {num.cols = modelled.nc$dim$COL$len}
         if (length(which(dim.names =="ROW"))==0) {missing.dims = paste(missing.dims,"ROW")} else {num.rows = modelled.nc$dim$ROW$len}
         if (length(which(dim.names =="LAY"))==0) {missing.dims = paste(missing.dims,"LAY")} else {num.layers = modelled.nc$dim$LAY$len}
         if (nchar(missing.dims) != 0) {
            stop(paste("CMAQ data is not in the expected format; the following dimensions are missing:",missing.dims))
            quit("no")
         }
         
         #Calculate the CMAQ grid coordinates (in map projection)
         grid.xcoords = seq(from=(XORIG+0.5*XCELL),to=(XORIG+0.5*XCELL+(num.cols-1)*XCELL),by=XCELL)
         grid.ycoords = seq(from=(YORIG+0.5*YCELL),to=(YORIG+0.5*YCELL+(num.rows-1)*YCELL),by=YCELL)
         #Set the map projection (uses GEOmap package), currently assuming lambert cc (type=3)
         map.proj = setPROJ(type = 3, LAT0 = YCENT, LON0 = XCENT, LAT1 = P_ALP, LAT2 = P_BET)
         #Convert station coordinates from lat long to x y (m)
         xy = cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
         stations$map.pos.x = xy$x
         stations$map.pos.y = xy$y   
         allstations.location = matrix(data=c(xy$x,xy$y),nrow=length(stations$name), ncol=2)
         regular.grid = TRUE
         
         #Now do something with the CMAQ timestep data to get a date/time column
         #Convert start time in HHMMSS format to hours, mins and secs and get a start date
         start.time = convert.time(STIME,"string")
         SDATETIME = paste(SDATE,start.time)
         start.date = as.POSIXct(SDATETIME,format="%Y%j %H %M %S",tz="GMT")
         #Convert timestep in HHMMSS format to seconds
         timestep = convert.time(TSTEP,"seconds")
         #Set up the sequence of date/times in the CMAQ file
         timestep.datetime = seq(from = start.date, to = (start.date+(num.timesteps-1)*timestep),by=timestep)              
         
         #Need the surface layer. Which layer is that in this dataset?
         layer.required = num.layers
         
         if (btesting) {
            #Set axis labels for testing output
            xaxis.label = "X (m)"
            yaxis.label = "Y (m)"
         }
         
      }
      if (any(mod.ncdf.format==c("AIRSHEDS","CAMS Ensemble","CAMS Ensemble CERC archive","Prevair CERC archive"))) {
  
         #Get all the required dimensions
         missing.dims = ""
         dim.names = names(modelled.nc$dim)
         if (mod.ncdf.format=="Prevair CERC archive") {
            if (length(which(dim.names =="west_east"))==0) {missing.dims = paste(missing.dims,"west_east")} else {num.cols = modelled.nc$dim$west_east$len}
            if (length(which(dim.names =="south_north"))==0) {missing.dims = paste(missing.dims,"south_north")} else {num.rows = modelled.nc$dim$south_north$len}
            if (length(which(dim.names =="bottom_top"))==0) {missing.dims = paste(missing.dims,"bottom_top")} else {num.layers = modelled.nc$dim$bottom_top$len}
            num.timesteps=24
         } else {
            if (length(which(dim.names =="longitude"))==0) {missing.dims = paste(missing.dims,"longitude")} else {num.cols = modelled.nc$dim$longitude$len}
            if (length(which(dim.names =="latitude"))==0) {missing.dims = paste(missing.dims,"latitude")} else {num.rows = modelled.nc$dim$latitude$len}
            if (length(which(dim.names =="level"))==0) {missing.dims = paste(missing.dims,"level")} else {num.layers = modelled.nc$dim$level$len}
            if (length(which(dim.names =="time"))==0) {missing.dims = paste(missing.dims,"time")} else {num.timesteps = modelled.nc$dim$time$len}
         }
         if (nchar(missing.dims) != 0) {
            stop(paste(mod.ncdf.format,"data is not in the expected format; the following dimensions are missing:",missing.dims))
            quit("no")
         }
         
         #Get all the variable names in order to find the latitude and longitude variables
         var.names = rep("",modelled.nc$nvars)
         for (i in 1: modelled.nc$nvars) {
            var.names[i] = modelled.nc$var[[i]]$name
         }
         var.lon = which(var.names=="lon")
         if (length(var.lon)!=0) {
            longitude.var.name = "lon"
         } else {
            longitude.var.name = "longitude"
         }
         var.lat = which(var.names=="lat")
         if (length(var.lat)!=0) {
            latitude.var.name = "lat"
         } else {
            latitude.var.name = "latitude"
         }

         #Get the station coordinates in the right format for interpolation function
         allstations.location = matrix(data=c(stations$longitude,stations$latitude),nrow=length(stations$name), ncol=2)
         stations$map.pos.x = stations$longitude
         stations$map.pos.y = stations$latitude     
         
         #Get the AIRSHEDS grid
         grid.xcoords = get.var.ncdf(modelled.nc,longitude.var.name)
         grid.ycoords = get.var.ncdf(modelled.nc,latitude.var.name)
         #Get all the longitude values in the range [-180,180)
         grid.xcoords[which(grid.xcoords >= 180)] = grid.xcoords[which(grid.xcoords >= 180)] - 360 
         grid.xcoords[which(grid.xcoords < -180)] = grid.xcoords[which(grid.xcoords < -180)] + 360 
         #Work out whether the grid is regular or irregular
         num.dimensions.x = length(dim(grid.xcoords))
         num.dimensions.y = length(dim(grid.ycoords))
         if (num.dimensions.x==num.dimensions.y) {
            num.dimensions = num.dimensions.x
         } else {
            stop(paste("netCDF variables '",longitude.var.name,"' and '",latitude.var.name,"' do not have the same number of dimensions",sep=""))
            quit("no")
         }         
         if (num.dimensions==1) {
            regular.grid = TRUE
         } else if(num.dimensions==2) {
            regular.grid = FALSE
         } else {
            stop(paste("netCDF variables '",longitude.var.name,"' and '",latitude.var.name,"' have ",num.dimensions," dimensions but only 1 or 2 are allowed",sep=""))
            quit("no")
         }
         
         if (mod.ncdf.format == "Prevair CERC archive") {
            time = get.var.ncdf(modelled.nc,"Times",start=c(1,first.timestep+1),count=c(19,num.timesteps))
            timestep.datetime = as.POSIXct(time,format="%Y-%m-%d_%H:%M:%S",tz="GMT") 
         } else {
            time = get.var.ncdf(modelled.nc,"time",start=1,count=num.timesteps)
            secs.in.hour = 60*60
            secs.in.day = 24*60*60
            if (mod.ncdf.format == "AIRSHEDS") {
            #Date/time is given in 'time' variable as number of seconds/hours/days since time.origin
               time.units = att.get.ncdf(modelled.nc,"time","units")$value
               if (time.units == 0) {
                  #This means the 'units' attribute is not present
                  stop("'time' variable is missing the 'units' attribute")
                  quit("no")
               }
               btime.seconds = length(grep("seconds",time.units))>0
               btime.hours = length(grep("hours",time.units))>0
               btime.days = length(grep("days",time.units))>0
               if (btime.seconds) {
                  time.origin = sub("seconds since ","",time.units)
                  #Need to round to the nearest hour because accuracy is not sufficient
                  time.secs = round(time/secs.in.hour)*secs.in.hour
               } else if (btime.hours) {
                  time.origin = sub("hours since ","",time.units)
                  #Round to the nearest hour just for completeness
                  time.secs = round(time)*secs.in.hour
               } else if (btime.days) {
                  time.origin = sub("days since ","",time.units)
                  #Need to round to the nearest hour because accuracy is not sufficient
                  time.secs = (floor(time)+(1/24)*round(24*(time-floor(time))))*secs.in.day
               } else {
                  stop("Time units must be seconds, hours or days")
                  quit("no")
               }
            } else if (any(mod.ncdf.format == c("CAMS Ensemble","CAMS Ensemble CERC archive"))) {
               #Date/time is given in 'time' variable as number of hours since time.origin
               time.units = att.get.ncdf(modelled.nc,"time","long_name")$value
               if (time.units == 0) {
                  #This means the 'long_name' attribute is not present
                  stop("'time' variable is missing the 'long_name' attribute")
                  quit("no")
               }
               time.origin = sub ("FORECAST time from ","",time.units)
               time.origin = as.POSIXct(time.origin,format="%Y%m%d%H",tz="GMT")
               time.secs = time*secs.in.hour
            }
            timestep.datetime = as.POSIXct(time.secs,origin=time.origin,tz="GMT")
         }
         #Need the surface layer. Which layer is that in this dataset?
         layer.required = 1
         
         if (mod.ncdf.format == "AIRSHEDS") {
            #Get the global attributes for dumping in graphical output
            missing.atts = ""
            airsheds.author.att = att.get.ncdf(modelled.nc,0,"author")
            if (airsheds.author.att$hasatt) {airsheds.author = airsheds.author.att$value} else {missing.atts = paste(missing.atts,"author")}
            airsheds.institution.att = att.get.ncdf(modelled.nc,0,"institution")
            if (airsheds.institution.att$hasatt) {airsheds.institution = airsheds.institution.att$value} else {missing.atts = paste(missing.atts,"institution")}
            airsheds.version.att = att.get.ncdf(modelled.nc,0,"version")
            if (airsheds.version.att$hasatt) {airsheds.version = airsheds.version.att$value} else {missing.atts = paste(missing.atts,"version")}
            airsheds.model.att = att.get.ncdf(modelled.nc,0,"model")
            if (airsheds.model.att$hasatt) {airsheds.model = airsheds.model.att$value} else {missing.atts = paste(missing.atts,"model")}
            airsheds.expid.att = att.get.ncdf(modelled.nc,0,"expid")
            if (airsheds.expid.att$hasatt) {airsheds.expid = airsheds.expid.att$value} else {missing.atts = paste(missing.atts,"expid")}
            #Stop if any expected attributes are missing
            if (nchar(missing.atts)!=0) {
               stop(paste("AIRSHEDS data is not in the expected format; the following attributes are missing:",missing.atts))
               quit("no")
            }
            airsheds.text = paste("AIRSHEDS data attributes:: Author:",airsheds.author,"; Institution:",airsheds.institution,"; Version:",airsheds.version,"; Model:",airsheds.model,"; Experiment ID:",airsheds.expid,sep="")
            #Check the files all have the same global data attributes
            if (f>1) {
               if (airsheds.text != airsheds.text.old) {
                  stop(paste("AIRSHEDS data attributes for this file are different to those for the previous file."))
                  quit("no")
               }
            }
            airsheds.text.old = airsheds.text
         }
         
         if (btesting) {
            #Set axis labels for testing output
            xaxis.label = "Longitude"
            yaxis.label = "Latitude"
         }
         
      }
      if (any(mod.ncdf.format==c("CMAQ","AIRSHEDS","CAMS Ensemble","CAMS Ensemble CERC archive","Prevair CERC archive"))) {   

         #Find out the pollutants that are in this file and save the variable names
         pollutants$bmod = FALSE
         pollutants$var.name = NA
         for (j in 1:length(pollutants$pollutant)) {
            for (i in 1: modelled.nc$nvars) {
               vi = modelled.nc$var[[i]]
               if (pollutants$mod.alias[j] == tolower(vi$name)) {
                  pollutants$bmod[j] = TRUE
                  pollutants$var.name[j] = vi$name
                  break
               }
            }
         }
         
         #Loop over all required pollutants and extract data for each station
         bDoneFirstPol = FALSE
         for (j in 1:length(pollutants$pollutant)) {
            #Skip this pollutant if it is not in the netCDF file
            if (!pollutants$bmod[j]) next

            #Extract all the surface layer data for the current pollutant
            pollutant.data = get.var.ncdf(modelled.nc,pollutants$var.name[j],start=c(1,1,layer.required,first.timestep+1),count=c(num.cols,num.rows,1,num.timesteps))
            
            if (btesting) {        
               #Set up the colour scale for this pollutant
               data.range = convert.units(range(pollutant.data),pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
               num.bands = 20
               interval = (data.range[2] - data.range[1]) / num.bands
               colour.scale.boundaries = signif(seq(from=data.range[1],to=data.range[2],by=interval),digits=3)
               colour.scale.colours = rainbow(num.bands,start=1/12)
               #Open the PDF file for testing graphical output
               pdf(paste(output.file.path.root,"_",f,"_",pollutants$pollutant[j],".pdf",sep=""),paper="a4",width=0,height=0,title=pdf.title)
            }           
                 
            #Loop through every timestep in the modelled data file
            for (k in 1:num.timesteps) {
               if (regular.grid) {
                  #Get the gridded data for the current timestep in the right format
                  if (num.timesteps==1) {
                     pollutant.data.obj = list(x=grid.xcoords,y=grid.ycoords,z=pollutant.data)
                  } else {
                     pollutant.data.obj = list(x=grid.xcoords,y=grid.ycoords,z=pollutant.data[,,k])
                  }
                  #Interpolate the gridded data to the station locations
                  interp.conc = interp.surface( pollutant.data.obj, allstations.location)
               } else {
                  #Interpolate from an irregular grid
                  if (num.timesteps==1) {
                     interp.conc = interpp(as.vector(grid.xcoords),as.vector(grid.ycoords),as.vector(pollutant.data),stations$longitude,stations$latitude)$z
                  } else {
                     interp.conc = interpp(as.vector(grid.xcoords),as.vector(grid.ycoords),as.vector(pollutant.data[,,k]),stations$longitude,stations$latitude)$z
                  }
               }
               onepol.allstations.onetimestep = data.frame(date=timestep.datetime[k],station=stations$name,mod=interp.conc)
               if (k==1) {
                  onepol.allstations.alltimesteps = onepol.allstations.onetimestep
               } else {
                  onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
               }
               if (btesting) {
                  #Testing only: maps of data with interpolated values marked on
                  #Convert units to output units
                  if (num.timesteps==1) {
                     pollutant.data.output.units = convert.units(pollutant.data,pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
                  } else {
                     pollutant.data.output.units = convert.units(pollutant.data[,,k],pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
                  }
                  #Work out colouring for each point
                  onepol.allstations.onetimestep$mod.output.units = convert.units(onepol.allstations.onetimestep$mod,pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
                  onepol.allstations.onetimestep$x=stations$map.pos.x
                  onepol.allstations.onetimestep$y=stations$map.pos.y
                  onepol.allstations.onetimestep$colour="white"
                  for (a in 1:length(onepol.allstations.onetimestep$colour)) {
                     for (b in 1:(length(colour.scale.boundaries)-1)) {
                        if (!is.na(onepol.allstations.onetimestep$mod.output.units[a])) {
                           if ((onepol.allstations.onetimestep$mod.output.units[a] >= colour.scale.boundaries[b]) &
                              (onepol.allstations.onetimestep$mod.output.units[a] < colour.scale.boundaries[b+1])) {
                              onepol.allstations.onetimestep$colour[a] = colour.scale.colours[b]
                              break
                           }
                        }
                     }
                  }
                  if (regular.grid) {      
                     #Map of whole modelled region
                     mymap=map('world',xlim=range(grid.xcoords),ylim=range(grid.ycoords),plot=FALSE)
                     plot1 = levelplot(
                        pollutant.data.output.units,
                        row.values=grid.xcoords,column.values=grid.ycoords,
                        at=colour.scale.boundaries,col.regions=colour.scale.colours,
                        colorkey=list(labels=as.character(colour.scale.boundaries)),
                        xlab=xaxis.label,ylab=yaxis.label,
                        main=list(paste("Checking interpolation of gridded netcdf data: gridded data and interpolated data;",
                                  toupper(pollutants$pollutant[j]),pollutants$output.units[j],timestep.datetime[k]),cex=0.7),
                        mymap=mymap,
                        panel=function(...,mymap,header.text) {
                           panel.levelplot(...)
                           panel.lines(mymap$x,mymap$y)
                        }
                     )
                     #Map of region where stations are
                     mymap=map('world',xlim=range(stations$map.pos.x),ylim=range(stations$map.pos.y),plot=FALSE)
                     plot2 = levelplot(
                        pollutant.data.output.units,
                        row.values=grid.xcoords,column.values=grid.ycoords,
                        at=colour.scale.boundaries,col.regions=colour.scale.colours,
                        colorkey=list(labels=as.character(colour.scale.boundaries)),
                        xlab=xaxis.label,ylab=yaxis.label,
                        xlim=range(stations$map.pos.x),ylim=range(stations$map.pos.y),
                        main=list(paste("Checking interpolation of gridded netcdf data: gridded data and interpolated data;",
                                  toupper(pollutants$pollutant[j]),pollutants$output.units[j],timestep.datetime[k]),cex=0.7),
                        sub=list(airsheds.text,cex=0.7),
                        mymap=mymap,
                        onepol.allstations.onetimestep=onepol.allstations.onetimestep,
                        panel=function(...,mymap,onepol.allstations.onetimestep) {
                           panel.levelplot(...)
                           panel.lines(mymap$x,mymap$y)
                           panel.points(x=onepol.allstations.onetimestep$x,
                                       y=onepol.allstations.onetimestep$y,
                                       pch=21,
                                       cex=1.5,
                                       col="black",
                                       fill=onepol.allstations.onetimestep$colour)
                        }
                     )
                     print(plot1,position=c(0,0,1,1),split=c(1,1,1,2),more=T)
                     print(plot2,position=c(0,0,1,1),split=c(1,2,1,2),more=F)
                  } else {
                     #irregular grid
                     #Map of whole modelled region
                     quilt.plot(
                        as.vector(grid.xcoords),
                        as.vector(grid.ycoords),
                        as.vector(pollutant.data.output.units),
                        nx=dim(grid.xcoords)[1],ny=dim(grid.xcoords)[2],
                        xlab=xaxis.label,ylab=yaxis.label,
                        col=colour.scale.colours,zlim=data.range,
                        main=list(paste("Checking interpolation of gridded netcdf data: gridded data and interpolated data;",
                                  toupper(pollutants$pollutant[j]),pollutants$output.units[j],timestep.datetime[k]),cex=0.7),
                        sub=list(airsheds.text,cex=0.7)
                     )
                     mymap=map('world',xlim=range(grid.xcoords),ylim=range(grid.ycoords),add=TRUE)
                     points(x=onepol.allstations.onetimestep$x,
                                    y=onepol.allstations.onetimestep$y,
                                    pch=21,
                                    cex=1.5,
                                    col="black",
                                    bg=onepol.allstations.onetimestep$colour,add=TRUE)
                  }
               }
            }
            if (btesting) {graphics.off()}
            #Rename the 'mod' column with the pollutant name
            names(onepol.allstations.alltimesteps)[3] = pollutants$pollutant[j]
            #Append to the data frame for all pollutants, timesteps and stations
            if (!bDoneFirstPol) {
               modelled.data.file = onepol.allstations.alltimesteps
               bDoneFirstPol = TRUE
            } else {
               modelled.data.file = merge(modelled.data.file,onepol.allstations.alltimesteps,all=TRUE)
            }
         } #End of loop over pollutants
         if (!bDoneFirstPol) next
         if (mod.ncdf.format == "CAMS Ensemble CERC archive") {
           if (f <= length(archive.dirs)*length(archive.timesteps)) {
             if (!bDoneFirstFile.a) {
               modelled.data.a = modelled.data.file
               bDoneFirstFile.a = TRUE
             } else {
               modelled.data.a = merge(modelled.data.a,modelled.data.file,all=TRUE)
             }
           } else {
             if (!bDoneFirstFile.b) {
               modelled.data.b = modelled.data.file
               bDoneFirstFile.b = TRUE
             } else {
               modelled.data.b = merge(modelled.data.b,modelled.data.file,all=TRUE)
             }
           }
         } else {
           if (!bDoneFirstFile) {
             modelled.data = modelled.data.file
             bDoneFirstFile = TRUE
           } else {
             modelled.data = merge(modelled.data,modelled.data.file,all=TRUE)
           }
         }
      } #End of test on if ncdf format is CMAQ or AIRSHEDS
      close.ncdf(modelled.nc)
   } #End of loop over files
   if (mod.ncdf.format == "CAMS Ensemble CERC archive") {
     if (!bDoneFirstFile.a & !bDoneFirstFile.b) {
       bDoneFirstFile = FALSE
     } else if (bDoneFirstFile.a & bDoneFirstFile.b) {
       modelled.data = merge(modelled.data.a,modelled.data.b,all=TRUE)
       modelled.data$nox = modelled.data$no+modelled.data$no2
       bDoneFirstFile = TRUE
     } else if (bDoneFirstFile.a) {
       modelled.data = modelled.data.a
       bDoneFirstFile = TRUE
     } else if (bDoneFirstFile.b) {
       modelled.data = modelled.data.b
       bDoneFirstFile = TRUE
     }
   }
   if (mod.ncdf.format != "DELTA") {
      if (!bDoneFirstFile) {
         stop("No data for any required pollutants in the netCDF data")
         quit("no")
      }
      #Check that there are data for each pollutant
      bstop = FALSE
      for (j in 1:length(pollutants$pollutant)) {
        if (!any(names(modelled.data)==pollutants$pollutant[j])) {
           bstop = TRUE
           cat(paste("Pollutant",toupper(pollutants$pollutant[j]),"is missing from the netCDF data. Please remove this pollutant from the pollutant definitions file.\n"))
        }
      }
      if (bstop) {
         stop("Missing pollutant(s) in netCDF data")
         quit("no")
      }
   }
} else if (bmod.csv) {
   #CSV file(s)
   mod.missing.data = mod.csv.missing.data
   if (mod.csv.separator=="comma") {
      mod.csv.separator=","
   } else if (mod.csv.separator=="semicolon") {
      mod.csv.separator=";"
   }
   if (bmod.csv.select.file) {
      #Read the specified file
      modelled.data = read.csv(mod.csv.file,header=TRUE,sep=mod.csv.separator,check.names=FALSE)
      #Check for duplicate or invalid column headers
      if (!all(names(modelled.data) == make.names(names(modelled.data),unique=TRUE))) {
         stop(paste("Invalid or duplicate column headers in model CSV file ", mod.csv.file))
         quit("no")
      }
      if (length(which(tolower(names(modelled.data))=="station"))==0) {
         #There is no 'station' column, so set the station name to be the filename
         modelled.data$station = strip.name(mod.csv.file)
      }
      
   } else if (bmod.csv.select.dir) {
      #Read all files in the specified directory
      modelled.data = lapply(list.files(path=mod.csv.dir,pattern=".csv",full.names=TRUE),
                             open.csv.and.add.station,csv.separator=mod.csv.separator)
      modelled.data = do.call(rbind.fill,modelled.data)
   }
   #Make sure all the column headers are in lower case
   names(modelled.data) = tolower(names(modelled.data))
   #Check all required date columns are given
   date.expected.col.names = c("year", "month", "day", "hour")
   for (i in 1:length(date.expected.col.names)) {
      if (length(which(names(modelled.data)==date.expected.col.names[i]))==0) {
         stop(paste("Missing",date.expected.col.names[i],"column in modelled .csv data file."))
         quit("no")
      }
   }
   #Turn the year, month, day, hour columns into a 'date' column that R can recognise
   #(Hours are sometimes given as 1-24 which is not handled by ISOdate function)
   modelled.data$date = ISOdate(modelled.data$year,modelled.data$month,modelled.data$day,modelled.data$hour,min=0,sec=0)
   modelled.data$date[modelled.data$hour==24] = ISOdate(modelled.data$year[modelled.data$hour==24],
                                                        modelled.data$month[modelled.data$hour==24],
                                                        modelled.data$day[modelled.data$hour==24],
                                                        23,min=0,sec=0)+3600
   #Check for erroneous dates
   if (any(is.na(modelled.data$date))) {
      stop(paste("Date in the modelled data is not in the expected format."))
      quit("no")
   }
   #Look for missing data and replace with NA
   #The data read in from csv is tested to see if it is numeric or contains text
   #If numeric, the missing values can be set to NA
   #If the data contains text, read.csv sets the column class as 'factor', which we need to process to coerce it to numeric
   # and deal with any additional text that wasn't expected.  By setting the missing value to -999
   # and applying as.numeric to the data, any NA values that are returned are unexpected and 
   # an error can be given. In order to do this, -999 needs to be added as a factor 'level'.
   # When coercing the data to numeric, the data must first be set to character strings 
   #   in case other text confuses the as.numeric function (text like #NaN can lead to 
   #   numbers being changed unexpectedly by the as.numeric function because the data doesn't look like numbers)
   # Once the data is numeric with -999 as missing values, then set -999 to NA and continue
   #Need to apply to each pollutant set individually becasue the read.csv will have set the class for each column
   pollutant.data.names = names(modelled.data)[names(modelled.data)!="date" & 
                            names(modelled.data)!="year" & 
                            names(modelled.data)!="month" & 
                            names(modelled.data)!="day" & names(modelled.data)!="hour" 
                            & names(modelled.data)!="station" ]
   for (pol in 1:length(pollutant.data.names)) { 
     #Find the missing values, as defined by the user
     #Note that RGGRunner passes NA as text whereas RScript reads in NA as special value
     # and if mod.missing.data is NA we can't use logical "==", so have to split test  
     if (is.na(mod.missing.data) || mod.missing.data=="NA") {
       missing.val.ids = is.na(modelled.data[,pollutant.data.names[pol]])
     } else {
       missing.val.ids = modelled.data[,pollutant.data.names[pol]]==mod.missing.data 
     }
     #If data is numeric (or special values) then just need to set values
     if (is.numeric(modelled.data[,pollutant.data.names[pol]])) {
       modelled.data[missing.val.ids,pollutant.data.names[pol]] = NA
     } else if (is.factor(modelled.data[,pollutant.data.names[pol]])) { 
       #Otherwise, read.csv will have read in "factor" class with mixture of numeric and text
       #Set missing value identifier to -999 (adding to the levels first) 
       levels(modelled.data[ ,pollutant.data.names[pol]]) = c(levels(modelled.data[ ,pollutant.data.names[pol]]), "-999")
       modelled.data[missing.val.ids,pollutant.data.names[pol]] = -999
       #coerce to numeric using 'try' silent to avoid warning message to user
       modelled.data[ ,pollutant.data.names[pol]] = try(as.numeric(as.character(modelled.data[ ,pollutant.data.names[pol]])), silent=TRUE)
       #If there are any NAs they have been converted from non-numeric and not the missing value so stop with error
       if (any(is.na(modelled.data[ ,pollutant.data.names[pol]]))) {
         stop(paste("Invalid pollutant data in model .csv data file."))
         quit("no")   
       }
       #otherwise continue, setting the missing values to NA
       modelled.data[missing.val.ids,pollutant.data.names[pol]] = NA
     } else { #Unidentified pollutant values - stop with error
       stop(paste("Internal error: model .csv file not read in as numeric or factor class."))
       quit("no")         
     }
   }
   #Check that there is a column for each required pollutant
   for (j in 1:length(pollutants$pollutant)) {
      this.pol.col = which(names(modelled.data)==pollutants$mod.alias[j])
      num.cols = length(this.pol.col)
      #Can only be 0 or 1 since duplicates have already been checked for
      if (num.cols==0) {
         stop(paste("Modelled CSV data is missing pollutant",toupper(pollutants$mod.alias[j])))
         quit("no")
      }
      #Rename with the main pollutant name
      names(modelled.data)[this.pol.col]=pollutants$pollutant[j]
   }

} else if (bmod.adms.pst) {
   #ADMS PST file(s)
   mod.missing.data = mod.pst.missing.data
   if (bmod.pst.select.file) {
      #Read the specified file
      modelled.data = importADMS(mod.pst.file)
   }
   if (bmod.pst.select.dir) {
      #Read all files in the specified directory
      modelled.data = lapply(list.files(path=mod.pst.dir,pattern=".pst",full.names=TRUE),
                         function(.file) importADMS(.file))
      modelled.data = do.call(rbind.fill,modelled.data)
   }
   #Change 'receptor.name' column header to 'station'
   names(modelled.data)[2]="station"
   #Remove trailing spaces from the receptor (station) names
   modelled.data$station = sub(" +$", "", modelled.data$station)
   #Check the pollutant column headers and assign the mod.alias
   for (j in 1:length(pollutants$pollutant)) {
      #Find any columns where the header is the pollutant name
      this.pol.col = which(names(modelled.data)==pollutants$mod.alias[j])
      num.cols = length(this.pol.col)
      #Check that there is one and only one column for each pollutants
      if (num.cols>1) {
         stop(paste("There must be only one column of",toupper(pollutants$mod.alias[j]),"in the .pst file(s)"))
         quit("no")      
      } else if (num.cols==0) {
         #Check whether a column header includes the pollutant name at the beginning or end
         partial.pol.col = grep(paste("^",pollutants$mod.alias[j],"\\.|\\.",pollutants$mod.alias[j],"$",sep=""), names(modelled.data))
         if (length(partial.pol.col)>1) {
            stop(paste("There must be only one column of",toupper(pollutants$mod.alias[j]),"in the .pst file(s)"))
            quit("no")          
         } else if (length(partial.pol.col)==0) {
            stop(paste(toupper(pollutants$mod.alias[j]),"is missing from the .pst file(s)"))
            quit("no") 
         } 
         this.pol.col = partial.pol.col
      } 
      #Rename the pollutant column with just the pollutant name
      names(modelled.data)[this.pol.col]=pollutants$pollutant[j]
   }
   #Look for missing data and replace with NA
   temp.modelled.data = modelled.data
   temp.modelled.data$date=NULL #remove date column because otherwise get errors in missing data search
   temp.modelled.data[temp.modelled.data==mod.missing.data]=NA
   temp.modelled.data$date=modelled.data$date #add date column back in
   modelled.data = temp.modelled.data
}
if (!(bmod.ncdf & mod.ncdf.format=="DELTA")) {
  modelled.data$station = tolower(as.character(modelled.data$station))
  #Concentrations are only valid if they lie in the valid data range
  for (j in 1:length(pollutants$pollutant)) {
    modelled.data[which(modelled.data[[pollutants$pollutant[j]]]<pollutants$min.allowed[j]),pollutants$pollutant[j]]=NA
    modelled.data[which(modelled.data[[pollutants$pollutant[j]]]>pollutants$max.allowed[j]),pollutants$pollutant[j]]=NA
  }
}

if (bNewData) {
   #OBSERVED DATA
   #Observed data can either be all in one file or in all files in one directory.
   cat("Importing observed data...\n")
   if (bobs.csv) {
      if (obs.csv.separator=="comma") {
         obs.csv.separator=","
      } else if (obs.csv.separator=="semicolon") {
         obs.csv.separator=";"
      }
      if (bobs.csv.select.file) {
         #Read the specified file
         all.obs.data = read.csv(obs.csv.file,header=TRUE,sep=obs.csv.separator,check.names=FALSE)
         #Check for duplicate or invalid column headers
         if (!all(names(all.obs.data) == make.names(names(all.obs.data),unique=TRUE))) {
            stop(paste("Invalid or duplicate column headers in observed CSV file ", obs.csv.file))
            quit("no")
         }
         if (length(which(tolower(names(all.obs.data))=="station"))==0) {
            #There is no 'station' column, so set the station name to be the filename
            all.obs.data$station = strip.name(obs.csv.file)
         }
      } else if (bobs.csv.select.dir) {
         #Read all files in the specified directory
         all.obs.data = lapply(list.files(path=obs.csv.dir,pattern=".csv",full.names=TRUE),
                                open.csv.and.add.station,csv.separator=obs.csv.separator)
         all.obs.data = do.call(rbind.fill,all.obs.data)
      }
      #Make sure all the column headers are in lower case
      names(all.obs.data) = tolower(names(all.obs.data))
      #Check all required date columns are given
      date.expected.col.names = c("year", "month", "day", "hour")
      for (i in 1:length(date.expected.col.names)) {
         if (length(which(names(all.obs.data)==date.expected.col.names[i]))==0) {
            stop(paste("Missing",date.expected.col.names[i],"column in observed .csv data file."))
            quit("no")
         }
      }
      #Turn the year, month, day, hour columns into a 'date' column that R can recognise
      #(Hours are sometimes given as 1-24 which is not handled by ISOdate function)
      all.obs.data$date = ISOdate(all.obs.data$year,all.obs.data$month,all.obs.data$day,all.obs.data$hour,min=0,sec=0)
      all.obs.data$date[all.obs.data$hour==24] = ISOdate(all.obs.data$year[all.obs.data$hour==24],
                                                  all.obs.data$month[all.obs.data$hour==24],
                                                  all.obs.data$day[all.obs.data$hour==24],
                                                  23,min=0,sec=0)+3600
      #Check for erroneous dates
      if (any(is.na(all.obs.data$date))) {
         stop(paste("Date in the observed data is not in the expected format."))
         quit("no")
      }
      #Look for missing data and replace with NA
      #The data read in from csv is tested to see if it is numeric or contains text
      #If numeric, the missing values can be set to NA
      #If the data contains text, read.csv sets the column class as 'factor', which we need to process to coerce it to numeric
      # and deal with any additional text that wasn't expected.  By setting the missing value to -999
      # and applying as.numeric to the data, any NA values that are returned are unexpected and 
      # an error can be given. In order to do this, -999 needs to be added as a factor 'level'.
      # When coercing the data to numeric, the data must first be set to character strings 
      #   in case other text confuses the as.numeric function (text like #NaN can lead to 
      #   numbers being changed unexpectedly by the as.numeric function because the data doesn't look like numbers)
      # Once the data is numeric with -999 as missing values, then set -999 to NA and continue
      #Need to apply to each pollutant set individually becasue the read.csv will have set the class for each column
      pollutant.data.names = names(all.obs.data)[names(all.obs.data)!="date" & 
                            names(all.obs.data)!="year" &  
                            names(all.obs.data)!="month" &  
                            names(all.obs.data)!="day" &  names(all.obs.data)!="hour" 
                            & names(all.obs.data)!="station" ]
      for (pol in 1:length(pollutant.data.names)) {    
        #Find the missing values, as defined by the user
        #Note that RGGRunner passes NA as text whereas RScript reads in NA as special value
        # and if mod.missing.data is NA we can't use logical "==", so have to split test  
        if (is.na(obs.csv.missing.data) || obs.csv.missing.data=="NA") {
          missing.val.ids = is.na(all.obs.data[,pollutant.data.names[pol]])
        } else {
          missing.val.ids = all.obs.data[,pollutant.data.names[pol]]==obs.csv.missing.data 
        }
        #If data is numeric (or special values) then just need to set values
        if (is.numeric(all.obs.data[,pollutant.data.names[pol]])) {
          all.obs.data[missing.val.ids,pollutant.data.names[pol]] = NA
        } else if (is.factor(all.obs.data[,pollutant.data.names[pol]])) {        
          #Otherwise, read.csv will have read in "factor" class with mixture of numeric and text
          #Set missing value identifier to -999 (adding to the levels first) 
          levels(all.obs.data[ ,pollutant.data.names[pol]]) = c(levels(all.obs.data[ ,pollutant.data.names[pol]]), "-999")
          all.obs.data[missing.val.ids,pollutant.data.names[pol]] = -999
          #coerce to numeric using 'try' silent to avoid warning message to user
          all.obs.data[ ,pollutant.data.names[pol]] = try(as.numeric(as.character(all.obs.data[ ,pollutant.data.names[pol]])), silent=TRUE)
          #If there are any NAs they have been converted from non-numeric and not the missing value so stop with error
          if (any(is.na(all.obs.data[ ,pollutant.data.names[pol]]))) {
            stop(paste("Invalid pollutant data in observed .csv data file."))
            quit("no")   
          }
          #otherwise continue, setting the missing values to NA
          all.obs.data[missing.val.ids,pollutant.data.names[pol]] = NA
        } else { #Unidentified pollutant values - stop with error
          stop(paste("Internal error: observed .csv file not read in as numeric or factor class."))
          quit("no")         
        }
      }
      #Make sure all station names are in lower case
      all.obs.data$station = tolower(as.character(all.obs.data$station))
      #Change pollutant column headers to the main pollutant name
      for (j in 1:length(pollutants$pollutant)) {
         this.pol.col = which(names(all.obs.data)==pollutants$obs.alias[j])
         num.cols = length(this.pol.col)
         #Can only be 0 or 1 since duplicates have already been checked for
         if (num.cols==1) {
            names(all.obs.data)[this.pol.col]=pollutants$pollutant[j]
         }
         #Concentrations are only valid if they lie in the valid data range
         all.obs.data[which(all.obs.data[[pollutants$pollutant[j]]]<pollutants$min.allowed[j]),pollutants$pollutant[j]]=NA
         all.obs.data[which(all.obs.data[[pollutants$pollutant[j]]]>pollutants$max.allowed[j]),pollutants$pollutant[j]]=NA         
      }
   }
   #Note that if KCL or AURN data are to be read in this is done later inside the loop over stations.
} #end if bNewData

#Check that the date ranges of the observed and modelled data intersect
#Don't need to do this for DELTA modelled data, because it has no date info
if (!(bmod.ncdf & mod.ncdf.format=="DELTA")) {
   mod.start.date = range(modelled.data$date)[1]
   mod.end.date = range(modelled.data$date)[2]
   if (bNewData) {
      if (bobs.kcl | bobs.aurn) {
         #Find the KCL year required (bearing in mind that kcl data is hour beginning, 0:23)
         mod.start.year = as.numeric(format(mod.start.date-3600,"%Y"))
         mod.end.year = as.numeric(format(mod.end.date-3600,"%Y"))
         obs.year = seq(from=mod.start.year,to=mod.end.year,by=1)
         if (bobs.kcl) {
            obs.type = "London KCL"
         } else if (bobs.aurn) {
            obs.type = "UK AURN"
         }
      } else {
         obs.start.date = range(all.obs.data$date)[1]
         obs.end.date = range(all.obs.data$date)[2]
         if (!branges.intersect(r1start=obs.start.date,r1end=obs.end.date,
                                 r2start=mod.start.date,r2end=mod.end.date)) {                     
            cat(paste("Observed data date range:",obs.start.date,"to",obs.end.date,"\n"))
            cat(paste("Modelled data date range:",mod.start.date,"to",mod.end.date,"\n"))
            stop("Date ranges of modelled and observed datasets do not overlap")
            quit("no")   
         }
      } 
   }
}
if (bNewData) {
   if (bobs.kcl | bobs.aurn) {
      years = paste(obs.year,collapse=", ")
      cat(paste(obs.type,"data will be downloaded for the following year(s):",years,"\n"))
   }
}

if (bNewData) {
   #Assign start and end date of modelled data, or observed data if DELTA data is being used, to new variables
   #These will be output as the date range on each graph.
   if (bmod.ncdf & mod.ncdf.format=="DELTA") { 
      data.start.date.for.graphs = range(all.obs.data$date)[1]
      data.end.date.for.graphs = range(all.obs.data$date)[2]
   } else {
      data.start.date.for.graphs = mod.start.date
      data.end.date.for.graphs = mod.end.date
   }
   
   #Get the start and end dates in a format suitable for output to graphs
   data.start.date.for.graphs = format(data.start.date.for.graphs,"%d/%m/%y %H:%M")
   data.end.date.for.graphs = format(data.end.date.for.graphs, "%d/%m/%y %H:%M")
}  # end if bNewData

#Loop over all stations and all pollutants
cat("Processing data...\n")
bDoneFirstStation = FALSE
if (bNewData) {
   for (i in 1:length(stations$name)) {
      #Output progress to screen
      cat(paste("Processing station ",i," of ",length(stations$name),": ",toupper(stations$name[i]),"...\n",sep=""))
      bDoneFirstPol = FALSE
      #Cycle on to the next station if this station was not modelled
      if (!(bmod.ncdf & mod.ncdf.format=="DELTA")) {
         mod.data.station = subset(modelled.data,station==stations$name[i])
         if (length(mod.data.station$date)==0) {
            cat(paste("Warning: no modelled data for station",toupper(stations$name[i]),"- skipping this station\n"))
            next
         }
      }
      
      #GET OBSERVED DATA
      if (bobs.csv) {
         #Get the data from the global data frame
         station.obs.data = subset(all.obs.data,station==stations$name[i])
         if (length(station.obs.data$date)==0) {
            cat(paste("Warning: no observed data for station",toupper(stations$name[i]),"- skipping this station\n"))
            next
         }
      } else if (bobs.kcl | bobs.aurn) {
         cat(paste("Importing ",obs.type," data for station ",toupper(stations$name[i]),"...\n",sep=""))
         #Loop over years required
         bDoneFirstYear = FALSE
         for (k in 1:length(obs.year)) {
            if (bobs.kcl) {
               station.obs.data.oneyear = tryCatch(importKCL(site=stations$name[i],year=obs.year[k]),
                   error = function(e) return("error"))
            } else if (bobs.aurn) {
               station.obs.data.oneyear = tryCatch(importAURN(site=stations$name[i],year=obs.year[k],hc=TRUE),
                   error = function(e) return("error"))
            }
            #If there's a problem extracting the data for this year, skip to the next year
            if (class(station.obs.data.oneyear)=="character") {
              if (station.obs.data.oneyear=="error") {
                next
              }
            }
            if (!bDoneFirstYear) {
               station.obs.data = station.obs.data.oneyear
               bDoneFirstYear = TRUE
            } else {
               station.obs.data = rbind.fill(station.obs.data,station.obs.data.oneyear)
            }
         }
         #If no data has been downloaded for this station, skip to next station
         if (!bDoneFirstYear) {
            cat(paste("Warning: no observed data for station",toupper(stations$name[i]),"- skipping this station\n"))
            next
         }
         #If London KCL or AURN downloaded through openair tool, 
         # add 1 hour to dates to adjust from hour beginning to hour ending
         # and truncate to dates corresponding to date range of modelled data
         if (bobs.kcl | bobs.aurn) {
            station.obs.data$date = station.obs.data$date + 3600
            station.obs.data = subset(station.obs.data, (station.obs.data$date>=mod.start.date & station.obs.data$date<=mod.end.date))
            if (length(station.obs.data$date)==0) {
              cat(paste("Warning: no observed data for station",toupper(stations$name[i]),"- skipping this station\n"))
              next
            }
         }
         for (j in 1:length(pollutants$pollutant)) {
            #Change pollutant column headers to the main pollutant name
            this.pol.col = which(names(station.obs.data)==pollutants$obs.alias[j])
            num.cols = length(this.pol.col)
            if (num.cols==1) {
               names(station.obs.data)[this.pol.col]=pollutants$pollutant[j]
            } else if (num.cols>1) {
               stop(paste("There must be no more than one column of",toupper(pollutants$obs.alias[j]),"in the observed data"))
               quit("no")      
            }
            #Concentrations are only valid if they lie in the valid data range
            station.obs.data[which(station.obs.data[[pollutants$pollutant[j]]]<pollutants$min.allowed[j]),pollutants$pollutant[j]]=NA
            station.obs.data[which(station.obs.data[[pollutants$pollutant[j]]]>pollutants$max.allowed[j]),pollutants$pollutant[j]]=NA
         }
      }
   
      for (j in 1:length(pollutants$pollutant)) {
         #Output progress to screen
         cat(paste("   Processing pollutant ",toupper(pollutants$pollutant[j]),"...\n",sep=""))
   
         #OBS DATA
         if (length(which(names(station.obs.data)==pollutants$pollutant[j]))==0) {
            cat(paste("   Warning: no observed data for pollutant",pollutants$pollutant[j],"- skipping this pollutant\n"))
            next
         }
         #Proceed with this pollutant for this station 
         obs.data = subset(station.obs.data,select=c("date",pollutants$pollutant[j]))
         #Check length of record against time averaging
         if (length(obs.data[!is.na(obs.data[,c(pollutants$pollutant[j])]), 1]) <= 1 & !(pollutants$mod.avg.time.hours[j] == pollutants$obs.avg.time.hours[j] & pollutants$mod.statistic[j] == pollutants$obs.statistic[j])) {
            cat(paste("Warning: insufficient observed data for pollutant",toupper(pollutants$pollutant[j]),"for station",toupper(stations$name[i]),"- skipping this pollutant\n"))
            next         
         }
         #Rename pollutant column as "obs"
         names(obs.data)[2]="obs"
         #Look for negative concentrations and replace with NA
         obs.data$obs[obs.data$obs < 0] = NA
         #Convert units to output units
         obs.data$obs = convert.units(obs.data$obs,pollutants$obs.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
         obs.statistic = pollutants$obs.statistic[j]
   
         #MODELLED DATA
         if (pollutants$mod.avg.time.hours[j]==24) {
            mod.avg.time = "day"
         } else {
            mod.avg.time = paste(pollutants$mod.avg.time.hours[j],"hour")
         }
         mod.statistic = pollutants$mod.statistic[j]
         if (bmod.ncdf & mod.ncdf.format=="DELTA") {
            #Get one of the fields from the netcdf data
            dataset.name = paste(toupper(stations$name[i]),"_",toupper(pollutants$mod.alias[j]),sep="")
            netcdf.data.station.pol = as.numeric(get.var.ncdf(modelled.nc,dataset.name))
            #Look for missing data and replace with NA
            netcdf.data.station.pol[netcdf.data.station.pol==mod.missing.data]=NA
            #Convert units to output units
            netcdf.data.station.pol = convert.units(netcdf.data.station.pol,pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
            #Now bind the modelled data to the observed data
            #(note we 'bind' rather than 'merge' here because the netcdf file of modelled data contains no dates)
            all.data.onestation.onepol = cbind(obs.data,netcdf.data.station.pol)
            names(all.data.onestation.onepol)[names(all.data.onestation.onepol)=="netcdf.data.station.pol"]="mod"
            
         } else {
            #Get the modelled data for this pol and station from the global data frame
            mod.data = subset(mod.data.station,select=c("date","station",pollutants$pollutant[j]))
            names(mod.data)[3] = "mod"
            #Convert units to output units
            mod.data$mod = convert.units(mod.data$mod,pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
            #Before combining with the observed data, get them both to have the same averaging time and statistic
            if (mod.statistic=="rolling mean") {
               if (obs.statistic!="rolling mean") {
                  obs.data = rollingMean(obs.data,pollutant="obs",width=pollutants$mod.avg.time.hours[j],new.name="rolling",align="right")
                  obs.data = subset(obs.data,select=c("date","rolling"))  
                  names(obs.data)[2]="obs"
               }
            } else {
               if (!(pollutants$mod.avg.time.hours[j] == pollutants$obs.avg.time.hours[j] & mod.statistic == obs.statistic)) {
                  #Need to do this to modelled data as well to get averaged values at the right hours
                  mod.data = timeAverage(mod.data,avg.time=mod.avg.time,statistic=mod.statistic)
                  obs.data = timeAverage(obs.data,avg.time=mod.avg.time,statistic=mod.statistic)
               }
            } 
            #Now merge the observed and modelled datasets (concentrations all in output units)
            all.data.onestation.onepol = merge(obs.data,mod.data,all.x=TRUE)
            names(all.data.onestation.onepol)[names(all.data.onestation.onepol)=="mod"] = import.model.name 
            rm(obs.data,mod.data) # to save virtual memory
         }
   
         #Only continue with this pollutant if there are any valid observed values to compare!         
         num.valid.obs.values = length(all.data.onestation.onepol$date[!is.na(all.data.onestation.onepol$obs)])
         if (num.valid.obs.values==0) {
            cat("   Warning: no valid observed values for comparison - skipping this pollutant\n")
            next
         }
         #Give a warning if there are no valid modelled values for this model
         num.valid.mod.values = length(all.data.onestation.onepol$date[!is.na(all.data.onestation.onepol[,c(import.model.name)])])
         if (num.valid.mod.values==0) {
            cat("   Warning: no valid modelled values for comparison\n")
         }         
                     
         #Collect the raw data (in output units) for output to CSV
         all.data.onestation.onepol$pollutant = pollutants$pollutant[j]
         all.data.onestation.onepol$station = stations$name[i]
         all.data.onestation.onepol$type = stations$type[i]
               
         if (bDoneFirstPol) {
            all.data.onestation = rbind(all.data.onestation,all.data.onestation.onepol)
         } else {
            all.data.onestation = all.data.onestation.onepol
            bDoneFirstPol = TRUE
         }
      } #End of loop over pollutants
      
      if (bDoneFirstPol) {
         if (bDoneFirstStation) {
            all.data = rbind(all.data,all.data.onestation)
         } else {
            all.data = all.data.onestation
            bDoneFirstStation = TRUE
         }
         rm(all.data.onestation) # to save virtual memory
      }
   
   } #End of loop over stations

   if (!bDoneFirstStation) { 
      stop("No valid data for any stations - stopping")  #no valid observed data
      quit("no")
   }

   #Check for valid data by pollutant
   for (j in 1:length(pollutants$pollutant)) {
      #Observed
      if (length(subset(all.data, pollutant==pollutants$pollutant[j])[,1])==0) { 
            cat("Warning: no valid values for pollutant",toupper(pollutants$pollutant[j]),"- check the obs.alias in your pollutants file\n")
            next
      }   
      #Modelled
      if (length(which(!is.na(subset(all.data, pollutant==pollutants$pollutant[j])[,import.model.name])))==0) { 
            cat("Warning: no valid modelled values for pollutant",toupper(pollutants$pollutant[j]),"- check the mod.alias in your pollutants file\n")
            next
      }    
   }
   
} else if (bAddData) {
   #loop over stations in case importing DELTA format modelled data
   for (i in 1:length(stations$name)) {
      #Output progress to screen
      cat(paste("Processing station ",i," of ",length(stations$name),": ",toupper(stations$name[i]),"...\n",sep=""))
      #Cycle on to the next station if this station was not modelled
      if (!(bmod.ncdf & mod.ncdf.format=="DELTA")) {
         mod.data.station = subset(modelled.data,station==stations$name[i])
         if (length(mod.data.station$date)==0) {
            cat(paste("   Warning: no modelled data for station",toupper(stations$name[i]),"- skipping this station\n"))
            next
         }
      }
      #Check the global data frame for observed data
      station.obs.data = subset(all.data,all.data$station==stations$name[i])
      if (length(station.obs.data$date)==0) {
         cat(paste("   Warning: no stored observed data for station",toupper(stations$name[i]),"- skipping this station\n"))
         next
      }
      
      bDoneFirstPol = FALSE
      for (j in 1:length(pollutants$pollutant)) {
         #Output progress to screen
         cat(paste("   Processing pollutant ",toupper(pollutants$pollutant[j]),"...\n",sep=""))

         #check for existing data for pollutant at station
         station.pollutant.obs.data = subset(all.data,all.data$station==stations$name[i] & all.data$pollutant==pollutants$pollutant[j])
         if (length(station.pollutant.obs.data$date)==0) {
            cat(paste("   Warning: no stored observed data for",toupper(pollutants$pollutant[j]),"at station",toupper(stations$name[i]),"- skipping this pollutant\n"))
            next
         }

         #MODELLED DATA
         if (pollutants$mod.avg.time.hours[j]==24) {
            mod.avg.time = "day"
         } else {
            mod.avg.time = paste(pollutants$mod.avg.time.hours[j],"hour")
         }
         mod.statistic = pollutants$mod.statistic[j]
         if (bmod.ncdf & mod.ncdf.format=="DELTA") {
            #EXTRACT DATE, OBS AND OTHER MODELLED DATA from workspace
            # merge obs from all.data by date to pick out data stored in the loaded workspace 
            obs.data.onestation.onepol = subset(all.data[ ,c("date","type","obs",existing.modelled.data.names)], all.data$station==stations$name[i] & all.data$pollutant==pollutants$pollutant[j])        
            #LOAD MODELLED DATA
            #Get one of the fields from the netcdf data
            dataset.name = paste(toupper(stations$name[i]),"_",toupper(pollutants$mod.alias[j]),sep="")
            netcdf.data.station.pol = as.numeric(get.var.ncdf(modelled.nc,dataset.name))
            #Look for missing data and replace with NA
            netcdf.data.station.pol[netcdf.data.station.pol==mod.missing.data]=NA
            #Convert units to output units
            netcdf.data.station.pol = convert.units(netcdf.data.station.pol,pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
            #Now bind the modelled data to the timeline
            #(note we 'bind' rather than 'merge' here because the netcdf file of modelled data contains no dates)
            if (length(timeline) != length(netcdf.data.station.pol)) {
               stop("The length of the DELTA format modelled data does not match the observed data in the loaded workspace.")
               quit("no")
            }
            # recreate timeline from range and time averaging
            timeline = seq(from=range(all.data$date)[1], to=range(all.data$date)[2], by=paste(pollutants$mod.avg.time.hours[j],"hour"))
            all.data.onestation.onepol = as.data.frame(timeline)
            names(all.data.onestation.onepol)="date" 
            all.data.onestation.onepol$mod = netcdf.data.station.pol
            all.data.onestation.onepol = merge(obs.data.onestation.onepol,all.data.onestation.onepol,by="date",all.x=TRUE)
            rm(obs.data.onestation) # to save virtual memory
         } else {
            #EXTRACT DATE and OBS from workspace
            all.data.onestation.onepol = subset(all.data[ ,c("date","station","type","obs",existing.modelled.data.names)], all.data$station==stations$name[i] & all.data$pollutant==pollutants$pollutant[j])
            #LOAD MODELLED DATA 
            #Get the new modelled data for this pol and station from the global data frame
            mod.data = subset(mod.data.station,select=c("date",pollutants$pollutant[j]))
            names(mod.data)[names(mod.data)==pollutants$pollutant[j]] = "mod"
            #Convert units to output units
            mod.data$mod = convert.units(mod.data$mod,pollutants$mod.units[j],pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
            #convert the modelled data to the time averaging given for previously loaded data
            if (!(pollutants$mod.avg.time.hours[j] == pollutants$obs.avg.time.hours[j] & mod.statistic == obs.statistic)) {
               #Need to do this to modelled data to get averaged values at the right hours
               mod.data = timeAverage(mod.data,avg.time=mod.avg.time,statistic=mod.statistic)
            }            
            #Now merge the observed and modelled datasets (concentrations all in output units)
            all.data.onestation.onepol = merge(all.data.onestation.onepol,mod.data,by="date",all.x=TRUE)  
            rm(mod.data) # to save virtual memory        
         }
         
         #Only continue with this pollutant if there are any valid observed values to compare!         
         num.valid.obs.values = length(all.data.onestation.onepol$date[!is.na(all.data.onestation.onepol$obs)])
         if (num.valid.obs.values==0) {
            cat("   Warning: no valid observed values for comparison - skipping this pollutant\n")
            next
         }
         #Give a warning if there are no valid modelled values for this model
         num.valid.mod.values = length(all.data.onestation.onepol$date[!is.na(all.data.onestation.onepol$mod)])
         if (num.valid.mod.values==0) {
            cat("   Warning: no valid modelled values for comparison\n")
         }

         # now tidy all.data.onestation.onepol to format of one data
         names(all.data.onestation.onepol)[names(all.data.onestation.onepol)=="mod"] = import.model.name
         all.data.onestation.onepol$pollutant = pollutants$pollutant[j]
      
         if (bDoneFirstPol) {
            all.data.onestation = rbind(all.data.onestation,all.data.onestation.onepol)
         } else {
            all.data.onestation = all.data.onestation.onepol
            bDoneFirstPol = TRUE
         }
                                       
      } #End of loop over pollutants
  
      if (bDoneFirstPol) {  
         if (bDoneFirstStation) {
            all.data.new = rbind(all.data.new,all.data.onestation)
         } else {
            all.data.new = all.data.onestation
            bDoneFirstStation = TRUE
         }
         rm(all.data.onestation) # to save virtual memory
      }

   } #End of loop over stations

   #The all.data.new data frame combines the previously stored observations and 'old' model
   # data, with 'new' added/imported model data. If there is valid observed data
   # then bDoneFirstStation will be TRUE but we need to check there is some valid new model data 
   num.valid.mod.values = length(all.data.new$date[!is.na(all.data.new[,import.model.name])])
   if (num.valid.mod.values==0) {
      # There must have been an error in importing the new model data 
      stop("No valid new modelled data for any stations - stopping") 
      quit("no")
   } else {
      rm(all.data)
      all.data = all.data.new
      rm(all.data.new)
   }

}

if (boutput.stats.as.csv) {
   cat("Exporting raw data to CSV...\n")
   write.csv(all.data, file=paste(output.file.path.root,"_raw_data.csv", sep=""), row.names=FALSE)
}

#Save this workspace for importing later into the evaluation and diagnostics tools
cat("Saving workspace...\n")
save.image(file=paste(output.file.path.root,".RData",sep=""))

cat("\nFinished!\n")
