R/plot.daily_output.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Plot method for daily dendrometer statistics
#'
#' @description
#' Unified S3 plotting method for objects returned by \code{daily.data()}.
#' Missing grouping labels are removed before grouped boxplots are drawn so
#' \code{NA} does not appear on the x-axis.
#'
#' @param x Object of class \code{"daily_output"}.
#' @param y Unused.
#' @param Year Optional numeric year or vector of years for subsetting.
#' @param DOY Optional numeric vector of length 2 giving the start and end
#'   day-of-year for subsetting.
#' @param type Plot type. One of
#'   \code{"summary"}, \code{"amplitude"}, \code{"minmax_ribbon"},
#'   \code{"timing"}, \code{"timing_violin"}, \code{"lag"}, \code{"change"},
#'   \code{"boxplot"}, or \code{"heatmap"}.
#' @param stat Statistic used for \code{type = "boxplot"} or
#'   \code{type = "heatmap"}. One of
#'   \code{"amplitude"}, \code{"Max_diff"}, \code{"mean"},
#'   \code{"median"}, \code{"Max"}, \code{"Min"},
#'   \code{"lag_h"}, \code{"Time_min_h"}, \code{"Time_max_h"}.
#' @param by Grouping for \code{type = "boxplot"}:
#'   \code{"month"}, \code{"month_of_year"}, or \code{"year"}.
#'   \code{"month"} uses a continuous calendar-date axis spanning the selected
#'   data window and supports multi-year plotting, for example 2022--2024.
#'   Empty calendar months may appear as gaps, but months with available data
#'   are not clipped at the start or end of the axis.
#' @param box_style Style for \code{type = "boxplot"}:
#'   \code{"boxplot"}, \code{"violin"}, or \code{"both"}.
#' @param status_cols Colors for \code{Day_status}.
#' @param timing_segment_cols Named vector of two colors for the timing connector:
#'   one for \code{"max after min"} and one for \code{"max before/equal min"}.
#' @param ... Unused.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @details
#' Plot types:
#' \itemize{
#'   \item \code{"summary"}: daily Min, mean, median, and Max.
#'   \item \code{"amplitude"}: daily amplitude through time.
#'   \item \code{"minmax_ribbon"}: ribbon between daily Min and Max.
#'   \item \code{"timing"}: per-day line connecting daily minimum and maximum timing.
#'   \item \code{"timing_violin"}: violin distribution of \code{Time_min_h}
#'         and \code{Time_max_h} for the selected window.
#'   \item \code{"lag"}: signed lag between daily maximum and minimum time.
#'   \item \code{"change"}: day-to-day change in daily maximum
#'         \code{Max_diff}.
#'   \item \code{"boxplot"}: grouped distributions of a chosen daily statistic.
#'         For \code{by = "month"}, monthly groups are plotted on a continuous
#'         date axis. The axis limits are extended beyond the first and last
#'         plotted month so edge months are not clipped.
#'   \item \code{"heatmap"}: year x day-of-year heatmap of a chosen daily statistic.
#' }
#'
#' For \code{type = "timing"}, the red and blue points are not connected across
#' days. Instead, each day has its own segment connecting minimum and maximum
#' time.
#'
#' @examples
#' \donttest{
#' data(nepa17)
#' dd <- daily.data(df = nepa17[1:1000, ], TreeNum = 1)
#'
#' plot(dd)
#' plot(dd, type = "timing")
#' plot(dd, type = "timing", Year = 2017, DOY = c(1, 6))
#' plot(dd, type = "timing_violin", Year = 2017, DOY = c(1, 6))
#' plot(dd, type = "boxplot", stat = "amplitude", by = "month")
#' plot(dd, type = "boxplot", stat = "amplitude", by = "month_of_year",
#'      box_style = "both")
#' plot(dd, type = "change")
#' }
#'
#' @method plot daily_output
#' @importFrom dplyr mutate filter %>%
#' @importFrom tidyr pivot_longer
#' @importFrom lubridate yday year month floor_date
#' @importFrom ggplot2 ggplot aes geom_line geom_ribbon geom_col geom_point
#'   geom_segment geom_boxplot geom_violin geom_tile geom_hline
#'   theme_minimal theme element_text labs scale_color_manual
#'   scale_fill_manual scale_y_continuous scale_x_continuous scale_x_discrete
#'   scale_x_date position_jitter expansion
#' @export
plot.daily_output <- function(x,
                              y = NULL,
                              Year = NULL,
                              DOY = NULL,
                              type = c("summary", "amplitude", "minmax_ribbon",
                                       "timing", "timing_violin", "lag", "change",
                                       "boxplot", "heatmap"),
                              stat = c("amplitude", "Max_diff", "mean",
                                       "median", "Max", "Min",
                                       "lag_h", "Time_min_h", "Time_max_h"),
                              by = c("month", "month_of_year", "year"),
                              box_style = c("boxplot", "violin", "both"),
                              status_cols = c(
                                growing = "forestgreen",
                                shrinking = "firebrick",
                                stable = "grey60"
                              ),
                              timing_segment_cols = c(
                                "max after min" = "grey65",
                                "max before/equal min" = "black"
                              ),
                              ...) {

  DATE <- value <- Statistic <- Hour <- group_label <- Day_status <- DOYv <- Yearv <- NULL
  Time_min_h <- Time_max_h <- order_flag <- NULL
  Min <- Max <- amplitude <- Remarks <- lag_h <- Max_diff <- NULL
  group_date <- month_index <- NULL

  if (!inherits(x, "daily_output")) {
    stop("Input must be an object of class 'daily_output'.")
  }

  type <- match.arg(type)
  stat <- match.arg(stat)
  by <- match.arg(by)
  box_style <- match.arg(box_style)

  .subset_daily <- function(dat, Year = NULL, DOY = NULL) {
    out <- dat

    if (!is.null(Year)) {
      out <- out[lubridate::year(out$DATE) %in% Year, , drop = FALSE]
    }

    if (!is.null(DOY)) {
      if (length(DOY) != 2) {
        stop("'DOY' must be a numeric vector of length 2.")
      }

      doy_min <- min(DOY)
      doy_max <- max(DOY)

      out <- out[
        lubridate::yday(out$DATE) >= doy_min &
          lubridate::yday(out$DATE) <= doy_max,
        ,
        drop = FALSE
      ]
    }

    out <- out[order(out$DATE), , drop = FALSE]

    if (nrow(out) == 0) {
      stop("No data available for the selected Year/DOY window.")
    }

    out
  }

  dat <- .subset_daily(x, Year = Year, DOY = DOY)

  if (type == "summary") {
    pdat <- tidyr::pivot_longer(
      dat,
      cols = c(Min, mean, median, Max),
      names_to = "Statistic",
      values_to = "value"
    )

    p <- ggplot2::ggplot(
      pdat,
      ggplot2::aes(x = DATE, y = value, color = Statistic)
    ) +
      ggplot2::geom_line(linewidth = 0.6) +
      ggplot2::scale_color_manual(
        values = c(
          Min = "firebrick",
          mean = "black",
          median = "grey40",
          Max = "steelblue"
        )
      ) +
      ggplot2::labs(
        x = "Date",
        y = "Daily statistic",
        color = "Statistic"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "amplitude") {
    p <- ggplot2::ggplot(dat, ggplot2::aes(x = DATE, y = amplitude)) +
      ggplot2::geom_col(fill = "steelblue", alpha = 0.8) +
      ggplot2::geom_line(color = "black", linewidth = 0.35) +
      ggplot2::labs(
        x = "Date",
        y = "Daily amplitude"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "minmax_ribbon") {
    p <- ggplot2::ggplot(dat, ggplot2::aes(x = DATE)) +
      ggplot2::geom_ribbon(
        ggplot2::aes(ymin = Min, ymax = Max),
        fill = "grey80",
        alpha = 0.8
      ) +
      ggplot2::geom_line(
        ggplot2::aes(y = Min),
        color = "firebrick",
        linewidth = 0.5
      ) +
      ggplot2::geom_line(
        ggplot2::aes(y = Max),
        color = "steelblue",
        linewidth = 0.5
      ) +
      ggplot2::geom_line(
        ggplot2::aes(y = mean),
        color = "black",
        linewidth = 0.6
      ) +
      ggplot2::labs(
        x = "Date",
        y = "Daily Min-Max envelope"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "timing") {
    tdat <- dat %>%
      dplyr::mutate(
        order_flag = ifelse(
          Remarks == "*",
          "max after min",
          "max before/equal min"
        )
      )

    p <- ggplot2::ggplot(tdat, ggplot2::aes(x = DATE)) +
      ggplot2::geom_segment(
        ggplot2::aes(
          xend = DATE,
          y = Time_min_h,
          yend = Time_max_h,
          color = order_flag
        ),
        linewidth = 0.55
      ) +
      ggplot2::geom_point(
        ggplot2::aes(y = Time_min_h),
        color = "firebrick",
        size = 1.9
      ) +
      ggplot2::geom_point(
        ggplot2::aes(y = Time_max_h),
        color = "steelblue",
        size = 1.9
      ) +
      ggplot2::scale_color_manual(
        values = timing_segment_cols,
        name = "Order"
      ) +
      ggplot2::scale_y_continuous(
        limits = c(0, 24),
        breaks = seq(0, 24, 4)
      ) +
      ggplot2::labs(
        x = "Date",
        y = "Hour of day"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        legend.position = "top",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "timing_violin") {
    pdat <- tidyr::pivot_longer(
      dat,
      cols = c(Time_min_h, Time_max_h),
      names_to = "Statistic",
      values_to = "Hour"
    ) %>%
      dplyr::filter(!is.na(Hour))

    if (nrow(pdat) == 0) {
      stop("No timing values available for the selected Year/DOY window.")
    }

    p <- ggplot2::ggplot(
      pdat,
      ggplot2::aes(x = Statistic, y = Hour, fill = Statistic)
    ) +
      ggplot2::geom_violin(alpha = 0.7, trim = FALSE) +
      ggplot2::geom_point(
        position = ggplot2::position_jitter(width = 0.08, height = 0),
        size = 1,
        alpha = 0.35
      ) +
      ggplot2::scale_fill_manual(
        values = c(
          Time_min_h = "firebrick",
          Time_max_h = "steelblue"
        ),
        labels = c(
          Time_min_h = "Time of minimum",
          Time_max_h = "Time of maximum"
        )
      ) +
      ggplot2::scale_y_continuous(
        limits = c(0, 24),
        breaks = seq(0, 24, 4)
      ) +
      ggplot2::labs(
        x = NULL,
        y = "Hour of day",
        fill = "Timing"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        legend.position = "top",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "lag") {
    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = DATE, y = lag_h, fill = lag_h >= 0)
    ) +
      ggplot2::geom_col(show.legend = FALSE) +
      ggplot2::geom_hline(yintercept = 0, linetype = 2) +
      ggplot2::scale_fill_manual(
        values = c(`TRUE` = "steelblue", `FALSE` = "firebrick")
      ) +
      ggplot2::labs(
        x = "Date",
        y = "Lag (hours): Time_max - Time_min"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "change") {
    dat$Day_status <- factor(
      dat$Day_status,
      levels = c("shrinking", "stable", "growing")
    )

    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = DATE, y = Max_diff, fill = Day_status)
    ) +
      ggplot2::geom_col(na.rm = TRUE) +
      ggplot2::geom_hline(yintercept = 0, linetype = 2) +
      ggplot2::scale_fill_manual(
        values = status_cols,
        na.value = "grey80",
        drop = FALSE
      ) +
      ggplot2::labs(
        x = "Date",
        y = "Change in daily maximum vs previous day",
        fill = "Day status"
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "boxplot") {
    bdat <- dat
    use_continuous_time <- FALSE

    if (by == "month") {
      use_continuous_time <- TRUE

      bdat$group_date <- as.Date(lubridate::floor_date(bdat$DATE, "month"))

      month_seq <- seq(
        min(bdat$group_date, na.rm = TRUE),
        max(bdat$group_date, na.rm = TRUE),
        by = "month"
      )

      n_months <- length(month_seq)
      n_years <- length(unique(lubridate::year(month_seq)))

      x_lab <- "Date"

      box_width <- 20
      point_jitter_width <- 5

      date_breaks <- if (n_months <= 24) {
        "1 month"
      } else if (n_months <= 60) {
        "3 months"
      } else if (n_months <= 120) {
        "6 months"
      } else {
        "1 year"
      }

      date_labels <- if (n_years > 1) {
        "%Y-%m"
      } else {
        "%b"
      }

      x_limits <- c(
        min(month_seq, na.rm = TRUE) - 15,
        max(month_seq, na.rm = TRUE) + 35
      )

    } else if (by == "month_of_year") {
      bdat$group_label <- factor(
        as.character(lubridate::month(bdat$DATE, label = TRUE, abbr = TRUE)),
        levels = month.abb
      )

      x_lab <- "Month of year"

    } else {
      full_levels <- as.character(sort(unique(lubridate::year(bdat$DATE))))

      bdat$group_label <- factor(
        as.character(lubridate::year(bdat$DATE)),
        levels = full_levels
      )

      x_lab <- "Year"
    }

    bdat$value <- bdat[[stat]]

    if (use_continuous_time) {
      bdat <- bdat %>%
        dplyr::filter(!is.na(value), !is.na(group_date))
    } else {
      bdat <- bdat %>%
        dplyr::filter(!is.na(value), !is.na(group_label))
    }

    if (nrow(bdat) == 0) {
      stop("No non-missing values available for the selected statistic.")
    }

    if (use_continuous_time) {
      counts <- table(format(bdat$group_date, "%Y-%m"))
    } else {
      counts <- table(bdat$group_label)
    }

    if (all(counts <= 1)) {
      warning("Each group has only one value; distributions collapse to lines/points.")
    }

    if (use_continuous_time) {
      p <- ggplot2::ggplot(
        bdat,
        ggplot2::aes(x = group_date, y = value, group = group_date)
      )

      if (box_style %in% c("violin", "both")) {
        p <- p +
          ggplot2::geom_violin(
            fill = "grey80",
            alpha = 0.8,
            trim = FALSE
          )
      }

      if (box_style %in% c("boxplot", "both")) {
        p <- p +
          ggplot2::geom_boxplot(
            width = if (box_style == "both") 8 else box_width,
            outlier.shape = NA,
            fill = if (box_style == "both") "white" else "grey80"
          )
      }

      p <- p +
        ggplot2::geom_point(
          position = ggplot2::position_jitter(
            width = point_jitter_width,
            height = 0
          ),
          size = 1.2,
          alpha = 0.45,
          color = "grey25"
        ) +
        ggplot2::scale_x_date(
          limits = x_limits,
          date_breaks = date_breaks,
          date_labels = date_labels,
          expand = ggplot2::expansion(mult = c(0, 0))
        ) +
        ggplot2::labs(
          x = x_lab,
          y = stat
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme(
          axis.text.x = ggplot2::element_text(
            angle = 90,
            hjust = 1,
            vjust = 0.5
          ),
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 11)
        )

    } else {
      p <- ggplot2::ggplot(
        bdat,
        ggplot2::aes(x = group_label, y = value)
      )

      if (box_style %in% c("violin", "both")) {
        p <- p +
          ggplot2::geom_violin(
            fill = "grey80",
            alpha = 0.8,
            trim = FALSE
          )
      }

      if (box_style %in% c("boxplot", "both")) {
        p <- p +
          ggplot2::geom_boxplot(
            width = if (box_style == "both") 0.18 else 0.65,
            outlier.shape = NA,
            fill = if (box_style == "both") "white" else "grey80"
          )
      }

      p <- p +
        ggplot2::geom_point(
          position = ggplot2::position_jitter(width = 0.12, height = 0),
          size = 1.2,
          alpha = 0.45,
          color = "grey25"
        ) +
        ggplot2::scale_x_discrete(drop = FALSE) +
        ggplot2::labs(
          x = x_lab,
          y = stat
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme(
          axis.text.x = ggplot2::element_text(
            angle = 90,
            hjust = 1,
            vjust = 0.5
          ),
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 11)
        )
    }

    print(p)
    return(invisible(p))
  }

  if (type == "heatmap") {
    hdat <- dat

    hdat$Yearv <- lubridate::year(hdat$DATE)
    hdat$DOYv <- lubridate::yday(hdat$DATE)
    hdat$value <- hdat[[stat]]

    hdat <- hdat %>%
      dplyr::filter(!is.na(value), !is.na(Yearv), !is.na(DOYv))

    if (nrow(hdat) == 0) {
      stop("No non-missing values available for the selected statistic.")
    }

    p <- ggplot2::ggplot(
      hdat,
      ggplot2::aes(x = DOYv, y = factor(Yearv), fill = value)
    ) +
      ggplot2::geom_tile() +
      ggplot2::labs(
        x = "Day of year",
        y = "Year",
        fill = stat
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      )

    print(p)
    return(invisible(p))
  }
}