#' @title Calculate NowCast values
#'
#' @description
#' This function calculates the NowCast version of the incoming data.
#'
#' @param x Vctor of ordered PM2.5 data in units of ug/m3.
#' @param version Character identity specifying the type of nowcast algorithm to
#' be used. For details see \code{AirMonitor::monitor_nowcast()}.
#' @param includeShortTerm Logical specifying whether to calculate preliminary
#' NowCast values starting with the 2nd hour
#'
#' @return Vector of unitless NowCast values of the same length as \code{x}.
#'
#' @export
.nowcast <- function(
x,
version = "pm",
includeShortTerm = FALSE
) {
# Set parameters based on version
if (version == 'pm') {
numHrs <- 12
weightFactorMin <- 0.5
digits <- 1
} else if (version == 'pmAsian') {
numHrs <- 3
weightFactorMin <- 0.1
digits <- 1
} else if (version == 'ozone') {
numHrs <- 8
weightFactorMin <- NA # negative values adjusted up to 0 in .weightFactor()
digits <- 3 # NOTE: digits=3 assumes Ozone values given in ppm; update to 0 if values given in ppb
} else if (version == 'identity') {
NULL
}
if ( includeShortTerm ) {
firstHr <- 2
} else {
firstHr <- numHrs
}
# Start at the end end of the data (most recent hour) and work backwards
# The oldest hour for which we can calculate nowcast is numHrs, unless
# includeShortTerm=TRUE in which case we can go back to the 2nd hour.
for ( i in length(x):firstHr ) {
# Apply nowcast algorithm to numHrs data points in order with more recent first
concByHour <- x[i:max(1, i - numHrs + 1)]
if ( sum( is.na(concByHour[1:3]) ) >= 2 ) {
# If two or more of the most recent 3 hours are missing, no valid Nowcast will be reported
x[i] <- NA
} else if ( is.na(concByHour[1]) ) {
# If the current hour is missing, no valid Nowcast will be reported
# NOTE: This conflicts with the algorithm as described here:
# NOTE: https://www3.epa.gov/airnow/ani/pm25_aqi_reporting_nowcast_overview.pdf
# NOTE:
# NOTE: But experience shows that NowCast replacements for missing values for missing
# NOTE: PM2.5 values are very problematic.
# NOTE:
# NOTE: The Wikipedia page: https://en.wikipedia.org/wiki/NowCast_(air_quality_index)
# NOTE: has the following statement without citation:
# NOTE: "Because the most recent hours of data are weighted so heavily in the NowCast when
# NOTE: PM levels are changing, EPA does not report the NowCast when data is missing for c1 or c2."
# NOTE:
# NOTE: We take a compromise approach and only invalidate NowCast when data is missing for c1.
x[i] <- NA
} else {
# Calculate the weight factor according to the type of air quality data
weightFactor <- .weightFactor(concByHour, weightFactorMin)
# NOTE: We need to create vectors so that we can sum at the end with na.rm=TRUE
weightedConcs <- rep(as.numeric(NA),numHrs)
weightFactors <- rep(as.numeric(NA),numHrs)
# Loop over hours to get individual elements
for (j in 1:numHrs) {
if ( !is.na( concByHour[j]) ) {
weightedConcs[j] <- concByHour[j] * weightFactor^(j-1)
weightFactors[j] <- weightFactor^(j-1)
}
}
x[i] <- sum(weightedConcs, na.rm = TRUE) / sum(weightFactors, na.rm = TRUE)
}
}
# Set missing values when there are not enough preceding hours
x[1:(firstHr-1)] <- NA
return(x)
}
# Calculate the weight factor ('w' in the nowcast formula)
# concByHour: vector of hourly concentration values
# weightFactorMin (optional): wight factor minimum
# Assumes concByHour has at least one valid value to calculate min & max. In fact, .nowcast won't even call
# this function if more than one of the three most recent hours is invalid.
.weightFactor <- function(concByHour, weightFactorMin) {
min <- min(concByHour, na.rm = TRUE)
max <- max(concByHour, na.rm = TRUE)
# Calculate weight factor
# NOTE: https://forum.airnowtech.org/t/the-nowcast-for-ozone-and-pm/172 says that there is "no minimum
# NOTE: weight factor" for ozone; however, we limit the value to zero since otherwise it would be possible
# NOTE: to get negative weights, even as large as -Inf (i.e. if min<0 & max=0).
# NOTE: Otherwise, we don't worry about negatives, per the following:
# NOTE: https://forum.airnowtech.org/t/how-does-airnow-handle-negative-hourly-concentrations/143
weightFactor <- 1 - (max - min)/max
weightFactor <- min(weightFactor, 1, na.rm = TRUE)
weightFactor <- max(weightFactor, weightFactorMin, 0, na.rm = TRUE)
return(weightFactor)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.