R/plot.SC_output.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Plotting output of stem cycle approach
#'
#' @description
#' S3 plot method for output from \code{phase.sc()}.
#' The function supports old \code{SC_phase} tables
#' (\code{TIME}, \code{dm}, \code{Phases}) and newer tables with tree-level
#' information (\code{Tree}, \code{TIME}, \code{dm}, \code{Phases}).
#'
#' @param x An object of class \code{"SC_output"}, returned by \code{phase.sc()}.
#' @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{"points"}, \code{"ribbon"},
#'   \code{"both"}, \code{"shrinkage"}, \code{"expansion"}, or \code{"increment"}.
#' @param cols Character vector of three colours, one for each phase.
#' @param phNames Character vector with three phase names.
#' @param point_size Numeric point size used for phase points.
#' @param ribbon_alpha Numeric transparency for phase ribbons.
#' @param show_tree Logical. If \code{TRUE}, tree information is shown in the
#'   plot subtitle when available.
#' @param stat Numeric cycle column for boxplots; in ZG, an alias for metric.
#' @param phase Phase names to include in boxplots, or all.
#' @param temporal Calendar grouping for boxplots: raw, daily, weekly, monthly or yearly.
#' @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.
#'
#' @examples
#' library(dendRoAnalyst)
#' data(gf_nepa17)
#' sc.phase <- phase.sc(df = gf_nepa17, TreeNum = 1, smoothing = 12)
#' plot(sc.phase, DOY = c(50, 60), Year = 2017)
#' plot(sc.phase, DOY = c(50, 60), Year = 2017, type = "ribbon")
#' plot(sc.phase, DOY = c(50, 60), Year = 2017, type = "increment")
#'
#' @method plot SC_output
#' @export
plot.SC_output <- function(x,
                           DOY = NULL,
                           Year = NULL,
                           type = c("points", "ribbon", "both",
                                    "shrinkage", "expansion", "increment", "boxplot"),
                           cols = c("#fee8c8", "#fdbb84", "#e34a33"),
                           phNames = c("Shrinkage", "Expansion", "Increment"),
                           point_size = 3,
                           ribbon_alpha = 0.35,
                           show_tree = TRUE,
                           stat = "Magnitude",
                           temporal = "raw",
                           phase = "all",
                           ...) {

  SC_output <- x

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

  type <- match.arg(type)
  dots <- list(...)
  if(length(dots)) stop("Unsupported plot arguments: ",paste(names(dots),collapse=", "),call.=FALSE)
  if(type == "boxplot") return(.dm_cycle_distribution(x,"SC",stat,temporal,phase,Year,DOY))
  if(!identical(phase,"all") || !identical(temporal,"raw") || !missing(stat)) stop("stat, temporal and phase apply to boxplots; choose a phase-specific type for series plots.",call.=FALSE)


  if (!is.character(phNames) || length(phNames) != 3) {
    stop("'phNames' must be a character vector of length 3.", call. = FALSE)
  }

  if (!is.character(cols) || length(cols) != 3) {
    stop("'cols' must be a character vector of length 3.", call. = FALSE)
  }

  if (!is.numeric(ribbon_alpha) || length(ribbon_alpha) != 1 ||
      is.na(ribbon_alpha) || ribbon_alpha < 0 || ribbon_alpha > 1) {
    stop("'ribbon_alpha' must be a numeric value between 0 and 1.", call. = FALSE)
  }

  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]
  }

  sc_phase <- SC_output$SC_phase
  if (is.null(sc_phase)) {
    sc_phase <- SC_output[[2]]
  }

  if (is.null(sc_phase) || !is.data.frame(sc_phase)) {
    stop("SC_output must contain a valid 'SC_phase' data frame.", call. = FALSE)
  }

  tree_col <- .first_existing(sc_phase, c("Tree", "tree", "TREE"))
  time_col <- .first_existing(sc_phase, c("TIME", "Time", "time", "Date", "date", "timestamp", "Timestamp"))
  dm_col <- .first_existing(sc_phase, c("dm", "DM", "Dm", "value", "Value"))
  phase_col <- .first_existing(sc_phase, c("Phases", "Phase", "phase", "PHASES", "PHASE"))

  if (is.na(time_col)) {
    time_candidates <- which(vapply(sc_phase, function(z) {
      inherits(z, "Date") || inherits(z, "POSIXct") || inherits(z, "POSIXt")
    }, logical(1)))

    if (length(time_candidates) > 0) {
      time_col <- names(sc_phase)[time_candidates[1]]
    } else {
      time_col <- if (!is.na(tree_col) && ncol(sc_phase) >= 2) names(sc_phase)[2] else names(sc_phase)[1]
    }
  }

  if (is.na(dm_col)) {
    numeric_candidates <- names(sc_phase)[vapply(sc_phase, is.numeric, logical(1))]
    numeric_candidates <- setdiff(numeric_candidates, c(tree_col, phase_col))
    if (length(numeric_candidates) > 0) {
      dm_col <- numeric_candidates[1]
    } else {
      dm_col <- if (!is.na(tree_col) && ncol(sc_phase) >= 3) names(sc_phase)[3] else names(sc_phase)[2]
    }
  }

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

  td <- tibble::tibble(
    TIME = sc_phase[[time_col]],
    dm = sc_phase[[dm_col]],
    Phases = sc_phase[[phase_col]]
  )

  if (!is.na(tree_col)) {
    td$Tree <- sc_phase[[tree_col]]
  } else {
    sc_cycle <- SC_output$SC_cycle
    if (is.null(sc_cycle)) {
      sc_cycle <- SC_output[[1]]
    }

    if (!is.null(sc_cycle) && is.data.frame(sc_cycle) && "Tree" %in% names(sc_cycle)) {
      tree_value <- unique(stats::na.omit(sc_cycle$Tree))
      td$Tree <- if (length(tree_value) > 0) tree_value[1] else NA_character_
    } else {
      td$Tree <- NA_character_
    }
  }

  td <- td[, c("TIME", "dm", "Phases", "Tree")]

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

    if (!all(c("TIME", "dm", "Phases") %in% names(td))) {
      names(td)[seq_len(min(3, ncol(td)))] <- c("TIME", "dm", "Phases")[seq_len(min(3, ncol(td)))]
    }
    if (!"Tree" %in% names(td)) {
      td$Tree <- NA_character_
    }
  }

  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") ~ phNames[1],
    phase_chr %in% c("2", "2L") ~ phNames[2],
    phase_chr %in% c("3", "3L") ~ phNames[3],
    phase_chr %in% phNames ~ phase_chr,
    TRUE ~ phase_chr
  )

  td$Phases <- factor(td$Phases, levels = phNames)

  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 = ", "))
  }

  phase_specific_types <- c("shrinkage", "expansion", "increment")
  selected_phase <- NULL
  if (type %in% phase_specific_types) {
    selected_phase <- phNames[match(type, phase_specific_types)]
  }

  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))

  if (!is.null(selected_phase)) {
    ribbon_data <- dplyr::filter(ribbon_data, .data$Phases == selected_phase)
  }

  if (nrow(ribbon_data) > 0) {
    resolution_sec <- suppressWarnings(stats::median(diff(as.numeric(as.POSIXct(td$TIME))), na.rm = TRUE))
    if (!is.finite(resolution_sec) || is.na(resolution_sec)) {
      resolution_sec <- 0
    }
    zero_width <- as.POSIXct(ribbon_data$xmin) == as.POSIXct(ribbon_data$xmax)
    ribbon_data$xmax[zero_width] <- as.POSIXct(ribbon_data$xmax[zero_width]) + resolution_sec
  }

  point_data <- td
  if (!is.null(selected_phase)) {
    point_data <- dplyr::filter(td, .data$Phases == selected_phase)
  }

  p <- ggplot2::ggplot(td, ggplot2::aes(x = as.POSIXct(.data$TIME), y = .data$dm))

  if (type %in% c("ribbon", "both") || !is.null(selected_phase)) {
    p <- p +
      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
      )
  }

  p <- p + ggplot2::geom_line(color = "grey35")

  if (type %in% c("points", "both") || !is.null(selected_phase)) {
    p <- p +
      ggplot2::geom_point(
        data = point_data,
        ggplot2::aes(color = .data$Phases),
        size = point_size
      )
  }

  if (type %in% c("points", "both") || !is.null(selected_phase)) {
    p <- p +
      ggplot2::scale_color_manual(
        values = cols,
        na.value = "grey70",
        name = "Phases",
        limits = phNames,
        drop = FALSE
      )
  }

  if (type %in% c("ribbon", "both") || !is.null(selected_phase)) {
    p <- p +
      ggplot2::scale_fill_manual(
        values = cols,
        na.value = "grey70",
        name = "Phases",
        limits = phNames,
        drop = FALSE
      )
  }

  title_text <- if (!is.null(selected_phase)) {
    paste0("Stem-cycle phase: ", selected_phase)
  } else {
    "Stem-cycle phases"
  }

  p +
    ggplot2::theme_minimal() +
    ggplot2::labs(
      title = title_text,
      subtitle = tree_subtitle,
      x = "Time",
      y = "Stem size variation"
    ) +
    ggplot2::theme(
      axis.title = ggplot2::element_text(size = 14),
      axis.text = ggplot2::element_text(size = 12)
    )
}