#Toolkit version number
r.tool.name = "Data Input"
r.data.toolkit.version = "5.4"
tki.version = "tkifileversion15"
vsi.version="VSIfileversion15"
svn.version.data = "$ProjectRevision: SVN 13255 (2026-06-03) $"
svn.prj.name.data = "$ProjectName: /cerc_cerc/Products/Utilities/ModelEvaluationToolkit/Tags/ModelEvaluationToolkit_5_4_1 $"

#Suppress warnings that we don't want to worry the user with
options(warn=-1)
#Flag to enable extra output for development testing; should always be false in checked in versions
btesting = FALSE

#Find the install folder
argv = commandArgs(trailingOnly = FALSE)
install_dir = dirname(substring(argv[grep("--file=", argv)], 8))

#Set up input and output
settings.data.file = commandArgs(trailingOnly=TRUE)[[1]]

# Determine if settings file is TKI or VSI
settings.file.type=tolower(substr(settings.data.file, nchar(settings.data.file)-3, nchar(settings.data.file)))
# Check settings file type is valid and get output.file.path.root 
output.file.path.root = substr(settings.data.file,0,nchar(settings.data.file)-4)
if (settings.file.type ==".vsi"){
   settings.file.type="VSI"
} else if (settings.file.type==".tki"){
   settings.file.type="TKI"
} else {
   stop("The settings file is not a recognised file type. Settings files for the Data Input Tool must either be a TKI file or a VSI file.")
   quit("no")
}

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)) {
   if(settings.file.type=="VSI"){
      stop(paste("The .VSI file specified does not exist: ",settings.data.file, sep=""))
      quit("no")   
   } else if (settings.file.type=="TKI"){
      stop(paste("The .TKI file specified does not exist: ",settings.data.file, sep=""))
      quit("no")
   }
}

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

#If settings file is a VSI file, remove commented out lines
if(settings.file.type=="VSI"){
   VSI.data=readLines(settings.data.file)
   VSI.data=VSI.data[-which(startsWith(VSI.data,"#"))]  
   tki.data=VSI.data
} else {
   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 TKI/VSI file is the correct version 
if (settings.file.type=="TKI"){
   tki.version.in = as.character(tki.data[1])
   if (tki.version.in != tki.version) {
      stop(paste("The TKI settings file version (",tki.version.in,") does not match the expected version (",tki.version,")", sep=""))
      quit("no")   
   }   
} else if (settings.file.type =="VSI"){
   vsi.version.in = as.character(tki.data[1])
   if (vsi.version.in != vsi.version) {
      stop(paste("The VSI settings file version (",vsi.version.in,") does not match the expected version (",vsi.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")
   }
}

#MODELLED DATA
netCDF.formats = c("ADMS output netCDF","IOAPI","WRF-Chem","EMEP","CAMS Regional ADS","CAMS Regional Reanalysis","Generic netCDF","CHIMERE")
nMod = which(tki.data=="Concentration data")
import.model.name = tki.data[nMod+2]
pollutants.mod.file = tki.data[nMod+3]
bSpeciesMapping=tki.data[nMod+4]=="Y"
species.mapping.file=tki.data[nMod+5]
mod.format=tki.data[nMod+6]
mod.time.diff=as.numeric(tki.data[nMod+7])
cmaq.mcip.file=tki.data[nMod+8]
chimere.cco.file=tki.data[nMod+9]
bfile.mod=tki.data[nMod+10]=="file"
bdir.mod=tki.data[nMod+10]=="dir"
stations.interp=tki.data[nMod+11]=="Interpolation"
file.path.mod=tki.data[nMod+12]
dir.path.mod=tki.data[nMod+13]
bUseTemplatePath = tki.data[nMod+14]=="Y"
template.path=tki.data[nMod+15]
start.datetime=tki.data[nMod+16]
finish.datetime=tki.data[nMod+17]
data.separator.mod=tki.data[nMod+18]
missing.data.mod=tki.data[nMod+19]
#NETCDF option
if (mod.format=="CMAQ" | mod.format=="CAMX"){ #If modelled data format is CAMx or CAMQ, change mod.format to I/O API for processing
   mod.format="IOAPI"
}
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
#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
#Defra Background Maps option
bmod.dbm = mod.format=="Defra Background Maps"
bmod.dbm.select.file=bfile.mod
mod.dbm.file=file.path.mod
bmod.dbm.select.dir=bdir.mod
mod.dbm.dir=dir.path.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 (tolower(import.model.name) == "mod"){
   stop("Modelled data label must not be 'mod'. Please change the modelled data label.")
   quit("no") 
}
if(grepl("EnsembleMedian",import.model.name,ignore.case=TRUE)){
   stop("Modelled data label must not be 'EnsembleMedian'.")
   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 folder selected.")
   quit("no")
}
if(bSpeciesMapping){
   if (species.mapping.file==""){
      stop("Species mapping has been selected, but no species mapping file has been provided.")
      quit("no")
   }
   if (!file.exists(species.mapping.file)){
      stop("Species mapping file cannot be found. Check file path provided.")
      quit("no")
   }
}
if (bUseTemplatePath) {
   if (template.path==""){
      stop("Template file path has been selected, but no template file path has been provided.")
      quit("no")
   }
   if (start.datetime==""){
      stop("A template file path is being used, but no start date for the modelled dataset has been provided.")
      quit("no")
   }
   if (finish.datetime==""){
      stop("A template file path is being used, but no finish date for the modelled dataset has been provided.")
      quit("no")
   }
   if (nchar(start.datetime)!=10){
      stop("The start date for the modelled dataset has been provided in the wrong format; the correct format is YYYYMMDDHH.")
      quit("no")
   }
   if (nchar(finish.datetime)!=10){
      stop("The finish date for the modelled dataset has been provided in the wrong format; the correct format is YYYYMMDDHH.")
      quit("no")
   }
}
if (mod.format=="CHIMERE" & chimere.cco.file==""){
   stop("You have selected CHIMERE netCDF as the modelled data format but have not supplied a cco file with the required projection parameters.\n")
   quit("no")
}
if (mod.time.diff > 26 | mod.time.diff < -26){
   stop("Modelled data time difference must be between -26 and 26 \n")
   quit("no")
}

#OBSERVED DATA
nObs=which(tki.data=="Observed data")
benter.obs.data = tki.data[nObs+1]=="Y"
bobs.csv= tki.data[nObs+2]=="Y"
pollutants.obs.file = tki.data[nObs+3]
bobs.csv.select.file=tki.data[nObs+4]=="file"
bobs.csv.select.dir =tki.data[nObs+4]=="dir"
obs.csv.file=tki.data[nObs+5]
obs.csv.dir =tki.data[nObs+6]
obs.csv.separator=tki.data[nObs+7]
obs.csv.missing.data =tki.data[nObs+8]

##Meteorological data
nMet=which(tki.data=="Meteorological data")
bMetData=tki.data[nMet+1]=="Y"
MetDataFormat=tki.data[nMet+2]
MetDataFilePath=tki.data[nMet+3]
bMetOutputFile=(MetDataFormat %in% c("ADMSmetnc","wrf")) & tki.data[nMet+4]=="Y"

#Error handling:
if (bNewData) {
   if (!benter.obs.data) {
      stop("Observed data must be entered when a new workspace is being created.")
      quit("no")
   }
   if (bobs.csv) {
      if (pollutants.obs.file=="") {
         stop("No observed data pollutants file selected.")
         quit("no")
      }
      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 folder selected.")
         quit("no")
      }
   }
   if (bMetData) {
      if (MetDataFilePath==""){
         stop("Met data has been selected, but no met file path has been provided.")
         quit("no")
      }
   } else {
      cat("Warning: No met data has been selected, so plots requiring met data (e.g. polar plots) will not be available.\n")
   }
} 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")
   }
   if (bMetData) {
      stop("No met 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) {
  #Function to return the part of an averaging time
  avg_time_num = function (avg.time) {
    #Split into number and string
    split.avg.time=strsplit(avg.time," ")[[1]]
    if (length(split.avg.time)==1) {
      avg.time.num = 1
    } else if (length(split.avg.time)==2) {
      avg.time.num = as.numeric(split.avg.time[1])
    } else {
      stop(paste("Invalid averaging time (",avg.time,").",sep=""), call. = FALSE)
      quit("no")
    }
    #Only continue if the numeric part is actually numeric
    if (is.na(avg.time.num)) {
      stop(paste("Invalid averaging time (",avg.time,"). The optional numerical specifier must be numeric, not text.",sep=""), call. = FALSE)
      quit("no")
    }
    #Only continue if the number is an integer (the rest of the code assumes this)
    if (avg.time.num %% 1 != 0) {
      stop(paste("Invalid averaging time (",avg.time,"). The optional numerical specifier must be a whole number.",sep=""), call. = FALSE)
      quit("no")
    }
    return(avg.time.num)
  }
  
  #Function to return the string part of an averaging time
  avg_time_string = function (avg.time) {
    #Split into number and string
    split.avg.time=strsplit(avg.time," ")[[1]]
    if (length(split.avg.time)==1) {
      avg.time.string = split.avg.time
    } else if (length(split.avg.time)==2) {
      avg.time.string = split.avg.time[2]
    } else {
      stop(paste("Invalid averaging time (",avg.time,").",sep=""), call. = FALSE)
      quit("no")
    }
    #If averaging time ends in 's', remove the s so it will work with timeAverage()
    if (substring(avg.time.string,nchar(avg.time.string),nchar(avg.time.string)) == "s"){
      avg.time.string=substring(avg.time.string,1,nchar(avg.time.string)-1)
    }
    return(avg.time.string)
  }
  
  #Function to return the averaging time in terms of hours
  avg_time_hours = function (avg.time) {
    #Split into number and string
    avg.time.num = avg_time_num(avg.time)
    avg.time.string = avg_time_string(avg.time)
    #Find the number of hours in the string period
    if (avg.time.string == "hour") {
      avg_time_hours = 1
    } else if (avg.time.string == "day") {
      avg_time_hours = 24
    } else if (avg.time.string == "week") {
      avg_time_hours = 7*24
    } else if (avg.time.string == "month") {
      avg_time_hours = NA
    } else if (avg.time.string == "year") {
      avg_time_hours = 365*24
    } else {
      cat("Invalid averaging time. \n")
      return(NA)
    }
    #Calculate total number of hours in the averaging period
    avg_time_hours = avg.time.num*avg_time_hours
    return(avg_time_hours)
  }
  
  #Function to check for valid averaging times and statistics.
  #Used for checking 4 files: pollutants file, obs pollutants file, network pollutants file and mod pollutants file
  check_and_reformat_averaging_times_and_statistics = function(avg.time, statistic, polnames, string.descriptor) {
    #Say which file is being checked, in case of un-labelled errors
    cat(gsub("\\s{2,}", " ", paste("Checking averaging time and statistic settings in the ",string.descriptor," pollutants definition file...\n")))
    #First a sanity check to make sure that all vectors have the same length
    if (length(avg.time) != length(statistic) | length(avg.time) != length(polnames)) {
      #Coding error - shouldn't happen
      stop("Internal error in check_and_reformat_averaging_times_and_statistics", call. = FALSE)
      quit("no")
    }
    #Make sure averaging time and statistic are lower case character vectors with extraneous spaces removed
    avg.time = tolower(trimws(as.character(avg.time)))
    statistic = tolower(trimws(as.character(statistic)))
    #Check averaging times and statistics are valid
    for (i in 1:length(avg.time)) {
      #First check it's not numeric (i.e. old format)
      if(!is.na(as.numeric(avg.time[i]))){
        stop(gsub("\\s{2,}", " ", paste("The averaging time setting for pollutant ", polnames[i],
                              " in the ",string.descriptor," pollutants definition file is numeric (",avg.time[i],
                              ").\nThe setting must be one word (options are ",paste(allowed.avg.time.values,collapse = ", "),
                              ") with an optional numerical specifier e.g. '1 hour','7 days'.\n", sep="")), call. = FALSE)
        quit("no")
      }
      #Next check that one of the allowed averaging time strings has been used
      if (!(avg_time_string(avg.time[i])) %in% allowed.avg.time.values) {
        stop(gsub("\\s{2,}", " ", paste("The averaging time setting for pollutant ", polnames[i],
                              " in the ",string.descriptor," pollutants definition file is not an allowed setting (",
                              avg.time[i],").\nThe setting must be one word (options are ",paste(allowed.avg.time.values,collapse = ", "),
                              ") with an optional numerical specifier e.g. '1 hour','7 days'.\n", sep="")), call. = FALSE)
        quit("no")
      }
      #Finally, check that one of the allowed averaging statistics has been used
      if (!(statistic[i] %in% allowed.statistics)) {
        stop(gsub("\\s{2,}", " ", paste("The averaging statistic for pollutant ", polnames[i],
                              " in the ",string.descriptor," pollutants definition file is not an allowed setting (",
                              statistic[i],").\nThe setting must contain one of: ", paste(allowed.statistics, collapse=", "),".\n", sep="")), call. = FALSE)
        quit("no")
      }
      #Re-format the averaging time to remove trailing 's'
      avg.time[i]=paste(avg_time_num(avg.time[i]),avg_time_string(avg.time[i]))
      #Special case: replace 24-hour averaging time with daily averaging time
      if (avg.time[i]=="24 hour") avg.time[i]="1 day"
      #Special case: change 'week' to '7 day', because timeAverage() treats these differently, and 7 days is less likely to lead to inconsistencies
      if (avg_time_string(avg.time[i]) == "week") {
        new.avg.time = 7*avg_time_num(avg.time[i])
        avg.time[i]=paste(new.avg.time,"day")
      }
    }
    return(data.frame(avg.time,statistic))
  }
  
  #For an input averaging time and statistic to be compatible with the output averaging time and statistic, they must be:
  #a.	The same; or
  #b.	The input avg time is hourly mean; or
  #c.	Both means, with the input averaging time a factor of the output averaging time.
  #   (month is handled as a special case, because it's only compatible with month or year)
  avg_compatibility_check = function (pollutant.name,in.descriptor,in.avg.time,in.statistic,
                                           out.descriptor,out.avg.time,out.statistic) {
    #Checking compatibility of averaging times and statistics 
    bCompatible = TRUE
    #a. Are they the same?
    if ((out.avg.time==in.avg.time) & (out.statistic==in.statistic)) {
      #Yes, so that's fine
    } else {
      #No, they are not the same
      #b. Is the input hourly mean?
      if ((in.avg.time=="1 hour") & (in.statistic=="mean")) {
        #Yes, so that's fine, but don't allow rolling means to be monthly, would cause a crash in rollingMean function
        if ((out.statistic=="rolling mean") & (avg_time_string(out.avg.time)=="month")) {
          bCompatible = FALSE
        }
      } else {
        #No, it is not hourly mean, but might still be compatible, depending...
        bCompatible = FALSE
        #c1. Are they both means?
        if (!((out.statistic=="mean") & (in.statistic=="mean"))) {
          #No, they are not both means.
          #The compatibility criteria are not met
        } else if (avg_time_string(out.avg.time)=="month" | avg_time_string(in.avg.time)=="month") {
          #c2. They are both means, but is the input averaging time a factor of the output averaging time?
          #Month is a special case, because the number of hours can't be calculated
          #Hours, weeks and days are not compatible with month
          #Month is only compatible with month and year
          if (avg_time_string(out.avg.time)==avg_time_string(in.avg.time)) {
            if ((avg_time_num(out.avg.time) %% avg_time_num(in.avg.time))==0) {
              #Both month, and the input number of months is a factor of the output number of months
              bCompatible = TRUE
            }
          } else if (avg_time_string(out.avg.time)=="year") {
            #Output averaging time is years, input avg time is month, which is ok as long as the input number of months is a factor of the output number of months
            if (((12*avg_time_num(out.avg.time)) %% avg_time_num(in.avg.time))==0) {
              bCompatible = TRUE
            }
          } else if (in.avg.time=="1 day") {
            #Input is 1 day mean, output is n-monthly mean, so that's fine
            bCompatible = TRUE
          }
        } else if ((avg_time_hours(out.avg.time) %% avg_time_hours(in.avg.time))==0) {
          bCompatible = TRUE
        }
      }
    }
    if (!bCompatible) {
      #The compatibility criteria are not met, therefore stop with an error
      cat(paste0(toupper(pollutant.name)," ",in.descriptor," averaging time and statistic (",in.avg.time," ",in.statistic,") is incompatible ",
                 "with the ",toupper(pollutant.name)," ",out.descriptor," averaging time and statistic (",out.avg.time," ",out.statistic,").\n"))
    }
    return(bCompatible)
  }
  
  #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))
   }
   
   # Function cerc.OSGB.to.WGS84
   # This function will convert from OSGB to WGS84 coordinate systems 
   # Requires packages:proj4, rnrfa, sf
   # Requires arguments: OSGB x coordinates as a vector (x) and OSGB y coordinates as a vector (y)
   # Vectors x and y must be the same length i.e. they must be every x and y combination not defining grid values
   # Use the expand.grid function if they are defining grid values to get all points 
   cerc.OSGB.to.WGS84 = function (x,y) {
      if (length(x)!=length(y)) {
         stop(paste("X and Y must be the same length.")) # Check that x and y are the same length
      }
      OSGB.xy=as.data.frame(cbind(x,y))  # create dataframe of OSGB coordiantes
      names(OSGB.xy)=c("x.coords","y.coords")
      # st_as_sf is used to convert the coordinates to an sf (special features) object and takes the 
      # argument crs = 27700 to reference OSGB system
      OSGB.xy=st_as_sf(OSGB.xy,coords=c("x.coords","y.coords"), crs=27700)
      # st_transform to convert from OSGB to WGS84, crs=4326 for WGS84
      WGS84.ll=st_transform(OSGB.xy,crs=4326)
      # extract coordinates 
      WGS84.ll=as.data.frame(st_coordinates(WGS84.ll))
      return(WGS84.ll)
   }
   
   #Function WGS84.to.OSGB 
   # This function will convert from WGS84 to OSGB coordinate systems
   # Requires packages:proj4, rnrfa, sf
   # Requires arguments: Latitudes as a vector (lat) and longitudes as a vector (lon)
   # lat and lon must be the same length, use expand.grid function if required 
   cerc.WGS84.to.OSGB = function (lat,lon){
      if (length(lat)!=length(lon)) {
         stop(paste("Latitude and Longitude must be the same length.")) # check lat and lon are same length
      }
      WGS84.ll=as.data.frame(cbind(lat,lon)) # Create dataframe of latitude longitude data 
      names(WGS84.ll)=c("latitude","longitude")
      # st_as_sf is used to convert the coordinates to an sf (special features) object and takes the 
      # argument crs = 4326 to reference WGS84 system
      WGS84.ll=st_as_sf(WGS84.ll,coords=c("longitude","latitude"),crs=4326)
      # st_transform to convert from WGS84 to OSGB. crs=27700 for OSGB
      OSGB.xy=st_transform(WGS84.ll,crs=27700)
      # extract coordinates 
      OSGB.xy=as.data.frame(st_coordinates(OSGB.xy))
      return(OSGB.xy)
   }
   
   #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)
      #Change all column headings to lower case
      names(file.data) = tolower(names(file.data))
      #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(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)
   }
   # Open and read ADMS output netCDF file 
   read.ADMSoutput.nc=function(mod.ncdf.file){
      
      #Open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
      
      #Check that nc file has opened and if not return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("ADMS output netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
      
      #check that the file is an ADMS-Urban or ADMS-Local netCDF
      model.att=ncatt_get(modelled.nc,0,"Model_Name") #check if Model_Name attribute exists
      if(model.att$hasatt){model=model.att$value} else {
         stop(paste("The netCDF file is not in the expected format for an ADMS-Urban or ADMS-Local netCDF file, cannot find the global attribute Model_Name."))
         quit("no")
      }
      if (!((model=="ADMS-Urban")|(model=="ADMS-Local"))){ #check Model_Name attribute is ADMS-Urban
         stop(paste("The ADMS output netCDF file Model_Name global attribute is not ADMS-Urban or ADMS-Local."))
         quit("no")
      }
      
      #Check that all the pollutants have the save averaging time
      if (length(unique(pollutants$mod.avg.time)) != 1){
         stop(paste("Check modelled data pollutant definitions file.\n mod.avg.time must be the same for all pollutants when using ADMS output netCDF data format"))
         quit("no")
      }
      #Converting modelled averaging time into hours
      for (j in 1:length(pollutants$pollutant)){
         pollutants$mod.avg.time.hours[j]=avg_time_hours(pollutants$mod.avg.time[j])
         #Checking the number of hours is valid
         if (is.na(pollutants$mod.avg.time.hours[j])){
            stop(paste("Unable to calculate modelled data averaging time in hours. Check averaging time in modelled data pollutant definitions file.\n"))
            quit("no")
         }
      }
      #Defining dataset variables
      if (pollutants$mod.avg.time.hours[1]==1){ #Only checking the first element because we have already confirmed that all averaging times are the same
         data.names = "DatasetNames"
         data.set="Dataset"
      } else{
         data.names = "Ave_DatasetNames"
         data.set="Ave_Dataset"
      }
      
      #Check required variables are present in NC file
      var.names=names(modelled.nc$var) #extract names of all variables in nc file
      missing.var=""
      if (!("PointX_XYZ" %in% var.names)) {missing.var = paste(missing.var,"PointX_XYZ")} else {grid.xcoords = as.vector(ncvar_get(modelled.nc,"PointX_XYZ"))}
      if (!("PointY_XYZ" %in% var.names)) {missing.var = paste(missing.var,"PointY_XYZ")} else {grid.ycoords = as.vector(ncvar_get(modelled.nc,"PointY_XYZ"))}
      if (!("PointZ_XYZ" %in% var.names)) {missing.var = paste(missing.var,"PointZ_XYZ")} else {grid.zcoords = as.vector(ncvar_get(modelled.nc, "PointZ_XYZ"))}
      if (!("PointName_XYZ" %in% var.names)) {missing.var = paste(missing.var,"PointName_XYZ")} else{point.name=as.vector(ncvar_get(modelled.nc,"PointName_XYZ"))} 
      if (!("Met_Line" %in% var.names)){missing.var = paste(missing.var,"Met_Line")}else{met.line=ncvar_get(modelled.nc,"Met_Line")}
      if (!("Output_Points_Type" %in% var.names)){missing.var = paste(missing.var,"Output_Points_Type")} else{ output.point.types=ncvar_get(modelled.nc,"Output_Points_Type")}
      if (!("Number_Of_Output_Points_Of_Each_Type" %in% var.names)){missing.var = paste(missing.var,"Number_Of_Output_Points_Of_Each_Type")} else{ number.output.points.by.type=ncvar_get(modelled.nc,"Number_Of_Output_Points_Of_Each_Type")}
      if (!(data.names %in% var.names)) {missing.var = paste(missing.var,data.names)} else{dataset.names =as.vector(ncvar_get(modelled.nc,data.names))}
      #Check if any of the variables are missing and if so write them to the log 
      if (nchar(missing.var)!=0) {
         stop(paste("ADMS output netCDF data is not in the expected format; the following variables are missing:",missing.var))
         quit("no")
      }
      
      #check required dimensions are present in NC file
      dim.names=names(modelled.nc$dim)
      missing.dims=""
      if (!("nPoints_XYZ" %in% dim.names)) {missing.dims = paste(missing.dims,"nPoints_XYZ")} else {num.points=modelled.nc$dim$nPoints_XYZ$len}
      if (!("nMetLines" %in% dim.names)) {missing.dims = paste(missing.dims,"nMetLines")} else {num.timesteps=modelled.nc$dim$nMetLines$len}
      if (nchar(missing.dims)!=0){
         stop(paste("ADMS output netCDF is not in the expected format; the following dimensions are missing:",missing.dims))
         quit("no")
      }
      
      if (pollutants$mod.avg.time.hours[1]!=1){ #Only checking the first element because we have already confirmed that all averaging times are the same
         num.timesteps=1 #Setting num.timesteps to 1 because period averaged data only has one time step of data per station
      }
      
      #Checks on data output point types
      #Check some specified points are available 
      if ((number.output.points.by.type[which(output.point.types=="Specified Points")])==0){
         stop(paste("ADMS output netCDF file is not in the format expected, and does not contain specified point output point types."))
         quit("no")
      }
      #If other types of points are available, send warning message to log that only specified points are used 
      if (any(!((number.output.points.by.type[-(which(output.point.types=="Specified Points"))]))==0)){
         cat(paste("Warning: ADMS output netCDF file contains several output point types. Only specified point types will be used.\n"))
      }
      
      
      #Check there are some pollutant datasets in the nc file
      if (length(dataset.names)==0){
         stop(paste("ADMS output netCDF file does not contain pollutant datasets"))
         quit("no")
      }
    
      nc.pol.name=c()
      nc.pol.unit=c()
      nc.pol.avg.time=c()
      for (i in (1:length(dataset.names))){
         pol.var.name=paste(data.set,i,sep="")
         nc.pol.name[i]=ncatt_get(modelled.nc,pol.var.name,attname="Pollutant_Name")$value
         nc.pol.unit[i]=ncatt_get(modelled.nc,pol.var.name,attname="Units")$value
         nc.pol.avg.time[i]=ncatt_get(modelled.nc,pol.var.name,attname="Av_Time_s")$value
         if (nchar(nc.pol.name[i])==0){
            stop(paste("Pollutant name is missing in ADMS output netCDF file for ", dataset.names[i],sep=""))
            quit("no")
         }
         if (nchar(nc.pol.unit[i])==0){
            stop(paste("Pollutant unit is missing in ADMS output netCDF file for ",  dataset.names[i],sep=""))
            quit("no")
         }
         if (nchar(nc.pol.avg.time[i])==0){
            stop(paste("Pollutant average time is missing in ADMS output netCDF file for ",  dataset.names[i],sep=""))
            quit("no")
         }
         #Checking the dataset averaging time is the same as the modelled data averaging time
         nc.pol.avg.time.hours=nc.pol.avg.time[i]/3600 # Converting from seconds to hours
         if (nc.pol.avg.time.hours != pollutants$mod.avg.time.hours[1]){ #Checking the averaging times are the same. Only need to compare to the first element in modelled avearaging time because they are all the same
            stop(paste("The dataset",dataset.names[i],"does not match the averaging time in the modelled data pollutants definitions file.\nDataset: ",nc.pol.avg.time.hours,"hour, Modelled data pollutant definitions file: ",pollutants$mod.avg.time.hours[1],"hour"))
            quit("no")
         }
      }
      nc.pol.avg.time.hour=as.integer(as.numeric(nc.pol.avg.time)/3600)
      
      #extract the dataset names from the variables e.g. Dataset1
      dataset.var=var.names[which(str_detect(var.names,data.set))]
      dataset.var=dataset.var[-(which(dataset.var==data.names))] #remove DatasetNames from the list 
      
      #Identify which pollutants in the pollutant definition file are found in the netCDF file
      pollutants$bmod = FALSE #add bmod col to pollutans 
      pollutants$dataset.name = NA #add dataset name col to pollutants 
      pollutants$var.name=NA #add var.name col to pollutants 
      for (j in 1:length(pollutants$pollutant)) {
         for (i in 1: length(nc.pol.name)) {
            vi = nc.pol.name[i]
            if (tolower(pollutants$mod.alias[j]) == tolower(vi)) { #if pol exists in nc dataframe 
               pollutants$bmod[j] = TRUE #set bmod to TRUE if pol exists in nc file
               pollutants$dataset.name[j] = dataset.names[i] #record corresponding dataset name
               pollutants$var.name[j] = dataset.var[i] #record corresponding variable name 
               break
            }
         }
      }
      
      #stop if no defined pollutants exist in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in ADMS-Urban/AMDS-Local netCDF."))
         quit("no")
      }
      
      #remove the "Receptor" string from the station names 
      sp.point.name=gsub("Receptor, ","",point.name)
      #check that there are some receptors that are specified points, and not just gridded or intelligent points
      valid.points=which(!(startsWith((sp.point.name),"|I|") | startsWith((sp.point.name),"|G|")))
      if (length(valid.points)==0){
         stop(paste("No specified points found in netCDF file, named in expected format 'Receptor, [point_name]"))
      }
      
      
      #get time profile for the file, tz is assumed as "GMT" for MAQS-Health but caution should be taken for other ADMS output netCDF files
      #Currently no attribute for tz in ADMS comprehensive output netcdf 
      date.time=as.POSIXct(met.line,format="%Y_%j_%H",tz="GMT")
      
      bDoneFirstPol=FALSE
      for(j in 1:length(pollutants$pollutant)){
         #skip pollutant if not in netCDF file 
         if (!pollutants$bmod[j]) next
         for (s in 1:num.points){
           pollutant.data=ncvar_get(modelled.nc, pollutants$var.name[j], start=c(s,1,1),count=c(1,1,num.timesteps))
           onepol.onestation.alltimesteps=data.frame(date=date.time, station=sp.point.name[s], mod=pollutant.data)
           if (s==1){
             onepol.allstations.alltimesteps=onepol.onestation.alltimesteps
           } else {
             onepol.allstations.alltimesteps=rbind(onepol.allstations.alltimesteps, onepol.onestation.alltimesteps)
           }
         }
         names(onepol.allstations.alltimesteps)[3]=pollutants$pollutant[j]
         if(!bDoneFirstPol){
            nc.data=onepol.allstations.alltimesteps
            bDoneFirstPol=TRUE
         }else{
            onepol.allstations.alltimesteps=subset(onepol.allstations.alltimesteps, select=pollutants$pollutant[j])
            nc.data=cbind(nc.data,onepol.allstations.alltimesteps)
         }
      }
      
      #convert station from factor to character
      nc.data$station=as.character(nc.data$station)
      #remove gridded and intelligent gridded points
      if(any(((startsWith(nc.data$station,"|I|"))|startsWith(nc.data$station,"|G|")))){
         nc.data=nc.data[-(which((startsWith(nc.data$station,"|I|"))|startsWith(nc.data$station,"|G|"))),]
      }
      
      #Replace -999 values with na
      missing.ADMSoutput.data=-999
      nc.data.temp=subset(nc.data, select=-c(date,station)) #remove date and station columns to leave just pollutant data columns
      nc.data.temp.cols=names(nc.data.temp)  
      #set all  -999 values to NA
      for (i in 1:length(nc.data.temp.cols)){
         nc.data.temp[[nc.data.temp.cols[i]]][which(nc.data.temp[[nc.data.temp.cols[i]]]==missing.ADMSoutput.data)]=NA #set all  -999 values to NA  
      }
      nc.data=cbind(subset(nc.data, select=c("date","station")),nc.data.temp) #join pollutant data columns and date and station columns back together 
      
      return(nc.data)
      nc_close(modelled.nc)
      
   }
   
   #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 folder
   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))
   }
   
   #Install and load all required R packages
   install.load.packages = function(packages.required) {
      #Get the CRAN mirror URL
      cat("Getting the CRAN mirror from CRANMirror.txt...\n")
      mirror.file = file.path(install_dir,"CRANMirror.txt")
      mirror=readLines(mirror.file)[1]
      mirror.data = tryCatch(url(mirror),error=function(e) return(NULL))
      if (is.null(mirror.data)) {
         stop("Invalid CRAN mirror web address. Please refer to the 'R Support' menu for a complete list of available CRAN mirrors.")
         quit("no")
      }
     
      cat("Setting the CRAN mirror...\n")
      r = getOption("repos")
      r["CRAN"] = mirror
      options(repos=r)
     
      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"))
         }
         #Windows-only code to make a write-able libraries folder if necessary.
         #Detect operating system
         if(Sys.info()[['sysname']]=="Windows") {
            cat("Checking library folder write access before installing packages...\n")
            #Find out the default library path
            lib.path = .libPaths()[1]
            #Do we have write access?
            if (file.access(lib.path,2)==0) {
               #We have write access. Happy Days! Don't need to do anything else.
            } else {
               #We don't have write access, so make a new libraries folder somewhere else
               new.lib.path=Sys.getenv("R_LIBS_USER")
               #Check if it already exists
               if (dir.exists(new.lib.path)) {
                  #It does exist. Do we have write access?
                  if (file.access(new.lib.path,2)==0) {
                     #We have write access. Happy Days! Don't need to do anything else.
                  } else {
                     #We don't have write access. Stop with an error.
                     stop(paste("Can't install R packages because the R libraries folder is read-only:",new.lib.path,
                          "\nWe recommend you uninstall R and re-install R in a folder where you have write access."))
                     quit("no")
                  }
               } else {
                  #It doesn't exist, so create the new folder
                  dir.create(new.lib.path,recursive = TRUE)
                  #Add the new path to the library tree, as the first in the list and therefore the default
                  .libPaths(new.lib.path)
               }
            }
         }
         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)}
      }
     
   }
   
   #function for merging modeled data from a ncdf folder
   #requires modelled.data (already created from first file) and nc.data, the data extracted from the next file 
   merge.nc.data.modelled.data=function(nc.data, modelled.data){
      pol.names = names(nc.data)[3:length(names(nc.data))]
      for (j in 1:length(pol.names)) {
         nc.data.onepol= subset(nc.data,select=c("date","station",pol.names[j]))
         if (any(names(modelled.data)==pol.names[j])) {
            #Pull out the existing pollutant column
            modelled.data.onepol = subset(modelled.data,select=c("date","station",pol.names[j]))
            #Remove NAs
            modelled.data.onepol = modelled.data.onepol[!is.na(modelled.data.onepol[,3]),]
            #Merge with the new data for this pollutant
            nc.data.onepol = merge(modelled.data.onepol,nc.data.onepol,all=TRUE)
            #Strip this pollutant out of the master data frame
            modelled.data = modelled.data[ , !(names(modelled.data) == pol.names[j])]
         }
         #Merge this pollutant with the master data frame
         modelled.data = merge(modelled.data,nc.data.onepol, all=TRUE)
      }
      return(modelled.data)
   }
   
   #functions required for netCDF modelled data formats: 
   #checks which pollutants from pollutats dataframe (using mod.alias) are present in ncdf file
   #adds two new columns to pol df, bmod (logical for if it is present in ncdf file) and var.name that stores the name in the ncdf file if present
   check.ncdf.pol=function(pollutants, modelled.nc){
      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 (tolower(pollutants$mod.alias[j]) == tolower(vi$name)) {
               pollutants$bmod[j] = TRUE
               pollutants$var.name[j] = vi$name
               break
            }
         }
      }
      return(pollutants)
   }
   
   #function for gridded data to find the grid pollutant value for the grid the station is located in 
   #requires packages: RANN, dplyr
   find.met.grid.value= function(stations,grid.xcoords, grid.ycoords,met.parameter.grid,num.timesteps,time.date,k,regular.grid){
      stations.grid.values=subset(stations, select=c("map.pos.x","map.pos.y"))
      #get a grid combination of x and y points
      #for reg grid need to expand all combinations of x and y
      #for ireg gid, can just bind x and y columns 
      if (regular.grid){
         grid.values=expand.grid(as.vector(grid.xcoords), as.vector(grid.ycoords))
      } else {
         grid.values=cbind(as.vector(grid.xcoords),as.vector(grid.ycoords))
      }
      #find the nearest grid point to each station. This gives an index of the nearest grid point
      near.points.index=nn2(grid.values,stations.grid.values,k=1)[[1]]
      stations.nn=cbind(stations,near.points.index)
      stations.nn$nearest.point.x=grid.values[stations.nn$near.points.index,1]#add col to stations with nearest grid x val
      stations.nn$nearest.point.y=grid.values[stations.nn$near.points.index,2]#add col to stations with nearest grid y val
      #For reg and ireg grid can assume that the x and y grid (expansion) follows the same structure as the pollutant conc
      #each x value is cycled through for a single y value, then the next y value is selected and each x cycled through again
      #Therefore can assume that the point index of the nearest point is the same as the pol conc index 
      stations.nn$met.value=(as.vector(met.parameter.grid))[stations.nn$near.points.index]
      met.allstations.onetimestep = data.frame(date= time.date,station=stations.nn$name,met.value=stations.nn$met.value)
      #Identify stations outside the modelled area and set pollutant concentration value to NA
      stations.outside=which(!(between(stations$map.pos.x,range(grid.xcoords)[1], range(grid.xcoords)[2])) |!(between(stations$map.pos.y,range(grid.ycoords)[1], range(grid.ycoords)[2])))
      if (length(stations.outside)>0){
         met.allstations.onetimestep$met.value[stations.outside]=NA
         if(all(is.na(met.allstations.onetimestep$met.value))){
            stop("Valid met data can't be calculated for your monitoring stations. Please check file domain and format.\n")
            quit("no")
         }
      }
      return(met.allstations.onetimestep)
   }
   
   #function for gridded data to find the grid pollutant value for the grid the station is located in 
   #requires packages: RANN, dplyr
   find.pol.grid.value= function(stations,grid.xcoords, grid.ycoords,pollutant.data,num.timesteps,time.date,k,regular.grid){
      stations.grid.values=subset(stations, select=c("map.pos.x","map.pos.y"))
      #get a grid combination of x and y points
      #for reg grid need to expand all combinations of x and y
      #for ireg gid, can just bind x and y columns 
      if (regular.grid){
         grid.values=expand.grid(as.vector(grid.xcoords), as.vector(grid.ycoords))
      } else {
         grid.values=cbind(as.vector(grid.xcoords),as.vector(grid.ycoords))
      }
      #find the nearest grid point to each station. This gives an index of the nearest grid point
      near.points.index=nn2(grid.values,stations.grid.values,k=1)[[1]]
      stations.nn=cbind(stations,near.points.index)
      stations.nn$nearest.point.x=grid.values[stations.nn$near.points.index,1]#add col to stations with nearest grid x val
      stations.nn$nearest.point.y=grid.values[stations.nn$near.points.index,2]#add col to stations with nearest grid y val
      #For reg and ireg grid can assume that the x and y grid (expansion) follows the same structure as the pollutant conc
      #each x value is cycled through for a single y value, then the next y value is selected and each x cycled through again
      #Thefore can assume that the point index of the nearest point is the same as the pol conc index 
      if (num.timesteps==1){
         stations.nn$pol.conc=(as.vector(pollutant.data))[stations.nn$near.points.index]
         
      } else {
         stations.nn$pol.conc=(as.vector(pollutant.data[,,k]))[stations.nn$near.points.index]
      }
      onepol.allstations.onetimestep = data.frame(date= time.date[k],station=stations.nn$name,mod=stations.nn$pol.conc)
      #Identify stations outside the modelled area and set pollutant concentration value to NA
      stations.outside=which(!(between(stations$map.pos.x,range(grid.xcoords)[1], range(grid.xcoords)[2])) |!(between(stations$map.pos.y,range(grid.ycoords)[1], range(grid.ycoords)[2])))
      if (length(stations.outside)>0){
         onepol.allstations.onetimestep$mod[stations.outside]=NA
      }
      return(onepol.allstations.onetimestep)
   }
   
 	#function for finding the station concentration via interpolation from gridded values to station location
 	#using akima function interpp for irregular grids, and field function interp.surface for regular grids
	interp.station.value=function (grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid, num.timesteps,time.date,k){
      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$map.pos.x,stations$map.pos.y)$z
         } else {
            interp.conc = interpp(as.vector(grid.xcoords),as.vector(grid.ycoords),as.vector(pollutant.data[,,k]),stations$map.pos.x,stations$map.pos.y)$z
         }
      }
      onepol.allstations.onetimestep = data.frame(date=time.date[k],station=stations$name,mod=interp.conc)
      return(onepol.allstations.onetimestep)
   }
   
	 
	#function for reading WRF-Chem netCDF files as modelled data format
	#function dependencies: stations data frame, pollutants data frame
	#stations.interp logical variable, cerc.lambert.cc.xy() function
	#packages required: ncdf4, GEOmap, fields, akima, dplyr, RANN
	read.WRFChem.nc=function(mod.ncdf.file){
	   
	   #Open nc file and write file name to log
	   modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
	   cat(paste("   Importing",mod.ncdf.file,"...\n"))
	   
	   #Check that nc file has opened and if not return error message to log 
	   if (is.null(modelled.nc)){
	      stop(paste("WRF-Chem data file cannot be opened, please check the file path and file format."))
	      quit("no")
	   }
	   
	   #Check to see if WRF-Chem projection is compatible, either LCC, PS or UTM
	   missing.atts =""
	   MAP_PROJ.att = ncatt_get(modelled.nc,0,attname="MAP_PROJ") #Extract map projection attribute 
	   if (MAP_PROJ.att$hasatt) {MAP_PROJ = MAP_PROJ.att$value} else {missing.atts = paste(missing.atts,"MAP_PROJ")}
	   #Stop if map proj attribute doesn't exist 
	   if (nchar(missing.atts)!=0) {
	      stop(paste("WRF-Chem data is not in the expected format; the following attributes are missing:",missing.atts))
	      quit("no")
	   }
	   #Check the map projection is either LCC (1), PS(2), or UTM(4)
	   if (!(MAP_PROJ %in% c(1,2,4)))  {
	      stop("WRF-Chem data grid projection must be in either Lambert Conformal Conic, Polar Stereographic or Universal Transverse Mercator.")
	      quit("no")
	   }
	   
	   #check required attributes are present and store values in parameters. 
	   #For PS STAND_LON TRUELAT1, POLE_LAT attributes are required. 
	   #For LCC TRUELAT1, TRUELAT2, STAND_LON, MOAD_CENT_LAT attributes are required.  
	   #For UTM STAND_LON is required which gives the UTM zone number.
	   if (MAP_PROJ == 1){
	      #check if LCC TRUELAT1, TRUELAT2, STAND_LON, MOAD_CEN_LAT are present
	      TRUELAT1.att = ncatt_get(modelled.nc,0,attname="TRUELAT1")
	      if (TRUELAT1.att$hasatt) {TRUELAT1 = TRUELAT1.att$value} else {missing.atts = paste(missing.atts,"TRUELAT1")}
	      TRUELAT2.att = ncatt_get(modelled.nc,0,attname="TRUELAT2")
	      if (TRUELAT2.att$hasatt) {TRUELAT2 = TRUELAT2.att$value} else {missing.atts = paste(missing.atts,"TRUELAT2")}
	      STAND_LON.att = ncatt_get(modelled.nc,0,attname="STAND_LON")
	      if (STAND_LON.att$hasatt) {STAND_LON = STAND_LON.att$value} else {missing.atts = paste(missing.atts,"STAND_LON")}
	      MOAD_CEN_LAT.att = ncatt_get(modelled.nc,0,attname="MOAD_CEN_LAT")
	      if (MOAD_CEN_LAT.att$hasatt) {MOAD_CEN_LAT = MOAD_CEN_LAT.att$value} else {missing.atts = paste(missing.atts,"MOAD_CEN_LAT")}
	      if (nchar(missing.atts)!=0) {
	         stop(paste("WRF-Chem data is not in the expected format; the following attributes are missing:",missing.atts))
	         quit("no")
	     }
	   } else if (MAP_PROJ ==2)  {
	     #check if PS, STAND_LON, TRUELAT1 and POLE_LAT present - needs testing
	      STAND_LON.att = ncatt_get(modelled.nc,0,attname="STAND_LON")
	      if (STAND_LON.att$hasatt) {STAND_LON = STAND_LON.att$value} else {missing.atts = paste(missing.atts,"STAND_LON")}
	      TRUELAT1.att = ncatt_get(modelled.nc,0,attname="TRUELAT1")
	      if (TRUELAT1.att$hasatt) {TRUELAT1 = TRUELAT1.att$value} else {missing.atts = paste(missing.atts,"TRUELAT1")}
	      POLE_LAT.att = ncatt_get(modelled.nc,0,attname="POLE_LAT")
	      if (POLE_LAT.att$hasatt) {POLE_LAT = POLE_LAT.att$value} else {missing.atts = paste(missing.atts,"POLE_LAT")}
	      if (nchar(missing.atts)!=0) {
	         stop(paste("WRF-Chem data is not in the expected format; the following attributes are missing:",missing.atts))
	         quit("no")
	      }
	   } else if (MAP_PROJ ==4){
	      #Check if UTM proj, present - needs testing
	      STAND_LON.att = ncatt_get(modelled.nc,0,attname="STAND_LON")
	      if (STAND_LON.att$hasatt) {STAND_LON = STAND_LON.att$value} else {missing.atts = paste(missing.atts,"STAND_LON")}
	      if (nchar(missing.atts)!=0) {
	         stop(paste("WRF-Chem data is not in the expected format; the following attributes are missing:",missing.atts))
	         quit("no")
	      }
	   }
	   
	   #Check for variables for grid and layers are present
	   #Extract and store XLAT XLONG (for creating grid) PH PHB and HGT (for calculating heights and identifying surface layer)
	   var.names=names(modelled.nc$var)
	   missing.vars=""
	   if (!("XLAT" %in% var.names)){missing.vars = paste(missing.vars,"XLAT")} else {XLAT=ncvar_get(modelled.nc, "XLAT")}
	   if (!("XLONG" %in% var.names)){missing.vars = paste(missing.vars,"XLONG")} else {XLONG=ncvar_get(modelled.nc, "XLONG")}
	   if (!("PH" %in% var.names)){missing.vars = paste(missing.vars,"PH")} else {PH=ncvar_get(modelled.nc, "PH")}
	   if (!("PHB" %in% var.names)){missing.vars = paste(missing.vars,"PHB")} else {PHB=ncvar_get(modelled.nc, "PHB")}
	   if (!("HGT" %in% var.names)){missing.vars = paste(missing.vars,"HGT")} else {HGT=ncvar_get(modelled.nc, "HGT")}
	   if (nchar(missing.vars) != 0) {
	      stop(paste("WRF-Chem data is not in the expected format; the following variables are missing:",missing.vars))
	      quit("no")
	   }
	   
	   #Check expected dimensions and set values of parameters
	   missing.dims = ""
	   dim.names = names(modelled.nc$dim)
	   if (length(which(dim.names =="Time"))==0){missing.dims = paste(missing.dims,"Time")} else {num.timesteps = modelled.nc$dim$Time$len}
	   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}
	   if (nchar(missing.dims) != 0) {
	      stop(paste("WRF-Chem data is not in the expected format; the following dimensions are missing:",missing.dims))
	      quit("no")
	   }
	   
	   
	   #Underlying grid for WRF-Chem is regular orthogonal so create map projection for conversion 
	   #How the map projection is created depends on whether the WRF-Chem file projection is LCC, PS or UTM
	   if (MAP_PROJ ==1 ){ #If map projection is LCC
	      #Create map projection using MOAD_CEN_LAT for latitude of origin, STAND_LON for longitude of origin and TRUELAT1 and TRUELAT2 for the standard parallels
	      #Type = 3 for LCC
	      map.proj=setPROJ(type=3, LAT0=MOAD_CEN_LAT,LON0=STAND_LON, LAT1 = TRUELAT1, LAT2 = TRUELAT2)
	      #Create the projected regular grid from the XLAT and XLONG variables
	      grid.xy = cerc.lambert.cc.xy(XLAT[,,1],XLONG[,,1],map.proj) #Extract 1 time [,,1](1 matrix slice) as the x and y grid points do not vary with time
	      #Extract x and y grid coordinates
	      grid.xcoords=grid.xy$x
	      grid.ycoords=grid.xy$y
	      #Project stations onto the projected grid using lat long from stations file and map projection created 
	      stations.xy=cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
	      #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
	      stations$map.pos.x = stations.xy$x
	      stations$map.pos.y = stations.xy$y  
	      allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
	      #Create plots to test the projection 
	      #test.plot.before projection 
	      #plot(x=as.vector(XLAT[,,1]), y=as.vector(XLONG[,,1]))
	      #points(x=stations$latitude, y=stations$longitude, col="red", cex=2)
	      #test.plot.after
	      #plot(x=as.vector(grid.xy$x),y=as.vector(grid.xy$y))
	      #points(x=stations.xy$x,y=stations.xy$y, col="red",cex=2)
	   } else if (MAP_PROJ ==2 ){ #If map projection is polar stereographic 
	      #Create map projection using STAND_LON for Longitude of Origin, LAT1 as the true parallel and POLE_LAT as Latitude of origin as it is polar stereographic 
	      #Type = 4 for PS
	      map.proj=setPROJ(type=4, LON0=STAND_LON, LAT1=TRUELAT1, LAT0 = POLE_LAT)
	      grid.xy = stereo.sphr.xy(XLAT[,,1],XLONG[,,1],map.proj)
	      #Extract x and y grid coordinates
	      grid.xcoords=grid.xy$x
	      grid.ycoords=grid.xy$y
	      #Project stations onto the projected grid using lat long from stations file and map projection created 
	      stations.xy=stereo.sphr.xy(stations$latitude,stations$longitude,map.proj)
	      #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
	      stations$map.pos.x = stations.xy$x
	      stations$map.pos.y = stations.xy$y  
	      allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
	      #Create plots to test the projection 
	      #test.plot.before projection 
	      #plot(x=as.vector(XLAT[,,1]), y=as.vector(XLONG[,,1]))
	      #points(x=stations$latitude, y=stations$longitude)
	      #test.plot.after
	      #plot(x=as.vector(grid.xy$x),y=as.vector(grid.xy$y))
	      #points(x=stations.xy$x,y=stations.xy$y, col="red")
	   } else if (MAP_PROJ == 4){ #If map projection is UTM
	      #Convert STAND_LON which gives UTM zone number to LON0 required for creating the projection
	      lon.0= 6.0 * (STAND_LON - 31.0) + 3.0
	      #Use Lon.0 as longitude of origin for projection
	      #Type = 7 for utm.wgs84
	      map.proj=setPROJ(type=7, LON0=lon.0)
	      #get projected grid xy coordinates from the latitude and longitude variables and map projection 
	      grid.xy = utm.wgs84.xy(XLAT[,,1],XLONG[,,1],map.proj)
	      #Extract x and y grid coordinates
	      grid.xcoords=grid.xy$x
	      grid.ycoords=grid.xy$y
	      #Project stations onto the projected grid using lat long from stations file and map projection created 
	      stations.xy=utm.wgs84.xy(stations$latitude,stations$longitude,map.proj)
	      #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
	      stations$map.pos.x = stations.xy$x
	      stations$map.pos.y = stations.xy$y  
	      allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
	     
	   }
	   
	   #Check for regular or irregular grid as this determines method of interpolation. 
	   num.dimensions.x = length(dim(grid.xcoords))
	   num.dimensions.y = length(dim(grid.ycoords))
	   #Check dimensions of x are equal to dimensions of y and stop if not
	   if (num.dimensions.x==num.dimensions.y) {
	      num.dimensions = num.dimensions.x
	   } else {
	      stop(paste("netCDF variables XLAT and XLONG do not have the same number of dimensions",sep=""))
	      quit("no")
	   } 
	   #If a regular grid, x and y have 1 dimension, if irregular x and y have 2 dimensions.
	   if (num.dimensions==1) {
	      regular.grid = TRUE
	   } else if(num.dimensions==2) {
	      regular.grid = FALSE
	   } else {
	      #If x and y have more than 2 dimensions, stop as this is not compatible 
	      stop(paste("netCDF variables XLAT and XLONG have ",num.dimensions," dimensions but only 1 or 2 are allowed",sep=""))
	      quit("no")
	   }
	   
	   #Get time profile for the file from the time variable and convert to POSIXct format
	   #GMT timezone assumed so compatible with the RML
	   time.date = as.POSIXct(ncvar_get(modelled.nc,"Times"), format="%Y-%m-%d_%H:%M:%S", tz="GMT")
	   
	   #Find the layer nearest the surface and check it is close enough to the surface to use  
	   PH_PHB=(PH[,,,1]+PHB[,,,1])/9.81 #only use one hour
	   heights=array(numeric(), dim(PH_PHB)) #create an empty array to use to store height values in
	   #Loop through and subtract HGT from PHB_PH and consider different dimensions of arrays
	   for (h in (1:dim(PH_PHB)[3])){
	      heights[,,h]=PH_PHB[,,h]-HGT[,,1]
	   }
	   #Calculate mean value for each height averaged over the modelled domain
	   mean.heights=colMeans(heights,dims=2)
	   surf.ref=which.min(mean.heights)
	   
	   #Check that min height is close to surface, 1m is used to be consistent with the MAQS-Health Coupled System 
	   if ((mean.heights[surf.ref]) > 1){
	      stop(paste("The lowest height in WRF-Chem netCDF file is",mean.heights[surf.ref],". This is not close enough to the surface to be considered as surface concentration.",sep=""))
	      quit("no")
	   }  
	   #If min layer is close enough to the surface, set surface layer as the layer required 
	   layer.required=surf.ref
	   
	   #Find out the pollutants that are in this nc file and save the variable names
	   pollutants=check.ncdf.pol(pollutants, modelled.nc)
	   
	   #Stop if no defined pollutants exist in netCDF file
	   if(!(any(pollutants$bmod == TRUE))){
	      stop(paste("No defined pollutants found in WRF-Chem netCDF file."))
	      quit("no")
	   }
	   
	   #Loop through the pollutants and extract the interpolated or grid value modeled data for each pollutant 
	   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
	      #Loop through every timestep in the modeled data file
	      for (k in 1:num.timesteps) {
	         #Extract all the surface layer data for the current pollutant
	         pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,layer.required,k),count=c(num.cols,num.rows,1,1))
	         #Calculating modelled station values by interpolation 
	         if(stations.interp){
	            onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
	         } 
	         #Calculating modelled station concentrations by grid value
	         if (!stations.interp){
	            onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
	         }
	         #create dataframe for that pollutant for all timesteps
	         if (k==1) { #If first timestep set as dataframe, 
	            onepol.allstations.alltimesteps = onepol.allstations.onetimestep
	         } else { #If later timestep merge with existing
	            onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
	         }
	      }
	     
	      #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) {
	         nc.data = onepol.allstations.alltimesteps
	         bDoneFirstPol = TRUE
	      } else {
	         nc.data = merge(nc.data,onepol.allstations.alltimesteps,all=TRUE)
	      }
	   } #end loop over pollutants
	   
	   #convert pollutant columns from character to numeric class
	   nc.data.headers=names(nc.data)[3:ncol(nc.data)]
	   for (i in (1:length(nc.data.headers))){
	      nc.data[,2+i]=as.numeric(nc.data[,2+i])
	   }
	   return(nc.data)
	   nc_close(modelled.nc)
	   
   }
   
   #function for reading EMEP netCDF files as modelled data format
   #function dependencies: stations data frame, pollutants data frame, mod.ncdf.file
   #stations.interp logical variable, check.ncdf.pol function, find.pol.grid.value function,interp.station.value function 
   #packages required: ncdf4, GEOmap, fileds, akima, dplyr, RANN
   read.EMEP.nc=function(mod.ncdf.file) {
      #open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
      
      #Check that nc file has opened and if not, stop and return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("EMEP netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
      
      #Check to ensure EMEP is a compatible projection (either PS or reg lat lon)
      #Check there is a global projection attribute 
      missing.atts = ""
      projection.att = ncatt_get(modelled.nc,0,attname="projection")
      if (projection.att$hasatt) {projection = projection.att$value} else {missing.atts = paste(missing.atts,"projection")}
      if (nchar(missing.atts)!=0) { #Stop if projection attribute cannot be found 
         stop(paste("EMEP data is not in the expected format; the following attributes are missing:",missing.atts))
         quit("no")
      }
      #Check that the projection is either Stereographic or "lon lat"
      if (!((projection == "Stereographic")|(projection == "lon lat"))) {
         stop("EMEP data must be in either regular latitude longitude or Polar Stereographic projection") # check value of projection for LCC may also want to add conical lambert conformal here too although this is not standard
         quit("no")
      }
      #Interpolation to station location is only allowed with regular gridded data
      #and if projection is polar stereographic then lat long grid will be irregular
      if ((projection == "Stereographic") & stations.interp) {
         stop(paste("Interpolation to station location is only permitted for EMEP 'lon lat' projection.\n",
                   "Select 'Grid value of station' instead for 'Stereographic' projection."))
         quit("no")
      }
      
      
      #Check the dimensions based on the projection 
      #For stereographic projections, time, j and i are required. lev and ilev not required as these are not used for pollutant conc dims 
      if (projection =="Stereographic") {
         missing.dims = ""
         dim.names = names(modelled.nc$dim)
         if (!("time" %in% dim.names)) {missing.dims = paste(missing.dims,"time")} else {num.timesteps = modelled.nc$dim$time$len}
         if (!("j"%in%dim.names)) {missing.dims = paste(missing.dims,"j")} else {num.cols = modelled.nc$dim$j$len}
         if (!("i"%in% dim.names)) {missing.dims = paste(missing.dims,"i")} else {num.rows = modelled.nc$dim$i$len}
         if (nchar(missing.dims) != 0) { #Stop if any dimensions missing 
            stop(paste("EMEP data is not in the expected format; the following dimensions are missing:",missing.dims))
            quit("no")
         }
      } else if (projection =="lon lat") { #for lon lat projections, time, lat and lon are required 
         missing.dims=""
         dim.names = names(modelled.nc$dim)
         if (!("time" %in% dim.names)) {missing.dims = paste(missing.dims,"time")} else {num.timesteps = modelled.nc$dim$time$len}
         if (!("lat" %in% dim.names)) {missing.dims = paste(missing.dims,"lat")} else {num.cols = modelled.nc$dim$lat$len}
         if (!("lon" %in% dim.names)) {missing.dims = paste(missing.dims,"lon")} else {num.rows = modelled.nc$dim$lon$len}
         if (nchar(missing.dims) != 0) { #Stop if any dimensions missing 
            stop(paste("EMEP data is not in the expected format; the following dimensions are missing:",missing.dims))
            quit("no")
         }
      }
      
      #Extract the grid and time values, have checked dimensions so variables should be present.
      grid.xcoords=ncvar_get(modelled.nc,"lon")
      grid.ycoords=ncvar_get(modelled.nc,"lat")
      time= ncvar_get(modelled.nc,"time")
      
      #Create allstations.location matrix/map.pos for interpolation/nearest point.
      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
      
      #check dimensions 
      num.dimensions.x = length(dim(grid.xcoords))
      num.dimensions.y = length(dim(grid.ycoords))
      if (num.dimensions.x==num.dimensions.y) { #check x and y have the same dimensions, if not stop. 
         num.dimensions = num.dimensions.x
      } else {
         stop(paste("netCDF variables for defining the grid (e.g. lat and lon or i and j) do not have the same number of dimensions",sep=""))
         quit("no")
      }         
      if (num.dimensions==1) { #If x/y only have 1 dimension, regular.grid=TRUE
         regular.grid = TRUE
      } else if(num.dimensions==2) {#If x/y have 2 dimensions, regular.grid=FALSE
         regular.grid = FALSE
      } else { #If x/y have more than 2 dimensions, stop and send error to log. Must have only 1 or 2 dims. 
         stop(paste("netCDF variables for defining the grid (e.g. lat and lon or i and j) have",num.dimensions," dimensions but only 1 or 2 are allowed",sep=""))
         quit("no")
      }

      #Get time profile for file, time has already been extracted but now needs to be converted to the correct format 
      time.unit.att=ncatt_get(modelled.nc, "time","units") #Extract the units for time to give origin date for jday
      if (time.unit.att$hasatt) {time.unit = time.unit.att$value} else {
         stop(paste("EMEP file does not contain origin date for the time to convert julian days to dates."))
      }
      orig.day=gsub("days since ","",time.unit) #remove string to get just the origin date
      time.date=as.POSIXct(86400*time,origin=orig.day,format="%Y-%m-%d %H:%M",tz="GMT") #convert to POSIXct format, using origin date for jday to cal date
      time.date=ceiling_date(time.date, unit="hours") #round to nearest hour, round up so hour ending 
      
      #Find out the pollutants that are in this nc file and save the variable names
      pollutants=check.ncdf.pol(pollutants,modelled.nc)
      #Check there are defined pollutants in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in EMEP netCDF file."))
         quit("no")
      }
      
      #see if EMEP file has a height dimension
      dim.names=names(modelled.nc$dim)
      if ("lev" %in% dim.names){
         #identify surface layer
         surf.ref=which.min(modelled.nc$dim$lev$vals)
      }
      
      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 data for the current pollutant
         for (k in 1:num.timesteps){ #Reading in one hour of data at a time from the EMEP file
            if (modelled.nc$var[[pollutants$var.name[j]]]$ndims==4){
               pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,surf.ref,k),count=c(num.rows,num.cols,1,1)) #Gets one hour of data at a time.
            } else if (modelled.nc$var[[pollutants$var.name[j]]]$ndims==3){
               pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,k),count=c(num.rows,num.cols,1))
            } else{
               stop(paste("Can only process EMEP files where variables have 3 or 4 dimensions."))
            }
            if (stations.interp){
               onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
            } else{ #If !stations.interp
               onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
            }
            if (k==1) {
               onepol.allstations.alltimesteps = onepol.allstations.onetimestep
            } else {
               onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
            }
         }
         
         #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) {
            nc.data = onepol.allstations.alltimesteps
            bDoneFirstPol = TRUE
         } else {
            onepol.allstations.alltimesteps=onepol.allstations.alltimesteps[,-c(1,2)]
            nc.data = cbind(nc.data,onepol.allstations.alltimesteps)
            names(nc.data)[which(names(nc.data)=="onepol.allstations.alltimesteps")]=pollutants$pollutant[j]
         }
      } #end loop over pollutants
      #Ensure date column is in the correct format
      nc.data$date=as.POSIXct(nc.data$date, format="%Y-%m-%d %H:%M", tz="GMT")#timezone GMT is assumed
      return(nc.data)
   }
   
    
   #requires custom functions: cerc.lambert.cc.xy, fmod, convert.time, check.ncdf.pol, find.pol.grid.value, interp.station.value
   #requires packages: ncdf4, GEOmap, fields, akima, dplyr, RANN
   read.IOAPI.nc=function(mod.ncdf.file){
       
      #Open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
       
      #Check that nc file has opened and if not return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("I/O API data file cannot be opened, please check the file path and file format."))
         quit("no")
      }
       
      #Check that the data are gridded and the projection is lambert conformal conic
      missing.atts = ""
      FTYPE.att = ncatt_get(modelled.nc,0,attname="FTYPE") #FTYPE attribute =1 if gridded data 
      if (FTYPE.att$hasatt) {FTYPE = FTYPE.att$value} else {missing.atts = paste(missing.atts,"FTYPE")} 
      GDTYP.att = ncatt_get(modelled.nc,0,attname="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) { #If FTYPE or GDTYP are missing, stop and send error to log
         stop(paste("I/O API file is not in the expected format; the following attributes are missing:",missing.atts))
         quit("no")
      }
      if (!(FTYPE == 1)) {
         stop("I/O API file data must be gridded")
         quit("no")
      }
      if (!(GDTYP %in% c(2,5,6))) { #GDTYP = 2 for LCC, 5 for UTM, 6 for PS
         stop("I/O API file data grid projection must be either Lambert Conformal Conic, Polar Stereographic or Universal Transverse Mercator")
         quit("no")
      }
       
      #Get all the required attributes can be found and use to store values
      # For LCC, P_ALP, P_BET, YCENT, XCENT are required
      # For PS YCENT and XCENT are required, P_ALP gives hemisphere
      # For UTM P_ALP gives UTM Zone 
      #XORIG, YORIG, XCELL, YCELL required for grid , SDATE, STIME, TSTEP for time profile
      missing.atts = ""
      XORIG.att = ncatt_get(modelled.nc,0,attname="XORIG") #used as base point of grid
      if (XORIG.att$hasatt) {XORIG = XORIG.att$value} else {missing.atts = paste(missing.atts,"XORIG")}
      YORIG.att = ncatt_get(modelled.nc,0,attname="YORIG") #used as base point of grid
      if (YORIG.att$hasatt) {YORIG = YORIG.att$value} else {missing.atts = paste(missing.atts,"YORIG")}
      XCELL.att = ncatt_get(modelled.nc,0,attname="XCELL") #used to define grid
      if (XCELL.att$hasatt) {XCELL = XCELL.att$value} else {missing.atts = paste(missing.atts,"XCELL")}
      YCELL.att = ncatt_get(modelled.nc,0,attname="YCELL") #use to define grid
      if (YCELL.att$hasatt) {YCELL = YCELL.att$value} else {missing.atts = paste(missing.atts,"YCELL")}
      Sigma.att = ncatt_get(modelled.nc,0,attname= "VGLVLS") #used to get sigma values to identify surface layer
      if (Sigma.att$hasatt) {sigma.values = Sigma.att$value} else {missing.atts = paste(missing.atts, "VGLVLS")}
      SDATE.att = ncatt_get(modelled.nc,0,attname="SDATE") #used for start date
      if (SDATE.att$hasatt) {SDATE = SDATE.att$value} else {missing.atts = paste(missing.atts,"SDATE")}
      STIME.att = ncatt_get(modelled.nc,0,attname="STIME") #used for start time
      if (STIME.att$hasatt) {STIME = STIME.att$value} else {missing.atts = paste(missing.atts,"STIME")}
      TSTEP.att = ncatt_get(modelled.nc,0,attname="TSTEP") #used to create time profile 
      if (TSTEP.att$hasatt) {TSTEP = TSTEP.att$value} else {missing.atts = paste(missing.atts,"TSTEP")}
      VGTYP.att = ncatt_get(modelled.nc,0,attname="VGTYP") #vertical co-ordinate type 
      if (VGTYP.att$hasatt) {VGTYP = VGTYP.att$value} else {missing.atts = paste(missing.atts,"VGTYP")}
      if (GDTYP ==2){ #LCC proj attributes
         XCENT.att = ncatt_get(modelled.nc,0,attname="XCENT") #used to define LON0
         if (XCENT.att$hasatt) {XCENT = XCENT.att$value} else {missing.atts = paste(missing.atts,"XCENT")}
         YCENT.att = ncatt_get(modelled.nc,0,attname="YCENT") #used to define LAT0
         if (YCENT.att$hasatt) {YCENT = YCENT.att$value} else {missing.atts = paste(missing.atts,"YCENT")}
         P_ALP.att = ncatt_get(modelled.nc,0,attname="P_ALP") #Used to define TRUELAT1
         if (P_ALP.att$hasatt) {P_ALP = P_ALP.att$value} else {missing.atts = paste(missing.atts,"P_ALP")}
         P_BET.att = ncatt_get(modelled.nc,0,attname="P_BET") #Used to define TRUELAT2
         if (P_BET.att$hasatt) {P_BET = P_BET.att$value} else {missing.atts = paste(missing.atts,"P_BET")}
      } else if (GDTYP ==6){ #PS proj attributes
         XCENT.att = ncatt_get(modelled.nc,0,attname="XCENT") #used to define LON0
         if (XCENT.att$hasatt) {XCENT = XCENT.att$value} else {missing.atts = paste(missing.atts,"XCENT")}
         YCENT.att = ncatt_get(modelled.nc,0,attname="YCENT") #used to define LAT0
         if (YCENT.att$hasatt) {YCENT = YCENT.att$value} else {missing.atts = paste(missing.atts,"YCENT")}
         P_ALP.att = ncatt_get(modelled.nc,0,attname="P_ALP") #Used to define TRUELAT1 (90 or -90)
         if (P_ALP.att$hasatt) {P_ALP = P_ALP.att$value} else {missing.atts = paste(missing.atts,"P_ALP")}
      } else if (GDTYP ==5){ #UTM proj attributes
         P_ALP.att = ncatt_get(modelled.nc,0,attname="P_ALP") #Used to define UTM Zone number 
         if (P_ALP.att$hasatt) {P_ALP = P_ALP.att$value} else {missing.atts = paste(missing.atts,"P_ALP")}
      }
      #Stop if any expected attributes are missing
      if (nchar(missing.atts)!=0) {
         stop(paste("I/O API 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("I/O API data is not in the expected format; the following dimensions are missing:",missing.dims))
         quit("no")
      }
       
      #Calculate the I/O API grid coordinates (in map projection)
      #Create x coords from XORIG using num.cols and XCELL
      grid.xcoords = seq(from=(XORIG+0.5*XCELL),to=(XORIG+0.5*XCELL+(num.cols-1)*XCELL),by=XCELL)
      #Create y coords from YORIG using num.rows and YCELL
      grid.ycoords = seq(from=(YORIG+0.5*YCELL),to=(YORIG+0.5*YCELL+(num.rows-1)*YCELL),by=YCELL)
      regular.grid = TRUE # regular.grid is always true for I/O API files when created as above
      #Set the map projection (uses GEOmap package), currently assuming lambert cc (type=3)
       
      #Project the stations onto the grid
      if (GDTYP == 2){
         #Set projection for LCC I/O API files, type = 3
         map.proj = setPROJ(type = 3, LAT0 = YCENT, LON0 = XCENT, LAT1 = P_ALP, LAT2 = P_BET)
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
      } else if (GDTYP == 5){
         #Set map projection for UTM I/O API files , type = 7
         #Convert P_ALP which gives UTM zone number to LON0 required for creating the projection
         lon.0= 6.0 * (P_ALP - 31.0) + 3.0
         map.proj = setPROJ(type = 7, LAT0 = YCENT, LON0 = XCENT, LAT1 = P_ALP)
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=utm.wgs84.xy(stations$latitude,stations$longitude,map.proj)
      } else if (GDTYP == 6){
         #Set map projection for PS I/O API files, type = 4  
         map.proj = setPROJ(type = 4, LAT0 = YCENT, LON0 = XCENT, LAT1 = P_ALP)
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=stereo.sphr.xy(stations$latitude,stations$longitude,map.proj)
      }
       
      #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
      stations$map.pos.x = stations.xy$x
      stations$map.pos.y = stations.xy$y  
      allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
       
      #Extract time profile from I/O API file 
      #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 I/O API file
      time.date = seq(from = start.date, to = (start.date+(num.timesteps-1)*timestep),by=timestep)              
      
      # Identify surface layer
      if (num.layers == 1) {
        surf.ref = 1
      } else if (VGTYP %in% c(1,2,7,8)) { # vertical co-ordinates are sigma-P values 
        surf.ref = which.max(sigma.values) # layer nearest ground has max sigma value
      } else if (VGTYP %in% c(3,4,5)) {
        stop(paste0("The vertical co-ordinate type (VGTYP) is ",VGTYP,". Only the following types are supported: 1, 2, 6, 7, 8 and -9999."))
        quit("no")
      } else if (VGTYP == 6) { # vertical co-ordinates are height above ground
        surf.ref = which.min(sigma.values) # layer nearest ground has min value
      } else if (VGTYP == -9999) { # CMAQ version >= 5.3 so MCIP file required
        # Check path to MCIP file
        if (cmaq.mcip.file == "") {
          stop("The MCIP file path has not been supplied.")
          quit("no")
        }
        cmaq.mcip.file.ext = tolower(substr(cmaq.mcip.file, nchar(cmaq.mcip.file)-2, nchar(cmaq.mcip.file))) 
        if (cmaq.mcip.file.ext != ".nc") {
          stop(paste0("The MCIP file extension is invalid. It should be \".nc\""))
          quit("no")
        }
        if (!file.exists(cmaq.mcip.file)) {
          stop("The MCIP file cannot be found. Check the file path provided.")
          quit("no")
        }
        
        # Read MCIP file
        mcip.data = tryCatch(nc_open(cmaq.mcip.file),error=function(e) return(NULL))
        cat(paste("   Importing",cmaq.mcip.file,"...\n"))
        
        #Check that nc file has opened and if not return error message to log 
        if (is.null(mcip.data)) {
          stop("MCIP file cannot be opened. Please check the file format.")
          quit("no")
        }
        
        # Check dimensions in MCIP file and CMAQ output file match
        missing.dims = ""
        mcip.dim.names = names(mcip.data$dim)
        if (length(which(mcip.dim.names =="TSTEP"))==0) {missing.dims = paste(missing.dims,"TSTEP")} else {num.timesteps.mcip = mcip.data$dim$TSTEP$len}
        if (length(which(mcip.dim.names =="COL"))==0) {missing.dims = paste(missing.dims,"COL")} else {num.cols.mcip = mcip.data$dim$COL$len}
        if (length(which(mcip.dim.names =="ROW"))==0) {missing.dims = paste(missing.dims,"ROW")} else {num.rows.mcip = mcip.data$dim$ROW$len}
        if (length(which(mcip.dim.names =="LAY"))==0) {missing.dims = paste(missing.dims,"LAY")} else {num.layers.mcip = mcip.data$dim$LAY$len}
        if (nchar(missing.dims) != 0) {
          stop(paste("MCIP data is not in the expected format; the following dimensions are missing:",missing.dims))
          quit("no")
        }
        if (num.timesteps.mcip != num.timesteps) {
          stop("The dimension TSTEP differs between the CMAQ concentration file and the MCIP file.")
          quit("no")
        }
        if (num.cols.mcip != num.cols) {
          stop("The dimension COL differs between the CMAQ concentration file and the MCIP file.")
          quit("no")
        }
        if (num.rows.mcip != num.rows) {
          stop("The dimension ROW differs between the CMAQ concentration file and the MCIP file.")
          quit("no")
        }
        if (num.layers.mcip != num.layers) {
          stop("The dimension LAY differs between the CMAQ concentration file and the MCIP file.")
          quit("no")
        }
        
        # Loop over timesteps
        surf.ref.vector <- c()
        for (t in 1:num.timesteps) { 
          # Get ZF values
          layer.heights.tstep <- ncvar_get(mcip.data, "ZF", start=c(1,1,1,t), count=c(-1,-1,-1,1))
          # average ZF values over x and y
          layer.heights.mean.tstep <- apply(layer.heights.tstep, 3, FUN=mean) # mean height per layer
          # find index of minimum ZF value for each timestep
          surf.ref.element <- which.min(layer.heights.mean.tstep)
          surf.ref.vector <- append(surf.ref.vector, surf.ref.element)
        }
        
        nc_close(mcip.data)
      }
      
      #check which pollutants from definition file are in IO API file
      pollutants=check.ncdf.pol(pollutants, modelled.nc)
      
      #Check there are defined pollutants in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in I/O API netCDF file."))
         quit("no")
      }
       
      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

         #Loop through every timestep in the modelled data file
         for (k in 1:num.timesteps) {
            #Extract all data for the current pollutant
            if (VGTYP == -9999 & num.layers > 1) {
              surf.ref = surf.ref.vector[num.timesteps]
            }
            pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,surf.ref,k),count=c(num.cols,num.rows,1,1)) #4 dimensions x, y lay and time
            if (stations.interp){
               #Get the gridded data for the current timestep in the right format
               onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
            }
            if (!stations.interp){
               onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
            }
             
            if (k==1) {
               onepol.allstations.alltimesteps = onepol.allstations.onetimestep
            } else {
               onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
            }
         }
          
         #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) {
            nc.data = onepol.allstations.alltimesteps
            bDoneFirstPol = TRUE
         } else {
            onepol.allstations.alltimesteps=onepol.allstations.alltimesteps[,-c(1,2)]
            nc.data = cbind(nc.data,onepol.allstations.alltimesteps)
            names(nc.data)[which(names(nc.data)=="onepol.allstations.alltimesteps")]=pollutants$pollutant[j]
         }
      } #end loop over pollutants
       
      #Ensure date column is in the correct format
      nc.data$date=as.POSIXct(nc.data$date, format="%Y-%m-%d %H:%M", tz="GMT")#timezone GMT is assumed
      return(nc.data)
   }
    
   #function for reading CAMS Regional ADS netCDF files as modelled data format
   #function dependencies: stations data frame, pollutants data frame, mod.ncdf.file
   #stations.interp logical variable, check.ncdf.pol function, find.pol.grid.value function,interp.station.value function 
   #packages required: ncdf4, GEOmap, fileds, akima, dplyr, RANN
   read.CAMSRegionalADS.nc=function(mod.ncdf.file) {
      #open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
       
      #Check that nc file has opened and if not, stop and return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("CAMS Regional ADS netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
       
      #Check the dimensions
      missing.dims=""
      dim.names = names(modelled.nc$dim)
      if (!("time" %in% dim.names)) {missing.dims = paste(missing.dims,"time")} else {num.timesteps = modelled.nc$dim$time$len}
      if (!("latitude" %in% dim.names)) {missing.dims = paste(missing.dims,"latitude")} else {num.cols = modelled.nc$dim$latitude$len}
      if (!("longitude" %in% dim.names)) {missing.dims = paste(missing.dims,"longitude")} else {num.rows = modelled.nc$dim$longitude$len}
      if (!("level" %in% dim.names)) {missing.dims = paste(missing.dims,"level")} else {num.levels = modelled.nc$dim$level$len}
      if (nchar(missing.dims) != 0) { #Stop if any dimensions missing 
         stop(paste("CAMS Regional ADS data is not in the expected format; the following dimensions are missing:",missing.dims))
         quit("no")
      }
       
      #Extract the grid values, have checked dimensions so variables should be present.
      grid.xcoords=ncvar_get(modelled.nc,"longitude")
      grid.ycoords=ncvar_get(modelled.nc,"latitude")
      time= ncvar_get(modelled.nc,"time")
      level=ncvar_get(modelled.nc,"level")
       
      #Create allstations.location matrix/map.pos for interpolation/nearest point.
      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
       
      #Can assume a regular grid for CAMS lat-long files
      regular.grid=TRUE
       
      #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 
       
      #Pull out only the region required
      xdim=(range(grid.xcoords)[2]-range(grid.xcoords)[1])/num.cols
      ydim=(range(grid.ycoords)[2]-range(grid.ycoords)[1])/num.rows
      min.x.req = min(stations$map.pos.x)-2*xdim
      max.x.req = max(stations$map.pos.x)+2*xdim
      min.y.req = min(stations$map.pos.y)-2*ydim
      max.y.req = max(stations$map.pos.y)+2*ydim
      ix.req = which(grid.xcoords>=min.x.req & grid.xcoords<=max.x.req)
      iy.req = which(grid.ycoords>=min.y.req & grid.ycoords<=max.y.req)
      ix.req.start = min(ix.req)
      ix.req.count = length(ix.req)
      iy.req.start = min(iy.req)
      iy.req.count = length(iy.req)
      grid.xcoords = grid.xcoords[ix.req.start:(ix.req.start+ix.req.count-1)]
      grid.ycoords = grid.ycoords[iy.req.start:(iy.req.start+iy.req.count-1)]
       
      #Get the date and time in the right format
      secs.in.hour = 60*60
      secs.in.day = 24*60*60
      time.secs = time*secs.in.hour
      time.unit.att=ncatt_get(modelled.nc, "time","long_name") #Extract the units for time to give origin date
      if (time.unit.att$hasatt) {time.unit = time.unit.att$value} else {
         stop(paste("'time' variable in the CAMS Regional ADS netCDF file does not contain the expected 'long_name' attribute."))
      }
      orig.day=gsub("FORECAST time from ","",time.unit) #remove string to get just the origin date
      orig.day = as.POSIXct(orig.day,format="%Y%m%d",tz="GMT")
      time.date = as.POSIXct(time.secs,origin=orig.day,tz="GMT")
    
      #check level units are metres
      level.unit.att=ncatt_get(modelled.nc, "level","units") 
      if (level.unit.att$hasatt) {level.unit = level.unit.att$value} else {
         stop(paste("'level' variable in the CAMS Regional ADS netCDF file does not contain the expected 'units' attribute."))
      }
      if (level.unit!="m") {
         stop(paste("'level' variable in the CAMS Regional ADS netCDF file has units of",level.unit,"; only m is supported."))
      }
      #identify surface layer
      surf.ref=which.min(level) #layer nearest ground has lowest level 
       
      #Find out the pollutants that are in this nc file and save the variable names
      pollutants=check.ncdf.pol(pollutants,modelled.nc)
      #Check there are defined pollutants in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in CAMS Regional ADS netCDF file."))
         quit("no")
      }
       
      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
         #Loop through every timestep in the modelled data file
         for (k in 1:num.timesteps) {
            #Extract all data for the current pollutant
               pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(ix.req.start,iy.req.start,surf.ref,k),count=c(ix.req.count,iy.req.count,1,1))
            if (stations.interp){
               #Get the gridded data for the current timestep in the right format
               onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
            }
            if (!stations.interp){
               onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
            }
             
            if (k==1) {
               onepol.allstations.alltimesteps = onepol.allstations.onetimestep
            } else {
               onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
            }
         }
          
         #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) {
            nc.data = onepol.allstations.alltimesteps
            bDoneFirstPol = TRUE
         } else {
            onepol.allstations.alltimesteps=onepol.allstations.alltimesteps[,-c(1,2)]
            nc.data = cbind(nc.data,onepol.allstations.alltimesteps)
            names(nc.data)[which(names(nc.data)=="onepol.allstations.alltimesteps")]=pollutants$pollutant[j]
         }
      } #end loop over pollutants
      #Ensure date column is in the correct format
      nc.data$date=as.POSIXct(nc.data$date, format="%Y-%m-%d %H:%M", tz="GMT")#timezone GMT is assumed
      return(nc.data)
   }
    
   #function for reading CAMS Regional Reanalysis netCDF files as modelled data format (2D surface files)
   #function dependencies: stations data frame, pollutants data frame, mod.ncdf.file
   #stations.interp logical variable, check.ncdf.pol function, find.pol.grid.value function,interp.station.value function 
   #packages required: ncdf4, GEOmap, fileds, akima, dplyr, RANN
   read.CAMSRegionalReanalysis.nc=function(mod.ncdf.file) {
      #open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
       
      #Check that nc file has opened and if not, stop and return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("CAMS Regional Reanalysis netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
       
      #Check the dimensions
      missing.dims=""
      dim.names = names(modelled.nc$dim)
      if (!("time" %in% dim.names)) {missing.dims = paste(missing.dims,"time")} else {num.timesteps = modelled.nc$dim$time$len}
      if (!("lat" %in% dim.names)) {missing.dims = paste(missing.dims,"lat")} else {num.cols = modelled.nc$dim$lat$len}
      if (!("lon" %in% dim.names)) {missing.dims = paste(missing.dims,"lon")} else {num.rows = modelled.nc$dim$lon$len}
      if (nchar(missing.dims) != 0) { #Stop if any dimensions missing 
         stop(paste("CAMS Regional Reanalysis data is not in the expected format; the following dimensions are missing:",missing.dims))
         quit("no")
      }
       
      #Extract the grid values, have checked dimensions so variables should be present.
      grid.xcoords=ncvar_get(modelled.nc,"lon")
      grid.ycoords=ncvar_get(modelled.nc,"lat")
      time= ncvar_get(modelled.nc,"time")
       
      #Create allstations.location matrix/map.pos for interpolation/nearest point.
      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
       
      #Can assume a regular grid for CAMS lat-long files
      regular.grid=TRUE
       
      #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 
       
      #Pull out only the region required
      xdim=(range(grid.xcoords)[2]-range(grid.xcoords)[1])/num.cols
      ydim=(range(grid.ycoords)[2]-range(grid.ycoords)[1])/num.rows
      min.x.req = min(stations$map.pos.x)-2*xdim
      max.x.req = max(stations$map.pos.x)+2*xdim
      min.y.req = min(stations$map.pos.y)-2*ydim
      max.y.req = max(stations$map.pos.y)+2*ydim
      ix.req = which(grid.xcoords>=min.x.req & grid.xcoords<=max.x.req)
      iy.req = which(grid.ycoords>=min.y.req & grid.ycoords<=max.y.req)
      ix.req.start = min(ix.req)
      ix.req.count = length(ix.req)
      iy.req.start = min(iy.req)
      iy.req.count = length(iy.req)
      grid.xcoords = grid.xcoords[ix.req.start:(ix.req.start+ix.req.count-1)]
      grid.ycoords = grid.ycoords[iy.req.start:(iy.req.start+iy.req.count-1)]
       
      #Get the date and time in the right format
      hour=formatC(24*(time-floor(time)),width=2,flag="0")
      day=format(floor(time),digits=6)
      time.date=as.POSIXct(paste(day,hour,sep=""),format="%Y%m%d%H",tz="GMT")
       
      #Find out the pollutants that are in this nc file and save the variable names
      pollutants=check.ncdf.pol(pollutants,modelled.nc)
      #Check there are defined pollutants in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in CAMS Regional Reanalysis netCDF file."))
         quit("no")
      }
       
      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
         #Loop through every timestep in the modelled data file
         for (k in 1:num.timesteps) {
            #Extract all data for the current pollutant
            pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(ix.req.start,iy.req.start,k),count=c(ix.req.count,iy.req.count,1))
            if (stations.interp){
               #Get the gridded data for the current timestep in the right format
               onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
            }
            if (!stations.interp){
               onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
            }
             
            if (k==1) {
               onepol.allstations.alltimesteps = onepol.allstations.onetimestep
            } else {
               onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
            }
         }
          
         #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) {
            nc.data = onepol.allstations.alltimesteps
            bDoneFirstPol = TRUE
         } else {
            onepol.allstations.alltimesteps=onepol.allstations.alltimesteps[,-c(1,2)]
            nc.data = cbind(nc.data,onepol.allstations.alltimesteps)
            names(nc.data)[which(names(nc.data)=="onepol.allstations.alltimesteps")]=pollutants$pollutant[j]
         }
      } #end loop over pollutants
      #Ensure date column is in the correct format
      nc.data$date=as.POSIXct(nc.data$date, format="%Y-%m-%d %H:%M", tz="GMT")#timezone GMT is assumed
      return(nc.data)
   }
  
   #function for reading Generic netCDF files as modelled data format
   #function dependencies: stations data frame, pollutants data frame, mod.ncdf.file
   #stations.interp logical variable, check.ncdf.pol function, find.pol.grid.value function,interp.station.value function 
   #packages required: ncdf4, GEOmap, fileds, akima, dplyr, RANN, sf
   read.Generic.nc=function(mod.ncdf.file){
       
      #open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
       
      #Check that nc file has opened and if not, stop and return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("Generic netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
       
      # Check the global attribute File_version to see if the generic netCDF is annual background or concentration data
      missing.atts=""
      genericnc.type.att.names=ncatt_get(modelled.nc,0)
      names(genericnc.type.att.names)=tolower(names(genericnc.type.att.names))
      genericnc.type.att=genericnc.type.att.names$file_version
      rm(genericnc.type.att.names) #No longer needed
      if (!is.null(genericnc.type.att)) {genericnc.type = genericnc.type.att} else {missing.atts = paste(missing.atts,"File_version")}
         # Stop if no File_version attribute is found 
      if (nchar(missing.atts)!=0){
         stop(paste("Generic netCDF data is not in the expected format; the following attributes are missing: ",paste(missing.atts,collapse=","), sep=""))
      }
      
      #Check the value of the title attribute matches ones of the following: ADMS 3D hourly concentration file,
      # MAQS-Health 3D Hourly concentration file, ADMS 2D Average Concentration file , MAQS-Health 2D Average Concentration file 
      if (genericnc.type == "ADMS_2DCF_v1.0"){
         genericnc.type="2D"
      } else if (genericnc.type == "ADMS_3DCF_v1.0"){
         genericnc.type="3D"
      } else {
         stop(paste("The File_version attribute in the Generic netCDF file: ", genericnc.type, " \n is not compatible with supported formats. \n The File_version attribute should be either 'ADMS_3DCF_v1.0' or 'ADMS_2DCF_v1.0'.", sep=""))
      }
       
      # Check the dimensions in the Generic netCDF file
      missing.dims=""
      dim.names=names(modelled.nc$dim)
      if (!("x" %in% dim.names)){missing.dims=paste(missing.dims,"x")} else {num.cols=modelled.nc$dim$x$len}
      if(!("y" %in% dim.names)){missing.dims=paste(missing.dims,"y")} else {num.rows=modelled.nc$dim$y$len}
      if (genericnc.type=="3D"){ #time and z should also be read for 3D hourly concentration files
         if(!("z" %in% dim.names)){missing.dims=paste(missing.dims,"z")} else {num.layers=modelled.nc$dim$z$len}
         if(!("time" %in% dim.names)){missing.dims=paste(missing.dims,"time")} else {num.timesteps=modelled.nc$dim$time$len}
      }
      if (nchar(missing.dims) != 0) { # Stop if any dimensions missing 
         stop(paste("Generic netCDF data is not in the expected format; the following dimensions are missing:",missing.dims))
         quit("no")
      }
      if (genericnc.type=="2D"){ # For 2D average files, only 1 timestep per file
         num.timesteps=1
      }
       
      # Extract variables for grid, for 2D files extract x and y, for 3D files extract x,y,z and time
      grid.xcoords=ncvar_get(modelled.nc,"x")
      grid.ycoords=ncvar_get(modelled.nc,"y")
      if (genericnc.type=="3D"){
         grid.layers=ncvar_get(modelled.nc,"z")
         time=ncvar_get(modelled.nc,"time")
      }
       
      # Extract name of variable containg projection information 
      projection.var.att=ncatt_get(modelled.nc,"x","grid_mapping")
      if (projection.var.att$hasatt){projection.var=projection.var.att$value} else {
         stop("Projection attribute 'grid_mapping' is missing on the x variable.")
         quit("no")
      }
       
      # Get projection information
      projection.type.att=ncatt_get(modelled.nc,projection.var,"grid_mapping_name")
      if(projection.type.att$hasatt){projection.type=projection.type.att$value} else{
         stop(paste("The required attribute 'grid_mapping_name' cannot be found on the projection variable: ", projection.var, sep=""))
         quit("no")
      }
       
      # Check the projection is allowed
      # For 3D files, LCC, PS and UTM are allowed 
      if (genericnc.type=="3D"){
         allowed.GenericnetCDF.proj=c("lambert_conformal_conic","transverse_mercator","polar_stereographic")
         if (!(projection.type %in% allowed.GenericnetCDF.proj)){
            stop(paste("The projection ",projection.type," is not supported for Generic netCDF formats. \n Only the following projections are supported for 3D Generic netCDF files: ",paste(allowed.GenericnetCDF.proj,collapse=", "), sep=""))
         }
         if (projection.type=="lambert_conformal_conic"){ # abbreviate lambert_conformal_conic to LCC
            projection.type="LCC"
         }
         if (projection.type=="transverse_mercator"){# abbreviate universial_transverse_mercator to UTM
            projection.type="UTM"
         } 
         if (projection.type=="polar_stereographic"){ # abbreviate polar_stereographic to PS
            projection.type="PS"
         }
      } else if (genericnc.type=="2D"){
         allowed.GenericnetCDF.proj=c("osgb")
         if (!(projection.type %in% allowed.GenericnetCDF.proj)){
            stop(paste("The projection ",projection.type," is not supported for Generic netCDF formats. \n Only the following projections are supported for 2D Generic netCDF files: ",paste(allowed.GenericnetCDF.proj,collapse=", "), sep=""))
         }
      }
       
      # Check the required projection parameters are present as attributes on the grid mapping variable 
      # Note no projection parameters are required for OSGB
      # For LCC, require LAT0, LON0, LAT1, LAT2, FE, FN
      missing.atts=""
      if (projection.type=="LCC"){
         lat0.att=ncatt_get(modelled.nc,projection.var,"latitude_of_projection_origin")
         if (lat0.att$hasatt){LAT0=lat0.att$value} else {missing.atts=paste(missing.atts, "latitude_of_projection_origin")}
         lon0.att=ncatt_get(modelled.nc,projection.var,"longitude_of_central_meridian")
         if (lon0.att$hasatt){LON0=lon0.att$value} else {missing.atts=paste(missing.atts, "longitude_of_central_meridian")}
         standardpar.att=ncatt_get(modelled.nc,projection.var,"standard_parallel")
         if (standardpar.att$hasatt){
            LAT1=standardpar.att$value[1]
            LAT2=standardpar.att$value[2]
         } else {missing.atts=paste(missing.atts,"standard_parallel")}
         FE.att=ncatt_get(modelled.nc,projection.var,"false_easting")
         if (FE.att$hasatt) {FE=FE.att$value}else{missing.atts=paste(missing.atts,"FE")}
         FN.att=ncatt_get(modelled.nc, projection.var, "false_northing")
         if (FN.att$hasatt) {FN=FN.att$value}else{missing.atts=paste(missing.atts,"FN")}
      } 
      # For PS, require LAT0, LON0, FE, FN 
      if (projection.type=="PS"){
         lat0.att=ncatt_get(modelled.nc,projection.var,"latitude_of_projection_origin")
         if (lat0.att$hasatt){LAT0=lat0.att$value} else {missing.atts=paste(missing.atts, "latitude_of_projection_origin")}
         lon0.att=ncatt_get(modelled.nc,projection.var,"straight_vertical_longitude_from_pole")
         if (lon0.att$hasatt){LON0=lon0.att$value} else {missing.atts=paste(missing.atts, "straight_vertical_longitude_from_pole")}
         FE.att=ncatt_get(modelled.nc,projection.var,"false_easting")
         if (FE.att$hasatt) {FE=FE.att$value}else{missing.atts=paste(missing.atts,"FE")}
         FN.att=ncatt_get(modelled.nc, projection.var, "false_northing")
         if (FN.att$hasatt) {FN=FN.att$value}else{missing.atts=paste(missing.atts,"FN")}
         sf.att=ncatt_get(modelled.nc,projection.var,"scale_factor_at_projection_origin")
         if (sf.att$hasatt){sf=sf.att$value} # sf may be supplied but not required
      }
       
      # For UTM, require longitude of central meridian
      if (projection.type=="UTM"){
         UTM.att=ncatt_get(modelled.nc,projection.var,"longitude_of_central_meridian")
         if (UTM.att$hasatt){UTM=UTM.att$value} else {missing.atts=paste(missing.atts, "longitude_of_central_meridian")}
      }
       
      # Stop if any of the required projection attributes are missing
      if (!nchar(missing.atts==0)){
         stop(paste("The following required projection attributes are missing: ", missing.atts, sep=""))
         quit("no")
      }
       
      # Create projections and project stations onto the projected grid 
      if (projection.type=="LCC"){
         # Set up the map projection: type = 3 for LCC
         map.proj=setPROJ(type=3, LAT0=LAT0,LON0=LON0, LAT1 = LAT1, LAT2 = LAT2, FE=FE, FN=FN)
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      }else if (projection.type=="PS"){
         #Set up the map projection: type = 4 for PS, depends on whether standard parallel or scale factor at projected origin has been supplied
         map.proj=setPROJ(type=4, LAT0 =LAT0, LON0 = LON0, FE=FE, FN=FN)
         #Get projected map positions of the stations
         stations.xy=stereo.sphr.xy(stations$latitude, stations$longitude, map.proj)
         # If a scale factor has been provided
         if (sf.att$hasatt){
            # Apply scale factor
            stations.xy$x=(stations.xy$x)*sf
            stations.xy$y=(stations.xy$y)*sf
         }
         #Store projected positions in stations table and in station locations matrix
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      }else if (projection.type=="UTM"){
         LON0=UTM
         # Set up the map projection: type = 7 for utm.wgs84
         map.proj=setPROJ(type=7, LON0=LON0)
         # Project stations onto the projected grid using lat long from stations file and map projection created
         stations.xy=utm.wgs84.xy(stations$latitude,stations$longitude,map.proj)
         # Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      }else if (projection.type=="osgb"){
         # Project stations onto the projected grid using lat long from stations file  
         stations.xy=cerc.WGS84.to.OSGB(stations$latitude, stations$longitude)
         # Add projected grid position columns to the station definition file with the projected x and y coordinates for each station 
         stations$map.pos.x = stations.xy$X
         stations$map.pos.y = stations.xy$Y  
         allstations.location = matrix(data=c(stations.xy$X,stations.xy$Y),nrow=length(stations$name), ncol=2)
      }
       
      # check dimensions and define whether it is a regular grid 
      num.dimensions.x = length(dim(grid.xcoords))
      num.dimensions.y = length(dim(grid.ycoords))
      if (num.dimensions.x==num.dimensions.y) { #check x and y have the same dimensions, if not stop. 
         num.dimensions = num.dimensions.x
      } else {
         stop(paste("netCDF variables for defining the grid (e.g. x and y) do not have the same number of dimensions.",sep=""))
         quit("no")
      }         
      # Generic netCDF files have to have a regular grid so set regular.grid=TRUE
      regular.grid = TRUE
       
      # For 3D files extract surface layer 
      if (genericnc.type=="3D"){
         layers=ncvar_get(modelled.nc,"z") # Extract all layers
         surf.ref=which.min(layers) # Identify surface layer 
      }
       
      # Get time profile for the file
      if (genericnc.type=="3D"){ # Convert the time from j days to calendar days 
         # Extract units for jday e.g days since
         origin.date.att = ncatt_get(modelled.nc,"time","units")
         if (origin.date.att$hasatt){origin.date=origin.date.att$value}else{
            stop("No origin date can be found for the time variable.")
         }
         origin.date=gsub("days since ","",origin.date)
         if (nchar(origin.date)==10){ # For the case when HH:MM:SS is not added
            origin.date=paste(origin.date,"00:00:00")
         }
         origin.date=as.POSIXct(origin.date, "%Y-%m-%d %H:%M:%S", tz="GMT") # Assumed tz = GMT for MAQS-Health
         time.date=as.POSIXct(86400*time,origin=origin.date,format="%Y-%m-%d %H:%M",tz="GMT") #convert to POSIXct format, using origin date for jday to cal date
         time.date=round_date(time.date, "hour")
      } 
      missing.atts=""
      if (genericnc.type=="2D"){ # extract the start and end date for the averaging period
         start.date.att=ncatt_get(modelled.nc, 0, "Period_start_date")
         if (start.date.att$hasatt){start.date=start.date.att$value}else{missing.atts=paste(missing.atts,"Period_start_date")}
         end.date.att=ncatt_get(modelled.nc,0,"Period_end_date")
         if (end.date.att$hasatt){end.date=end.date.att$value}else{missing.atts=paste(missing.atts, "Period_end_date")}
         if (!(nchar(missing.atts==0))){
            stop(paste("The following attributes cannot be found and are required for 2D Average Generic netCDF files: ",missing.atts,sep=""))
            quit("no")
         }
         start.date=as.POSIXct(start.date, "%Y-%m-%d",tz="GMT")
         time.date=start.date
      }
       
      #Find out the pollutants that are in this nc file and save the variable names
      pollutants=check.ncdf.pol(pollutants,modelled.nc)
      #Check there are defined pollutants in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in Generic netCDF file."))
         quit("no")
      }
       
      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
         #Loop through every timestep in the modelled data file
         for (k in 1:num.timesteps) {
            if (genericnc.type=="3D"){
              pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,surf.ref,k),count=c(num.cols,num.rows,1,1))
            } else if (genericnc.type=="2D"){
               # Assume x or y is never equal to 1
               pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1),count=c(num.cols,num.rows))
            }
            if (stations.interp){
               #Get the gridded data for the current timestep in the right format
               onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
            }
            if (!stations.interp){
               onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
            }
             
            if (k==1) {
               onepol.allstations.alltimesteps = onepol.allstations.onetimestep
            } else {
               onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
            }
         }
          
         #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) {
            nc.data = onepol.allstations.alltimesteps
            bDoneFirstPol = TRUE
         } else {
            onepol.allstations.alltimesteps=onepol.allstations.alltimesteps[,-c(1,2)]
            nc.data = cbind(nc.data,onepol.allstations.alltimesteps)
            names(nc.data)[which(names(nc.data)=="onepol.allstations.alltimesteps")]=pollutants$pollutant[j]
         }
      } #end loop over pollutants
       
      return(nc.data)
       
   }
    
   #function for reading CHIMERE netCDF files as modelled data format
   #function dependencies: stations data frame, pollutants data frame, mod.ncdf.file
   #stations.interp logical variable, check.ncdf.pol function, find.pol.grid.value function,interp.station.value function 
   #packages required: ncdf4, GEOmap, fileds, akima, dplyr, RANN, sf
   read.CHIMERE.nc = function (mod.ncdf.file, chimere.cco.file){
      #open nc file and write file name to log
      modelled.nc = tryCatch(nc_open(mod.ncdf.file),error=function(e) return(NULL))
      cat(paste("   Importing",mod.ncdf.file,"...\n"))
      
      #Check that nc file has opened and if not, stop and return error message to log 
      if (is.null(modelled.nc)){
         stop(paste("CHIMERE netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
      
      #Check to see if opened file is a CHIMERE file
      missing.atts =""
      title.att = ncatt_get(modelled.nc,0,attname="Title") #Extract Title global attribute
      if (title.att$hasatt) {Title = title.att$value} else {missing.atts = paste(missing.atts,"Title")}
      #Stop if map proj attribute doesn't exist 
      if (nchar(missing.atts)!=0) {
         stop(paste("CHIMERE netCDF data is not in the expected format; the following attributes are missing:",missing.atts))
         quit("no")
      }
      title=title.att$value
      if (title != "CHIMERE SUITE"){
         stop(paste("CHIMERE netCDF file has unexpected title: ", title,". Expected title is CHIMERE SUITE.\nPlease check file format."))
         quit("no")
      }
      
      #Read cco file 
      CHIMEREcco=readLines(chimere.cco.file)
      #find what projection the file is
      CHIMERE.proj=""
      if (length(which(CHIMEREcco=="UTM"))>0){CHIMERE.proj="UTM"}
      if (length(which(CHIMEREcco=="LCC"))>0){CHIMERE.proj="LCC"}
      if (length(which(CHIMEREcco=="PS"))>0){CHIMERE.proj="PS"}
      #If no projection is found, stop and send error message to log
      if(CHIMERE.proj==""){
         stop(paste("Cannot find valid projection in CHIMERE cco file. Projection must be UTM, LCC or PS"))
         quit("no")
      }
      
      #Extract proj params, depends on each type of proj
      if (CHIMERE.proj=="UTM"){ #For UTM zone number is required to be converted into LON0
         proj.line=(which(CHIMEREcco=="UTM")) 
         UTM_Zone=as.numeric(CHIMEREcco[proj.line+1])
         UTM_hemi=CHIMEREcco[proj.line+2]
      } else if (CHIMERE.proj=="LCC"){ #For LCC LAT1, LAT2, LAT0, LON0 required 
         proj.line=(which(CHIMEREcco=="LCC")) 
         Ref.Long=CHIMEREcco[proj.line+1]
         Ref.Lat=CHIMEREcco[proj.line+2]
         STDP1=CHIMEREcco[proj.line+3]
         STDP2=CHIMEREcco[proj.line+4]
      } else if (CHIMERE.proj=="PS"){ #For PS LAT0 and LON0 required
         proj.line=(which(CHIMEREcco=="PS"))
         Ref.Long=CHIMEREcco[proj.line+1]
         Ref.Lat=CHIMEREcco[proj.line+2]
      }
      
      #Check and extract dimensions 
      dim.names=names(modelled.nc$dim)
      missing.dims=""
      if (!("Time" %in% dim.names)){missing.dims = paste(missing.dims,"Time")} else {num.timesteps=modelled.nc$dim$Time$len}
      if (!("west_east" %in% dim.names)){missing.dims = paste(missing.dims,"west_east")} else {num.cols=modelled.nc$dim$west_east$len}
      if (!("south_north" %in% dim.names)){missing.dims = paste(missing.dims,"south_north")} else {num.rows=modelled.nc$dim$south_north$len}
      if (!("bottom_top" %in% dim.names)){missing.dims = paste(missing.dims,"bottom_top")} else {num.layers=modelled.nc$dim$bottom_top$len}
      if (nchar(missing.dims) != 0) {
         stop(paste("CHIMERE netCDF is not in the expected format; the following dimensions are missing:",missing.dims))
         quit("no")
      }
      
      #Check and extract required variables 
      var.names=names(modelled.nc$var)
      missing.vars=""
      if (!("lat" %in% var.names)){missing.vars = paste(missing.vars,"lat")} else {lat=ncvar_get(modelled.nc, "lat")}
      if (!("lon" %in% var.names)){missing.vars = paste(missing.vars,"lon")} else {lon=ncvar_get(modelled.nc, "lon")}
      if (!("Times" %in% var.names)){missing.vars = paste(missing.vars,"Times")} else {times=ncvar_get(modelled.nc, "Times")}
      if (!("hlay" %in% var.names)){missing.vars = paste(missing.vars,"hlay")} else {layers=ncvar_get(modelled.nc, "hlay")}
      if (nchar(missing.vars) != 0) {
         stop(paste("CHIMERE netCDF file is not in the expected format; the following variables are missing:",missing.vars))
         quit("no")
      }
      
      
      #Convert times vector to POSIXct format 
      time.date=as.POSIXct(times,format="%Y-%m-%d_%H:%M:%S",tz="GMT")
      
      #create the projection to project grid and stations. Depends on file projection
      if (CHIMERE.proj=="UTM") { #For UTM
         #Convert UTM_zone which gives UTM zone number to LON0 required for creating the projection
         lon.0= 6.0 * (UTM_Zone - 31.0) + 3.0
         #Use Lon.0 as longitude of origin for projection
         #Type = 7 for utm.wgs84
         map.proj=setPROJ(type=7, LON0=lon.0)
         #get projected grid xy coordinates from the latitude and longitude variables and map projection 
         grid.xy = utm.wgs84.xy(lat,lon,map.proj) #lat lon have only EW and NS dims
         #Extract x and y grid coordinates
         grid.xcoords=grid.xy$x
         grid.ycoords=grid.xy$y
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=utm.wgs84.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      } else if (CHIMERE.proj=="LCC") { #for LCC
         #Create map projection using LON0, LAT0, LAT1 and LAT2
         #type = 3 for LCC 
         map.proj=setPROJ(type=3, LAT0=,LON0=, LAT1 = , LAT2 =)
         #Create the projected regular grid from the XLAT and XLONG variables
         grid.xy = cerc.lambert.cc.xy(lat,lon,map.proj) 
         #Extract x and y grid coordinates
         grid.xcoords=grid.xy$x
         grid.ycoords=grid.xy$y
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      } else if (CHIMERE.proj=="PS") { #for PS
         #Create map projection using LON0 and LAT0 
         #Type = 4 for PS
         map.proj=setPROJ(type=4, LON0=STAND_LON, LAT0 = )
         grid.xy = stereo.sphr.xy(lat,lon,map.proj)
         #Extract x and y grid coordinates
         grid.xcoords=grid.xy$x
         grid.ycoords=grid.xy$y
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=stereo.sphr.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      }
      
      #Check for regular or irregular grid as this determines method of interpolation. 
      num.dimensions.x = length(dim(grid.xcoords))
      num.dimensions.y = length(dim(grid.ycoords))
      #Check dimensions of x are equal to dimensions of y and stop if not
      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 a regular grid, x and y have 1 dimension, if irregular x and y have 2 dimensions.
      if (num.dimensions==1) {
         regular.grid = TRUE
      } else if(num.dimensions==2) {
         regular.grid = FALSE
      } else {
         #If x and y have more than 2 dimensions, stop as this is not compatible 
         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")
      }
      
      
      #identify the surface layer 
      if((num.timesteps==1) & (num.layers==1)){
         mean.height=mean(as.matrix(layers))
      } else if (num.layers==1) {
         layers=layer[,,1]
         mean.height=mean(as.matrix(layers))
      } else if (num.timesteps==1){
         mean.heights=colMeans(layers,dims=2)
         surf.ref=which.min(mean.heights)
      } else {
         layers=layer[,,,1]
         mean.heights=colMeans(layers,dims=2)
         surf.ref=which.min(mean.heights)
      }
      
      
      #Find out the pollutants that are in this nc file and save the variable names
      pollutants=check.ncdf.pol(pollutants,modelled.nc)
      #Check there are defined pollutants in netCDF file
      if(!(any(pollutants$bmod == TRUE))){
         stop(paste("No defined pollutants found in CHIMERE netCDF file."))
         quit("no")
      }
      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 data for the current pollutant
         for (k in 1:num.timesteps){ #Reading in one hour of data at a time from the CHIMERE file
            if (modelled.nc$var[[pollutants$var.name[j]]]$ndims==4){
               pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,k,surf.ref),count=c(num.cols,num.rows,1,1)) #Gets one hour of data at a time.
            } else if (modelled.nc$var[[pollutants$var.name[j]]]$ndims==3){
               pollutant.data = ncvar_get(modelled.nc,pollutants$var.name[j],start=c(1,1,k),count=c(num.cols,num.rows,1))
            } else{
               stop(paste("Can only process EMEP files where variables have 3 or 4 dimensions."))
            }
            if (stations.interp){
               onepol.allstations.onetimestep=interp.station.value(grid.xcoords,grid.ycoords,pollutant.data, allstations.location,stations,regular.grid,1,time.date,k)
            } else{ #If !stations.interp
               onepol.allstations.onetimestep=find.pol.grid.value(stations,grid.xcoords, grid.ycoords, pollutant.data,1, time.date, k,regular.grid)
            }
            if (k==1) {
               onepol.allstations.alltimesteps = onepol.allstations.onetimestep
            } else {
               onepol.allstations.alltimesteps = rbind(onepol.allstations.alltimesteps,onepol.allstations.onetimestep)
            }
         }
        
         #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) {
            nc.data = onepol.allstations.alltimesteps
            bDoneFirstPol = TRUE
         } else {
            onepol.allstations.alltimesteps=onepol.allstations.alltimesteps[,-c(1,2)]
            nc.data = cbind(nc.data,onepol.allstations.alltimesteps)
            names(nc.data)[which(names(nc.data)=="onepol.allstations.alltimesteps")]=pollutants$pollutant[j]
         }
      } #end loop over pollutants
      #Ensure date column is in the correct format
      nc.data$date=as.POSIXct(nc.data$date, format="%Y-%m-%d %H:%M", tz="GMT")#timezone GMT is assumed
      return(nc.data)
      
   }
   
   # function for reading Defra Background Maps csv files as modelled data format
   # function dependencies: stations data frame, pollutants data frame, dbm.file, stations.interp logical variable,
   #                        cerc.WGS84.to.OSGB function
   # packages required: RANN, interp, data.table
   read.DefraBackgroundMaps.csv=function(dbm.file){
      # open and read Defra Background Maps csv file
      # skip first 5 lines which contain useful information but not the data
      dbm.data=tryCatch(fread(dbm.file, skip=5),error=function(e) return(NULL))
      #check the file could be opened and read, if not stop and send error to log
      if (is.null(dbm.data)){
         stop(paste("Defra Background Maps data file cannot be opened. \n Please check the correct file path and format has been supplied."))
         quit("no")
      }  
       
      #Check for missing headers 
      missing.headers = c()
      if (!("x" %in% names(dbm.data))) missing.headers = c(missing.headers,"'x'")
      if (!("y" %in% names(dbm.data))) missing.headers = c(missing.headers,"'y'")
      if (!("ukgridcode" %in% names(dbm.data) | "gridcode" %in% names(dbm.data))) {
         missing.headers = c(missing.headers,"'ukgridcode'/'gridcode'")
      }
      if (length(missing.headers)>0){
         stop(paste("Expected column header(s)",paste(missing.headers, collapse=", "),
                   "are missing from the Defra Background Maps file. Invalid format."))
      }
     
      #Extract the information from the first 4 lines of the dbm file 
      dbm.info=readLines(dbm.file, n=4)
      dbm.pollutant=tolower(trimws(gsub(",","",dbm.info[1])))
      dbm.year=as.numeric(trimws(gsub(",","",dbm.info[2])))
      dbm.statistic=tolower(trimws(gsub(",","",dbm.info[3])))
      dbm.units=tolower(trimws(gsub(",","",dbm.info[4])))
       
      #Check pollutant is in pollutants definitions
      if (!(dbm.pollutant %in% pollutants$mod.alias)){
         stop(paste("The pollutant in the Defra Background Maps CSV file, ", toupper(dbm.pollutant),
                    ", is not present in the modelled data pollutant definitions file.\n", sep=""))
         quit("no")
      } 
      
      #Check that the statistic is annual mean.
      #(All Defra Background Map datasets are annual mean except O3, which is the number of days where the 
      # daily max 8-hour rolling mean exceeds 120 ug/m3 - this dataset isn't supported)
      if (dbm.statistic != "annual mean"){
         stop(paste("The Defra Background Maps CSV file statistic is '", dbm.statistic, 
                   "', however only annual mean Defra Background Map datasets are supported.\n", sep=""))
         quit("no")
      }
       
      #Set units to recognized value
      #e.g. some files have "ug m-3 (NO2)" or "ug m-3 (gravitational units)"
      if (grepl("ug m-3", dbm.units)) dbm.units="ug/m3"
      if (grepl(" m-3", dbm.units)) dbm.units=gsub(" m-3", "/m3", dbm.units)
      mod.units = pollutants$mod.units[which(pollutants$mod.alias == dbm.pollutant)]
      if (dbm.units != mod.units){
         stop(paste(toupper(dbm.pollutant)," has units '",dbm.units,"' in the Defra Background Maps CSV file ",
                   "but units '",mod.units,"' in the modelled pollutants definitions file.\n",
                   "  Please correct the units in the modelled pollutants definitions file.",sep=""))
         quit("no")
      }
       
      #Set the date to the beginning of the year 
      dbm.date=ISOdatetime(dbm.year,1,1,0,0,0,tz="GMT")
       
      #Assume that the concentration data is in the 4th column - it always should be
      names(dbm.data)[4]="mod"
       
      #Convert the pollutant data values of "MISSING" to NA for consistency with rest of the Data Input Tool
      dbm.data$mod[dbm.data$mod == "MISSING"]=NA
       
      #Remove ukgridcode/gridcode column from the dataframe
      dbm.data=subset(dbm.data, select = c("x","y","mod"))
       
      #Remove NA values from the dbm data (these can't be used in the interpolation function) and convert to numeric
      dbm.data=subset(dbm.data,!is.na(mod))
      dbm.data$mod=as.numeric(dbm.data$mod)
    
      #Convert the stations lat lon WGS84 values to OSGB x y coordinates 
      stations.xy=cerc.WGS84.to.OSGB(stations$latitude, stations$longitude)
      stations$map.pos.x = stations.xy$X #store x values in stations df
      stations$map.pos.y = stations.xy$Y #store y values in stations df
      
      if (stations.interp){
         #Calculate concentration by interpolation (interp function gives NA to any station outside modelled domain)
         interp.conc = interp::interp(dbm.data$x,dbm.data$y,dbm.data$mod,
                             stations$map.pos.x,stations$map.pos.y,output="points")$z
         onepol.allstations=data.frame(date=rep(dbm.date,length(stations$map.pos.x)), 
                                      station=stations$name,mod=interp.conc) 
      } else {
         #Get concentration from nearest grid point
         stations.grid.values=subset(stations, select=c("map.pos.x","map.pos.y"))
         grid.values=cbind(dbm.data$x, dbm.data$y)
         near.points.index=nn2(grid.values,stations.grid.values,k=1)[[1]]
         stations.nn=cbind(stations,near.points.index)
         stations.nn$nearest.point.x=grid.values[stations.nn$near.points.index,1]#add col to stations with nearest grid x val
         stations.nn$nearest.point.y=grid.values[stations.nn$near.points.index,2]#add col to stations with nearest grid y val
         stations.nn$pol.conc=dbm.data$mod[stations.nn$near.points.index]
        
         onepol.allstations = data.frame(date=dbm.date,station=stations.nn$name,mod=stations.nn$pol.conc)
        
         #Identify stations outside the modelled area and set pollutant concentration value to NA
         stations.outside=which(!(between(stations$map.pos.x,range(dbm.data$x)[1], range(dbm.data$x)[2])) |!(between(stations$map.pos.y,range(dbm.data$y)[1], range(dbm.data$y)[2])))
         if (length(stations.outside)>0) onepol.allstations$mod[stations.outside]=NA
          
      }
      #Re-name concentration column with the pollutant name 
      names(onepol.allstations)[names(onepol.allstations)=="mod"]=dbm.pollutant
       
      return(onepol.allstations)
       
   }
    
   # function to read in met data in a ADMS .pro format
   # requires packages data.table 
   read.pro=function(pro.file.path){
      #Read the PRO file and return an error if file is unreadable
      cat(paste("   Importing PRO met data file: ",pro.file.path,"...\n", sep=""))
      pro.data=tryCatch(fread(pro.file.path),error=function(e) return(NULL))
      if (is.null(pro.data)) {
         stop(paste("PRO met data file cannot be opened or read, please check the file path and format.\n ",sep=""))
         quit("no")
      }
      #Make sure the PRO file contains the expected columns and stop with an error if not
      expected.cols=c("Year","Day","Hour","z(m)","U(z)(m/s)","Phi(z)(deg)","Temp(z)(C)")
      missing.cols = c() #empty list
      for (i in 1:length(expected.cols)) {
         if (!(expected.cols[i] %in% names(pro.data))) {
            missing.cols = c(missing.cols,expected.cols[i])
         }
      }
      if (length(missing.cols)>0) {
         stop(paste("PRO met data file is not in the correct format.\n",
                    "These expected columns are missing: ",paste(missing.cols,collapse = ", "),"\n",sep=""))
         quit("no")
      }
      #Check the PRO file only contains data for a single height
      nHeights=length(unique(pro.data$`z(m)`))
      if (nHeights>1) {
         stop(paste("PRO met data file contains data for more than one height.",
                    "Only single height PRO met data files are supported.",sep="\n"))
         quit("no")
      }
      #We have all the columns we need, so proceed to processing
      pro.data$date=as.POSIXct(paste(pro.data$Year, "-",pro.data$Day,"-",pro.data$Hour, sep=""), format="%Y-%j-%H", tz="GMT")
      pro.data=subset(pro.data, select=c("date","U(z)(m/s)","Phi(z)(deg)","Temp(z)(C)"))
      names(pro.data)=c("date","ws","wd","temp")
      return(pro.data)
   }
   
   # function to read in met data in a .mio format
   # requires packages data.table 
   read.mio=function(mio.file.path){
      #Read the MIO file and return an error if file is unreadable
      cat(paste("   Importing MIO met data file: ",mio.file.path,"...\n", sep=""))
      mio.data=tryCatch(fread(mio.file.path),error=function(e) return(NULL))
      if (is.null(mio.data)) {
         stop(paste("MIO met data file cannot be opened or read, please check the file path and format.\n ",sep=""))
         quit("no")
      }
      #Make sure the MIO file contains the expected columns and stop with an error if not
      expected.cols=c("YEAR","DAY","HOUR","PROCESSED_DATA:")
      missing.cols = c() #empty list
      for (i in 1:length(expected.cols)) {
         if (!(expected.cols[i] %in% names(mio.data))) {
            missing.cols = c(missing.cols,expected.cols[i])
         }
      }
      if (length(missing.cols)>0) {
         stop(paste("MIO met data file is not in the correct format.\n",
                     "These expected columns are missing: ",paste(missing.cols,collapse = ", "),"\n",sep=""))
         quit("no")
      }
      #Make sure the MIO file contains the expected met data columns to the right of "PROCESSED_DATA:" and stop with an error if not
      processed.ref=which(names(mio.data)=="PROCESSED_DATA:")
      met.names=names(mio.data)[(processed.ref+1):(length(names(mio.data)))]
      expected.cols=c("U","PHI","T0C","H/LMO")
      missing.cols = c() #empty list
      for (i in 1:length(expected.cols)) {
        if (!(expected.cols[i] %in% met.names)) {
          missing.cols = c(missing.cols,expected.cols[i])
        }
      }
      if (length(missing.cols)>0) {
        stop(paste("MIO met data file is not in the correct format.\n",
                   "These expected columns to the right of 'PROCESSED_DATA:' are missing: ",paste(missing.cols,collapse = ", "),"\n",sep=""))
        quit("no")
      }
      #We have all the columns we need, so proceed to processing
      mio.data$date=as.POSIXct(paste(mio.data$YEAR, "-",mio.data$DAY,"-",mio.data$HOUR, sep=""), format="%Y-%j-%H", tz="GMT")
      mio.data=mio.data[,(processed.ref+1):(ncol(mio.data))]
      mio.data=subset(mio.data, select=c("date","U","PHI","T0C","H/LMO"))
      names(mio.data)=c("date","ws","wd","temp","h.over.lmo")
      return(mio.data)
   }
    
   #function for extracting the required met data from WRF and ADMS generic gridded data
   #function dependencies: find.met.grid.value
   #packages required: 
   extract.met=function(grid.xcoords,grid.ycoords,regular.grid,times,wind.v,wind.u,stations,
                        bTemp,temp,bHeatFlux,heatflux,bBL,BL.height,bCloud,cloudcover,bPr,precipitation){ 
       
      #How many grid points in the met data?
      num.x.grids=length(grid.xcoords)
      num.y.grids=length(grid.ycoords)
      num.timesteps = length(times)
    
      #EXTRACT THE MET VALUES FOR THE REQUIRED GRIDPOINT
      bOneGridpoint = FALSE
      if (num.x.grids==1 & num.y.grids==1){
         bOneGridpoint = TRUE
         #Only one grid point, so use this for every station
         met.data=data.frame(date=times, wind.v.comp=wind.v, wind.u.comp=wind.u)
         if (bTemp){
           met.data$temp=temp
         }
         if (bHeatFlux){
           met.data$HEAT_FLUX=heatflux
         }
         if (bBL){
           met.data$BLAYER_HGT=BL.height
         }
         if (bCloud){
           met.data$CLOUD_COVER=cloudcover
         }
         if (bPr){
           met.data$Pr=precipitation
         }
         #Make an empty data frame with a timestep for every station
         stations.dates=data.frame(date=rep(times,length(stations$name)),station=rep(stations$name,length(times))[order(rep(stations$name,length(times)))])
         #Merge with the met data to get a data frame of met values for every station for every timestep 
         met.data = merge(stations.dates,met.data,by="date")
      } else {
         # loop through each time step and extract the met parameters for the required grid point 
         # wind v
         bDoneFirstTimestep=FALSE
         for (i in 1:num.timesteps){
            if (num.timesteps==1) {
               time.date = times
               wind.u.grid = wind.u
               wind.v.grid = wind.v
               if (bTemp) temp.grid = temp
               if (bHeatFlux) heatflux.grid = heatflux
               if (bBL) BL.height.grid = BL.height
               if (bCloud) cloud.grid = cloudcover
               if (bPr) precipitation.grid = precipitation
            } else {
               time.date = times[i]
               wind.u.grid = wind.u[,,i]
               wind.v.grid = wind.v[,,i]
               if (bTemp) temp.grid = temp[,,i]
               if (bHeatFlux) heatflux.grid = heatflux[,,i]
               if (bBL) BL.height.grid = BL.height[,,i]
               if (bCloud) cloud.grid = cloudcover[,,i]
               if (bPr) precipitation.grid = precipitation[,,i]
            }
            one.timestep.wind.u.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                            met.parameter.grid =  wind.u.grid, num.timesteps =  num.timesteps, time.date =  time.date, 
                                                            k=1, regular.grid = regular.grid )
            names(one.timestep.wind.u.data)[names(one.timestep.wind.u.data)=="met.value"]="wind.u.comp"
            one.timestep.wind.v.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                            met.parameter.grid =  wind.v.grid, num.timesteps =  num.timesteps, time.date = time.date, 
                                                            k=1, regular.grid = regular.grid )
            names(one.timestep.wind.v.data)[names(one.timestep.wind.v.data)=="met.value"]="wind.v.comp"
            one.timestep.data = merge(one.timestep.wind.u.data,one.timestep.wind.v.data,by=c("date","station"))
            if (bTemp) {
              one.timestep.temp.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                           met.parameter.grid =  temp.grid, num.timesteps =  num.timesteps, time.date =  time.date, 
                                                           k=1, regular.grid = regular.grid )
              names(one.timestep.temp.data)[names(one.timestep.temp.data)=="met.value"]="temp"
              one.timestep.data = merge(one.timestep.data,one.timestep.temp.data,by=c("date","station"))
            }
            if (bHeatFlux) {
              one.timestep.HF.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                         met.parameter.grid =  heatflux.grid, num.timesteps =  num.timesteps, time.date =  time.date, 
                                                         k=1, regular.grid = regular.grid )
              names(one.timestep.HF.data)[names(one.timestep.HF.data)=="met.value"]="HEAT_FLUX"
              one.timestep.data = merge(one.timestep.data,one.timestep.HF.data,by=c("date","station"))
            }
            if (bBL) {
              one.timestep.BL.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                         met.parameter.grid =  BL.height.grid, num.timesteps =  num.timesteps, time.date =  time.date, 
                                                         k=1, regular.grid= regular.grid )
              names(one.timestep.BL.data)[names(one.timestep.BL.data)=="met.value"]="BLAYER_HGT"
              one.timestep.data = merge(one.timestep.data,one.timestep.BL.data,by=c("date","station"))
            }
            if (bCloud) {
              one.timestep.Cl.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                         met.parameter.grid =  cloud.grid, num.timesteps =  num.timesteps, time.date =  time.date, 
                                                         k=1, regular.grid= regular.grid )
              names(one.timestep.Cl.data)[names(one.timestep.Cl.data)=="met.value"]="CLOUD_COVER"
              one.timestep.data = merge(one.timestep.data,one.timestep.Cl.data,by=c("date","station"))
            }
            if (bPr) {
              one.timestep.Pr.data = find.met.grid.value(stations=stations, grid.xcoords =  grid.xcoords, grid.ycoords =  grid.ycoords,
                                                         met.parameter.grid =  precipitation.grid, num.timesteps =  num.timesteps, time.date =  time.date, 
                                                         k=1, regular.grid= regular.grid )
              names(one.timestep.Pr.data)[names(one.timestep.Pr.data)=="met.value"]="Pr"
              one.timestep.data = merge(one.timestep.data,one.timestep.Pr.data,by=c("date","station"))
            }
            
            if (!bDoneFirstTimestep){
               met.data=one.timestep.data
               bDoneFirstTimestep=TRUE
            } else {
               met.data=rbind(met.data,one.timestep.data)
            }
         }
      }
      # from the u and v wind components, calculated wind speed and wind direction 
      met.data$ws=sqrt(met.data$wind.v.comp^2 + met.data$wind.u.comp^2)
      met.data$wd=270-atan2(met.data$wind.v.comp,met.data$wind.u.comp)*180/pi
      if (any(met.data$wd>360,na.rm=TRUE)) {
         #Subtract 360 degrees from any wind directions greater than 360
         met.data$wd[!is.na(met.data$wd) & met.data$wd>360]=met.data$wd[!is.na(met.data$wd) & met.data$wd>360]-360
      }
      if (any(met.data$wd<0,na.rm=TRUE)) {
         #Add 360 degrees to any wind directions less than zero
         met.data$wd[!is.na(met.data$wd) & met.data$wd<0]=met.data$wd[!is.na(met.data$wd) & met.data$wd<0]+360
      }
      # Convert temp from K to 0C
      if (bTemp) met.data$temp=met.data$temp-273.15
      return(list("met.data"=met.data,"bOneGridpoint"=bOneGridpoint))
   }
        
   #function for reading WRF netCDF files as met data format
   #function dependencies: extract.met(), cerc.lambert.cc.xy(), stations data frame
   #packages required: ncdf4, GEOmap, fields, akima, dplyr
   read.wrf.met=function(wrf.file.path){ 
      cat(paste("   Importing WRF met data file: ",wrf.file.path,"...\n", sep=""))
      #Open nc file and return an error if file is unreadable
      wrf.met.nc = tryCatch(nc_open(wrf.file.path),error=function(e) return(NULL))
      if (is.null(wrf.met.nc)) {
         stop(paste("WRF met data file cannot be opened, please check the file path and file format."))
         quit("no")
      }
      # Extract map projection attributes
      map.proj=ncatt_get(wrf.met.nc,0, "MAP_PROJ")
      if (!map.proj$hasatt){
         stop("Cannot find 'MAP_PROJ' global attribute in WRF met file.\nPlease check file format.")
         quit("no")
      } else {
         map.proj=map.proj$value
      }
      # Set projection code to projection name
      if (map.proj==1){
         WRF.proj="LCC"
      } else if (map.proj==2){
         WRf.proj="PS"
      } else if (map.proj==4){
         WRf.proj="UTM"
      } else {
         stop("Map projection of WRF met file is not supported. Only Lambert Conformal Conic, Polar Stereographic and Universal Transverse Mercator are supported.\nPlease check file format.")
         quit("no")
      }
      # Check other required attributes are present for map projection 
      if(WRF.proj=="LCC"){ #Lambert Conformal Conic
         lon0=ncatt_get(wrf.met.nc, 0, "STAND_LON")
         if (!lon0$hasatt){
            stop("Cannot find the attribute 'STAND_LON' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lon0=lon0$value
         }
         lat0=ncatt_get(wrf.met.nc, 0, "MOAD_CEN_LAT")
         if (!lat0$hasatt){
            stop("Cannot find the attribute 'MOAD_CEN_LAT' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lat0=lat0$value
         }
         lat1=ncatt_get(wrf.met.nc, 0, "TRUELAT1")
         if (!lat1$hasatt){
            stop("Cannot find the attribute 'TRUELAT1' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lat1=lat1$value
         }
         lat2=ncatt_get(wrf.met.nc, 0, "TRUELAT2")
         if (!lat2$hasatt){
            stop("Cannot find the attribute 'TRUELAT2' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lat2=lat2$value
         }
        
      } else if (WRF.proj=="PS"){ #Polar Stereographic
         lon0=ncatt_get(wrf.met.nc, 0, "STAND_LON")
         if (!lon0$hasatt){
            stop("Cannot find the attribute 'STAND_LON' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lon0=lon0$value
         }
         lat1=ncatt_get(wrf.met.nc, 0, "TRUELAT1")
         if (!lat1$hasatt){
            stop("Cannot find the attribute 'TRUELAT1' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lat1=lat1$value
         }
         lat0=ncatt_get(wrf.met.nc, 0, "POLE_LAT")
         if (!lat0$hasatt){
            stop("Cannot find the attribute 'POLE_LAT' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lat0=lat0$value
         }
      } else if (WRF.proj=="UTM"){ #Univeral Transverse Mercator
         lon0=ncatt_get(wrf.met.nc, 0, "STAND_LON")
         if (!lon0$hasatt){
            stop("Cannot find the attribute 'STAND_LON' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lon0=lon0$value
         }
         lat0=ncatt_get(wrf.met.nc, 0, "MOAD_CEN_LAT")
         if (!lat0$hasatt){
            stop("Cannot find the attribute 'MOAD_CEN_LAT' in the WRF met file required for setting the projection.\nPlease check file format.")
            quit("no")
         } else {
            lat0=lat0$value
         }
      }
      
      #Extract the required data from the WRF netCDF file
      var.names=names(wrf.met.nc$var)
      if (!"XLAT" %in% var.names){
         stop("The variable 'XLAT' cannot be found in the WRF met file.\nPlease check file format.")  
         quit("no")
      } else{
         XLAT=ncvar_get(wrf.met.nc, "XLAT")
      } 
      if (!"XLONG" %in% var.names){
         stop("The variable 'XLONG' cannot be found in the WRF met file.\nPlease check file format.")  
         quit("no")
      } else{
         XLONG=ncvar_get(wrf.met.nc, "XLONG")
      } 
      if (!"V10" %in% var.names){
         stop("The variable 'V10' cannot be found in the WRF met file.\nPlease check file format.")  
         quit("no")
      } else{
         wind.v=ncvar_get(wrf.met.nc, "V10")
      } 
      if (!"U10" %in% var.names){
         stop("The variable 'U10' cannot be found in the WRF met file.\nPlease check file format.")  
         quit("no")
      } else{
         wind.u=ncvar_get(wrf.met.nc, "U10")
      } 
      if (!"Times" %in% var.names){
         stop("The variable 'Times' cannot be found in the WRF met file.\nPlease check file format.")  
         quit("no")
      } else{
         times=ncvar_get(wrf.met.nc, "Times")
      } 
      # Check for temperature variable in WRF met file and extract if present
      if (!"T2" %in% var.names){
        temp=NULL
        bTemp=FALSE  
      } else{
        temp=ncvar_get(wrf.met.nc, "T2")
        bTemp=TRUE
      }
      # Check for heat flux variable in WRF met file and extract if present
      if (!"HFX" %in% var.names){
        heatflux=NULL
        bHeatFlux=FALSE  
      } else{
        heatflux=ncvar_get(wrf.met.nc, "HFX")
        bHeatFlux=TRUE
      }
      # Check for boundary layer height variable in WRF met file and extract if present
      if (!"PBLH" %in% var.names){
        BL.height=NULL
        bBL=FALSE  
      } else{
        BL.height=ncvar_get(wrf.met.nc, "PBLH")
        bBL=TRUE
      } 
      # convert time sequence to POSIXct class
      times=as.POSIXct(times, format="%Y-%m-%d_%H:%M:%S", tz="GMT")# assume tz of GMT
      num.timesteps=length(times)
      if (num.timesteps>1) {
         #Just pick the grid locations for the first timestep
         if (is.null(dim(XLAT))) {
            #Only one grid point, because XLAT is a list not a matrix
            xlat.for.proj=XLAT[1]
            xlon.for.proj=XLONG[1]
         } else if (length(dim(XLAT))==2) {
            #Must be only one row or one column because matrix only has two dimensions
            xlat.for.proj=XLAT[,1]
            xlon.for.proj=XLONG[,1]
         } else if (length(dim(XLAT))==3) {
            xlat.for.proj=XLAT[,,1]
            xlon.for.proj=XLONG[,,1]
         }
      } else {
         #Only one timestep so use all the grid location data
         xlat.for.proj=XLAT
         xlon.for.proj=XLONG
      }
      
      # set projection 
      #Underlying grid for WRF is regular orthogonal so create map projection for conversion 
      #How the map projection is created depends on whether the WRF file projection is LCC, PS or UTM
      if (WRF.proj=="LCC" ){ #If map projection is LCC
         #Create map projection using MOAD_CEN_LAT for latitude of origin, STAND_LON for longitude of origin and TRUELAT1 and TRUELAT2 for the standard parallels
         #Type = 3 for LCC
         map.proj=setPROJ(type=3, LAT0=lat0,LON0=lon0, LAT1 = lat1, LAT2 = lat2)
         #Create the projected regular grid from the XLAT and XLONG variables
         grid.xy = cerc.lambert.cc.xy(xlat.for.proj,xlon.for.proj,map.proj) #Extract 1 time [,,1](1 matrix slice) as the x and y grid points do not vary with time
         #Extract x and y grid coordinates
         grid.xcoords=grid.xy$x
         grid.ycoords=grid.xy$y
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      } else if (WRF.proj=="PS"){ #If map projection is polar stereographic 
         #Create map projection using STAND_LON for Longitude of Origin, LAT1 as the true parallel and POLE_LAT as Latitude of origin as it is polar stereographic 
         #Type = 4 for PS
         map.proj=setPROJ(type=4, LON0=lon0, LAT1=lat1, LAT0 =lat0)
         grid.xy = stereo.sphr.xy(XLAT[,,1],XLONG[,,1],map.proj)
         #Extract x and y grid coordinates
         grid.xcoords=grid.xy$x
         grid.ycoords=grid.xy$y
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=stereo.sphr.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      } else if (WRF.proj=="UTM"){ #If map projection is UTM
         #Convert STAND_LON which gives UTM zone number to LON0 required for creating the projection
         lon.0= 6.0 * (lon0 - 31.0) + 3.0
         #Use Lon.0 as longitude of origin for projection
         #Type = 7 for utm.wgs84
         map.proj=setPROJ(type=7, LON0=lon.0)
         #get projected grid xy coordinates from the latitude and longitude variables and map projection 
         grid.xy = utm.wgs84.xy(XLAT[,,1],XLONG[,,1],map.proj)
         #Extract x and y grid coordinates
         grid.xcoords=grid.xy$x
         grid.ycoords=grid.xy$y
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=utm.wgs84.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
         allstations.location = matrix(data=c(stations.xy$x,stations.xy$y),nrow=length(stations$name), ncol=2)
      }
      
      #Check for regular or irregular grid as this determines method of finding the grid value 
      num.dimensions.x = length(dim(grid.xcoords))
      num.dimensions.y = length(dim(grid.ycoords))
      #Check dimensions of x are equal to dimensions of y and stop if not
      if (num.dimensions.x==num.dimensions.y) {
         num.dimensions = num.dimensions.x
      } else {
         stop(paste("WRF met variables XLAT and XLONG do not have the same number of dimensions",sep=""))
         quit("no")
      } 
      #If a regular grid, x and y have 1 dimension, if irregular x and y have 2 dimensions. If only one grid point then num.dimensions=0.
      if (num.dimensions %in% c(0,1)) {
         regular.grid = TRUE
      } else if(num.dimensions==2) {
         regular.grid = FALSE
      } else {
         #If x and y have more than 2 dimensions, stop as this is not compatible 
         stop(paste("WRF met variables XLAT and XLONG have ",num.dimensions," dimensions but only 1 or 2 are allowed",sep=""))
         quit("no")
      }
      
      #Extract the met data for the required grid points
      extract.met.result=extract.met(grid.xcoords,grid.ycoords,regular.grid,times,wind.v,wind.u,stations,
                                     bTemp,temp,bHeatFlux,heatflux,bBL,BL.height,FALSE,NULL,FALSE,NULL)
      
      return(list("met.data"=extract.met.result$met.data,"stations"=stations,"bOneGridpoint"=extract.met.result$bOneGridpoint))
   }
    
   # Function to read in ADMS meteorological data
   # ADMS met data is gridded so a met value is retrieved via taking the grid value for each station for each hour
   # function dependencies: extract.met(), cerc.lambert.cc.xy(), cerc.WGS84.to.OSGB(), stations data frame
   # packages required: ncdf4, GEOmap, lubridate, dplyr
   read.ADMS.nc.met=function(ADMS.nc.met.path){
      cat(paste("   Importing ADMS generic gridded met file: ",ADMS.nc.met.path,"...\n", sep=""))
      #Open nc file
      met.nc = tryCatch(nc_open(ADMS.nc.met.path),error=function(e) return(NULL))
      #Check that nc file has opened and if not return error message to log 
      if (is.null(met.nc)){
         stop(paste("ADMS generic gridded met netCDF file cannot be opened, please check the file path and file format."))
         quit("no")
      }
       
      # Check the file version attribute
      version.att=ncatt_get(met.nc, 0, "File_version")
      if (version.att$hasat){
         ADMS.nc.met.version=version.att$value
      } else {
         stop("The global attribute 'File_version' cannot be found in the ADMS generic gridded met netCDF file and is required for processing.\nPlease check file format")
      }
      if(ADMS.nc.met.version != "ADMS_3DMF_v1.0"){
         stop(paste("The ADMS generic gridded met netCDF file contains an unexpected value for the File_version attribute: ",ADMS.nc.met.version, ". Only the following value is allowed: ADMS_3DMF_v1.0." ))
         quit("no")
      }
       
      # Check required variables present   
      var.names=names(met.nc$var)
      if(!"ADMS_grid_mapping" %in% var.names){
         stop("The variable 'ADMS_grid_mapping' cannot be found in the ADMS generic gridded met netCDF file.\nPlease check the file format.")
         quit("no")
      }
      mapping.type=ncatt_get(met.nc, "ADMS_grid_mapping","grid_mapping_name")
      # Check the file is in an allowed coordinate system
      if (! mapping.type$hasatt){
         stop("The attribute: grid_mapping_name on the variable ADMS_grid_mapping cannot be found.\nPlease check the file format.")
         quit("no")
      } else {
         mapping.type=mapping.type$value
      }
      allowed.ADMSmetnc.proj=c("lambert_conformal_conic","transverse_mercator","polar_stereographic","osgb")
      if (!mapping.type %in% allowed.ADMSmetnc.proj){
         stop(paste("The grid_mapping_name attribute in the ADMS generic gridded met netCDF file is not supported.\n",
                     "Only 'osgb', 'lambert_conformal_conic', 'polar_stereographic' and 'universal_transverse_mercator' are supported.\n",
                     "Please check the file format."))
         quit("no")
      }
      if (mapping.type %in% "osgb"){
         projection.type="OSGB"
      } else if (mapping.type %in% "lambert_conformal_conic"){
         projection.type="LCC"
      } else if (mapping.type %in% "transverse_mercator"){
         projection.type="UTM"
      } else if (mapping.type %in% "polar_stereographic"){
         projection.type="PS"
      }
       
      # Check the required projection parameters are present as attributes on the grid mapping variable 
      # Note no projection parameters are required for OSGB
      # For LCC, require LAT0, LON0, LAT1, LAT2, FE, FN
      missing.atts=""
      if (projection.type=="LCC"){
         lat0.att=ncatt_get(met.nc,"ADMS_grid_mapping","latitude_of_projection_origin")
         if (lat0.att$hasatt){LAT0=lat0.att$value} else {missing.atts=paste(missing.atts, "latitude_of_projection_origin")}
         lon0.att=ncatt_get(met.nc,"ADMS_grid_mapping","longitude_of_central_meridian")
         if (lon0.att$hasatt){LON0=lon0.att$value} else {missing.atts=paste(missing.atts, "longitude_of_central_meridian")}
         standardpar.att=ncatt_get(met.nc,"ADMS_grid_mapping","standard_parallel")
         if (standardpar.att$hasatt){
            LAT1=standardpar.att$value[1]
            LAT2=standardpar.att$value[2]
         } else {missing.atts=paste(missing.atts,"standard_parallel")}
         FE.att=ncatt_get(met.nc,"ADMS_grid_mapping","false_easting")
         if (FE.att$hasatt) {FE=FE.att$value}else{missing.atts=paste(missing.atts,"FE")}
         FN.att=ncatt_get(met.nc, "ADMS_grid_mapping", "false_northing")
         if (FN.att$hasatt) {FN=FN.att$value}else{missing.atts=paste(missing.atts,"FN")}
      } 
      # For PS, require LAT0, LON0, FE, FN 
      if (projection.type=="PS"){
         lat0.att=ncatt_get(met.nc,"ADMS_grid_mapping","latitude_of_projection_origin")
         if (lat0.att$hasatt){LAT0=lat0.att$value} else {missing.atts=paste(missing.atts, "latitude_of_projection_origin")}
         lon0.att=ncatt_get(met.nc,"ADMS_grid_mapping","straight_vertical_longitude_from_pole")
         if (lon0.att$hasatt){LON0=lon0.att$value} else {missing.atts=paste(missing.atts, "straight_vertical_longitude_from_pole")}
         FE.att=ncatt_get(met.nc,"ADMS_grid_mapping","false_easting")
         if (FE.att$hasatt) {FE=FE.att$value}else{missing.atts=paste(missing.atts,"FE")}
         FN.att=ncatt_get(met.nc, "ADMS_grid_mapping", "false_northing")
         if (FN.att$hasatt) {FN=FN.att$value}else{missing.atts=paste(missing.atts,"FN")}
         sf.att=ncatt_get(met.nc,"ADMS_grid_mapping","scale_factor_at_projection_origin")
         if (sf.att$hasatt){sf=sf.att$value} # sf may be supplied but not required
      }
       
      # For UTM, require longitude of central meridian
      if (projection.type=="UTM"){
         UTM.att=ncatt_get(met.nc,"ADMS_grid_mapping","longitude_of_central_meridian")
         if (UTM.att$hasatt){UTM=UTM.att$value} else {missing.atts=paste(missing.atts, "longitude_of_central_meridian")}
      }
       
      # Stop if any of the required projection attributes are missing
      if (!nchar(missing.atts==0)){
         stop(paste("The following required projection attributes are missing: ", missing.atts, sep=""))
         quit("no")
      }
    
      # Create projections and project stations onto the projected grid 
      if (projection.type=="LCC"){
         # Set up the map projection: type = 3 for LCC
         map.proj=setPROJ(type=3, LAT0=LAT0,LON0=LON0, LAT1 = LAT1, LAT2 = LAT2, FE=FE, FN=FN)
         #Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=cerc.lambert.cc.xy(stations$latitude,stations$longitude,map.proj)
         #Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
      } else if (projection.type=="PS"){
         #Set up the map projection: type = 4 for PS, depends on whether standard parallel or scale factor at projected origin has been supplied
         map.proj=setPROJ(type=4, LAT0 =LAT0, LON0 = LON0, FE=FE, FN=FN)
         #Get projected map positions of the stations
         stations.xy=stereo.sphr.xy(stations$latitude, stations$longitude, map.proj)
         # If a scale factor has been provided
         if (sf.att$hasatt){
            # Apply scale factor
            stations.xy$x=(stations.xy$x)*sf
            stations.xy$y=(stations.xy$y)*sf
         }
         #Store projected positions in stations table and in station locations matrix
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
      } else if (projection.type=="UTM"){
         LON0=UTM
         # Set up the map projection: type = 7 for utm.wgs84
         map.proj=setPROJ(type=7, LON0=LON0)
         # Project stations onto the projected grid using lat long from stations file and map projection created 
         stations.xy=utm.wgs84.xy(stations$latitude,stations$longitude,map.proj)
         # Add projected grid position columns to the station definition file with the projected x and y coordinates for each station
         stations$map.pos.x = stations.xy$x
         stations$map.pos.y = stations.xy$y  
      } else if (projection.type=="OSGB"){
         # Project stations onto the projected grid using lat long from stations file  
         stations.xy=cerc.WGS84.to.OSGB(stations$latitude, stations$longitude)
         # Add projected grid position columns to the station definition file with the projected x and y coordinates for each station 
         stations$map.pos.x = stations.xy$X
         stations$map.pos.y = stations.xy$Y  
      }
       
      # Check time x and y variables (hidden variables)
      # Extract time sequence
      time.steps=tryCatch(ncvar_get(met.nc, "time"),error=function(e) return(NULL))
      if (is.null(time.steps)){
         stop("Cannot find the time variable in the ADMS generic gridded met netCDF file.\nPlease check the file format.")
         quit("no")
      }
       
      # create time sequence from time step 
      time.unit=met.nc$dim$time$units
      origin.date=gsub("days since ", "", time.unit)
      #convert to POSIXct format, using origin date for jday to cal date
      time.date=as.POSIXct(86400*time.steps,origin=origin.date,format="%Y-%m-%d %H:%M",tz="GMT")
      time.date=round_date(time.date, unit="hours") #round to nearest hour, round up so hour ending 
      num.timesteps=length(time.date)
      
      # Extract x coordinates
      grid.x.coords=tryCatch(ncvar_get(met.nc, "x"),error=function(e) return(NULL))
      if (is.null(grid.x.coords)){
         stop(paste("Cannot find the x grid dimension variable in the ADMS generic gridded met netCDF file.\n",
               "Please check the file format."))
         quit("no")
      }
      num.x.grids=length(grid.x.coords)
       
      # Extract y coordinates
      grid.y.coords=tryCatch(ncvar_get(met.nc, "y"),error=function(e) return(NULL))
      if (is.null(grid.y.coords)){
         stop("Cannot find the y grid dimension variable in the ADMS generic gridded met netCDF file.\nPlease check the file format.")
         quit("no")
      }
      num.y.grids=length(grid.y.coords)
       
      # Extract u wind component
      if(!"WIND_U10" %in% var.names){
         stop("The variable WIND_U10 cannot be found in the ADMS generic gridded met netCDF file.\nPlease check the file format.")
         quit("no")
      }
      wind_u10=ncvar_get(met.nc, "WIND_U10")
       
      # Extract v wind component
      if(!"WIND_V10" %in% var.names){
         stop("The variable WIND_V10 cannot be found in the ADMS generic gridded met netCDF file.\nPlease check the file format.")
         quit("no")
      }
      wind_v10=ncvar_get(met.nc, "WIND_V10")
       
      # Extract temperature
      if ("TEMP_T2" %in% var.names){
        bTemp=TRUE
        temp=ncvar_get(met.nc,"TEMP_T2")
      } else {
        bTemp=FALSE
        temp=NULL
      }
      # Extract heat flux
      if ("HEAT_FLUX" %in% var.names){
        bHeatFlux=TRUE
        heatflux=ncvar_get(met.nc,"HEAT_FLUX")
      } else {
        bHeatFlux=FALSE
        heatflux=NULL
      }
      
      # Extract planetary boundary layer 
      if ("BLAYER_HGT" %in% var.names){
        bBL=TRUE
        BL.height=ncvar_get(met.nc,"BLAYER_HGT")
      } else {
        bBL=FALSE
        BL.height=NULL
      }
      
      # Extract total cloud cover 
      if ("CL_COVER" %in% var.names){
        bCloud=TRUE
        cloudcover=ncvar_get(met.nc,"CL_COVER")
        # Convert from fraction to oktas
        cloudcover=cloudcover*8
        if (max(cloudcover>8) | min(cloudcover<0)){
          cat("   Cloud cover values are invalid. Ensure that cloud fractions are between 0 and 1.\n")
          bCloud=FALSE
          cloudcover=NULL
        }
      } else {
        bCloud=FALSE
        cloudcover=NULL
      }
      
      # Extract precipitation flux 
      if ("PRECIPITATION" %in% var.names){
        bPr=TRUE
        precipitation=ncvar_get(met.nc,"PRECIPITATION")
      } else {
        bPr=FALSE
        precipitation=NULL
      }
      
      # Check for regular or irregular grid 
      num.dimensions.x = length(dim(grid.x.coords))
      num.dimensions.y = length(dim(grid.y.coords))
      #Check dimensions of x are equal to dimensions of y and stop if not
      if (num.dimensions.x==num.dimensions.y) {
         num.dimensions = num.dimensions.x
      } else {
         stop(paste("ADMS generic gridded met netCDF variables x and y do not have the same number of dimensions",sep=""))
         quit("no")
      } 
      #If a regular grid, x and y have 1 dimension, if irregular x and y have 2 dimensions.
      if (num.dimensions==1) {
         regular.grid = TRUE
      } else if(num.dimensions==2) {
         regular.grid = FALSE
      } else {
         #If x and y have more than 2 dimensions, stop as this is not compatible 
         stop(paste("ADMS generic gridded met netCDF met variables x and y ",num.dimensions," dimensions but only 1 or 2 are allowed",sep=""))
         quit("no")
      }
       
      #Extract the met data for the required grid points
      extract.met.result=extract.met(grid.x.coords,grid.y.coords,regular.grid,time.date,wind_v10,wind_u10,stations,
                                     bTemp,temp,bHeatFlux,heatflux,bBL,BL.height,bCloud,cloudcover,bPr,precipitation)
      
      return(list("met.data"=extract.met.result$met.data,"stations"=stations,"bOneGridpoint"=extract.met.result$bOneGridpoint))
   }
    
    
   #Function to make a list of file paths based on a fixed directory, a template path and a start and finish time.
   #If the template path is empty, then this function returns all files in the fixed dir with the extension file.ext.
   make.file.list = function(fixed.dir,bUseTemplatePath,template.path,file.ext,start.datetime,finish.datetime) {
       
      if (!bUseTemplatePath) {
          
         #No template path specified, so use all the files in the fixed directory
         path.list = list.files(path=fixed.dir,pattern=file.ext,full.names=TRUE)
          
         #Check that there is at least one file in the path.list
         if (length(path.list)==0) {
            stop("There are no files of the specified type in the specified folder.")
            quit("no")
         }
          
      } else {
          
         #Convert start and finish times to date-time format from YYYYMMDDHH
         start.datetime = as.POSIXct(start.datetime,format="%Y%m%d%H",tz="GMT")
         if (is.na(start.datetime)) {
            stop("The start date and time has been provided in the wrong format; the correct format is YYYYMMDDHH.")
            quit("no")
         }
         finish.datetime = as.POSIXct(finish.datetime,format="%Y%m%d%H",tz="GMT")
         if (is.na(finish.datetime)) {
            stop("The finish date and time has been provided in the wrong format; the correct format is YYYYMMDDHH.")
            quit("no")
         }
          
         #Check that the finish date is later than the start date
         if (finish.datetime < start.datetime) {
            stop("The finish date and time must be later than the start date and time.")
            quit("no")
         }
          
         #Look in the template path for year (%Y), month (%M), day (%D), julian day (%J), hour (%h), minute (%m) and seconds (%s) 
         bYear.in.template = grepl("%Y",template.path,fixed = TRUE)
         bMonth.in.template = grepl("%M",template.path,fixed = TRUE)
         bDay.in.template = grepl("%D",template.path,fixed = TRUE)
         bJulianDay.in.template = grepl("%J",template.path,fixed = TRUE)
         bHour.in.template = grepl("%h",template.path,fixed = TRUE)
          
         if (bHour.in.template) {
            #Construct hourly time series
            time.series = seq.POSIXt(from=start.datetime,to=finish.datetime,by="hour")
         } else if (bDay.in.template | bJulianDay.in.template) {
            #Construct daily time series
            time.series = seq.Date(from=as.Date(start.datetime),to=as.Date(finish.datetime),by="day")
         } else if (bMonth.in.template) {
            #Construct monthly time series
            time.series = seq.Date(from=as.Date(start.datetime),to=as.Date(finish.datetime),by="month")
         } else if (bYear.in.template) {
            #Construct yearly time series
            time.series = seq.Date(from=as.Date(start.datetime),to=as.Date(finish.datetime),by="year")
         } else {
            #Didn't find any of the template date-time tags
            stop("Template file path does not contain any recognised date-time tags. Use %Y for 4-digit year, %M for 2-digit month, %D for 2-digit day, %J for 3-digit julian day and %h for 2-digit hour.")
            quit("no")
         }
          
         #Convert to a data frame
         time.series = data.frame(date=time.series)
         #Initialise the filepath to the template path
         time.series$filepath = template.path
         #Search for date-time placeholders and replace with appropriate strings
         if (bYear.in.template) time.series = mutate(time.series,filepath = str_replace_all(filepath,"%Y",format(date,format="%Y")))
         if (bMonth.in.template) time.series = mutate(time.series,filepath = str_replace_all(filepath,"%M",format(date,format="%m")))
         if (bDay.in.template) time.series = mutate(time.series,filepath = str_replace_all(filepath,"%D",format(date,format="%d")))
         if (bJulianDay.in.template) time.series = mutate(time.series,filepath = str_replace_all(filepath,"%J",format(date,format="%j")))
         if (bHour.in.template) time.series = mutate(time.series,filepath = str_replace_all(filepath,"%h",format(date,format="%H")))
          
         time.series$filepath.full = file.path(fixed.dir,time.series$filepath)
         time.series$filepath.full.exists = file.exists(time.series$filepath.full)
         time.series = subset(time.series,filepath.full.exists)
          
         path.list = time.series$filepath.full
          
         #Check that there is at least one file in the path.list
         if (length(path.list)==0) {
            stop("There are no files in the specified folder that match the specified template file path.")
            quit("no")
         }
          
      }
       
      return(path.list)
   }
}
# 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
   update.packages.ADD = update.packages
   import.model.name.ADD = import.model.name
   pollutants.mod.file.ADD = pollutants.mod.file
   bSpeciesMapping.ADD=bSpeciesMapping
   species.mapping.file.ADD=species.mapping.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
   bUseTemplatePath.ADD=bUseTemplatePath
   template.path.ADD=template.path
   start.datetime.ADD=start.datetime
   finish.datetime.ADD=finish.datetime
   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
   stations.interp.ADD=stations.interp

   
  
   #load workspace
   cat("Loading workspace...\n")
   data.loaded=tryCatch(load(existing.workspace.file.ADD),error=function(e) return(NULL))
   if (is.null(data.loaded)) {
      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" & names(all.data)!="network" & names(all.data)!="auto.download"]
   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
   update.packages = update.packages.ADD
   import.model.name = import.model.name.ADD
   pollutants.mod.file = pollutants.mod.file.ADD
   bSpeciesMapping=bSpeciesMapping.ADD
   species.mapping.file=species.mapping.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
   bUseTemplatePath = bUseTemplatePath.ADD
   template.path=template.path.ADD
   start.datetime=start.datetime.ADD
   finish.datetime=finish.datetime.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
   stations.interp=stations.interp.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")
allowed.networks = c("kcl","aurn","aqe","waqn","saqn","ni","europe","aqmesh")
allowed.avg.time.values=c("hour","day","week","month","year")

gridded.modelled.dataset.formats = c("IOAPI","WRF-Chem","EMEP","CAMS Regional ADS","CAMS Regional Reanalysis","Generic netCDF") #These are the modelled data formats that need latitude and longitude data in the stations file
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","avg.time","statistic")
   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) {
         #Allow the column ozone (optional column)
         if (!tolower(names(pollutants)[j])=="ozone") {
            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")
   }
   #If Ozone column is present check that it is logical if not add column ozone
   if (any(tolower(names(pollutants)=="ozone"))){
      pollutants$ozone=tolower(as.character(pollutants$ozone))
      if (any(! pollutants$ozone %in% c("yes","no"))){
         stop("The pollutant definitions file must contain either 'yes' or 'no' for Ozone.")
      }
   } else {
      pollutants$ozone="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")
     } 
   }
   
   #Perform checks on the pollutants averaging times and statistics
   reformatted.avg.time.statistic = check_and_reformat_averaging_times_and_statistics(pollutants$avg.time,pollutants$statistic,pollutants$pollutant,"")
   #Update using the re-formatted averaging times and statistics
   pollutants$avg.time = reformatted.avg.time.statistic$avg.time
   pollutants$statistic = reformatted.avg.time.statistic$statistic

   #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==gridded.modelled.dataset.formats))|bmod.dbm) {
      stn.expected.col.names = c("name","type","latitude","longitude","network","auto.download")
   } else {
      stn.expected.col.names = c("name","type","network","auto.download")
   }
   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")
      }
   }
   station.col.names=names(stations)
   # Check that if gridded met data is being used, the stations dataframe has a latitude and longitude column 
   if (bMetData & MetDataFormat %in% c("ADMSmetnc","wrf")){
      if (!("latitude" %in% station.col.names) | !("longitude" %in% station.col.names)){
         stop("When using a gridded met data format (e.g. ADMS generic gridded met data netCDF or WRF), the monitoring station definitions file must contain latitude and longitude data.\n")
         quit("no")
      }
   }
   
   #Checks if the stations data contains latitude and longitude column and outputs a warning message if not.
   bStationsLatLong=TRUE
   if (!("latitude" %in% station.col.names) | !("longitude" %in% station.col.names)){
      cat("   Warning: The monitoring station definitions file does not contain latitude and longitude data, therefore only non-gridded modelled datasets can be added to this workspace in future\n")
      bStationsLatLong=FALSE
   }
   rm(station.col.names) # No longer needed
   
   #Save input station names for use on output graphs
   stations$station.output.name = trimws(as.character(stations$name))
   #Make sure station names for internal use are in lower case with no leading/trailing spaces
   stations$name = trimws(tolower(as.character(stations$name)))
   stations$type = tolower(as.character(stations$type))
   stations$network=trimws(tolower(as.character(stations$network)))
   stations$auto.download=trimws(tolower(as.character(stations$auto.download)))
   #check auto.download inputs are either yes or no
   if(!(all(unique(stations$auto.download)%in% c("yes","no")))){
      stop(paste("Autodownload column must be either: yes  no"))
      quit("no")
   }
   #convert to logical inputs
   stations$auto.download[which(stations$auto.download=="yes")]=TRUE
   stations$auto.download[which(stations$auto.download=="no")]=FALSE
   stations$auto.download=as.logical(stations$auto.download)
   #check any networks for autodownloaded stations are in list of allowed networks
   if(any(!(unique(stations$network[which(stations$auto.download==TRUE)]) %in% allowed.networks))){
      stop(paste("Unrecognised network for auto-download station. Allowed networks for autodownload are: ",paste(allowed.networks,collapse="  "),"\n"))
      quit("no")
   }
   #if any stations are to be downloaded automatically, pick up the built in meta data file
   #if all stations are autodownloaded, don't use the csv file
   bUseNetworkFile=FALSE
   bUseCSVFile=TRUE
   if(any(stations$auto.download)){
      bUseNetworkFile=TRUE
      pollutants.obs.network.file=paste(install_dir, "PollutantsNetworkMetaData.csv", sep="/")
      if(all(stations$auto.download)){
         bUseCSVFile=FALSE
         #if all are automatically downloaded, ignore the input CSV if there is one given
         if(bobs.csv){
            cat("   All stations will be automatically downloaded, the input observed pollutants data file will be ignored \n")
            bobs.csv=FALSE
         }
      }
   } else if (!bobs.csv){
      #if there is no observed data from the csv and no stations are auto-download then there is no data to download
      stop("No observed data has been supplied. Either change the auto-download column in the station file or select to import from a CSV file")
      quit("no")
   }
   
  
   #READ AND PROCESS OBS DATA POLLUTANT DEFINITION FILE (separator here is ",")
   #header = "pollutant,network, obs.alias, obs.units, obs.avg.time, obs.statistic, hour.ending
   #if the observed data comes from a CSv file, read in pollutants.obs file from input file
   if(bUseCSVFile){
      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","network","obs.avg.time","obs.statistic","hour.ending")
      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")
         }
      }
      #Checking that each network in the stations definitions file is in the observed pollutant definitions file
      stations.networks = tolower(unique(stations$network[stations$auto.download==FALSE]))
      invalid.networks = which(!(stations.networks %in% tolower(unique(pollutants.obs$network))))
      if (length(invalid.networks)!=0){ #If there are any invalid stations then output an error
        stop(paste("Monitoring station definitions file contains networks not found in the observed pollutant definitions file:",paste(stations.networks[invalid.networks], collapse = ', '),"\n"))
        quit("no")
      }
      rm(stations.networks) #Removing data that is no longer needed
      rm(invalid.networks)
    
      #subset only the networks that are in the stations file
      pollutants.obs$network=trimws(tolower(as.character(pollutants.obs$network)))
      pollutants.obs=subset(pollutants.obs,network %in% unique(stations$network))
      if(length(pollutants.obs$network)==0){
         stop("Pollutant definitions file does not contain information for the networks in the station file")
         quit("no")
      }
      
      #Perform checks on the observed pollutants averaging times and statistics
      reformatted.avg.time.statistic = check_and_reformat_averaging_times_and_statistics(pollutants.obs$obs.avg.time,pollutants.obs$obs.statistic,pollutants.obs$pollutant,"observed")
      #Update using the re-formatted averaging times and statistics
      pollutants.obs$obs.avg.time = reformatted.avg.time.statistic$avg.time
      pollutants.obs$obs.statistic = reformatted.avg.time.statistic$statistic
      
      #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(unique(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")
         } 
      }
      
      #check the obs.aliases are the same for any networks that are not automatically downloaded
      for (j in 1:length(unique(pollutants.obs$pollutant))) {
         if(length(unique(pollutants.obs$obs.alias[which(pollutants.obs$pollutant==unique(pollutants.obs$pollutant[j]))]))>1){
            stop(paste("Obs pollutant aliases must be the same for all networks supplied in the CSV input file" ))
            quit("no")
         }
      }
      #check the obs statistics are the same for any networks that are not automatically downloaded
      for (j in 1:length(unique(pollutants.obs$pollutant))) {
         if (length(unique(pollutants.obs$obs.statistic[which(pollutants.obs$pollutant==unique(pollutants.obs$pollutant)[j])]))>1) {
            stop(paste("Obs statistics for",toupper(pollutants.obs$pollutant[j]),"must be the same across all the networks that are included in the input pollutant obs file"))
            quit("no")
         } 
      }
      #check the obs averaging times are the same for any networks that are not automatically downloaded
      for (j in 1:length(unique(pollutants.obs$pollutant))) {
         if (length(unique(pollutants.obs$obs.avg.time[which(pollutants.obs$pollutant==unique(pollutants.obs$pollutant)[j])]))>1) {
            stop(paste("Obs averaging time for",toupper(pollutants.obs$pollutant[j]),"must be the same across all the networks that are included in the input pollutant obs file"))
            quit("no")
         } 
      }
      #make sure the hour ending is either yes or no
      pollutants.obs$hour.ending=trimws(tolower(as.character(pollutants.obs$hour.ending)))
      if(!(all(unique(pollutants.obs$hour.ending)%in% c("yes","no")))){
         stop(paste("Hour Ending column must be either: yes  no"))
         quit("no")
      }
      #make sure hour ending conventions are the same for all pollutants within each network
      available.obs.networks=unique(pollutants.obs$network)
      for(i in 1:length(available.obs.networks)){
         if(length(unique(pollutants.obs$hour.ending[which(pollutants.obs$network==available.obs.networks[i])]))>1){
            stop(paste("Hour Ending column must be same for all pollutants within each network. Check network: ",available.obs.networks[i]))
            quit("no")
         }
      }
      #convert to logical inputs
      pollutants.obs$hour.ending[which(pollutants.obs$hour.ending=="yes")]=TRUE
      pollutants.obs$hour.ending[which(pollutants.obs$hour.ending=="no")]=FALSE
      pollutants.obs$hour.ending=as.logical(pollutants.obs$hour.ending)
      
   }
   
   #if some of the data is automatically downloaded, read in fixed network file
   if (bUseNetworkFile){
      cat("Importing auto-download monitoring network pollutant details...\n")
      pollutants.obs.network=read.csv(pollutants.obs.network.file,header=TRUE,sep=",")
      names(pollutants.obs.network) = tolower(names(pollutants.obs.network))
      pol.obs.expected.col.names = c("pollutant","obs.alias","obs.units","network","obs.avg.time","obs.statistic","hour.ending")
      for (j in 1:length(pol.obs.expected.col.names)) {
         if (length(which(names(pollutants.obs.network)==pol.obs.expected.col.names[j]))==0) {
            stop(paste("Missing column in network pollutant definitions file: ",pol.obs.expected.col.names[j]))
            quit("no")
         }
      }
      for (j in 1:length(names(pollutants.obs.network))) {
         if (length(which(pol.obs.expected.col.names==names(pollutants.obs.network)[j]))==0) {
            stop(paste("Network pollutant definitions file contains unexpected column: ", names(pollutants.obs.network)[j]))
            quit("no")
         }
      }
      #subset only the networks that are in the stations file
      pollutants.obs.network$network=trimws(tolower(as.character(pollutants.obs.network$network)))
      pollutants.obs.network=subset(pollutants.obs.network,network %in% unique(stations$network))
      if(length(pollutants.obs.network$network)==0){
         stop("Network pollutant definitions file does not contain information for the networks in the station file")
         quit("no")
      }

      #Convert pollutant names from factors to strings
      pollutants.obs.network$pollutant=tolower(as.character(pollutants.obs.network$pollutant))
      pollutants.obs.network$obs.alias=tolower(as.character(pollutants.obs.network$obs.alias))
      
      #Error handling:
      #Make sure that the units entered are allowed
      pollutants.obs.network$obs.units = tolower(as.character(pollutants.obs.network$obs.units))
      for (j in 1:length(unique(pollutants.obs.network$pollutant))) {
         if (length(which(allowed.units==pollutants.obs.network$obs.units[j]))==0) {
            stop(paste("Pollutant data concentration units for",toupper(pollutants.obs.network$pollutant[j]),"in the network pollutants file must be either ppb, ppm, ng/m3, ug/m3, mg/m3 or g/m3"))
            quit("no")
         } 
      }
      
      #Perform checks on the network pollutants averaging times and statistics
      reformatted.avg.time.statistic = check_and_reformat_averaging_times_and_statistics(pollutants.obs.network$obs.avg.time,pollutants.obs.network$obs.statistic,pollutants.obs.network$pollutant,"network")
      #Update using the re-formatted averaging times and statistics
      pollutants.obs.network$obs.avg.time = reformatted.avg.time.statistic$avg.time
      pollutants.obs.network$obs.statistic = reformatted.avg.time.statistic$statistic
      
      #make sure the hour ending is either yes or no
      pollutants.obs.network$hour.ending=trimws(tolower(as.character(pollutants.obs.network$hour.ending)))
      #check auto.download inputs are either yes or no
      if(!(all(unique(pollutants.obs.network$hour.ending)%in% c("yes","no")))){
         stop(paste("Hour Ending column in the network pollutants file must be either: yes  no"))
         quit("no")
      }
      #make sure hour ending conventions are the same for all pollutants within each network
      available.obs.networks=unique(pollutants.obs.network$network)
      for(i in 1:length(available.obs.networks)){
         if(length(unique(pollutants.obs.network$hour.ending[which(pollutants.obs.network$network==available.obs.networks[i])]))>1){
            stop(paste("Hour Ending column in the network pollutants file must be same for all pollutants within each network. Check network: ",available.obs.networks[i]))
            quit("no")
         }
      }
    
      #convert to logical inputs
      pollutants.obs.network$hour.ending[which(pollutants.obs.network$hour.ending=="yes")]=TRUE
      pollutants.obs.network$hour.ending[which(pollutants.obs.network$hour.ending=="no")]=FALSE
      pollutants.obs.network$hour.ending=as.logical(pollutants.obs.network$hour.ending)
      
      #if the CSV file is input also, combine the two into a larger dataframe. 
      #This dataframe can have multiple rows per pollutant
      if(bUseCSVFile){
         pollutants.obs=rbind(pollutants.obs,pollutants.obs.network)
      } else{
         pollutants.obs=pollutants.obs.network
      }
    
   }
   pollutants.obs=subset(pollutants.obs,pollutant %in% pollutants$pollutant)
   #create a networks dataframe containing the network and the hour-ending convention 
   networks=unique(subset(pollutants.obs,select=c("network","hour.ending")))
   networks.download=unique(subset(stations,select=c("network","auto.download")))
   networks=merge(networks,networks.download)
   #check the pollutants in the observed file are in the pollutants file
   if (length(pollutants.obs$pollutant)==0) {
      stop("The obs pollutant data file does not contain any of the same pollutants as the main pollutants data file")
      quit("no")      
   }
   #Create a column with the pollutant name with any '.' removed (i.e. pm25 not pm2.5)
   pollutants$pollutant.noperiod = gsub("[.]","",pollutants$pollutant)
   
} else {
  #Removed the old modelled data from pollutants data
  pollutants = subset(pollutants,select=c("pollutant", "output.units",  "conv.ugm3.ppb", "min.allowed", "max.allowed","ozone","pollutant.noperiod", "avg.time","statistic"))
} # end if bNewData

#Checking if latitude and longitude data is required in station definitions file
if ((!bStationsLatLong) & any(mod.ncdf.format==gridded.modelled.dataset.formats)){
   stop("You cannot add gridded modelled datasets to this workspace because latitude and longitude were missing from the station definitions file when the workspace was created. To include this modelled dataset you will need to re-run the Data Input Tool to re-create the workspace using a station definitions file that includes latitude and longitude data.\n")
   quit("no")
}

#Load the required packages
packages.required = c("latticeExtra","plyr","nlme","openair","tibble","reshape2","lubridate","data.table","dplyr","stringr")
if (bmod.ncdf|bmod.dbm|(bMetData & MetDataFormat %in% c("ADMSmetnc","wrf"))) packages.required = c(packages.required,"ncdf4","GEOmap","fields","maps","akima","RANN","dplyr","stringr","sf")
install.load.packages(packages.required)

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

#Check stations dataframe for duplicates stations
if (!nrow(stations)==nrow(unique(stations))){
   stations.freq=as.data.table(table(stations$name))
   stations.not.unique=subset(stations.freq, !(N == 1))
   stations.not.unique=paste(stations.not.unique$V1, collapse=", ")
   stop(paste("The stations definitions file contains duplicated stations for stations: ",stations.not.unique,".\nPlease removed duplicated stations. \n", sep=""))
   quit("no")
}

#READ AND PROCESS MOD DATA POLLUTANT DEFINITION FILE (separator here is ",")
#header = "pollutant, mod.alias, mod.units, mod.avg.time, 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))
#If no species mapping, modelled pollutant definition file requires pollutant, mod.alias mod.units, mod.avg.time and mod.statistic columns
if(!bSpeciesMapping){
   pol.mod.expected.col.names = c("pollutant","mod.alias","mod.units","mod.avg.time","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")
      }
   }
} else if (bSpeciesMapping){ #If species mapping is happening, no mod.units is required 
   pol.mod.expected.col.names = c("pollutant","mod.avg.time","mod.statistic","mod.alias") # check for required column headers
   if(!(length(pol.mod.expected.col.names[which(!(pol.mod.expected.col.names %in% names(pollutants.mod)))])==0)){
      stop(paste("Missing column in mod data pollutant definitions file: ", pol.mod.expected.col.names[which(!(pol.mod.expected.col.names %in% names(pollutants.mod)))], sep=""))
      quit("no")
   }
   
   #send warning message to log if modelled pollutants definition file contains mod.units column to say it won't be used
   if ("mod.units" %in% names(pollutants.mod)){
      warning("The mod data pollutant definitions file contains mod.units. These will not be used as the conversion factors in the species mapping file will be used.\n")
   }
   
   #set mod.units to be ug/m3 for the pollutants.mod as all units will be converted to ug/m3 after mapping (also note column is required for checking unexpected columns)
   pollutants.mod$mod.units="ug/m3"
   
   #add mod.units to expected pol mod definition column names 
   pol.mod.expected.col.names=c(pol.mod.expected.col.names,"mod.units")
}

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

#Perform checks on the modelled pollutants averaging times and statistics
reformatted.avg.time.statistic = check_and_reformat_averaging_times_and_statistics(pollutants.mod$mod.avg.time,pollutants.mod$mod.statistic,pollutants.mod$pollutant,"modelled")
#Update using the re-formatted averaging times and statistics
pollutants.mod$mod.avg.time = reformatted.avg.time.statistic$avg.time
pollutants.mod$mod.statistic = reformatted.avg.time.statistic$statistic

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

# If the modelled data format is Defra Background maps, ensure that for all pollutants:
# 1. the mod.avg.time is year 
# 2. the the mod.statistic is mean. 
if (bmod.dbm){
   if (!all(pollutants.mod$mod.avg.time=="year" | pollutants.mod$mod.avg.time=="1 year")){
      stop("For Defra Background Maps modelled data, only annual mean datasets are supported, however the mod.avg.time settings in the modelled pollutants definitions file are not all 'year'.\n")
      quit("no")
   }
   if (!all(pollutants.mod$mod.statistic=="mean")){
      stop("For Defra Background Maps modelled data, only annual mean datasets are supported, however the mod.statistic settings in the modelled pollutants definitions file are not all 'mean'.\n")
      quit("no")
   }
}

#Process species mapping file
if (bSpeciesMapping){
   #Read in species mapping file 
   species.mapping.data=readLines(species.mapping.file)
   #Remove trailing commas from the end of lines 
   species.mapping.data[1]=gsub(",","",species.mapping.data[1])
   #Check that species mapping file is species map version 2 (compatible format)
   if (!(species.mapping.data[1]=="SPECIESMAPVERSION2")){
      stop(paste("The species mapping file version ",species.mapping.data[1]," is unsupported. SPECIESMAPVERSION2 is required for the species mapping data file.",sep=""))
      quit("no")
   }
   #Remove trailing commas from the end of lines 
   species.mapping.data[2]=gsub(",","",species.mapping.data[2])
   if (!(species.mapping.data[2]=="ROWSTOCOLUMNS")){
     stop(paste("The species mapping file format ",species.mapping.data[2]," is unsupported. ROWSTOCOLUMNS is required for the species mapping data file.",sep=""))
     quit("no")
   }
   #Identify which line the data starts on, this is marked by the term "variables" at the beginning of the line
   mapping.data.start.line=which(grepl("variables",species.mapping.data,ignore.case=TRUE))
   #Re-read in species mapping file as a CSV, starting at the line the data starts on and convert to data frame 
   species.mapping.data=read.csv(species.mapping.file, skip=mapping.data.start.line-1)
   species.mapping.data=as.data.frame(species.mapping.data)
   #Extract the first column which is the names of the species mapping components
   component.species=species.mapping.data[,1]
   # Fix NAs for sodium interpreted as NA, assume there can be no actual NA's in species mapping file so change to "Na"
   # If any NAs in mapping file exist, change any NAs in definition file to "Na"
   if(any(is.na(component.species))){
      component.species[which(is.na(component.species))]="Na"
      if(any(is.na(pollutants.mod$mod.alias))){
         pollutants.mod$mod.alias[which(is.na(pollutants.mod$mod.alias))]="Na"
      }
      if(any(is.na(pollutants.mod$pollutant))){
         pollutants.mod$pollutant[which(is.na(pollutants.mod$pollutant))]="Na"
      }
   }
   #convert column and row names to lower case and remove the variable column and give rows their component names
   rownames(species.mapping.data)=tolower(component.species)
   names(species.mapping.data)=tolower(names(species.mapping.data))
   #Remove the column with the component names on now that these have been given to the row names 
   species.mapping.data=species.mapping.data[,-which(names(species.mapping.data)=="variables")]
   
   # Check that all pollutants exist in pollutant mapping file
   if (any(!(pollutants$pollutant %in% names(species.mapping.data)))){
      missing.pollutants=pollutants$pollutant[which(!(pollutants$pollutant %in% names(species.mapping.data)))]
      stop(paste("The following output pollutants are missing from the species mapping file: ",paste(missing.pollutants, collapse=", "), sep=""))
      quit("no")
   }
   
   #Check that species mapping contains all required pollutants/species and that definition files contain all required pollutants/species
   #subset species mapping file for just pollutants required in the output
   species.mapping.data=subset(species.mapping.data, select=(pollutants$pollutant))
   # Remove rows where all columns are 0 (un-mapped components so not required) from species mapping data frame 
   # Add a temporary row to sum all conversion factors across the row 
   species.mapping.data$sum=apply(species.mapping.data,1,sum) 
   # Check if any have summed to NA which will cause errors later 
   if (any(is.na(species.mapping.data$sum))){
      stop("The species mapping file must not contain NA for conversion factors.")
   }
   species.mapping.data=species.mapping.data[which(species.mapping.data$sum!=0),]# Remove the rows where the sum of conversion factors is 0
   if (nrow(species.mapping.data)==0){
      stop("Cannot evaluate species mapping file. Please ensure valid mappings are entered.")
      quit("no")
   }
   species.mapping.data=select(species.mapping.data, selct=-c(sum)) # Remove the sum column 
      
   
   #Check that all species required for output pollutants exist in pollutant mod definition file
   if (any(!(row.names(species.mapping.data) %in% pollutants.mod$mod.alias))){
      missing.species=row.names(species.mapping.data)[which(!(row.names(species.mapping.data)%in% pollutants.mod$mod.alias))]
      stop(paste("The following component species used in the species mapping file are missing in the mod pollutants definition file: ",paste(missing.species, collapse = ", "), sep=""))
      quit("no")
   }
   
   #subset pollutants.mod to only include required components 
   pollutants.mod=subset(pollutants.mod, mod.alias %in% row.names(species.mapping.data))
   
   #check that the output statistic and averaging time of the components are the same 
   for (i in 1:length(pollutants$pollutant)){#Loop through output pollutants
      output.pol=pollutants$pollutant[i] 
      one.output.pol.mapping=subset(species.mapping.data, select=output.pol)#Subset species.mapping.data for that output pollutant 
      #one.output.pol.mapping=one.output.pol.mapping[-which(one.output.pol.mapping[[output.pol]]==0),]
      one.output.pol.species.data=subset(pollutants.mod, mod.alias %in% row.names(one.output.pol.mapping)[which(!one.output.pol.mapping[[output.pol]]==0)])
      if (!length(unique(one.output.pol.species.data$mod.avg.time))==1|!length(unique(one.output.pol.species.data$mod.statistic)==1)){
         stop(paste("The species components for an output pollutant must have the same averaging time and statistic. \n The species components for pollutant ", output.pol, " have different output averaging times/statistics in the modelled pollutants definition file.", sep=""))
         quit("no")
      } else {
         #Set the mod.avg.time and mod.statistic for the output pollutants
         pollutants$mod.avg.time[i]=unique(one.output.pol.species.data$mod.avg.time)
         pollutants$mod.statistic[i]=unique(one.output.pol.species.data$mod.statistic)
      }
   }
   #set units for pollutants$mod.units to be "ug/m3" which will be the case after applying the conversion factor in the species mapping file
   pollutants$mod.units="ug/m3"
} else{
   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")   
}

#Check that the observed pollutants metadata contains all required networks for each required pollutant before proceeding and issue a warning if not.
#Not all networks may offer all pollutants.
networks.required = unique(stations$network)
networks.required.nondownload = unique(stations$network[!stations$auto.download])
for (j in 1:length(pollutants$pollutant)) {
  for (i in 1:length(networks.required)) {
    if (!(networks.required[i] %in% pollutants.obs$network[pollutants.obs$pollutant==pollutants$pollutant[j]])) {
      if (networks.required[i] %in% networks.required.nondownload) {
        #Any networks being imported from CSV file need metadata for all pollutants, because all pollutants are imported for all stations.
        stop(paste("Missing observed metadata for network",networks.required[i],"for pollutant",pollutants$pollutant[j],"in the observed pollutants data file.\n"))
        quit("no")
      } else {
        #Any networks being downloaded are allowed to have missing metadata for one or more pollutants, because stations
        #are downloaded individually and different networks might measure different pollutants.
        cat(paste("   Warning: Missing observed metadata for network",networks.required[i],"for pollutant",pollutants$pollutant[j],"\n"))
      }
    }
  }  
}

#Check that obs and mod averaging times and statistics are compatible with the required output averaging time and statistic
#For a mod/obs averaging time and statistic to be compatible with the output averaging time and statistic, they must be:
#a.	The same; or
#b.	The mod/obs avg time is hourly mean; or
#c.	Both means, with the mod/obs averaging time a factor of the main averaging time.
for (j in 1:length(pollutants$pollutant)){
  bOK = avg_compatibility_check(pollutants$pollutant[j],"modelled",pollutants$mod.avg.time[j],pollutants$mod.statistic[j],
                                     "output",pollutants$avg.time[j],pollutants$statistic[j])
  if (!bOK) {
    stop(paste0("Averaging time compatibility checks failed between modelled and output data for pollutant ",toupper(pollutants$pollutant[j]),"."))
    quit("no")
  }
  #Checking observed averaging times and statistics for each required network
  for(k in 1:length(networks$network)) {
    pollutants.onenetwork = subset(pollutants.obs, pollutant==pollutants$pollutant[j] & network==networks$network[k])
    if (length(pollutants.onenetwork$pollutant)==0) next #not all networks may provide all pollutants so skip to next network
    if (length(pollutants.onenetwork$pollutant)>1) {
      #Shouldn't happen, but rest of code assumes length 1 so better put in a sanity check to make sure!
      stop(paste0("More than one row in the observed pollutants data for network ",networks$network[k]," for pollutant ",pollutants$pollutant[j],"."))
      quit("no")
    }
    bOK.network = avg_compatibility_check(pollutants$pollutant[j],
                                               "observed",pollutants.onenetwork$obs.avg.time[1],pollutants.onenetwork$obs.statistic[1],
                                               "output",pollutants$avg.time[j],pollutants$statistic[j])
    if (!bOK.network) {
      stop(paste0("Averaging time compatibility checks failed between observed and output data for pollutant ",toupper(pollutants$pollutant[j])," and network ",networks$network[k],"."))
      quit("no")
    }
  }
}

if (bSpeciesMapping){
   pollutants.out=pollutants
   pollutants=pollutants.mod
}

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

if (bmod.ncdf) {
   if (mod.ncdf.format=="WRF-Chem"){
      #read in modelled data if format is WRF Chem
      if (bmod.ncdf.select.file){
         modelled.data=read.WRFChem.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if(bmod.ncdf.select.dir){
         WRFChem.files=make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                        start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in(1:length(WRFChem.files))){
            if(!bDoneFirstFile){
               modelled.data=read.WRFChem.nc(WRFChem.files[i])
               bDoneFirstFile=TRUE
            } else {
               nc.data=read.WRFChem.nc(WRFChem.files[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
      
   } else if (mod.ncdf.format=="EMEP"){
      #read in modelled data if format is EMEP
      if (bmod.ncdf.select.file) {
         modelled.data=read.EMEP.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if(bmod.ncdf.select.dir){
         EMEP.files=make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                      start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in(1:length(EMEP.files))){
            if(!bDoneFirstFile){
               modelled.data=read.EMEP.nc(EMEP.files[i])
               bDoneFirstFile=TRUE
            } else {
               nc.data=read.EMEP.nc(EMEP.files[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
      
   } else if (mod.ncdf.format=="IOAPI"){
      #read in the modelled data if format is IOAPI
      if (bmod.ncdf.select.file) {
         modelled.data=read.IOAPI.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if(bmod.ncdf.select.dir){
         IOAPI.files=make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                   start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in(1:length(IOAPI.files))){
            if(!bDoneFirstFile){
               modelled.data=read.IOAPI.nc(IOAPI.files[i])
               bDoneFirstFile=TRUE
            } else {
               nc.data=read.IOAPI.nc(IOAPI.files[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
   
   } else if (mod.ncdf.format=="ADMS output netCDF"){
      #read in modelled data if modelled data format is ADMS output ncdf
      if (bmod.ncdf.select.file) {
         modelled.data=read.ADMSoutput.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if (bmod.ncdf.select.dir){
         ADMS.files = make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                     start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in (1:length(ADMS.files))){
            if (!bDoneFirstFile) {
               modelled.data=read.ADMSoutput.nc(ADMS.files[i])
               bDoneFirstFile=TRUE
            } else{
               nc.data=read.ADMSoutput.nc(ADMS.files[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
   
   } else if (mod.ncdf.format=="Generic netCDF"){
       #read in the modelled data if format is Generic netCDF
      if (bmod.ncdf.select.file) {
         modelled.data=read.Generic.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if(bmod.ncdf.select.dir){
         Genericnc.files = make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                     start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in(1:length(Genericnc.files))){
            if(!bDoneFirstFile){
               modelled.data=read.Generic.nc(Genericnc.files[i])
               bDoneFirstFile=TRUE
            } else {
               nc.data=read.Generic.nc(Genericnc.files[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
   } else if (mod.ncdf.format=="CHIMERE"){
     #read in the modelled data if format is Generic netCDF
     if (bmod.ncdf.select.file) {
         modelled.data=read.CHIMERE.nc(mod.ncdf.file, chimere.cco.file)
         bDoneFirstFile=TRUE
     }
     if(bmod.ncdf.select.dir){
         CHIMERE.files = make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                        start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in(1:length(CHIMERE.files))){
            if(!bDoneFirstFile){
               modelled.data=read.CHIMERE.nc(CHIMERE.files[i], chimere.cco.file)
               bDoneFirstFile=TRUE
            } else {
               nc.data=read.CHIMERE.nc(CHIMERE.files[i],chimere.cco.file)
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
   } else if (mod.ncdf.format=="CAMS Regional ADS"){
      #read in modelled data if modelled data format is CAMS Regional ncdf
      if (bmod.ncdf.select.file) {
         modelled.data=read.CAMSRegionalADS.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if (bmod.ncdf.select.dir){
         CAMSRegionalADSFiles = make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                     start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in (1:length(CAMSRegionalADSFiles))){
            if (!bDoneFirstFile) {
               modelled.data=read.CAMSRegionalADS.nc(CAMSRegionalADSFiles[i])
               bDoneFirstFile=TRUE
            } else{
               nc.data=read.CAMSRegionalADS.nc(CAMSRegionalADSFiles[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
   
   } else if (mod.ncdf.format=="CAMS Regional Reanalysis"){
      #read in modelled data if modelled data format is CAMS Regional Reanalysis ncdf
      if (bmod.ncdf.select.file) {
         modelled.data=read.CAMSRegionalReanalysis.nc(mod.ncdf.file)
         bDoneFirstFile=TRUE
      }
      if (bmod.ncdf.select.dir){
         CAMSRegionalReanalysisFiles = make.file.list(fixed.dir = mod.ncdf.dir, bUseTemplatePath, template.path, file.ext = ".nc",
                                               start.datetime, finish.datetime)
         bDoneFirstFile=FALSE
         for (i in (1:length(CAMSRegionalReanalysisFiles))){
            if (!bDoneFirstFile) {
               modelled.data=read.CAMSRegionalReanalysis.nc(CAMSRegionalReanalysisFiles[i])
               bDoneFirstFile=TRUE
            } else{
               nc.data=read.CAMSRegionalReanalysis.nc(CAMSRegionalReanalysisFiles[i])
               modelled.data=merge.nc.data.modelled.data(nc.data,modelled.data)
            }
         }
      }
   }
  
   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.dbm){
  #Defra Background Maps 
   if (bmod.dbm.select.file ) {
      modelled.data=read.DefraBackgroundMaps.csv(mod.dbm.file)
      bDoneFirstFile=TRUE
   }
   if(bmod.dbm.select.dir ){
      dbm.csv.files = make.file.list(fixed.dir = mod.dbm.dir, bUseTemplatePath, template.path, file.ext = ".csv",
                                   start.datetime, finish.datetime)
      bDoneFirstFile=FALSE
      for (i in(1:length(dbm.csv.files))){
         if(!bDoneFirstFile){
            modelled.data=read.DefraBackgroundMaps.csv(dbm.csv.files[i])
            bDoneFirstFile=TRUE
         } else {
            dbm.data=read.DefraBackgroundMaps.csv(dbm.csv.files[i])
            modelled.data=merge.nc.data.modelled.data(dbm.data,modelled.data)
         }
      }
   }
   if (!bDoneFirstFile) {
      stop("No data for any required pollutants in the Defra Background Maps 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 Defra Background Maps data.",
                "Please remove this pollutant from the pollutant definitions file.\n"))
      }
   }
   if (bstop) {
      stop("Missing pollutant(s) in Defra Background Maps 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)
      #Make sure all the column headers are in lower case
      names(modelled.data) = tolower(names(modelled.data))
      #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(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 folder
      modelled.data = lapply(make.file.list(fixed.dir = mod.csv.dir, bUseTemplatePath, template.path, file.ext = ".csv",
                                            start.datetime, finish.datetime),
                             open.csv.and.add.station,csv.separator=mod.csv.separator)
      modelled.data = do.call(rbind.fill,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 = which(is.na(modelled.data[,pollutant.data.names[pol]]))
      } else {
         missing.val.ids = which(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
         #Replace 'missing value' text/number with NA
         modelled.data[missing.val.ids,pollutant.data.names[pol]] = NA
         #Find indices of all NA values (this will include missing pollutants)
         na.ids = which(is.na(modelled.data[,pollutant.data.names[pol]]))
         #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[na.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[na.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 folder
      pst.files = make.file.list(fixed.dir = mod.pst.dir, bUseTemplatePath, template.path, file.ext = ".pst",
                                start.datetime, finish.datetime)
      pst.files.nonzero = pst.files[file.info(pst.files)$size > 0] # Remove PST files with zero size from the list
      modelled.data = lapply(pst.files.nonzero,
                         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
}
#Make sure station names are in lower case with no leading/trailing spaces  
modelled.data$station = trimws(tolower(as.character(modelled.data$station)))

#Check no duplicates in modelled data of station and date combinations
modelled.data.unique.check=subset(modelled.data, select=c("date","station"))
if (!(nrow(modelled.data.unique.check)==nrow(unique(modelled.data.unique.check)))){
   modelled.data.unique.check$date_station=paste(modelled.data.unique.check$date," for station: ",modelled.data.unique.check$station,".", sep="")
   modelled.data.freq=as.data.table(table(modelled.data.unique.check$date_station))
   duplicate.ref=which(!(modelled.data.freq$N == 1))
   cat("The modelled data contains duplicated station date time data for the following station-date-time combinations: \n")
   for (i in 1:length(duplicate.ref)){
      duplicate.data=as.character(modelled.data.freq[duplicate.ref[i],1])
      cat(paste("   ", duplicate.data,"\n",sep=""))
   }
   stop("Correct the duplicated modelled data.\n")
   quit("no")
}

#Offsetting the modelled data based on the time difference between modelled and observed data
if (mod.time.diff != 0){ #Skipping the case when time difference is 0 to avoid unnecessary processing
   modelled.data$date=modelled.data$date+3600*mod.time.diff
   #Updating the year, month, day and hour columns for the new dates
   modelled.data$year=year(modelled.data$date)
   modelled.data$month=month(modelled.data$date)
   modelled.data$day=day(modelled.data$date)
   modelled.data$hour=hour(modelled.data$date)
}


# If species mapping is true, convert component species pollutants into output pollutants
if (bSpeciesMapping){
   cat("   Doing species mapping...\n")
   modelled.data.output.pol=subset(modelled.data, select=c("date","station"))
   for (i in 1:length(pollutants.out$pollutant)){ #loop through output pollutants
      one.pol.out=pollutants.out$pollutant[i]
      species.mapping.data.one.pol=subset(species.mapping.data, select=one.pol.out) # Subset species mapping data to get data for the one output pollutant
      component.species=row.names(species.mapping.data.one.pol)[which(!species.mapping.data.one.pol[[one.pol.out]]==0)] # Use the row names to get the names of the component species 
      component.species=pollutants$pollutant[which(pollutants$mod.alias %in% component.species)] # Convert mod.alias in species definition file to the pollutant name used in header
      cat(paste("      Mapping component species: ",paste(component.species, collapse = ", "), " to output pollutant: ",one.pol.out, "...\n",sep="" )) # This is using the pollutants name not mod.alias in log message
      modelled.data.component.species=subset(modelled.data, select=component.species) # Subset modeled data data frame to get the component pollutants 
      for (j in 1:ncol(modelled.data.component.species)){# Loop through the columns of the component pollutants and multiply by the conversion factor  
         single.component.species=names(modelled.data.component.species)[j]
         single.component.species.mod.alias=pollutants$mod.alias[which(pollutants$pollutant==single.component.species)]
         conversion.factor=species.mapping.data.one.pol[[one.pol.out]][which(row.names(species.mapping.data.one.pol)==single.component.species.mod.alias)]
         modelled.data.component.species[[component.species[j]]]=modelled.data.component.species[[component.species[j]]]*conversion.factor
      }
      modelled.data.component.species[[one.pol.out]]=apply(modelled.data.component.species,1,sum) # Sum the total of all components into a new column named with the output pollutant
      modelled.data.output.pol=cbind(modelled.data.output.pol,subset(modelled.data.component.species, select=one.pol.out))# Join summed output pollutant col to the stations/date/any already processed pollutants 
   }
   modelled.data=modelled.data.output.pol # Reset the modelled.data data frame to be the processed modelled data pollutant outputs  
   pollutants=pollutants.out
   rm(modelled.data.output.pol)
}
#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 & bobs.csv) {
   #OBSERVED DATA
   #Observed data can either be all in one file or in all files in one folder.
   cat("Importing observed data...\n")
   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)
      #Make sure all the column headers are in lower case
      names(all.obs.data) = tolower(names(all.obs.data))
      #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(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 folder
      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)
   }
   
   #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" ]
   #Removing any pollutant.data.names that aren't valid
   new.pol.data.names = c() #This will contain the valid pollutants
   for (i in 1:length(pollutant.data.names)){
      col.index=which(colnames(all.obs.data) == pollutant.data.names[i])
      if (!is.numeric(all.obs.data[,col.index])){ #All columns with valid pollutant data are numeric
         all.obs.data = all.obs.data[-c(col.index)] #Removes the column from observed data if it is not a valid pollutant
         cat(paste("  ",pollutant.data.names[i],"is not a valid pollutant column and will be ignored.\n"))
      }else{
         new.pol.data.names=c(new.pol.data.names,pollutant.data.names[i])
      }
   }
   pollutant.data.names=new.pol.data.names
   rm(new.pol.data.names) #Removing variable now it is no longer needed
   
   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 = which(is.na(all.obs.data[,pollutant.data.names[pol]]))
      } else {
         missing.val.ids = which(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
         #Replace 'missing value' text/number with NA
         all.obs.data[missing.val.ids,pollutant.data.names[pol]] = NA           
         #Find indices of all NA values (this will include missing pollutants)
         na.ids = which(is.na(all.obs.data[,pollutant.data.names[pol]]))
         #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[na.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[na.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 with no leading/trailing spaces
   all.obs.data$station = trimws(tolower(as.character(all.obs.data$station)))
   #Change pollutant column headers to the main pollutant name
   for (j in 1:length(pollutants$pollutant)) {
      #extract obs alias from pollutant obs data frame
      #the obs alias was checked that it was unique in the observed csv (despite multiple networks being permitted)
      obs.alias=unique(pollutants.obs$obs.alias[which(pollutants.obs$pollutant==pollutants$pollutant[j] & !(pollutants.obs$network %in% allowed.networks))])
      this.pol.col = which(names(all.obs.data)==obs.alias)
      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 data are to be read in from online networks this is done later inside the loop over stations.
   #Check no duplicates in observed data of station and date combinations
   obs.data.unique.check=subset(all.obs.data, select=c("year","month","day","hour","station"))
   if (!(nrow(obs.data.unique.check)==nrow(unique(obs.data.unique.check)))){
      obs.data.unique.check$date_station=paste(obs.data.unique.check$year,"-",obs.data.unique.check$month,"-",obs.data.unique.check$day," ",obs.data.unique.check$hour,":00 for station: ",obs.data.unique.check$station,".", sep="")
      obs.data.freq=as.data.table(table(obs.data.unique.check$date_station))
      duplicate.ref=which(!(obs.data.freq$N == 1))
      cat("The observed data contains duplicated station date time data for the following station-date-time combinations: \n")
      for (i in 1:length(duplicate.ref)){
         duplicate.data=as.character(obs.data.freq[duplicate.ref[i],1])
         cat(paste("   ",duplicate.data,"\n",sep=""))
      }
      stop("Correct the duplicated observed data.\n")
      quit("no")
   }
   
} #end if bNewData

#Check that the date ranges of the observed and modelled data intersect
mod.start.date = range(modelled.data$date)[1]
mod.end.date=c()
for (i in 1:length(pollutants$pollutant)){
   mod.dates = seq.POSIXt(from=range(modelled.data$date)[2],by=pollutants$mod.avg.time[i],length.out = 2)
   mod.end.date[i]=mod.dates[2]-3600
}
mod.end.date=as.POSIXct(max(mod.end.date),origin = "1970-01-01") #Choosing the latest date as the mod.end.date

networks$start.year=as.numeric(format(mod.start.date,"%Y"))
networks$end.year=as.numeric(format(mod.end.date,"%Y"))
if (bNewData) {
   if (any(stations$auto.download)) {
      for(i in 1:length(networks$network)){
         if(!networks$auto.download[i]) next
         if(!(networks$hour.ending[i])){
            #Find the year required for that network (accounting for the hour-ending convention of that network)
            networks$start.year[i] = as.numeric(format(mod.start.date-3600,"%Y"))
            networks$end.year[i] = as.numeric(format(mod.end.date-3600,"%Y"))
         }
         cat(paste("Observed data will be downloaded for network ",toupper(networks$network[i])," for ",networks$start.year[i]," to ",networks$end.year[i],"\n",sep=""))
      
      }
    
   }
   if(bobs.csv){
      cat("Checking the modelled and observed date ranges overlap...\n")
      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")   
      }
   } 
}


#Loop over all stations and all pollutants
cat("Processing data...\n")
bDoneFirstStation = FALSE
if (bNewData) {
   for (i in 1:length(stations$name)) {
      closeAllConnections()
      log.file = file(paste(output.file.path.root,".log",sep=""),open="at")
      sink(file=log.file,type="output",split=TRUE)
      sink(file=log.file,type="message")
      #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
      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 (!(stations$auto.download[i])) {
         #Get the data from the global data frame
         station.obs.data = subset(all.obs.data,station==stations$name[i])
         # If observed data is hour beginning 
         # add 1 hour to dates to adjust from hour beginning to hour ending
         if (!(networks$hour.ending[which(networks$network==stations$network[i])])) {
            station.obs.data$date = station.obs.data$date + 3600
         }
         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 download from the automatic networks
      } else {
         cat(paste("  Importing ",toupper(stations$network[i])," data for station ",toupper(stations$name[i]),"...\n",sep=""))
         obs.year=seq(from=networks$start.year[which(networks$network == stations$network[i])],to=networks$end.year[which(networks$network == stations$network[i])],by=1)
         #Loop over years required
         bDoneFirstYear = FALSE
         for (k in 1:length(obs.year)) {
            if (stations$network[i]=="kcl") {
               station.obs.data.oneyear = tryCatch(importKCL(site=stations$name[i],year=obs.year[k]),
                                              error = function(e) return(NULL))
            } else if (stations$network[i]=="aurn") {
               station.obs.data.oneyear = tryCatch(importAURN(site=stations$name[i],year=obs.year[k],hc=TRUE),
                                              error = function(e) return(NULL))
            } else if (stations$network[i]=="saqn") {
               station.obs.data.oneyear = tryCatch(importSAQN(site=stations$name[i],year=obs.year[k]),
                                              error = function(e) return(NULL))
            } else if (stations$network[i]=="waqn") {
               station.obs.data.oneyear = tryCatch(importWAQN(site=stations$name[i],year=obs.year[k]),
                                              error = function(e) return(NULL))
            } else if (stations$network[i]=="aqe") {
               station.obs.data.oneyear = tryCatch(importAQE(site=stations$name[i],year=obs.year[k]),
                                              error = function(e) return(NULL))
            } else if (stations$network[i]=="ni") {
               station.obs.data.oneyear = tryCatch(importNI(site=stations$name[i],year=obs.year[k]),
                                              error = function(e) return(NULL))
            } else if (stations$network[i]=="europe") {
               station.obs.data.oneyear = tryCatch(importEurope(site=stations$name[i],year=obs.year[k]),
                                              error = function(e) return(NULL))
            }
            #If there's a problem extracting the data for this year, skip to the next year
            if (is.null(station.obs.data.oneyear) | length(station.obs.data.oneyear$date)==0){
               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 available to download for station",toupper(stations$name[i]),"- skipping this station\n"))
            next
         }
        
         for (j in 1:length(pollutants$pollutant)) {
            #check the pollutant is one provided by the network
            if (!(pollutants$pollutant[j] %in% pollutants.obs$pollutant[which(pollutants.obs$network==stations$network[i])])) {
               cat("  No meta data is available for pollutant ", toupper(pollutants$pollutant[j]), " for network ",stations$network[i], ", therefore no observed data are available for this station for pollutant ", toupper(pollutants$pollutant[j]),"\n",sep="")
               #if data for this pollutant has been downloaded, remove this column
               if(pollutants$pollutant[j] %in% names(station.obs.data)){
                  station.obs.data=subset(station.obs.data,select=names(station.obs.data)[which(names(station.obs.data) != pollutants$pollutant[j])])
               }
               next
            }
            #Change pollutant column headers to the main pollutant name
            nrow.obs.alias=which(pollutants.obs$pollutant==pollutants$pollutant[j] & pollutants.obs$network == stations$network[i])
            this.pol.col = which(names(station.obs.data)==pollutants.obs$obs.alias[nrow.obs.alias])
            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$obs.alias[nrow.obs.alias]),"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
         }
         # If observed data is hour beginning 
         # 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 (!(networks$hour.ending[which(networks$network==stations$network[i])])) {
            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)) {
         #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 there is valid data for the pollutant
         non.missing.obs=obs.data[!is.na(obs.data[,c(pollutants$pollutant[j])]),]
         if (length(non.missing.obs$date) == 0) {
            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"
         #Convert units to output units
         nrow.obs.meta=which(pollutants.obs$pollutant==pollutants$pollutant[j] & pollutants.obs$network==stations$network[i])
         obs.data$obs = convert.units(obs.data$obs,pollutants.obs$obs.units[nrow.obs.meta],
                                      pollutants$output.units[j],pollutants$conv.ugm3.ppb[j])
         #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 the modelled and observed data, get them both to have the output averaging time and statistic
         obs.statistic = pollutants.obs$obs.statistic[nrow.obs.meta]
         mod.statistic = pollutants$mod.statistic[j]
         avg.time = pollutants$avg.time[j]
         statistic = pollutants$statistic[j]
         if (statistic=="rolling mean") {
           if (obs.statistic!="rolling mean") {
             #Can assume that the obs data is hourly mean, otherwise earlier checks would have failed
             obs.data = rollingMean(obs.data,pollutant="obs",width=avg_time_hours(avg.time),new.name="rolling",align="right")
             obs.data = subset(obs.data,select=c("date","rolling"))  
             names(obs.data)[2]="obs"
           }
           if (mod.statistic!="rolling mean") {
             #Can assume that the mod data is hourly mean, otherwise earlier checks would have failed
             mod.data = rollingMean(mod.data,pollutant="mod",width=avg_time_hours(avg.time),new.name="rolling",align="right")
             mod.data = subset(mod.data,select=c("date","rolling"))  
             names(mod.data)[2]="mod"
           }
         } else {
            #Can assume that the obs and mod data are either hourly mean, or the same averaging time and statistic as the output averaging time and statistic,
            #or mean with an averaging time a factor of the output averaging time and the output statistic is mean, or else earlier checks would have failed
            obs.data = timeAverage(obs.data,avg.time=avg.time,statistic=statistic)
            mod.data = timeAverage(mod.data,avg.time=avg.time,statistic=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$station.output.name = stations$station.output.name[i]
         all.data.onestation.onepol$type = stations$type[i]
         all.data.onestation.onepol$network = stations$network[i]
         all.data.onestation.onepol$auto.download = stations$auto.download[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
   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
      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
         }
         #Extract existing obs and mod data from the workspace
         all.data.onestation.onepol = subset(all.data[ ,c("date","station","type","network","auto.download","obs",existing.modelled.data.names)], all.data$station==stations$name[i] & all.data$pollutant==pollutants$pollutant[j])
         #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])
         #Before combining the modelled and existing data, get the modelled data in the output averaging time and statistic
         mod.statistic = pollutants$mod.statistic[j]
         avg.time = pollutants$avg.time[j]
         statistic = pollutants$statistic[j]
         if (statistic=="rolling mean") {
           if (mod.statistic!="rolling mean") {
             #Can assume that the mod data is hourly mean, otherwise earlier checks would have failed
             mod.data = rollingMean(mod.data,pollutant="mod",width=avg_time_hours(avg.time),new.name="rolling",align="right")
             mod.data = subset(mod.data,select=c("date","rolling"))  
             names(mod.data)[2]="mod"
           }
         } else {
           #Can assume that the obs and mod data are either hourly mean, or the same averaging time and statistic as the output averaging time and statistic,
           #or mean with an averaging time a factor of the output averaging time and the output statistic is mean, or else earlier checks would have failed
           mod.data = timeAverage(mod.data,avg.time=avg.time,statistic=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")
   all.data.out = subset(all.data, select=-station.output.name)
   write.csv(all.data.out, file=paste(output.file.path.root,"_raw_data.csv", sep=""), row.names=FALSE)
   rm(all.data.out)
}

if (bMetData){
   cat("Processing meteorological data...\n")
   ## Check extension of met data file is compatible with specified met format 
   if (MetDataFormat == "mop"){ 
      met.extension=tolower(substr(MetDataFilePath, nchar(MetDataFilePath)-3, nchar(MetDataFilePath)))
      if (met.extension != ".mop"){
         stop(paste("The meteorological data file has the extension: ", met.extension,
                    ", however in the interface the following meteorological data format has been specified: .MOP\n", sep=""))
         quit("no")
      }
      met.data = importADMS(MetDataFilePath)
   } else if (MetDataFormat == "met"){
      met.extension=tolower(substr(MetDataFilePath, nchar(MetDataFilePath)-3, nchar(MetDataFilePath)))
      if (met.extension != ".met"){
         stop(paste("The meteorological data file has the extension: ", met.extension,
                    ", however in the interface the following meteorological data format has been specified: .MET\n", sep=""))
         quit("no")
      }
      met.data=importADMS(MetDataFilePath)
   } else if (MetDataFormat == "pro"){
      met.extension=tolower(substr(MetDataFilePath, nchar(MetDataFilePath)-3, nchar(MetDataFilePath)))
      if (met.extension != ".pro"){
         stop(paste("The meteorological data file has the extension: ", met.extension,
                    ", however in the interface the following meteorological data format has been specified: .PRO\n", sep=""))
         quit("no")
      }
      met.data=read.pro(MetDataFilePath)
   } else if (MetDataFormat == "mio"){
      met.extension=tolower(substr(MetDataFilePath, nchar(MetDataFilePath)-3, nchar(MetDataFilePath)))
      if (met.extension != ".mio"){
         stop(paste("The meteorological data file has the extension: ", met.extension,
                    ", however in the interface the following meteorological data format has been specified: .MIO\n", sep=""))
         quit("no")
      }
      met.data=read.mio(MetDataFilePath)
   } else if (MetDataFormat == "wrf"){
      met.extension=tolower(substr(MetDataFilePath, nchar(MetDataFilePath)-2, nchar(MetDataFilePath)))
      if (met.extension != ".nc"){
         stop(paste("The meteorological data file has the extension: ", met.extension,
                    ", however in the interface the following meteorological data format has been specified: .nc (WRF met)\n", sep=""))
         quit("no")
      }
      read.wrf.met.output=read.wrf.met(MetDataFilePath)
      met.data=read.wrf.met.output$met.data
      stations.met=read.wrf.met.output$stations
      bOneGridpoint = read.wrf.met.output$bOneGridpoint
   } else if (MetDataFormat == "ADMSmetnc"){
      met.extension=tolower(substr(MetDataFilePath, nchar(MetDataFilePath)-2, nchar(MetDataFilePath)))
      if (met.extension != ".nc"){
         stop(paste("The meteorological data file has the extension: ", met.extension,
                    ", however in the interface the following meteorological data format has been specified: .nc (ADMS generic gridded netCDF)\n", sep=""))
         quit("no")
      }
      read.ADMS.nc.met.output=read.ADMS.nc.met(MetDataFilePath)
      met.data=read.ADMS.nc.met.output$met.data
      stations.met=read.ADMS.nc.met.output$stations
      bOneGridpoint = read.ADMS.nc.met.output$bOneGridpoint
   }

   #Check met data date range and all.data.for.plotting date range overlap
   cat("   Checking the meteorological data and Modelled/Observed date ranges overlap...\n")
   met.data.start.date = range(met.data$date, na.rm=TRUE)[1]
   met.data.end.date = range(met.data$date, na.rm=TRUE)[2]
   plot.data.start.date=range(all.data$date, na.rm=TRUE)[1]
   plot.data.end.date=range(all.data$date, na.rm=TRUE)[2]
   if (!branges.intersect(r1start=met.data.start.date, r1end=met.data.end.date,
                          r2start=plot.data.start.date,r2end=plot.data.end.date)) {                     
      cat(paste("Meteorological data date range:",met.data.start.date,"to",met.data.end.date,"\n"))
      cat(paste("Modelled/Observed data date range:",plot.data.start.date,"to",plot.data.end.date,"\n"))
      stop("Date ranges of meteorological data and Modelled/Observed data do not overlap")
      quit("no")   
   }
   
   # write out met data as a ADMS .met file if required
   if (bMetOutputFile){
     cat("   Writing out extracted meteorological parameters...\n")
     met.data.to.write = data.frame(STATION=met.data$station,
                                    YEAR=year(met.data$date),
                                    TDAY=yday(met.data$date),
                                    THOUR=hour(met.data$date),
                                    U=met.data$ws,
                                    PHI=met.data$wd)
     if("temp" %in% names(met.data)) met.data.to.write$T0C = met.data$temp
     if("HEAT_FLUX" %in% names(met.data)) met.data.to.write$FTHETA0 = met.data$HEAT_FLUX
     if("BLAYER_HGT" %in% names(met.data)) met.data.to.write$H = met.data$BLAYER_HGT
     if("CLOUD_COVER" %in% names(met.data)) met.data.to.write$CL = met.data$CLOUD_COVER
     if("Pr" %in% names(met.data)) met.data.to.write$P = met.data$Pr
     
     #replace NA variables with -999
     met.data.to.write$U[which(is.na(met.data.to.write$U))]=-999
     met.data.to.write$PHI[which(is.na(met.data.to.write$PHI))]=-999
     if("T0C" %in% names(met.data.to.write)) met.data.to.write$T0C[which(is.na(met.data.to.write$T0C))] = -999
     if("FTHETA0" %in% names(met.data.to.write)) met.data.to.write$FTHETA0[which(is.na(met.data.to.write$FTHETA0))] = -999
     if("H" %in% names(met.data.to.write)) met.data.to.write$H[which(is.na(met.data.to.write$H))] = -999
     if("CL" %in% names(met.data.to.write)) met.data.to.write$CL[which(is.na(met.data.to.write$CL))] = -999
     if("P" %in% names(met.data.to.write)) met.data.to.write$P[which(is.na(met.data.to.write$P))] = -999
     
     comment.lines=c("=========================================================",
                     "Met parameters extracted from gridded met file:",MetDataFilePath,"",
                     "Station locations:","name, latitude, longitude, projectedX, projectedY",
                     str_replace_all(unlist(strsplit(paste(stations.met$name, stations.met$latitude, 
                                                           stations.met$longitude, stations.met$map.pos.x, 
                                                           stations.met$map.pos.y, collapse = ","),","))," ",", "))
     if (bOneGridpoint) comment.lines=c(comment.lines,"Values for a single gridpoint used for all stations")
     comment.lines=c(comment.lines,
                     "=========================================================")
     num.vars = as.character(ncol(met.data.to.write))
     vars = as.vector(colnames(met.data.to.write))
     header.lines = c("VARIABLES:",num.vars,vars,comment.lines,"","DATA:")
     
     writeLines(header.lines, file(paste(output.file.path.root,"_extracted.met",sep="")))
     fwrite(met.data.to.write, paste(output.file.path.root,"_extracted.met",sep=""), append=TRUE)
   }
   
   
}

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

