Nothing
#' @export
#'
#' @title Leaflet interactive map of monitor locations
#'
#' @param monitor \emph{mts_monitor} object.
#' @param slice Either a formatted time string, a time index, the name of a
#' (potentially user defined) function used to collapse the time axis.
# @param breaks set of breaks used to assign colors
# @param colors a set of colors for different levels of air quality data
# determined by \code{breaks}
# @param labels a set of text labels, one for each color
# @param legendTitle legend title
#' @param radius radius of monitor circles
#' @param opacity opacity of monitor circles
#' @param maptype optional name of leaflet ProviderTiles to use, e.g. "terrain"
# @param popupInfo a vector of column names from monitor$meta to be shown in
# a popup window
#' @param extraVars Character vector of additional column names from \code{monitor$meta}
#' to be shown in leaflet popups.
#' @param jitter Amount to use to slightly adjust locations so that multiple
#' monitors at the same location can be seen. Use zero or \code{NA} to see
#' precise locations.
#' @param NAAQS Version of NAAQS levels to use. See Note.
#' @param ... Additional arguments passed to \code{leaflet::addCircleMarker()}.
#'
#' @description
#' This function creates interactive maps that will be displayed in RStudio's
#' 'Viewer' tab. The \code{slice} argument is used to collapse a
#' \emph{mts_monitor} timeseries into a single value. If \code{slice} is an
#' integer, that row index will be selected from the \code{monitor$data}
#' dataframe. If \code{slice} is a function (unquoted), that function will be
#' applied to the timeseries with the argument \code{na.rm=TRUE} (e.g.
#' \code{max(..., na.rm=TRUE)}).
#'
#' If \code{slice} is a user defined function, it will be used with argument
#' \code{na.rm=TRUE} to collapse the time dimension. Thus, user defined
#' functions must accept \code{na.rm} as an argument.
#'
#' @details
#' The \code{maptype} argument is mapped onto leaflet "ProviderTile" names.
#' Current map types include:
#'
#' \describe{
#' \item{"roadmap"}{ -- "OpenStreetMap"}
#' \item{"satellite"}{ -- "Esri.WorldImagery"}
#' \item{"terrain"}{ -- "Esri.WorldTopoMap"}
#' \item{"toner"}{ -- "Stamen.Toner"}
#' }
#'
#' If a character string not listed above is provided, it will be used as the
#' underlying map tile if available. See
#' \url{https://leaflet-extras.github.io/leaflet-providers/} for a list of
#' "provider tiles" to use as the background map.
#'
#' @note
#' On February 7, 2024, EPA strengthened the National Ambient Air Quality
#' Standards for Particulate Matter (PM NAAQS) to protect millions of Americans
#' from harmful and costly health impacts, such as heart attacks and premature
#' death. Particle or soot pollution is one of the most dangerous forms of air
#' pollution, and an extensive body of science links it to a range of serious
#' and sometimes deadly illnesses. EPA is setting the level of the primary
#' (health-based) annual PM2.5 standard at 9.0 micrograms per cubic meter to
#' provide increased public health protection, consistent with the available
#' health science.
#' See \href{https://www.epa.gov/pm-pollution/final-reconsideration-national-ambient-air-quality-standards-particulate-matter-pm}{PM NAAQS update}.
#'
#' @return Invisibly returns a leaflet map of class "leaflet".
#'
#' @examples
#' \dontrun{
#' library(AirMonitor)
#' # Fail gracefully if any resources are not available
#' try({
#'
#' # Maximum AQI category at each site
#' monitor_loadLatest() %>%
#' monitor_filter(stateCode %in% CONUS) %>%
#' monitor_leaflet()
#'
#' # Mean AQI category at each site
#' monitor_loadLatest() %>%
#' monitor_filter(stateCode %in% CONUS) %>%
#' monitor_leaflet(
#' slice = "mean"
#' )
#'
#' # Mean AQI category at each site using the updated NAAQS
#' monitor_loadLatest() %>%
#' monitor_filter(stateCode %in% CONUS) %>%
#' monitor_leaflet(
#' slice = "mean",
#' NAAQS = "PM2.5_2024"
#' )
#'
#' }, silent = FALSE)
#' }
monitor_leaflet <- function(
monitor,
slice = "max",
radius = 10,
opacity = 0.7,
maptype = "terrain",
extraVars = NULL,
jitter = 5e-4,
NAAQS = c("PM2.5", "PM2.5_2024"),
...
) {
# Configurable values
na.color = "#cccccc"
# ----- Validate parameters --------------------------------------------------
monitor_check(monitor)
if ( monitor_isEmpty(monitor) )
stop("monitor object has no data")
slice <- MazamaCoreUtils::setIfNull(slice, "max")
maptype <- MazamaCoreUtils::setIfNull(maptype, "terrain")
if ( !is.null(extraVars) ) {
unrecognizedVars <- setdiff(extraVars, names(monitor$meta))
if ( length(unrecognizedVars) > 0 ) {
stop("variables in 'extraVars' not found in 'monitor$meta'")
}
}
if ( is.null(jitter) || is.na(jitter) ) {
jitter <- 0
}
NAAQS = match.arg(NAAQS)
# ----- Pollutant dependent AQI ----------------------------------------------
# See: https://www.airnow.gov/sites/default/files/2020-05/aqi-technical-assistance-document-sept2018.pdf
pollutant <- toupper(unique(monitor$meta$pollutant))
if ( length(pollutant) > 1 ) {
pollutantString <- paste0(pollutant, collapse = ", ")
stop(sprintf("multiple pollutants found: %s", pollutantString))
}
if ( pollutant == "CO" ) {
legendTitle <- "CO AQI Level"
units <- "ppm"
digits <- 1
} else if ( pollutant == "OZONE" ) {
legendTitle <- "Ozone AQI Level"
units <- "ppm"
digits <- 2
} else if ( pollutant == "PM2.5" ) {
legendTitle <- "PM2.5 AQI Level"
units <- "\U00B5g/m3"
digits <- 0
} else if ( pollutant == "PM10" ) {
legendTitle <- "PM10 AQI Level"
units <- "\U00B5g/m3"
digits <- 0
}
AQI_breaks <- US_AQI[[paste0("breaks_", pollutant)]]
AQI_colors <- US_AQI[[paste0("colors_", "EPA")]]
AQI_names <- US_AQI$names_eng
# Handle the added NAAQS argument
if ( pollutant == "PM2.5" && NAAQS == "PM2.5_2024" ) {
AQI_breaks <- US_AQI$breaks_PM2.5_2024
}
breaks <- AQI_breaks
colors <- AQI_colors
labels <- AQI_names
# ----- Create the 'slice' values --------------------------------------------
if ( is.numeric(slice) ) {
# * slice = number -----
if ( slice < 1 || slice > nrow(monitor$data) )
stop(sprintf("slice = %d is outside the range 1:%d", slice, nrow(monitor$data)))
slice <- round(slice)
popupValue <- as.numeric(dplyr::slice(monitor$data, slice)[-1])
popupWhen <-
strftime(monitor$data$datetime[slice], "on %B %d, %Y at %H:00", tz = "UTC", usetz = TRUE)
} else if ( is.character(slice) ) {
# * slice = character -----
if ( exists(slice) && is.function(get(slice)) ) {
# ** slice = function -----
FUN <- get(slice)
# NOTE: min/max will warn and return Inf/-Inf when all data are missing
# NOTE: while mean returns NaN so we need to suppress warnings and replace
# NOTE: all those non-finite values with NA.
suppressWarnings({
popupValue <- base::apply(
dplyr::select(monitor$data, -1),
2,
FUN,
na.rm = TRUE
)
})
popupValue[!is.finite(popupValue)] <- NA
popupWhen <- ""
legendTitle <- sprintf("%s %s", stringr::str_to_title(slice), legendTitle)
if ( slice == "max" ) {
# Can't find a good dplyr way to get the first occurance of each value so we roll our own
suppressWarnings({
dataBrick <- dplyr::select(monitor$data, -1)
sliceValueBrick <- matrix(rep(popupValue, nrow(dataBrick)), nrow = nrow(dataBrick), byrow = TRUE)
logicalBrick <- dataBrick == sliceValueBrick
logicalBrick[is.na(logicalBrick)] <- FALSE
firstRowAtMax <- base::apply(logicalBrick, 2, function(x) { min(which(x), na.rm = TRUE) })
firstRowAtMax[!is.finite(firstRowAtMax)] <- NA
firstTimeAtMax <- monitor$data$datetime[firstRowAtMax]
})
popupWhen <-
strftime(firstTimeAtMax, "on %B %d, %Y at %H:00", tz = "UTC", usetz = TRUE)
# TODO: Idea for local time not working yet.
# popupWhen <- vector("character", length(firstTimeAtMax))
# for ( i in seq_len(firstTimeAtMax) ) {
# popupWhen[i] <-
# strftime(firstTimeAtMax[i], "on %B %d, %Y at %H:00", tz = monitor$meta$timezone[i], usetz = TRUE)
# }
popupWhen[is.na(popupWhen)] <- ""
} # END of slice == "max"
} else {
# ** slice = datetime -----
result <- try({
sliceTime <- MazamaCoreUtils::parseDatetime(slice, timezone = "UTC")
}, silent = TRUE)
if ( "try-error" %in% class(result) ) {
stop("improper use of slice parameter")
} else {
# Now proceed as if slice were an integer
slice <- which(monitor$data$datetime == sliceTime)
popupValue <- as.numeric(dplyr::slice(monitor$data, slice)[-1])
popupWhen <-
strftime(monitor$data$datetime[slice], "on %B %d, %Y at %H:00", tz = "UTC", usetz = TRUE)
}
} # END slice = character
} else {
# * slice = neither -----
stop("improper use of slice parameter")
}
# ----- Order by popupValue --------------------------------------------------
# NOTE: This step is required if you want to have higher valued locations
# NOTE: plotted on top so that they aren't hidden by lower valued locations.
orderedIndices <- order(popupValue)
popupValue <- popupValue[orderedIndices]
popupWhen <- popupWhen[orderedIndices]
monitor$meta <- monitor$meta[orderedIndices,]
# monitor$data is no longer needed
# ----- Create colors and legend labels --------------------------------------
# If the user didn't use custom breaks then use AQI names and colors
if ( all.equal(breaks, AQI_breaks) && all.equal(colors, AQI_colors) ) {
# Ignore warnings from RColorBrewer as leaflet::colorBin does the right thing
suppressWarnings({
colorFunc <- leaflet::colorBin(
AQI_colors,
bins = AQI_breaks,
na.color = na.color
)
cols <- colorFunc(popupValue)
colors <- AQI_colors
labels <- AQI_names
})
} else {
if ( length(breaks) <= 2) {
stop("Please specify the correct vector of breaks")
}
if (! (length(breaks) - 1 == length(colors)) ) {
stop("The number of colors provided should be one less than the number of breaks")
}
if ( missing(labels) ){
labels <- paste(sprintf("%.1f", breaks[-length(breaks)]),
"--",
sprintf("%.1f", breaks[-1]))
} else if ( length(labels) != length(colors) ) {
stop("The number of labels should be equal to the number of colors")
}
# Create levels and use them to create a color mask
levels <- .bincode(popupValue, breaks, include.lowest = TRUE)
if ( !all(!is.na(levels)) ) {
print("NOTE that there are data points outside of your specified breaks, non-requested color(s) might be displayed on your map.")
}
cols <- colors[levels]
}
# ----- Create popup text ----------------------------------------------------
popupText <- paste0(
"<b>", monitor$meta$locationName, "</b><br>",
"<b>", round(popupValue, digits), " ", units, "</b> ", popupWhen, "<br><br>",
"<b>", monitor$meta$deviceDeploymentID, "</b><br>",
"locationID = ", monitor$meta$locationID, "<br>",
"deviceID = ", monitor$meta$deviceID, "<br><br>",
monitor$meta$countyName, " County, ", monitor$meta$stateCode, "<br>",
"timezone = ", monitor$meta$timezone, "<br>",
"longitude = ", monitor$meta$longitude, ", ", "latitude = ", monitor$meta$latitude, "<br>"
)
# Add extra vars
for ( i in seq_along(popupText) ) {
extraText <- vector("character", length(extraVars))
for ( j in seq_along(extraVars) ) {
var <- extraVars[j]
extraText[j] <- paste0(var, " = ", monitor$meta[i, var], "<br>")
}
extraText <- paste0(extraText, collapse = "")
popupText[i] <- paste0(popupText[i], "<hr>", extraText)
}
monitor$meta$popupText <- popupText
# ----- Create leaflet map ---------------------------------------------------
# Filter out missing location data
monitor$meta <-
monitor$meta %>%
dplyr::filter(!is.na(.data$longitude)) %>%
dplyr::filter(!is.na(.data$latitude))
# Spread out locations if requested
if ( jitter > 0 ) {
monitor$meta <-
monitor$meta %>%
dplyr::mutate(
longitude = jitter(.data$longitude, amount = 5e-4),
latitude = jitter(.data$latitude, amount = 5e-4)
)
}
# Convert maptype to a character string that addProviderTiles can read
if ( missing(maptype) || maptype == "terrain") {
providerTiles <- "Esri.WorldTopoMap"
} else if ( maptype == "roadmap" ) {
providerTiles <- "OpenStreetMap"
} else if ( maptype == "toner" ) {
providerTiles <- "Stamen.Toner"
} else if (maptype == "satellite" ) {
providerTiles <- "Esri.WorldImagery"
} else {
providerTiles <- maptype
}
# Determine appropriate zoom level
lonRange <- range(monitor$meta$longitude, na.rm = TRUE)
latRange <- range(monitor$meta$latitude, na.rm = TRUE)
# Create leaflet map
leafletMap <-
leaflet::leaflet(monitor$meta) %>%
leaflet::fitBounds(lonRange[1], latRange[1], lonRange[2], latRange[2]) %>%
leaflet::addProviderTiles(providerTiles) %>%
leaflet::addCircleMarkers(
lat = ~latitude,
lng = ~longitude,
radius = radius,
fillColor = cols,
fillOpacity = opacity,
stroke = FALSE,
popup = monitor$meta$popupText,
...) %>%
leaflet::addLegend(
position = "bottomright",
colors = rev(colors), # show low levels at the bottom
labels = rev(labels), # show low levels at the bottom
opacity = 1,
title = legendTitle
)
return(leafletMap)
}
# ===== DEBUGGING ==============================================================
if ( FALSE ) {
library(AirMonitor)
archiveBaseUrl <- "https://airfire-data-exports.s3.us-west-2.amazonaws.com/monitoring/v2"
monitor <-
airnow_loadDaily(archiveBaseUrl = archiveBaseUrl) %>%
##monitor_select(c("089a067f92712ad1_530750003", "d121a99bc6c2ac7f_160570005"))
monitor_filter(countryCode %in% c("US", "CA", "MX"))
slice = "max"
radius = 10
opacity = 0.7
maptype = "terrain"
extraVars = NULL
jitter = 5e-4
monitor_leaflet(
monitor,
slice = "max",
radius = 10,
opacity = 0.7,
maptype = "terrain",
extraVars = NULL,
jitter = 5e-4,
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.