R/block_average.R

In robwschlegel/heatwaveR: Detect Heatwaves and Cold-Spells

#' Calculate yearly means for event metrics.
#'
#' @param data Accepts the data returned by the \code{\link{detect_event}} function.
#' @param x This column is expected to contain a vector of dates as per the
#' specification of \code{ts2clm}. If a column headed \code{t} is present in
#' the dataframe, this argument may be omitted; otherwise, specify the name of
#' the column with dates here.
#' @param y This is a column containing the measurement variable. If the column
#' name differs from the default (i.e. \code{temp}), specify the name here.
#' @param report Specify either \code{full} or \code{partial}. Selecting \code{full} causes
#' the report to contain NAs for any years in which no events were detected
#' (except for \code{count}, which will be zero in those years), while \code{partial}
#' reports only the years wherein events were detected. The default is \code{full}.
#' @param returnDF The default (\code{TRUE}) tells the function to return the results as
#' type \code{data.frame}. \code{FALSE} will return the results as a \code{data.table}.
#'
#' @details
#' This function needs to be provided with the full output from the \code{detect_event}
#' or \code{exceedance} functions. Note that the yearly averages are calculated only for
#' complete years (i.e. years that start/end part-way through the year at the beginning
#' or end of the original time series are removed from the calculations).
#'
#' This function differs from the python implementation of the function of the
#' same name (i.e., \code{blockAverage}, see \url{https://github.com/ecjoliver/marineHeatWaves})
#' in that we only provide the ability to calculate the average (or aggregate)
#' event metrics in 'blocks' of one year, while the python version allows
#' arbitrary (integer) block sizes.
#'
#' Note that if this function is used on the output of \code{exceedance}, all of the metrics
#' (see below) with \code{relThresh} in the name will be returned as \code{NA} values.
#'
#' @return The function will return a data frame of the averaged (or aggregate)
#' metrics. It includes the following:
#'   \item{year}{The year over which the metrics were averaged.}
#'   \item{count}{The number of events per year.}
#'   \item{duration}{The average duration of events per year [days].}
#'   \item{duration_max}{The maximum duration of an event in each year [days].}
#'   \item{intensity_mean}{The average event "mean intensity" in each year [deg. C].}
#'   \item{intensity_max}{The average event "maximum (peak) intensity" in each year
#'   [deg. C].}
#'   \item{intensity_max_max}{The maximum event "maximum (peak) intensity" in
#'   each year [deg. C].}
#'   \item{intensity_var}{The average event "intensity variability" in each year
#'   [deg. C].}
#'   \item{intensity_cumulative}{The average event "cumulative intensity" in each year
#'   [deg. C x days].}
#'   \item{rate_onset}{Average event onset rate in each year [deg. C / days].}
#'   \item{rate_decline}{Average event decline rate in each year [deg. C / days].}
#'   \item{total_days}{Total number of events days in each year [days].}
#'   \item{total_icum}{Total cumulative intensity over all events in each year [deg. C x days].}
#'
#' \code{intensity_max_relThresh}, \code{intensity_mean_relThresh},
#' \code{intensity_var_relThresh}, and \code{intensity_cumulative_relThresh}
#' are as above except relative to the threshold (e.g., 90th percentile) rather
#' than the seasonal climatology.
#'
#' \code{intensity_max_abs}, \code{intensity_mean_abs}, \code{intensity_var_abs}, and
#' \code{intensity_cumulative_abs} are as above except as absolute magnitudes
#' rather than relative to the seasonal climatology or threshold.
#'
#' @author Albertus J. Smit, Eric C. J. Oliver, Robert W. Schlegel
#'
#' @references Hobday, A.J. et al. (2016), A hierarchical approach to defining
#' marine heatwaves, Progress in Oceanography, 141, pp. 227-238,
#' doi: 10.1016/j.pocean.2015.12.014
#'
#' @export
#'
#' @examples
#' ts <- ts2clm(sst_WA, climatologyPeriod = c("1983-01-01", "2012-12-31"))
#' res <- detect_event(ts)
#' out <- block_average(res)
#' summary(glm(count ~ year, out, family = "poisson"))
#'
#' library(ggplot2)
#'
#' ggplot(data = out, aes(x = year, y = count)) +
#'   geom_point(colour = "salmon") +
#'   geom_line() +
#'   labs(x = NULL, y = "Number of events")
#'
block_average <- function(data,
                          x = t,
                          y = temp,
                          report = "full",
                          returnDF = TRUE) {

  if("climatology" %in% names(data)) {
    clim <- data$climatology
    ts_x <- eval(substitute(x), data$climatology)
    ts_y <- eval(substitute(y), data$climatology)
    event <- data$event
  } else if("threshold" %in% names(data)) {
    clim <- data$threshold
    ts_x <- eval(substitute(x), data$threshold)
    ts_y <- eval(substitute(y), data$threshold)
    event <- data$exceedance
  } else {
    stop("Please ensure you are providing this function with the output of either detect_event() or exceedance().")
  }

  clim$t <- ts_x
  clim$temp <- ts_y
  event$year <- as.numeric(format(event$date_start, "%Y"))

  year <- temp <- date_start <- temp_mean <- temp_min <- temp_max <- NULL

  temp_yr <- data.frame(year = as.numeric(base::unique(format(clim$t, "%Y"))))

  duration <- duration_max <- count <- intensity_mean <- intensity_max <- intensity_var <- intensity_cumulative <-
    intensity_mean_relThresh <- intensity_max_relThresh <- intensity_var_relThresh <-
    intensity_cumulative_relThresh <- intensity_mean_abs <- intensity_max_abs <- intensity_var_abs <-
    intensity_cumulative_abs <- rate_onset <- rate_decline <- total_days <- total_icum <-
    duration_mean <- intensity_max_mean <- intensity_cumulative_mean <- intensity_cum_abs <-  NULL

  event_block <- data.table::setDT(event)
  if ("intensity_cumulative_relThresh" %in% names(event_block)) {
    event_block <- event_block[, list(
      count = length(duration),
      duration_mean = mean(duration, na.rm = T),
      duration_max = max(duration, na.rm = T),
      intensity_mean = mean(intensity_mean, na.rm = T),
      intensity_max_mean = mean(intensity_max, na.rm = T),
      intensity_max_max = max(intensity_max, na.rm = T),
      intensity_var = mean(intensity_var, na.rm = T),
      intensity_cumulative_mean = mean(intensity_cumulative, na.rm = T),
      intensity_mean_relThresh = mean(intensity_mean_relThresh, na.rm = T),
      intensity_max_relThresh = mean(intensity_max_relThresh, na.rm = T),
      intensity_var_relThresh = mean(intensity_var_relThresh, na.rm = T),
      intensity_cumulative_relThresh = mean(intensity_cumulative_relThresh, na.rm = T),
      intensity_mean_abs = mean(intensity_mean_abs, na.rm = T),
      intensity_max_abs = mean(intensity_max_abs, na.rm = T),
      intensity_var_abs = mean(intensity_var_abs, na.rm = T),
      intensity_cumulative_abs = mean(intensity_cumulative_abs, na.rm = T),
      rate_onset = mean(rate_onset, na.rm = T),
      rate_decline = mean(rate_decline, na.rm = T),
      total_days = sum(duration, na.rm = T),
      total_icum = sum(intensity_cumulative, na.rm = T)
    ), by = list(year)]

  } else {
    event_block <-  event_block[, list(
      count = length(duration),
      duration_mean = mean(duration, na.rm = T),
      duration_max = max(duration, na.rm = T),
      intensity_mean = mean(intensity_mean, na.rm = T),
      intensity_max_mean = mean(intensity_max, na.rm = T),
      intensity_max_max = max(intensity_max, na.rm = T),
      intensity_var = mean(intensity_var, na.rm = T),
      intensity_cumulative_mean = mean(intensity_cumulative, na.rm = T),
      intensity_mean_abs = mean(intensity_mean_abs, na.rm = T),
      intensity_max_abs = mean(intensity_max_abs, na.rm = T),
      intensity_var_abs = mean(intensity_var_abs, na.rm = T),
      intensity_cum_abs = mean(intensity_cum_abs, na.rm = T),
      rate_onset = mean(rate_onset, na.rm = T),
      rate_decline = mean(rate_decline, na.rm = T),
      total_days = sum(duration, na.rm = T),
      total_icum = sum(intensity_cumulative, na.rm = T)
    ), by = list(year)]
  }

  colnames(event_block)[c(3, 6, 9)] <- c("duration", "intensity_max", "intensity_cumulative")

  if (report == "full") {
    event_block <- base::merge(temp_yr, event_block, by = "year", all.x = TRUE)
  } else if (report == "partial") {
    event_block <- base::merge(temp_yr, event_block, by = "year", all.y = TRUE)
  } else stop("Oops, 'report' must be either 'full' or 'partial'!")

  event_block$count[is.na(event_block$count)] <- 0

  if(returnDF) {
    data.table::setDF(event_block)
  } else {
    data.table::setDT(event_block)
  }
  return(event_block)
}