R/plot.ZG_output.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Plotting output of zero-growth approach
#'
#' @description
#' S3 plot method for output from \code{phase.zg()}.
#' It supports tree-level \code{ZG_phase}/\code{ZG_cycle} outputs and older
#' vignette calls such as \code{plot(zg_out, temporal = "daily")}.
#'
#' @param x An object of class \code{"ZG_output"}, returned by \code{phase.zg()}.
#' @param DOY Numeric vector with initial and final day of year for plotting,
#'   e.g. \code{c(50, 60)}. If \code{NULL}, the full available period is plotted.
#' @param Year Numeric value indicating the calendar year for plotting.
#'   Required when \code{DOY} is supplied.
#' @param type Character. Plot type, including boxplot for cycle distributions. One of \code{"gro_twd"}, \code{"twd"},
#'   \code{"phase_ribbon"}, \code{"abr"}, \code{"phase_summary"},
#'   \code{"slope_diagnostics"}, or \code{"event_timeline"}.
#' @param temporal Character. Temporal resolution for phase-level plots.
#' @param metric Character. Event-level metric used for event plots.
#' @param phase_names Character vector with two phase labels.
#' @param cols Named character vector of colours for \code{GRO} and \code{TWD}.
#' @param point_size Numeric point size.
#' @param ribbon_alpha Numeric transparency for phase ribbons.
#' @param show_tree Logical. If \code{TRUE}, tree information is shown in the subtitle.
#' @param stat Numeric cycle column for boxplots; in ZG, an alias for metric.
#' @param phase Phase names to include in boxplots, or all.
#' @param source Table for scatter plots: ZG_cycle (events) or ZG_drought (days).
#' @param xvar,yvar Names of numeric columns in source for scatter axes.
#' @param fit Scatter fit: none or linear, with a 95 percent confidence band.
#' @param ... Unsupported arguments raise an error.
#'
#' @details For type = boxplot, each cycle/event contributes one observation.
#' Calendar groups use event start dates, retain the year, and do not aggregate
#' values before calculating quartiles. Non-finite measurements are excluded.
#' @return A \code{ggplot} object. Scatter plots carry \code{scatter_stats} and
#'   \code{scatter_model} attributes. Linear fit statistics assume independent
#'   observations and are pooled across the selected rows. Non-finite pairs are excluded.
#'
#' @examples
#' library(dendRoAnalyst)
#' data(gf_nepa17)
#' zg.phase <- phase.zg(df = gf_nepa17, TreeNum = 1)
#' plot(zg.phase, DOY = c(50, 51), Year = 2017)
#' plot(zg.phase, temporal = "daily")
#'
#' @method plot ZG_output
#' @export
plot.ZG_output <- function(x,
                           DOY = NULL,
                           Year = NULL,
                           type = c("gro_twd", "twd", "phase_ribbon",
                                    "abr", "phase_summary",
                                    "slope_diagnostics", "event_timeline", "boxplot", "scatter"),
                           temporal = c("raw", "subdaily", "daily", "weekly",
                                        "monthly", "yearly", "event", "cycle"),
                           metric = "ABr.value",
                           phase_names = c("TWD", "GRO"),
                           cols = c("GRO" = "#2c7fb8", "TWD" = "#d95f0e"),
                           point_size = 2.5,
                           ribbon_alpha = 0.35,
                           show_tree = TRUE,
                           stat = NULL,
                           phase = "all",
                           source = c("ZG_cycle", "ZG_drought"),
                           xvar = NULL, yvar = NULL,
                           fit = c("none", "linear"),
                           ...) {

  ZG_output <- x

  if (!inherits(ZG_output, "ZG_output")) {
    stop("The input list's class must be 'ZG_output'.", call. = FALSE)
  }

  type <- match.arg(type)
  dots <- list(...)
  if(length(dots)) stop("Unsupported plot arguments: ",paste(names(dots),collapse=", "),call.=FALSE)
  if (type == "scatter") return(.dm_zg_scatter(x, match.arg(source), xvar, yvar,
                                              match.arg(fit), Year, DOY, point_size))
  if(!is.null(stat)) {
    if(!missing(metric) && !identical(stat,metric)) stop("Supply either stat or metric, or matching values.",call.=FALSE)
    metric <- stat
  }
  if(type == "boxplot") return(.dm_cycle_distribution(x,"ZG",metric,temporal,phase,Year,DOY))
  if(!identical(phase,"all")) stop("phase filtering is supported for boxplots only.",call.=FALSE)


  if (missing(temporal) || is.null(temporal)) {
    temporal <- "raw"
  } else {
    temporal <- as.character(temporal)[1]
    allowed_temporal <- c("raw", "subdaily", "daily", "weekly",
                          "monthly", "yearly", "event", "cycle")
    if (!temporal %in% allowed_temporal) {
      temporal <- "raw"
    }
  }

  if (!is.null(DOY) && is.null(Year)) {
    stop("'Year' must be supplied when 'DOY' is supplied.", call. = FALSE)
  }

  if (is.null(DOY) && !is.null(Year)) {
    stop("'DOY' must be supplied when 'Year' is supplied.", call. = FALSE)
  }

  .first_existing <- function(data, candidates) {
    hit <- candidates[candidates %in% names(data)]
    if (length(hit) == 0) NA_character_ else hit[1]
  }

  .safe_max <- function(x) {
    x <- x[is.finite(x)]
    if (length(x) == 0) NA_real_ else max(x, na.rm = TRUE)
  }

  .safe_mean <- function(x) {
    x <- x[is.finite(x)]
    if (length(x) == 0) NA_real_ else mean(x, na.rm = TRUE)
  }

  .modal_value <- function(x) {
    x <- x[!is.na(x)]
    if (length(x) == 0) NA else names(sort(table(x), decreasing = TRUE))[1]
  }

  zg_phase <- ZG_output$ZG_phase
  if (is.null(zg_phase)) {
    zg_phase <- ZG_output[[2]]
  }

  if (is.null(zg_phase) || !is.data.frame(zg_phase)) {
    stop("ZG_output must contain a valid 'ZG_phase' data frame.", call. = FALSE)
  }

  tree_col <- .first_existing(zg_phase, c("Tree", "tree", "TREE"))
  time_col <- .first_existing(zg_phase, c("TIME", "Time", "time", "Date", "date", "timestamp", "Timestamp"))
  dm_col <- .first_existing(zg_phase, c("dm", "DM", "Dm", "value", "Value"))
  phase_col <- .first_existing(zg_phase, c("Phases", "Phase", "phase", "PHASES", "PHASE"))
  gro_col <- .first_existing(zg_phase, c("GRO", "gro", "Gro"))
  twd_col <- .first_existing(zg_phase, c("TWD", "twd", "Twd"))

  if (is.na(time_col)) {
    time_candidates <- which(vapply(zg_phase, function(z) {
      inherits(z, "Date") || inherits(z, "POSIXct") || inherits(z, "POSIXt")
    }, logical(1)))
    if (length(time_candidates) > 0) {
      time_col <- names(zg_phase)[time_candidates[1]]
    } else {
      time_col <- if (!is.na(tree_col) && ncol(zg_phase) >= 2) names(zg_phase)[2] else names(zg_phase)[1]
    }
  }

  if (is.na(dm_col)) {
    numeric_candidates <- names(zg_phase)[vapply(zg_phase, is.numeric, logical(1))]
    numeric_candidates <- setdiff(numeric_candidates, c(tree_col, phase_col, gro_col, twd_col))
    dm_col <- if (length(numeric_candidates) > 0) numeric_candidates[1] else names(zg_phase)[2]
  }

  if (is.na(phase_col)) {
    phase_col <- if (!is.na(tree_col) && ncol(zg_phase) >= 4) names(zg_phase)[4] else names(zg_phase)[3]
  }

  td <- tibble::tibble(
    TIME = zg_phase[[time_col]],
    dm = zg_phase[[dm_col]],
    Phases = zg_phase[[phase_col]],
    GRO = if (!is.na(gro_col)) zg_phase[[gro_col]] else NA_real_,
    TWD = if (!is.na(twd_col)) zg_phase[[twd_col]] else NA_real_,
    Tree = if (!is.na(tree_col)) zg_phase[[tree_col]] else NA_character_
  )

  if (!is.null(DOY) && !is.null(Year)) {
    td <- dendro.truncate(td, CalYear = Year, DOY = DOY)
    td <- as.data.frame(td)
  }

  td <- dplyr::arrange(td, .data$TIME)

  phase_chr <- as.character(td$Phases)
  td$Phases <- dplyr::case_when(
    is.na(phase_chr) ~ NA_character_,
    phase_chr %in% c("1", "1L") ~ phase_names[1],
    phase_chr %in% c("2", "2L") ~ phase_names[2],
    phase_chr %in% phase_names ~ phase_chr,
    TRUE ~ phase_chr
  )
  td$Phases <- factor(td$Phases, levels = phase_names)

  if (temporal %in% c("daily", "weekly", "monthly", "yearly")) {
    td$Period <- switch(
      temporal,
      daily = as.Date(td$TIME),
      weekly = as.Date(cut(as.Date(td$TIME), breaks = "week")),
      monthly = as.Date(cut(as.Date(td$TIME), breaks = "month")),
      yearly = as.Date(cut(as.Date(td$TIME), breaks = "year"))
    )

    td <- td |>
      dplyr::group_by(.data$Tree, .data$Period) |>
      dplyr::summarise(
        TIME = as.POSIXct(.data$Period[1]),
        dm = .safe_mean(.data$dm),
        GRO = .safe_max(.data$GRO),
        TWD = .safe_max(.data$TWD),
        Phases = .modal_value(as.character(.data$Phases)),
        .groups = "drop"
      )

    td$Phases <- factor(td$Phases, levels = phase_names)
  }

  tree_values <- unique(stats::na.omit(as.character(td$Tree)))
  tree_subtitle <- NULL
  if (isTRUE(show_tree) && length(tree_values) > 0) {
    tree_subtitle <- paste0("Tree: ", paste(tree_values, collapse = ", "))
  }

  temporal_label <- if (temporal %in% c("raw", "subdaily")) "raw" else temporal

  if (type == "gro_twd") {
    data_long <- tidyr::pivot_longer(
      td,
      cols = c("GRO", "TWD"),
      names_to = "Variable",
      values_to = "Value"
    )

    data_long$Value <- dplyr::if_else(
      data_long$Variable == "TWD",
      data_long$Value * -1,
      data_long$Value
    )

    return(
      ggplot2::ggplot(data_long, ggplot2::aes(x = as.POSIXct(.data$TIME), y = .data$Value)) +
        ggplot2::geom_line(
          data = dplyr::filter(data_long, .data$Variable == "GRO"),
          ggplot2::aes(color = .data$Variable)
        ) +
        ggplot2::geom_area(
          data = dplyr::filter(data_long, .data$Variable == "TWD"),
          ggplot2::aes(fill = .data$Variable, group = 1),
          alpha = 0.5
        ) +
        ggplot2::scale_fill_manual(values = cols) +
        ggplot2::scale_color_manual(values = cols) +
        ggplot2::facet_wrap(~ Variable, scales = "free_y", ncol = 1) +
        ggplot2::labs(
          title = "Zero-growth output: GRO and TWD",
          subtitle = paste(stats::na.omit(c(tree_subtitle, paste0("Temporal: ", temporal_label))), collapse = " | "),
          x = "Time",
          y = "Value"
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme(
          legend.position = "top",
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 12)
        )
    )
  }

  if (type == "twd") {
    return(
      ggplot2::ggplot(td, ggplot2::aes(x = as.POSIXct(.data$TIME), y = .data$TWD)) +
        ggplot2::geom_area(fill = cols["TWD"], alpha = 0.45) +
        ggplot2::geom_line(color = cols["TWD"]) +
        ggplot2::labs(
          title = "Tree water deficit",
          subtitle = paste(stats::na.omit(c(tree_subtitle, paste0("Temporal: ", temporal_label))), collapse = " | "),
          x = "Time",
          y = "TWD"
        ) +
        ggplot2::theme_minimal()
    )
  }

  if (type == "phase_ribbon") {
    if (nrow(td) > 1) {
      td$.phase_run <- cumsum(c(TRUE, as.character(td$Phases[-1]) != as.character(td$Phases[-nrow(td)])))
    } else {
      td$.phase_run <- 1L
    }

    ribbon_data <- td |>
      dplyr::group_by(.data$.phase_run, .data$Phases) |>
      dplyr::summarise(
        xmin = min(.data$TIME, na.rm = TRUE),
        xmax = max(.data$TIME, na.rm = TRUE),
        .groups = "drop"
      ) |>
      dplyr::filter(!is.na(.data$Phases))

    return(
      ggplot2::ggplot(td, ggplot2::aes(x = as.POSIXct(.data$TIME), y = .data$dm)) +
        ggplot2::geom_rect(
          data = ribbon_data,
          ggplot2::aes(
            xmin = as.POSIXct(.data$xmin),
            xmax = as.POSIXct(.data$xmax),
            ymin = -Inf,
            ymax = Inf,
            fill = .data$Phases
          ),
          inherit.aes = FALSE,
          alpha = ribbon_alpha
        ) +
        ggplot2::geom_line(color = "grey35") +
        ggplot2::scale_fill_manual(values = stats::setNames(unname(cols[c("TWD", "GRO")]), phase_names)) +
        ggplot2::labs(title = "Zero-growth phases", subtitle = tree_subtitle, x = "Time", y = "Stem size variation") +
        ggplot2::theme_minimal()
    )
  }

  zg_cycle <- ZG_output$ZG_cycle
  if (is.null(zg_cycle)) {
    zg_cycle <- ZG_output[[1]]
  }

  if (is.null(zg_cycle) || !is.data.frame(zg_cycle)) {
    stop("ZG_output must contain a valid 'ZG_cycle' data frame for this plot type.", call. = FALSE)
  }

  cycle <- as.data.frame(zg_cycle)
  cycle <- .dm_subset_window(cycle, time_col = "Start", Year = Year, DOY = DOY)
  if(type != "phase_summary" && !metric %in% names(cycle)) stop("Unknown event metric: ",metric,call.=FALSE)
  if(type != "phase_summary") cycle$metric_plot <- cycle[[metric]]

  if (type == "abr") {
    return(
      ggplot2::ggplot(cycle, ggplot2::aes(x = as.POSIXct(.data$Start), y = .data$metric_plot)) +
        ggplot2::geom_col(alpha = 0.8) +
        ggplot2::labs(title = paste0("Event-level ", metric), x = "Event start", y = metric) +
        ggplot2::theme_minimal()
    )
  }

  if (type == "phase_summary") {
    phase_col <- .first_existing(cycle, c("Phases", "Phase", "phase"))
    if (is.na(phase_col) || !"Duration_h" %in% names(cycle)) {
      stop("'Phases' and 'Duration_h' are required for type = 'phase_summary'.", call. = FALSE)
    }

    cycle$Phases_plot <- as.character(cycle[[phase_col]])
    cycle$Phases_plot <- dplyr::case_when(
      cycle$Phases_plot %in% c("1", "1L") ~ phase_names[1],
      cycle$Phases_plot %in% c("2", "2L") ~ phase_names[2],
      TRUE ~ cycle$Phases_plot
    )

    phase_summary <- cycle |>
      dplyr::group_by(.data$Phases_plot) |>
      dplyr::summarise(Duration_h = sum(.data$Duration_h, na.rm = TRUE), .groups = "drop")

    return(
      ggplot2::ggplot(phase_summary, ggplot2::aes(x = .data$Phases_plot, y = .data$Duration_h, fill = .data$Phases_plot)) +
        ggplot2::geom_col(alpha = 0.85) +
        ggplot2::labs(title = "Zero-growth phase duration", x = "Phase", y = "Total duration (h)") +
        ggplot2::theme_minimal()
    )
  }

  if (type == "slope_diagnostics") {
    if (!all(c("slope.load.twd", "recovery.slowness.std") %in% names(cycle))) {
      stop("Columns 'slope.load.twd' and 'recovery.slowness.std' are required for slope diagnostics.", call. = FALSE)
    }

    return(
      ggplot2::ggplot(cycle, ggplot2::aes(x = .data$`slope.load.twd`, y = .data$recovery.slowness.std, color = .data$metric_plot)) +
        ggplot2::geom_point(size = point_size, alpha = 0.75) +
        ggplot2::labs(title = "ABr slope diagnostics", x = "Loading slope of TWD", y = "Recovery slowness") +
        ggplot2::theme_minimal()
    )
  }

  if (type == "event_timeline") {
    return(
      ggplot2::ggplot(cycle) +
        ggplot2::geom_segment(
          ggplot2::aes(
            x = as.POSIXct(.data$Start),
            xend = as.POSIXct(.data$End),
            y = .data$metric_plot,
            yend = .data$metric_plot
          ),
          linewidth = 1.2,
          alpha = 0.8
        ) +
        ggplot2::labs(title = paste0("Event timeline: ", metric), x = "Time", y = metric) +
        ggplot2::theme_minimal()
    )
  }
}