R/monitor_rollingMeanPlot.R

In PWFSLSmoke: Utilities for Working with Air Quality Monitoring Data

#' @title Create Rolling Mean Plot
#'
#' @description
#' Creates a plot of individual (e.g. hourly) and rolling mean PM2.5 values for
#' a specific monitor.
#'
#' @param ws_monitor \emph{ws_monitor} object
#' @param monitorID Monitor ID for a specific monitor in the ws_monitor object
#'   (optional if only one monitor in the ws_monitor object).
#' @param width Number of periods to average (e.g. for hourly data, \code{width
#'   = 24} plots 24-hour rolling means).
#' @param align Alignment of averaging window relative to point being
#'   calculated; one of \code{"left|center|right"}.
#' @param data.thresh Minimum number of valid observations required as a percent
#'   of \code{width}; NA is returned if insufficicnet valid data to calculate.
#'   mean
#' @param ylim y limits for the plot.
#' @param tlim Optional vector with start and end times (integer or character
#'   representing YYYYMMDD[HH]).
#' @param localTime Logical specifying whether \code{tlim} is in local time or
#'   UTC.
#' @param shadedNight Add nighttime shading.
#' @param aqiLines Horizontal lines indicating AQI levels.
#' @param gridHorizontal Add dashed horizontal grid lines.
#' @param grid24hr Add dashed grid lines at day boundaries.
#' @param grid3hr Add dashed grid lines every 3 hours.
#' @param showLegend Include legend in top left.
#'
#' @details
#' \itemize{
#'   \item{\code{align = "left"}: Forward roll, using hour of interest and the
#'     (\code{width}-1) subsequent hours (e.g. 3-hr left-aligned roll for Hr 5
#'     will consist of average of Hrs 5, 6 and 7)
#'   }
#'   \item{\code{align = "right"}: Backwards roll, using hour of interest and
#'     the (\code{width}-1) prior hours (e.g. 3-hr right-aligned roll for Hr 5
#'     will consist of average of Hrs 3, 4 and 5)
#'   }
#'   \item{\code{align = "center"} for odd \code{width}: Average of hour of
#'     interest and (\code{width}-1)/2 on either side (e.g. 3-hr center-aligned
#'     roll for Hr 5 will consist of average of Hrs 4, 5 and 6)
#'   }
#'   \item{\code{align = "center"} for even \code{width}: Average of hour of
#'     interest and (\code{width}/2)-1 hours prior and \code{width}/2 hours
#'     after (e.g. 4-hr center-aligned roll for Hr 5 will consist of average of
#'     Hrs 4, 5, 6 and 7)
#'   }
#' }
#'
#' @note
#' This function attempts to provide a 'publication ready' rolling mean plot.
#'
#' @import graphics
#' @export
#' @keywords ws_monitor
#'
#' @examples
#' library(PWFSLSmoke)
#'
#' N_M <- Northwest_Megafires
#' Roseburg <- monitor_subset(N_M, tlim = c(20150821, 20150831),
#'                            monitorIDs = c("410190002_01"))
#' monitor_rollingMeanPlot(Roseburg, shadedNight = TRUE)

monitor_rollingMeanPlot <- function(
  ws_monitor,
  monitorID = NULL,
  width = 3,
  align = "center",
  data.thresh = 75,
  tlim = NULL,
  ylim = NULL,
  localTime = TRUE,
  shadedNight = FALSE,
  aqiLines = TRUE,
  gridHorizontal = FALSE,
  grid24hr = FALSE,
  grid3hr = FALSE,
  showLegend = TRUE
) {

  # Sanity check
  if ( monitor_isEmpty(ws_monitor) ) {
    stop("ws_monitor object contains zero monitors")
  }

  # ----- Style ---------------------------------------------------------------

  # points
  col_points <- adjustcolor("black", 0.2)
  pch_points <- 16
  cex_points <- .8

  # rolling mean
  col_mean <- "black"
  lwd_mean <- 1
  lty_mean <- "solid"

  # Grid lines
  lty_grid <- "dotted"
  col_grid <- "gray80"
  lwd_grid <- 1.6
  lwd_gridLight <- 0.8 # for 3 hour intervals

  # AQI lines
  lwd_aqi <- 6
  col_aqi <- adjustcolor(AQI$colors[2:6], 0.6)

  # ----- Data Preparaion ------------------------

  # Roll direction for legend
  if ( align == "left" ) {
    direction <- "Forward-"
  } else if ( align == "right" ) {
    direction <- "Backward-"
  } else {
    direction <- "Centered "
  }

  # Allow single monitor objects to be used without specifying monitorID
  if ( is.null(monitorID) ) {
    if ( nrow(ws_monitor$meta) == 1 ) {
      monitorID <- ws_monitor$meta$monitorID[1]
    } else {
      stop(
        paste0(
          "ws_monitor object contains data for >1 monitor. Please specify a monitorID from: '",
          paste(ws_monitor$meta$monitorID, collapse = "', '"),
          "'"
        )
      )
    }
  }

  # Subset ws_monitor object to monitor of interest and pull out meta & data
  ws_monitor <- PWFSLSmoke::monitor_subset(ws_monitor, monitorIDs = monitorID)
  meta <- ws_monitor$meta
  data <- ws_monitor$data

  # Pull out hour averages and create rolling means for dataframe prior to subsetting by tlim
  hourAvgs <- data[[monitorID]]
  rollingMeans <-
    PWFSLSmoke::monitor_rollingMean(
      ws_monitor,
      width = width,
      data.thresh = data.thresh,
      align = align
    )$data[[monitorID]]

  # Assign timeStamp based on localTime setting
  if ( localTime ) {
    timeStamp <- lubridate::with_tz(data$datetime, tzone = meta$timezone)
    tzLabel <- "(local)"
  } else {
    timeStamp <- data$datetime
    tzLabel <- "(UTC)"
  }

  # put the pieces of the new dataframe together
  df <- data.frame(timeStamp, hourAvgs, rollingMeans)

  # Time limit application
  if ( !is.null(tlim) ) {
    # TODO: Warn if no data for any dates within tlim?
    # TODO: add logic to check for tlim format
    timeMask <-
      timeStamp >= lubridate::ymd(tlim[1], tz = "UTC") &
      timeStamp < lubridate::ymd(tlim[2], tz = "UTC") + lubridate::days(1)

    if ( sum(timeMask) == 0 ) {
      monitor_noDataPlot(ws_monitor)
      stop("No data contained within specified time limits, please try again.")
    }
    df <- df[timeMask, ]
  }

  # Set y Limits
  if ( is.null(ylim) ) {
    # ylim <- c(min(0, min(df$hourAvgs, na.rm = TRUE)), max(df$hourAvgs, na.rm = TRUE))
    ylim <- c(0, max(c(df$hourAvgs, 275), na.rm = TRUE))
  }

  # Date ranges for grid line locatioins and chart title
  minTime <- df$timeStamp[1]
  maxTime <- utils::tail(df$timeStamp, 1)

  # Find first full day for grid lines and time axis labels
  min24Hour <- timeStamp[which(lubridate::hour(timeStamp) == 0)[1]]

  # Find start of first 3-hour period (e.g. Hr 0, 3, 6...)
  threeHrRemainder <- lubridate::hour(minTime) / 3 - floor(lubridate::hour(minTime) / 3)
  if ( threeHrRemainder == 0 ) { # first time stamp is right on a 3-hr marker
    min3Hour <- df$timeStamp[1]
  } else if ( threeHrRemainder < .5 ) { # first time stamp is one hour after a 3-hr marker
    min3Hour <- df$timeStamp[1] - lubridate::hours(1)
  } else {
    min3Hour <- df$timeStamp[1] - lubridate::hours(2)
  }

  # Define grid line locations
  xGrid24Hour <- seq(min24Hour, maxTime + lubridate::hours(1), "day")
  xGrid3Hour <- seq(min3Hour, maxTime + lubridate::hours(1), "3 hour")
  if ( grid24hr == TRUE) {
    xGrid3Hour <- xGrid3Hour[!(xGrid3Hour %in% xGrid24Hour)]
  }

  # Assign time axis label
  if ( nrow(df) <= 24 ) {
    xlab <- "Hour"
  } else {
    xlab <- "Date"
  }
  xlab <- paste(xlab, tzLabel)

  # ----- Plotting --------------------------------

  # Create blank plot
  plot(df$timeStamp, df$hourAvgs,
       type = "n",
       axes = FALSE,
       xlab = xlab,
       ylab = "",
       ylim = ylim
  )

  # Shaded Night; breaks if >1 time zone
  if ( shadedNight ) {
    if ( length(unique(meta$timezone)) > 1) {
      stop("Can't do shaded night for more than one time zone!!")
    } else {
      lon <- meta$longitude
      lat <- meta$latitude
      timezone <- meta$timezone
      if ( localTime ) {
        timeInfo <- PWFSLSmoke::timeInfo(timeStamp, lon, lat, timezone)
        PWFSLSmoke::addShadedNight(timeInfo)
      } else {
        timeInfo <- PWFSLSmoke::timeInfo(data$datetime, lon, lat, timezone)
        PWFSLSmoke::addShadedNight(timeInfo)
      }
    }
  }

  # AQI Lines
  if ( aqiLines ) {
    abline(h = AQI$breaks_24[2:6], col = col_aqi, lwd = lwd_aqi)
  }

  # Add grid lines
  if ( gridHorizontal ) {  # Horizontal
    #horizontal lines
    grid(nx = NA, ny = NULL, col = col_grid, lwd = lwd_grid, lty = lty_grid)
  }
  if ( grid3hr ) {  # 3-hour interval
    # Vertical lines - 3 Hr intervals
    abline(v = xGrid3Hour, col = col_grid, lwd = lwd_gridLight, lty = lty_grid)
  }
  if ( grid24hr ) {   # 24-hour interval
    abline(v = xGrid24Hour, col = col_grid, lwd = lwd_grid, lty = lty_grid)
  }

  # Box the plot area
  box()

  # Plot light grey circles for actual PM2.5 data
  points(
    df$timeStamp,
    df$hourAvgs,
    col = col_points,
    pch = pch_points,
    cex = cex_points
  )

  # Add rolling means
  lines(df$timeStamp, df$rollingMeans, col = "black", lwd = lwd_mean, lty = lty_mean)

  # ----- Annotations ---------------------

  ## Horizontal axis (see "Assign timeStamp based on localTime setting" above
  #  for horizontal axis label)
  if (nrow(df) <= 24) {
    axis.POSIXct(
      1,
      at = seq(minTime, maxTime + lubridate::hours(1), by = "3 hour"),
      format = "%H"
    )
  } else {
    axis.POSIXct(
      1,
      at = seq(min24Hour, maxTime + lubridate::hours(1), by = "1 day"),
      format = "%b %d"
    )
  }

  # Vertical axis
  axis(2, las = 1)
  mtext(expression(paste("PM"[2.5] * " (", mu, "g/m"^3 * ")")), 2, line = 2.5)

  # Title
  mtext(expression(paste("Hourly and Rolling Mean PM"[2.5])), line = 2, cex = 1.5)
  if ( lubridate::date(minTime) == lubridate::date(maxTime) ) {
    mtext(format(minTime, "%b. %d %Y"), line = .7, cex = 1.5)
  } else {
    mtext(paste(format(minTime, "%b. %d - "), format(maxTime, "%b. %d %Y")), line = .7, cex = 1.5)
  }

  # Legend
  if ( showLegend ) {
    legend(
      "topleft",
      legend = c("Hourly Averages", paste0(width,  "-hour ", direction, "Rolling Mean")),
      cex = par("cex.lab"),
      col = c(col_points, col_mean),
      lwd = c(NA, lwd_mean),
      lty = c(NA, lty_mean),
      pch = c(pch_points, NA))
  }

}