diff --git a/NEWS.md b/NEWS.md index 3fff1068..cc7f7399 100644 --- a/NEWS.md +++ b/NEWS.md @@ -22,6 +22,8 @@ from the US EPA, AirNow, AIRSIS, WRCC and others. * added language specific versions of `AQI`: `AQI_en` and `AQI_es` * added `monitor_isMonitor()` to valide the structure of a *ws_monitor* object * updated example code in all functions + * turning off grid lines with `~Lwd=0` in `monitorPlot_timeseries()` now works on Windows + * `monitorPlot_dailyBarplot()` accepts `tlim` argument of class `POSIXct` ### PWFSLSmoke 1.0.18 diff --git a/R/addShadedNight.R b/R/addShadedNight.R index 4e59208b..548dfe54 100644 --- a/R/addShadedNight.R +++ b/R/addShadedNight.R @@ -22,7 +22,7 @@ addShadedNight <- function(timeInfo, col=adjustcolor('black',0.1)) { if ( localTime[1] < sunrise[1] ) { rect(par('usr')[1], ybottom=par('usr')[3], xright=sunrise[1], ytop=par('usr')[4], - col=col, lwd=0) + col=col, border=NA) } # Complete nights @@ -30,14 +30,14 @@ addShadedNight <- function(timeInfo, col=adjustcolor('black',0.1)) { for (i in seq(length(sunset)-1)) { rect(xleft=sunset[i], ybottom=par('usr')[3], xright=sunrise[i+1], ytop=par('usr')[4], - col=col, lwd=0) + col=col, border=NA) } } # Last sunset to right edge rect(xleft=sunset[length(sunset)], ybottom=par('usr')[3], xright=par('usr')[2], ytop=par('usr')[4], - col=col, lwd=0) + col=col, border=NA) } diff --git a/R/monitorPlot_dailyBarplot.R b/R/monitorPlot_dailyBarplot.R index 85b6c1dd..d74813d2 100644 --- a/R/monitorPlot_dailyBarplot.R +++ b/R/monitorPlot_dailyBarplot.R @@ -5,7 +5,7 @@ #' @param ws_monitor \emph{ws_monitor} object #' @param monitorID monitor ID for a specific monitor in \code{ws_monitor} (optional #' if \code{ws_monitor} only has one monitor) -#' @param tlim optional vector with start and end times (integer or character representing YYYYMMDD[HH]) +#' @param tlim optional vector with start and end times (integer or character representing YYYYMMDD[HH] or \code{POSIXct}) #' @param minHours minimum number of valid data hours required to calculate each daily average #' @param gridPos position of grid lines either 'over', 'under' ('' for no grid lines) #' @param gridCol color of grid lines (see graphical parameter 'col') @@ -28,7 +28,7 @@ #' \dontrun{ #' N_M <- monitor_subset(Northwest_Megafires, tlim=c(20150715,20150930)) #' main <- "Daily Average PM2.5 for Omak, WA" -#' monitorPlot_dailyBarplot(N_M, monitorID="530470013", main=main, +#' monitorPlot_dailyBarplot(N_M, monitorID="530470013_01", main=main, #' labels_x_nudge=1) #' addAQILegend(fill=rev(AQI$colors), pch=NULL) #' } @@ -62,16 +62,21 @@ monitorPlot_dailyBarplot <- function(ws_monitor, } } - # TODO: tlim in monitorPlot_dailyBarplot should accept POSIXct # When tlim is specified in whole days we should add hours to get the requsted full days if ( !is.null(tlim) ) { - tlimStrings <- as.character(tlim) + if ( 'POSIXct' %in% class(tlim) ) { + tlimStrings <- strftime(tlim, "%Y%m%d", tz='UTC') + } else { + tlimStrings <- as.character(tlim) + } if ( stringr::str_length(tlimStrings)[1] == 8 ) { - tlim[1] <- paste0(tlim[1],'00') + tlimStrings[1] <- paste0(tlimStrings[1],'00') } if ( stringr::str_length(tlimStrings)[2] == 8 ) { - tlim[2] <- paste0(tlim[2],'23') + tlimStrings[2] <- paste0(tlimStrings[2],'23') } + # Recreate tlim + tlim <- tlimStrings } # Subset to a single monitor diff --git a/R/monitorPlot_timeOfDaySpaghetti.R b/R/monitorPlot_timeOfDaySpaghetti.R index ddbb832b..4295a3b4 100644 --- a/R/monitorPlot_timeOfDaySpaghetti.R +++ b/R/monitorPlot_timeOfDaySpaghetti.R @@ -149,17 +149,17 @@ monitorPlot_timeOfDaySpaghetti <- function(ws_monitor, # Left edge to sunrise rect(par('usr')[1], ybottom=par('usr')[3], xright=sunrise, ytop=par('usr')[4], - col=col_shadedNight, lwd=0) + col=col_shadedNight, border=NA) # Sunset to right edge rect(xleft=sunset, ybottom=par('usr')[3], xright=par('usr')[2], ytop=par('usr')[4], - col=col_shadedNight, lwd=0) + col=col_shadedNight, border=NA) } # AQI Lines - if ( aqiLines ) { + if ( aqiLines && lwd_aqi > 0 ) { abline(h=AQI$breaks_24[2:6], col=col_aqi, lwd=lwd_aqi) } diff --git a/R/monitorPlot_timeseries.R b/R/monitorPlot_timeseries.R index 2b82f095..43916e56 100644 --- a/R/monitorPlot_timeseries.R +++ b/R/monitorPlot_timeseries.R @@ -191,11 +191,16 @@ monitorPlot_timeseries <- function(ws_monitor, # Add vertical lines to denote days and/or hour breaks hour_indices <- which(as.numeric(strftime(times,format="%H",tz=timezone)) %% hourInterval == 0) day_indices <- which(as.numeric(strftime(times,format="%H",tz=timezone)) %% 24 == 0) - abline(v=times[hour_indices], lwd=hourLwd) # at beginning of hour - abline(v=times[day_indices], lwd=dayLwd) # at beginning of day + # NOTE: Windows doesn't support lwd=0 + if ( hourLwd > 0 ) { + abline(v=times[hour_indices], lwd=hourLwd) # at beginning of hour + } + if ( dayLwd > 0 ) { + abline(v=times[day_indices], lwd=dayLwd) # at beginning of day + } # Add horizontal grid lines (before points if grid=='under') - if ( gridPos == 'under' ) { + if ( gridPos == 'under' && gridLwd > 0 ) { abline(h=axTicks(2)[-1], col=gridCol, lwd=gridLwd, lty=gridLty) } @@ -274,7 +279,7 @@ monitorPlot_timeseries <- function(ws_monitor, } # Add horizontal grid lines (on top of points if grid=='over') - if ( gridPos == 'over' ) { + if ( gridPos == 'over' && gridLwd > 0 ) { abline(h=axTicks(2)[-1], col=gridCol, lwd=gridLwd, lty=gridLty) } diff --git a/R/rawPlot_timeOfDaySpaghetti.R b/R/rawPlot_timeOfDaySpaghetti.R index e0a5ed22..93f0806a 100644 --- a/R/rawPlot_timeOfDaySpaghetti.R +++ b/R/rawPlot_timeOfDaySpaghetti.R @@ -134,12 +134,12 @@ rawPlot_timeOfDaySpaghetti <- function(df, # Left edge to sunrise rect(par('usr')[1], ybottom=par('usr')[3], xright=sunrise, ytop=par('usr')[4], - col=col_shadedNight, lwd=0) + col=col_shadedNight, border=NA) # Sunset to right edge rect(xleft=sunset, ybottom=par('usr')[3], xright=par('usr')[2], ytop=par('usr')[4], - col=col_shadedNight, lwd=0) + col=col_shadedNight, border=NA) } diff --git a/R/rawPlot_timeseries.R b/R/rawPlot_timeseries.R index 65ec1dc5..ea410fe4 100644 --- a/R/rawPlot_timeseries.R +++ b/R/rawPlot_timeseries.R @@ -226,11 +226,16 @@ rawPlot_timeseries <- function(df, # Add vertical lines to denote days and/or hour breaks hour_indices <- which(as.numeric(strftime(times,format="%H",tz=timezone)) %% hourInterval == 0) day_indices <- which(as.numeric(strftime(times,format="%H",tz=timezone)) %% 24 == 0) - abline(v=times[hour_indices], lwd=hourLwd) # at beginning of hour - abline(v=times[day_indices], lwd=dayLwd) # at beginning of day + # NOTE: Windows doesn't support lwd=0 + if ( hourLwd > 0 ) { + abline(v=times[hour_indices], lwd=hourLwd) # at beginning of hour + } + if ( dayLwd > 0 ) { + abline(v=times[day_indices], lwd=dayLwd) # at beginning of day + } # Add horizontal grid lines (before points if grid=='under') - if ( gridPos == 'under' ) { + if ( gridPos == 'under' && gridLwd > 0 ) { abline(h=axTicks(2)[-1], col=gridCol, lwd=gridLwd, lty=gridLty) } @@ -252,11 +257,9 @@ rawPlot_timeseries <- function(df, # Add lines do.call(lines,argsList) - if ( gridPos == 'over' ) { - + if ( gridPos == 'over' && gridLwd > 0 ) { # Horizontal lines abline(h=axTicks(2)[-1], col=gridCol, lwd=gridLwd, lty=gridLty) - } } diff --git a/localExamples/Boise_SmokeMtg_2018.R b/localExamples/Boise_SmokeMtg_2018.R index 2674fe3c..cd9bd6c3 100644 --- a/localExamples/Boise_SmokeMtg_2018.R +++ b/localExamples/Boise_SmokeMtg_2018.R @@ -28,7 +28,7 @@ png('oridwa_hourly.png', width=800, height=600) par(mar=c(5.1,4.1,9.1,8.1)) monitorPlot_timeseries(jas, style='gnats', xlab='2017', ylim=c(0,400), xpd=NA) -addAQIStackedBar(pos='left', width=0.01, labels=FALSE, title=FALSE) +addAQIStackedBar(pos='left', width=0.01) addAQILines(lwd=2) text(par('usr')[2], AQI$breaks_24[2:6], AQI$names[2:6], pos=4, xpd=NA) mtext('Hourly PM2.5', side=3, line=6, adj=0, cex=2.0, font=2) @@ -51,7 +51,7 @@ png('oridwa_daily.png', width=800, height=600) par(mar=c(5.1,4.1,9.1,8.1)) monitorPlot_timeseries(dailyMean, style='gnats', xlab='2017', ylim=c(0,200), xpd=NA) -addAQIStackedBar(pos='left', width=0.01, labels=FALSE, title=FALSE) +addAQIStackedBar(pos='left', width=0.01) addAQILines(lwd=2) text(par('usr')[2], AQI$breaks_24[2:5], AQI$names[2:5], pos=4, xpd=NA) mtext('Daily Mean PM2.5', side=3, line=6, adj=0, cex=2.0, font=2) @@ -80,7 +80,7 @@ par(oldPar) # ----------------------------------------------------------------------------- # Nez Perce monitors -as <- monitor_subset(jas, tlim=c(20170801,20170915), tz="America/Los_Angeles") +as <- monitor_subset(jas, tlim=c(20170801,20170915), timezone="America/Los_Angeles") Lewiston <- monitor_subset(as, monitorIDs='160690012_01') Clarkston <- monitor_subset(as, monitorIDs='530030004_01') diff --git a/localExamples/esriLocationMap.R b/localExamples/esriLocationMap.R index f2c3f88b..6a322f7a 100644 --- a/localExamples/esriLocationMap.R +++ b/localExamples/esriLocationMap.R @@ -46,7 +46,7 @@ downloadEsriLocationMap <- function(latitude, longitude, maptype = "World_Street # Map for Wenatchee airnow <- airnow_loadLatest() -wenatchee <- monitor_subset(airnow, monitorIDs = "530070011") +wenatchee <- monitor_subset(airnow, monitorIDs = "530070011_01") downloadEsriLocationMap(wenatchee$meta$latitude, wenatchee$meta$longitude, zoom = 9, size = 250, maptype = "World_Street_Map", destination = "~/Desktop/esrimaps/worldstreet.png") diff --git a/man/monitorDygraph.Rd b/man/monitorDygraph.Rd index b46f052e..e6fc76fe 100644 --- a/man/monitorDygraph.Rd +++ b/man/monitorDygraph.Rd @@ -28,6 +28,7 @@ This function creates interactive graphs that will be displayed in RStudio's 'Vi } \examples{ \dontrun{ +# Napa Fires -- October, 2017 ca <- airnow_load(2017) \%>\% monitor_subset(tlim=c(20171001,20171101), stateCodes='CA') Vallejo <- monitor_subset(ca, monitorIDs='060950004_01') diff --git a/man/monitorLeaflet.Rd b/man/monitorLeaflet.Rd index bff8dd81..4f910937 100644 --- a/man/monitorLeaflet.Rd +++ b/man/monitorLeaflet.Rd @@ -58,6 +58,7 @@ to use as the background map. } \examples{ \dontrun{ +# Napa Fires -- October, 2017 ca <- airnow_load(2017) \%>\% monitor_subset(tlim=c(20171001,20171101), stateCodes='CA') v_low <- AQI$breaks_24[5] diff --git a/man/monitorMap_performance.Rd b/man/monitorMap_performance.Rd index f938b34f..565ee52f 100644 --- a/man/monitorMap_performance.Rd +++ b/man/monitorMap_performance.Rd @@ -60,11 +60,15 @@ the size/colors to remain constant. } \examples{ \dontrun{ -# Spokane summer of 2015 -wa <- airnow_load(20150701, 20150930, stateCodes='WA') -wa <- monitor_rollingMean(wa, width=3) -MonroeSt <- monitor_subset(wa, monitorIDs="530630047_01") -monitorMap_performance(wa, MonroeSt, cex=2) +# Napa Fires -- October, 2017 +ca <- airnow_load(2017) \%>\% + monitor_subset(tlim=c(20171001,20171101), stateCodes='CA') +Vallejo <- monitor_subset(ca, monitorIDs='060950004_01') +Napa_Fires <- monitor_subsetByDistance(ca, + longitude = Vallejo$meta$longitude, + latitude = Vallejo$meta$latitude, + radius = 50) +monitorMap_performance(ca, Vallejo, cex=2) title('Heidike Skill of monitors predicting another monitor.') } } diff --git a/man/monitorPlot_dailyBarplot.Rd b/man/monitorPlot_dailyBarplot.Rd index d49d7ee7..7663a7d9 100644 --- a/man/monitorPlot_dailyBarplot.Rd +++ b/man/monitorPlot_dailyBarplot.Rd @@ -14,7 +14,7 @@ monitorPlot_dailyBarplot(ws_monitor, monitorID = NULL, tlim = NULL, \item{monitorID}{monitor ID for a specific monitor in \code{ws_monitor} (optional if \code{ws_monitor} only has one monitor)} -\item{tlim}{optional vector with start and end times (integer or character representing YYYYMMDD[HH])} +\item{tlim}{optional vector with start and end times (integer or character representing YYYYMMDD[HH] or \code{POSIXct})} \item{minHours}{minimum number of valid data hours required to calculate each daily average} @@ -50,7 +50,7 @@ tweak the date labeling. Units used are the same as those in the plot. \dontrun{ N_M <- monitor_subset(Northwest_Megafires, tlim=c(20150715,20150930)) main <- "Daily Average PM2.5 for Omak, WA" -monitorPlot_dailyBarplot(N_M, monitorID="530470013", main=main, +monitorPlot_dailyBarplot(N_M, monitorID="530470013_01", main=main, labels_x_nudge=1) addAQILegend(fill=rev(AQI$colors), pch=NULL) } diff --git a/man/monitor_performance.Rd b/man/monitor_performance.Rd index 7db32ed8..87f9e7ac 100644 --- a/man/monitor_performance.Rd +++ b/man/monitor_performance.Rd @@ -36,16 +36,21 @@ all available metrics are returned. } \examples{ \dontrun{ -airnow <- airnow_load(20150101, 20151231) -airnow_dailyAvg <- monitor_dailyStatistic(airnow, mean) -airnow_dailyAvg <- monitor_subset(airnow_dailyAvg, stateCodes='WA') -airnow_m1 <- airnow_load(20141231, 20151230) -airnow_dailyAvg_m1 <- monitor_dailyStatistic(airnow_m1, mean) -airnow_dailyAvg_m1 <- monitor_subset(airnow_dailyAvg_m1, stateCodes='WA') -threshold <- AQI$breaks_24[3] -performanceMetrices <- monitor_performance(airnow_dailyAvg_m1, - airnow_dailyAve, - threshold, threshold) +# If daily avg data were the prediciton and Spokane were +# the observed, which WA State monitors had skill? +wa <- airnow_load(2017) \%>\% monitor_subset(stateCodes='WA') +wa_dailyAvg <- monitor_dailyStatistic(wa, mean) +Spokane_dailyAvg <- monitor_subset(wa_dailyAvg, monitorIDs='530630021_01') +threshold <- AQI$breaks_24[4] # Unhealthy +performanceMetrics <- monitor_performance(wa_dailyAvg, + Spokane_dailyAvg, + threshold, threshold) +monitorIDs <- rownames(performanceMetrics) +mask <- performanceMetrics$heidikeSkill & + !is.na(performanceMetrics$heidikeSkill) +skillfulIDs <- monitorIDs[mask] +skillful <- monitor_subset(wa_dailyAvg, monitorIDs=skillfulIDs) +monitorLeaflet(skillful) } } \seealso{ diff --git a/man/monitor_subsetByDistance.Rd b/man/monitor_subsetByDistance.Rd index 766957d7..e80d241c 100644 --- a/man/monitor_subsetByDistance.Rd +++ b/man/monitor_subsetByDistance.Rd @@ -34,9 +34,15 @@ of monitors (or grid cells) returned may be less than the specified \code{count} } \examples{ \dontrun{ -airnow <- airnow_load(20140913, 20141010) -KingFire <- monitor_subsetByDistance(airnow, longitude=-120.604, latitude=38.782, radius=50) -monitorLeaflet(KingFire) +# Napa Fires -- October, 2017 +ca <- airnow_load(2017) \%>\% + monitor_subset(tlim=c(20171001,20171101), stateCodes='CA') +Vallejo <- monitor_subset(ca, monitorIDs='060950004_01') +Napa_Fires <- monitor_subsetByDistance(ca, + longitude = Vallejo$meta$longitude, + latitude = Vallejo$meta$latitude, + radius = 50) +monitorLeaflet(Napa_Fires) } } \seealso{ diff --git a/man/skill_ROC.Rd b/man/skill_ROC.Rd index b9e4e295..521b8aa6 100644 --- a/man/skill_ROC.Rd +++ b/man/skill_ROC.Rd @@ -26,12 +26,13 @@ thresholds as well as the area under the ROC curve. } \examples{ \dontrun{ -# Spokane summer of 2015 -airnow <- airnow_load(20150701,20150930) -airnow <- monitor_rollingMean(airnow, width=3) -MonroeSt <- monitor_subset(airnow, monitorIDs="530630047_01") -EBroadway <- monitor_subset(airnow, monitorIDs="530639997_01") -rocList <- skill_ROC(EBroadway, MonroeSt, t1Range=c(0,100), t2=55) +# Napa Fires -- October, 2017 +ca <- airnow_load(2017) \%>\% + monitor_subset(tlim=c(20171001,20171101), stateCodes='CA') +Vallejo <- monitor_subset(ca, monitorIDs='060950004_01') +Napa <- monitor_subset(ca, monitorIDs='060550003_01') +t2 <- AQI$breaks_24[4] # 'Unhealthy' +rocList <- skill_ROC(Vallejo, Napa, t1Range=c(0,100), t2=t2) roc <- rocList$roc auc <- rocList$auc plot(roc$TPR ~ roc$FPR, type='S') diff --git a/man/skill_ROCPlot.Rd b/man/skill_ROCPlot.Rd index 95f5ee44..e4f7f6de 100644 --- a/man/skill_ROCPlot.Rd +++ b/man/skill_ROCPlot.Rd @@ -25,12 +25,12 @@ This function plots ROC curves for a variety of \code{observed} classification t } \examples{ \dontrun{ -# Spokane summer of 2015 -airnow <- airnow_load(20150701,20150930) -airnow <- monitor_rollingMean(airnow, width=3) -MonroeSt <- monitor_subset(airnow, monitorIDs="530630047_01") -EBroadway <- monitor_subset(airnow, monitorIDs="530639997_01") -skill_ROCPlot(EBroadway, MonroeSt) +# Napa Fires -- October, 2017 +ca <- airnow_load(2017) \%>\% + monitor_subset(tlim=c(20171001,20171101), stateCodes='CA') +Vallejo <- monitor_subset(ca, monitorIDs='060950004_01') +Napa <- monitor_subset(ca, monitorIDs='060550003_01') +skill_ROCPlot(Vallejo, Napa) } } \references{ diff --git a/vignettes/NowCast.Rmd b/vignettes/NowCast.Rmd index 8eefcb83..1d7d6a85 100644 --- a/vignettes/NowCast.Rmd +++ b/vignettes/NowCast.Rmd @@ -14,19 +14,19 @@ knitr::opts_chunk$set(fig.width=7, fig.height=5, comment=NA) options(width = 105) ``` -This vignette documents the basis and functionality of the `monitor_nowcast()` function, which converts a `ws_monitor` object's values to NowCast values. We provide details on the NowCast algorithm and our functional implementation thereof. We also provide examples to highlight specific attributes and potential points of confusion in the algorithm and/or our implementation. +This vignette documents the `monitor_nowcast()` function, which converts a `ws_monitor` object's values to NowCast values. We provide details on the NowCast algorithm and our implementation thereof. We also provide examples to highlight specific attributes and potential points of confusion in the algorithm. -This vignette also briefly covers the `monitor_aqi()` function, which uses the `monitor_nowcast()` function to convert raw PM~2.5~ data into NowCast-based AQI values. +This vignette also briefly covers the `monitor_aqi()` function, which uses the `monitor_nowcast()` function to convert raw PM2.5 data into NowCast-based AQI values. ## What is NowCast? NowCast is an air quality data smoothing algorithm that puts an emphasis on recent values when measurements are unstable, and approaches a long-term (*e.g.* 12-hour) average when measurements are stable. -The original algorithm, known as the Conroy method, was developed in 2003 to make real-time air quality measurements roughly comparable to established regulatory air quality health thresholds (*e.g.* 24‐hour PM~2.5~ standards). However, that method was shown to be slow to respond to rapidly changing air quality conditions, which, at best, reduced public confidence in disseminated AQI values, and at worst, had the potential to adversely affect the health of those in high impact areas (*e.g.* near wildfires). +The original algorithm, known as the Conroy method, was developed in 2003 to make real-time air quality measurements roughly comparable to established regulatory air quality health thresholds (*e.g.* 24‐hour PM2.5 standards). However, that method was shown to be slow to respond to rapidly changing air quality conditions, which, at best, reduced public confidence in disseminated AQI values, and at worst, had the potential to adversely affect the health of those in high impact areas (*e.g.* near wildfires). -In response, EPA developed a new method -- the Reff method -- in 2013 to be more responsive to rapidly changing air quality conditions. We find technical support for applying the new NowCast algorithm to hourly PM~2.5~, PM~10~, and O~3~ data, though theoretically it could be applied to regular-interval time series data of any type, including other criteria pollutants. +In response, EPA developed a new method -- the Reff method -- in 2013 to be more responsive to rapidly changing air quality conditions. We provide technical support for applying the new NowCast algorithm to hourly PM2.5, PM10, and O3 data, though theoretically it could be applied to regular-interval time series data of any type, including other criteria pollutants. -We provide algorithm specifics for PM~2.5~, PM~10~, and O~3~ below. +We provide algorithm specifics for PM2.5, PM10, and O3 below. Sources: https://www3.epa.gov/airnow/ani/pm25_aqi_reporting_nowcast_overview.pdf @@ -34,8 +34,6 @@ https://forum.airnowtech.org/t/the-nowcast-for-ozone-and-pm/172 ## The NowCast Algorithm (Reff Method) -The following sections outline the NowCast algorithm. - ### NowCast Equation The NowCast value for a given hour can be calculated as follows: @@ -54,9 +52,9 @@ https://www3.epa.gov/airnow/ani/pm25_aqi_reporting_nowcast_overview.pdf The NowCast algorithm uses hourly averages from the prior $N$ clock hours, where the value of $N$ depends on the pollutant being processed: - - **PM~2.5~**: $N = 12$ - - **PM~10~**: $N = 12$ - - **O~3~**: $N = 8$ + - **PM2.5**: $N = 12$ + - **PM10**: $N = 12$ + - **O3**: $N = 8$ The hourly averages are denoted below by $c_i$, where $i$ represents the number of hours before present. For example, $c_1$, $c_2$, $c_3$, $...$, $c_N$ represent the hourly averages for the most recent 1, 2, 3, $...$, $N$ hours. @@ -84,7 +82,7 @@ https://en.wikipedia.org/wiki/NowCast_(air_quality_index) ### Minimum Weight Factor -Before plugging into the NowCast equation, $w*$ is updated to $w$ for PM~2.5~ and PM~10~ as follows: +Before plugging into the NowCast equation, $w*$ is updated to $w$ for PM2.5 and PM10 as follows: $$ w = @@ -94,7 +92,7 @@ w^* & \text{if} & w^*>\frac{1}{2} \\ \end{cases} $$ -For O~3~, there is no minimum weight factor, so we define $w = w^*$ +For O3, there is no minimum weight factor, so we define $w = w^*$ Source: https://forum.airnowtech.org/t/the-nowcast-for-ozone-and-pm/172 @@ -103,16 +101,16 @@ https://forum.airnowtech.org/t/the-nowcast-for-ozone-and-pm/172 Final NowCast values are truncated based on the type of data being processed: -- **PM~2.5~**: 0.1 µg/m3 -- **PM~10~**: 1 µg/m3 -- **O~3~**: 1 ppb (0.001 ppm) +- **PM2.5**: 0.1 µg/m3 +- **PM10**: 1 µg/m3 +- **O3**: 1 ppb (0.001 ppm) Source: https://forum.airnowtech.org/t/the-nowcast-for-ozone-and-pm/172 ## Algorithm Expansion -For PM~2.5~, which uses 12 hours of data, the NowCast equation can be expanded as follows: +For PM2.5, which uses 12 hours of data, the NowCast equation can be expanded as follows: $$ NowCast = \frac @@ -133,15 +131,13 @@ $$ which is just a simple 12-hour arithmetic average. Incidentally, all 12 hourly averages and the 12-hour average itself would all be equivalent in this case. -In the case of highly variable PM~2.5~ data, $w$ would be set to the minimum value of $1/2$, and the most recent data would carry the majority of the weight in the equation above. +In the case of highly variable PM2.5 data, $w$ would be set to the minimum value of $1/2$, and the most recent data would carry the majority of the weight in the equation above. ## Implementation Details -The following sections expand on specifics of the NowCast algorithm and details regarding our functional implementation. - ### Missing Data -The NowCast algorithm ignores terms corresponding to hours for which a valid observation is not available. For example, suppose PM~2.5~ is invalid for all but the first three and last three hours of a 12-hour period. Then the PM~2.5~ NowCast equation takes the following form: +The NowCast algorithm ignores terms corresponding to hours for which a valid observation is not available. For example, suppose PM2.5 is invalid for all but the first three and last three hours of a 12-hour period. Then the PM2.5 NowCast equation takes the following form: $$ NowCast = \frac @@ -173,7 +169,7 @@ As mentioned above, the NowCast algorithm only requires two valid hours to calcu We assert that it would be inappropriate to do so, _usually_. -In most cases a user will have created a ws_monitor object using one of the `*_createMonitorObject()` or `*_load()`, functions, which return data for a specific period of time. Data before this period is not necessarily invalid, it was simply not retrieved. But the function itself has no way of knowing whether such earlier data exists, so it has no choice but to consider earlier hours "invalid". This means that, if followed by-the-book, the NowCast algorithm could return different values for a given hour depending on whether or not the earlier data had been retrieved. This is not a desirable behavior, so by default the `monitor_nowcast()` returns invalid NowCast values until the $N$^th^ hour of data. +In most cases a user will have created a ws_monitor object using one of the `~_createMonitorObject()` or `~_load()`, functions, which return data for a specific period of time. Data before this period is not necessarily invalid, it was simply not retrieved. But the function itself has no way of knowing whether such earlier data exists, so it has no choice but to consider earlier hours "invalid". This means that, if followed by-the-book, the NowCast algorithm could return different values for a given hour depending on whether or not the earlier data had been retrieved. This is not a desirable behavior, so by default the `monitor_nowcast()` returns invalid NowCast values until the $N$^th^ hour of data. However, we provide a manual override in `includeShortTerm=TRUE` which causes the function to return valid values as per the bare-bones data availability requirements described above, treating the hours before the beginning of the data as invalid. Thus, it can return "valid" NowCast values as early as the second (valid) hour in the data. @@ -187,7 +183,7 @@ Negative values are handled prior to converting NowCast or other values to AQI, ### Negative Weight Factors -While PM~2.5~ and PM~10~ have minimum weight factors of $1/2$, there is no minimum weight factor for O~3~. Does this mean O~3~ could have negative weight factors? We find no mention of this in the NowCast literature; however, we assume $w_{min}$ cannot be negative since otherwise the calculation could become overwhelmed by outrageous weight factors, even as large as negative infinity (e.g. if $c_{min}<0$ and $c_{max}=0$). This assumption is built into our `monitor_nowcast()` function. +While PM2.5 and PM10 have minimum weight factors of $1/2$, there is no minimum weight factor for O3. Does this mean O3 could have negative weight factors? We find no mention of this in the NowCast literature; however, we assume $w_{min}$ cannot be negative since otherwise the calculation could become overwhelmed by outrageous weight factors, even as large as negative infinity (e.g. if $c_{min}<0$ and $c_{max}=0$). This assumption is built into our `monitor_nowcast()` function. ### Argument: `version` @@ -208,12 +204,12 @@ The `version` argument sets defaults for the number of hours in the lookback $N$ The default setting is `version='pm25'` since this is the parameter most commonly stored in `ws_monitor` objects. -`version='ozone'` supports the O~3~ NowCast as described above. +`version='ozone'` supports the O3 NowCast as described above. `version='pmAsian'` supports an alternative shorter-term NowCast as proposed here: https://aqicn.org/faq/2015-03-15/air-quality-nowcast-a-beginners-guide/ -Although the NowCast algorithm itself supports PM~10~, we do not currently provide functionality for this parameter in the `monitor_nowcast()` function. +Although the NowCast algorithm itself supports PM10, we do not currently provide functionality for this parameter in the `monitor_nowcast()` function. In the future we may allow manual override of the settings described above to allow for custom NowCast-type algorithms. @@ -233,7 +229,7 @@ The following examples demonstrate the functionality of `monitor_nowcast()` and ### Setup -For the following examples we will use the Northwest Megafires data from the PWFSLSmoke package. In particular, we will look at PM~2.5~ data from Omak, WA, which was heavily impacted by smoke from wildfires during the second half of August, 2015: +For the following examples we will use the Northwest Megafires data from the PWFSLSmoke package. In particular, we will look at PM2.5 data from Omak, WA, which was heavily impacted by smoke from wildfires during the second half of August, 2015: ```{r} suppressPackageStartupMessages(library(PWFSLSmoke)) @@ -255,9 +251,9 @@ title("Hourly and Nowcast PM2.5 Values\nOmak, Washington; August, 2015") ### Example 1: Basic Formula Verification -In the code above we used the `monitor_nowcast()` function to calculate PM~2.5~ NowCast values for the Omak `ws_monitor` object. Let's see if we can verify the function's output for a single hour. +In the code above we used the `monitor_nowcast()` function to calculate PM2.5 NowCast values for the Omak `ws_monitor` object. Let's see if we can verify the function's output for a single hour. -Below is the Omak PM~2.5~ data for the first 12 hours of 8/21/15. Let's see if we can verify the accuracy of the NowCast value for the last hour in this series, 8/21/15 Hour 11. +Below is the Omak PM2.5 data for the first 12 hours of 8/21/15. Let's see if we can verify the accuracy of the NowCast value for the last hour in this series, 8/21/15 Hour 11. ```{r} Omak_2015_08_21 <- monitor_subset(Omak, tlim=c(2015082100, 2015082111)) @@ -276,7 +272,7 @@ Per the NowCast algorithm, we define $w*$ as $\frac{c_{min}}{c_{max}}$: (w_star <- min(example1_values)/max(example1_values)) ``` -We now define $w$ based on $w^*$ and the minimum weight factor; recall that $w_{min}=\frac{1}{2}$ for PM~2.5~: +We now define $w$ based on $w^*$ and the minimum weight factor; recall that $w_{min}=\frac{1}{2}$ for PM2.5: ```{r} (w <- max(1/2, w_star))