#Functions common to Model Evaluation and Model Diagnostics tools
r.common.functions.version = "5.4"
svn.version.common.functions = "$ProjectRevision: SVN 13255 (2026-06-03) $"
svn.prj.name.common.functions = "$ProjectName: /cerc_cerc/Products/Utilities/ModelEvaluationToolkit/Tags/ModelEvaluationToolkit_5_4_1 $"

# Function to read in user defined averaging time and statistic for Model Evaluation tool
# process output averaging and merge with pollutants data frame 
process.output.averaging = function(output.averaging.file, pol)  {
  # Firstly read in the data, checking its format
  output.averaging = read.csv(file=output.averaging.file)
  names(output.averaging) = tolower(names(output.averaging))
  # make sure each column exists
  if (is.null(output.averaging$pollutant) | is.null(output.averaging$output.avg.time) | is.null(output.averaging$output.statistic) | is.null(output.averaging$output.daily.max) ) {
    stop("You have chosen to perform concentration evaluation but you have not included all of the required information in the output averaging times file.")
    quit("no")
  }
  output.averaging$pollutant = tolower(trimws(as.character(output.averaging$pollutant)))
  output.averaging$output.daily.max = tolower(trimws(as.character(output.averaging$output.daily.max)))
  
  #Perform checks on the output averaging times and statistics
  reformatted.avg.time.statistic = check_and_reformat_averaging_times_and_statistics(output.averaging$output.avg.time,output.averaging$output.statistic,output.averaging$pollutant,"output")
  #Update using the re-formatted averaging times and statistics
  output.averaging$output.avg.time = reformatted.avg.time.statistic$avg.time
  output.averaging$output.statistic = reformatted.avg.time.statistic$statistic
  
  #Give a warning if ae.band not included
  if (is.null(output.averaging$ae.band)) {
    cat("WARNING: ae.band is missing from the output averaging times file, therefore AE_band will not be calculated.\n")
  }

  # check each pollutant in the pollutant data frame (which may be a single pollutant) is included in the user defined averaging file
  for (j in 1:length(pol$pollutant)) {   
    if (length(which(output.averaging$pollutant==pol$pollutant[j]))==0 ) {
      stop(paste("Pollutant",toupper(pol$pollutant[j]),"is missing from the output averaging times file."))
      quit("no")
    }
    if (length(which(output.averaging$pollutant==pol$pollutant[j]))>1 ) {
      stop(paste("Pollutant",toupper(pol$pollutant[j]),"is listed more than once in the output averaging times file."))
      quit("no")
    }
  }
  # now merge with pollutants file
  pol = merge(pol,output.averaging,by="pollutant")
  
  #Check that the input averaging times and statistics are compatible with the required output averaging time and statistic
  for (j in 1:length(pol$pollutant)){
    bOK = avg_compatibility_check(pol$pollutant[j],"input",pol$avg.time[j],pol$statistic[j],
                                  "output",pol$output.avg.time[j],pol$output.statistic[j])
    if (!bOK) {
      stop(paste0("Averaging time compatibility checks failed between input and output data for pollutant ",toupper(pol$pollutant[j]),"."))
      quit("no")
    }
  }

  #Make sure that if the daily max is to be calculated, then the averaging time is less than or equal to 
  #24 hours, unless statistic is rolling mean.
  for (j in 1:length(pol$pollutant)){
    pol$output.avg.time.hours[j]=avg_time_hours(pol$output.avg.time[j])
    
    #Check the daily max setting isn't missing and is either yes or no
    if (is.na(pol$output.daily.max[j])) {
      stop(paste("output.daily.max setting for",toupper(pol$pollutant[j]),"is missing"))
      quit("no")      
    } else if (!(pol$output.daily.max[j] %in% c("yes","no"))) {
      stop(paste("output.daily.max setting for",toupper(pol$pollutant[j]),"must be yes or no"))
    }
    #If daily max output of non-rolling averages is selected, then check modelled averaging time isn't greater than 24 hours
    if (pol$output.daily.max[j] == "yes" & pol$output.statistic[j] != "rolling mean") {
      ok.to.continue = TRUE
      if (is.na(pol$output.avg.time.hours[j])) {
        #This will be the case if the averaging time is in months
        ok.to.continue = FALSE
      } else if (pol$output.avg.time.hours[j] > 24) {
        ok.to.continue = FALSE
      }
      if (!ok.to.continue) {
        stop(paste("Averaging time for",toupper(pol$pollutant[j]),"must be equal to or less than a day for the daily max of a non-rolling statistic to be calculated"))
        quit("no")
      }
    }
  }

  return(pol)
}


#function to average data in Model Evaluation tool to user defined output time and statistic
average.data = function(original.dataframe, pol, j, stns, threshold) {
  # applying averaging to subset of data for pollutant and station ID
  if (pol$output.avg.time[j]==pol$avg.time[j] & pol$output.statistic[j]==pol$statistic[j]) {
    data.averaged = original.dataframe
  } else {
    if (pol$output.statistic[j]=="rolling mean") {
      #Can assume obs and mod are both hourly data - anything else would have caused an error during earlier checks
      data.averaged.obs = rollingMean(original.dataframe,pollutant="obs",width=pol$output.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right", type="site")
      data.averaged.obs = subset(data.averaged.obs,select=c("date","site","rolling"))
      names(data.averaged.obs)[names(data.averaged.obs)=="rolling"]="obs"
      data.averaged.mod = rollingMean(original.dataframe,pollutant="mod",width=pol$output.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right", type="site")
      data.averaged.mod = subset(data.averaged.mod,select=c("date","site","rolling"))
      names(data.averaged.mod)[names(data.averaged.mod)=="rolling"]="mod"
      data.averaged = merge(data.averaged.obs,data.averaged.mod)
    } else {
      data.averaged = timeAverage(original.dataframe,avg.time=pol$output.avg.time[j],statistic=pol$output.statistic[j],data.thresh=threshold,type="site",interval=pol$avg.time[j], progress=FALSE)
    } 
  }
  #Convert to daily maximum if required
  if (pol$output.daily.max[j]=="yes" & !(pol$output.avg.time[j] %in% c("day","1 day"))) {
    data.averaged = timeAverage(data.averaged,avg.time="day",statistic="max",data.thresh=threshold,type="site",interval=pol$output.avg.time[j], progress=FALSE)
  }
  # add station type back in by merge, to mimic all.data format
  names(stns)[names(stns) == "station"] = "site"
  data.averaged$site = as.character(data.averaged$site)
  data.averaged = merge(data.averaged, stns, by="site")
  return(data.averaged)  # because of placement in ModelEvaluation code, actually data averaged by pollutant
}

#Process index file and add to pollutants file
#Inputs index scales file and pollutants, outputs pollutants (merged with index scales)
process.index.scales = function(index.scales,pol) {
  #Read and process the index scales file
  indices = read.csv(file=index.scales)  
  #Make sure the index columns exist and there is at least one index band defined for each pollutant
  if (is.null(indices$pollutant) | is.null(indices$index.units) | is.null(indices$index.statistic) | is.null(indices$index.avg.time) | is.null(indices$index.daily.max) | is.null(indices$i1)) {
    stop("You have chosen to assess forecast indices but you have not included all the required forecast index information in the index scales file.")
    quit("no")
  }
  # convert strings to lower case for consistency within R
  indices$pollutant = tolower(trimws(as.character(indices$pollutant)))
  indices$index.units = tolower(trimws(as.character(indices$index.units)))
  indices$index.daily.max = tolower(trimws(as.character(indices$index.daily.max)))
  #Make sure that the units entered are allowed
  for (j in 1:length(indices$pollutant)) {
    if (length(which(allowed.units==indices$index.units[j]))==0) {
      stop(paste("Concentration units for",toupper(indices$pollutant[j]),"must be either ppb, ppm, ug/m3, mg/m3 or g/m3"))
      quit("no")
    } 
  }
  
  #Perform checks on the index averaging times and statistics
  reformatted.avg.time.statistic = check_and_reformat_averaging_times_and_statistics(indices$index.avg.time,indices$index.statistic,indices$pollutant,"index")
  #Update using the re-formatted averaging times and statistics
  indices$index.avg.time = reformatted.avg.time.statistic$avg.time
  indices$index.statistic = reformatted.avg.time.statistic$statistic
  
  # now merge the error-checked indices data with the pollutants data.frame
  pol = merge(pol,indices,by="pollutant")
  
  #Check that the input averaging times and statistics are compatible with the required index averaging time and statistic
  for (j in 1:length(pol$pollutant)){
    bOK = avg_compatibility_check(pol$pollutant[j],"input",pol$avg.time[j],pol$statistic[j],
                                  "index",pol$index.avg.time[j],pol$index.statistic[j])
    if (!bOK) {
      stop(paste0("Averaging time compatibility checks failed between input and index data for pollutant ",toupper(pol$pollutant[j]),"."))
      quit("no")
    }
  }
  
  #Make sure that if the daily max is to be calculated, then the averaging time is less than or equal to 
  #24 hours, unless statistic is rolling mean.
  for (j in 1:length(pol$pollutant)){
    pol$index.avg.time.hours[j]=avg_time_hours(pol$index.avg.time[j])
    #Check the daily max setting isn't missing and is either yes or no
    if (is.na(pol$index.daily.max[j])) {
      stop(paste("index.daily.max setting for",toupper(pol$pollutant[j]),"is missing"))
      quit("no")      
    } else if (!(pol$index.daily.max[j] %in% c("yes","no"))) {
      stop(paste("index.daily.max setting for",toupper(pol$pollutant[j]),"must be yes or no"))
    }
    #If daily max output of non-rolling averages is selected, then check modelled averaging time isn't greater than 24 hours
    if (pol$index.daily.max[j] == "yes" & pol$index.statistic[j] != "rolling mean") {
      ok.to.continue = TRUE
      if (is.na(pol$index.avg.time.hours[j])) {
        #This will be the case if the averaging time is in months - not less than 24 hours by definition
        ok.to.continue = FALSE
      } else if (pol$index.avg.time.hours[j] > 24) {
        ok.to.continue = FALSE
      }
      if (!ok.to.continue) {
        stop(paste("Averaging time for",toupper(pol$pollutant[j]),"must be equal to or less than a day for the daily max of a non-rolling statistic to be calculated"))
        quit("no")
      }
    }
  }

  return(pol)
}

#Obtain the lower bounds of each index
#Inputs pollutants, outputs lower bounds of each index
find.index.lower.bounds = function(pol) {
  
  #Find out the number of forecast index bands
  #We already know that there's an 'i1' column, so look for 'i2', 'i3' etc
  j=1
  while (length(which(names(pol)==paste("i",j,sep="")))!=0) {
    j=j+1
  }
  number.index.bands = j-1
  index.headers = paste("i",seq(from=1, to=number.index.bands, by=1),sep="")
  lower.bounds.of.index = pol[,which(names(pol)==index.headers[1])]
  if (number.index.bands > 1) {
    for (j in 2:number.index.bands) {
      start = (j-1)*length(pol$pollutant)+1
      finish = j*length(pol$pollutant)
      lower.bounds.of.index[start:finish] = pol[,which(names(pol)==index.headers[j])]
    }
  }
  lower.bounds.of.index = matrix(lower.bounds.of.index,nrow=length(pol$pollutant),ncol=number.index.bands)
  #Check that index lower bounds are given in ascending order
  for (j in 1:length(pol$pollutant)) {
    if (any(lower.bounds.of.index[j,]!=sort(lower.bounds.of.index[j,]))) {
      stop(paste("Index lower bounds are not in ascending order for pollutant",toupper(pol$pollutant[j])))
      quit("no")
    }
  }
  return(lower.bounds.of.index)
}

#Calculate indices of data
#Inputs all data, pollutants, number of pollutants and lower bounds and outputs data with indices applied
find.data.for.indexing = function(data,pol,j,lower.bounds.of.index,threshold) {
  
  data = subset(data, select=c("date","site","obs","mod")) 
  
  #GET OBS AND MOD TO SAME AVERAGING TIME AND STATISTIC AS INDEX
  if (pol$index.avg.time[j]==pol$avg.time[j] & pol$index.statistic[j]==pol$statistic[j]) {
    data.averaged = data
  } else {
    #Can assume obs and mod are both hourly data - anything else would have caused an error during earlier checks
    if (pol$index.statistic[j]=="rolling mean") {
      data.averaged.obs = rollingMean(data,pollutant="obs",width=pol$index.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right", type="site")
      data.averaged.obs = subset(data.averaged.obs,select=c("date","site","rolling"))
      names(data.averaged.obs)[names(data.averaged.obs)=="rolling"]="obs"
      data.averaged.mod = rollingMean(data,pollutant="mod",width=pol$index.avg.time.hours[j],new.name="rolling",data.thresh=threshold,align="right", type="site")
      data.averaged.mod = subset(data.averaged.mod,select=c("date","site","rolling"))
      names(data.averaged.mod)[names(data.averaged.mod)=="rolling"]="mod"
      data.averaged = merge(data.averaged.obs,data.averaged.mod)
    } else {
      data.averaged = timeAverage(data,avg.time=pol$index.avg.time[j],statistic=pol$index.statistic[j],data.thresh=threshold,type="site",interval=pol$avg.time[j], progress=FALSE)
    }
  }
  #Convert to daily maximum if required
  if (pol$index.daily.max[j]=="yes" & !(pol$index.avg.time[j] %in% c("day","1 day"))) {
    data.for.indexing = timeAverage(data.averaged,avg.time="day",statistic="max",data.thresh=threshold,type="site",interval=pol$index.avg.time[j], progress=FALSE)
  } else {
    data.for.indexing = data.averaged
  }
  #Convert concentrations to index units
  names(data.for.indexing)[which(names(data.for.indexing)=="obs")]="obs.output.units"
  names(data.for.indexing)[which(names(data.for.indexing)=="mod")]="mod.output.units"
  data.for.indexing$obs.index.units = convert.units(data.for.indexing$obs.output.units,pol$output.units[j],pol$index.units[j],pol$conv.ugm3.ppb[j])
  data.for.indexing$mod.index.units = convert.units(data.for.indexing$mod.output.units,pol$output.units[j],pol$index.units[j],pol$conv.ugm3.ppb[j])
  
  #APPLY INDEX BANDS TO CONCENTRATIONS
  data.for.indexing$index.obs = findInterval(data.for.indexing$obs.index.units,lower.bounds.of.index[j,])
  data.for.indexing$index.mod = findInterval(data.for.indexing$mod.index.units,lower.bounds.of.index[j,])
  
  #return data with site names as character strings for future manipulation 
  data.for.indexing$site = as.character(data.for.indexing$site)
  
  return(data.for.indexing)
}
