R/plot_climate_phase.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

# =========================================================
# Internal helpers for climate attachment and plotting
# =========================================================

.dm_subset_window <- function(dat, time_col, Year = NULL, DOY = NULL) {

  if (!time_col %in% names(dat)) {
    stop(sprintf("Time column '%s' was not found.", time_col))
  }

  tt <- dat[[time_col]]

  if (!inherits(tt, "Date") && !inherits(tt, "POSIXct")) {
    tt <- as.POSIXct(tt)
  }

  keep <- rep(TRUE, length(tt))

  if (!is.null(Year)) {
    keep <- keep & lubridate::year(tt) %in% Year
  }

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

    doy_min <- min(DOY, na.rm = TRUE)
    doy_max <- max(DOY, na.rm = TRUE)

    keep <- keep &
      lubridate::yday(tt) >= doy_min &
      lubridate::yday(tt) <= doy_max
  }

  out <- dat[keep, , drop = FALSE]

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

  out
}


.dm_label_phase <- function(ph, obj_class = c("ZG", "SC")) {
  obj_class <- match.arg(obj_class)

  if (obj_class == "ZG") {
    factor(ph, levels = c(1, 2), labels = c("TWD", "GRO"))
  } else {
    factor(ph, levels = c(1, 2, 3),
           labels = c("Shrinkage", "Expansion", "Increment"))
  }
}


.dm_detect_cycle_climate_cols <- function(dat, kind = c("ZG", "SC")) {
  kind <- match.arg(kind)

  core <- if (kind == "ZG") {
    c(
      "Phases", "Start", "End", "Duration_h", "Magnitude", "rate",
      "max.twd", "Max.twd.time", "ABr.value", "Avg.twd", "STD.twd",
      "AUC_load", "AUC_total", "DOY"
    )
  } else {
    c(
      "Phases", "Start", "End", "Duration_h", "Duration_m",
      "Magnitude", "rate", "DOY"
    )
  }

  setdiff(names(dat), core)
}


.dm_detect_point_climate_cols <- function(dat, kind = c("ZG", "SC")) {
  kind <- match.arg(kind)

  core <- if (kind == "ZG") {
    c("TIME", "dm", "Phases", "TWD", "GRO")
  } else {
    c("TIME", "dm", "Phases")
  }

  setdiff(names(dat), core)
}


# =========================================================
# Plot one climate variable
# =========================================================

#' Plot climate attached to dendrometer outputs
#'
#' @description
#' General plotting function for dendrometer outputs with attached climate
#' information. It works with objects of class \code{daily_output_clim},
#' \code{ZG_output_clim}, and \code{SC_output_clim}, usually produced by
#' \code{dm_add_climate()}, \code{dm_join_daily_clim()},
#' \code{dm_join_phase_clim()}, or \code{dm_join_subdaily_clim()}.
#'
#' The function can display one selected climate variable through time,
#' by phase or day status, as boxplots or violin plots, or as a
#' climate-response relationship.
#'
#' @param x A climate-augmented dendrometer object.
#' @param climate_var Character. Name of the climate variable to plot.
#' @param metric_var Optional character. Biological response variable used when
#'   \code{type = "relation"}. If \code{NULL}, a suitable default is selected
#'   when possible.
#' @param scale Character. Data level to use. One of \code{"auto"},
#'   \code{"daily"}, \code{"cycle"}, or \code{"point"}.
#' @param type Character. Plot type. One of \code{"timeseries"},
#'   \code{"boxplot"}, \code{"violin"}, \code{"both"}, or \code{"relation"}.
#' @param group_var Optional character. Grouping variable. If \code{NULL},
#'   \code{Day_status} is used for daily outputs when available, and
#'   \code{Phases} is used for cycle- or point-level phase outputs.
#' @param Year Optional numeric year or vector of years for subsetting.
#' @param DOY Optional numeric vector of length 2 giving the day-of-year range.
#' @param point_alpha Numeric. Alpha value for points.
#' @param add_smooth Logical. If \code{TRUE} and \code{type = "relation"},
#'   a linear trend line is added.
#' @param temporal Deprecated argument kept for compatibility with older
#'   examples. It is ignored.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @examples
#' \donttest{
#' data(gf_nepa17)
#' data(ktm_clim_hourly)
#'
#' zg <- phase.zg(df = gf_nepa17[1:800, ], TreeNum = 1)
#' zg_clim <- dm_add_climate(
#'   zg,
#'   ktm_clim_hourly,
#'   scale = "phase",
#'   mean_vars = c("temp", "VPD", "RH"),
#'   max_vars  = c("temp", "VPD"),
#'   sum_vars  = c("prec")
#' )
#'
#' dm_plot_climate(
#'   zg_clim,
#'   climate_var = "VPD_mean_phase",
#'   scale = "cycle",
#'   type = "boxplot"
#' )
#'
#' plot(
#'   zg_clim,
#'   climate_var = "VPD_mean_phase",
#'   scale = "cycle",
#'   type = "boxplot"
#' )
#' }
#'
#' @importFrom ggplot2 ggplot aes geom_line geom_point geom_boxplot geom_violin
#'   geom_smooth theme_minimal theme element_text labs position_jitter
#' @importFrom rlang .data
#' @importFrom tibble as_tibble
#' @importFrom lubridate year yday
#'
#' @export
dm_plot_climate <- function(x,
                            climate_var,
                            metric_var = NULL,
                            scale = c("auto", "daily", "cycle", "point"),
                            type = c("timeseries", "boxplot", "violin", "both", "relation"),
                            group_var = NULL,
                            Year = NULL,
                            DOY = NULL,
                            point_alpha = 0.6,
                            add_smooth = FALSE,
                            temporal = NULL) {

  scale <- match.arg(scale)
  type <- match.arg(type)

  dat <- NULL
  time_col <- NULL
  default_group <- NULL
  obj_kind <- NULL

  # -----------------------------------------------------
  # Choose data source
  # -----------------------------------------------------
  if (inherits(x, "daily_output_clim")) {

    dat <- tibble::as_tibble(x)
    time_col <- "DATE"
    obj_kind <- "daily"
    default_group <- if ("Day_status" %in% names(dat)) "Day_status" else NULL

  } else if (inherits(x, "ZG_output_clim") || inherits(x, "SC_output_clim")) {

    kind <- if (inherits(x, "ZG_output_clim") || inherits(x, "ZG_output")) {
      "ZG"
    } else {
      "SC"
    }

    cyc_name <- if (kind == "ZG") "ZG_cycle" else "SC_cycle"
    ph_name  <- if (kind == "ZG") "ZG_phase" else "SC_phase"

    cyc_dat <- tibble::as_tibble(x[[cyc_name]])
    ph_dat  <- tibble::as_tibble(x[[ph_name]])

    cyc_clim <- .dm_detect_cycle_climate_cols(cyc_dat, kind = kind)
    ph_clim  <- .dm_detect_point_climate_cols(ph_dat, kind = kind)

    if (scale == "auto") {
      if (climate_var %in% cyc_clim) {
        scale <- "cycle"
      } else if (climate_var %in% ph_clim) {
        scale <- "point"
      } else {
        stop(
          sprintf(
            "Climate variable '%s' was not found in cycle-level or point-level data.",
            climate_var
          )
        )
      }
    }

    if (scale == "daily") {
      stop("For ZG/SC outputs, use scale = 'cycle' or scale = 'point'.")
    }

    if (scale == "cycle") {
      dat <- cyc_dat
      time_col <- "Start"
    } else {
      dat <- ph_dat
      time_col <- "TIME"
    }

    obj_kind <- kind
    default_group <- "Phases"

    if ("Phases" %in% names(dat)) {
      dat$Phases <- .dm_label_phase(dat$Phases, obj_class = kind)
    }

  } else {
    stop(
      "'x' must be a climate-augmented object returned by dm_add_climate(), ",
      "dm_join_daily_clim(), dm_join_phase_clim(), or dm_join_subdaily_clim()."
    )
  }

  if (!climate_var %in% names(dat)) {
    stop(sprintf("Climate variable '%s' was not found in the selected data.", climate_var))
  }

  if (is.null(group_var)) {
    group_var <- default_group
  }

  if (!is.null(group_var) && !group_var %in% names(dat)) {
    stop(sprintf("Grouping variable '%s' was not found in the selected data.", group_var))
  }

  dat <- .dm_subset_window(dat, time_col = time_col, Year = Year, DOY = DOY)

  # -----------------------------------------------------
  # Time series
  # -----------------------------------------------------
  if (type == "timeseries") {

    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = .data[[time_col]], y = .data[[climate_var]])
    )

    if (!is.null(group_var)) {
      p <- p +
        ggplot2::geom_line(
          ggplot2::aes(color = .data[[group_var]]),
          linewidth = 0.5
        ) +
        ggplot2::geom_point(
          ggplot2::aes(color = .data[[group_var]]),
          size = 1.5,
          alpha = point_alpha
        ) +
        ggplot2::labs(color = group_var)
    } else {
      p <- p +
        ggplot2::geom_line(linewidth = 0.5) +
        ggplot2::geom_point(size = 1.5, alpha = point_alpha)
    }

    p <- p +
      ggplot2::labs(x = time_col, y = climate_var) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text  = ggplot2::element_text(size = 12)
      )

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

  # -----------------------------------------------------
  # Boxplot / violin / both
  # -----------------------------------------------------
  if (type %in% c("boxplot", "violin", "both")) {

    if (is.null(group_var)) {
      stop("A grouping variable is required for boxplot/violin style plots.")
    }

    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = .data[[group_var]], y = .data[[climate_var]])
    )

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

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

    p <- p +
      ggplot2::geom_point(
        position = ggplot2::position_jitter(width = 0.12, height = 0),
        size = 1.2,
        alpha = point_alpha,
        color = "grey25"
      ) +
      ggplot2::labs(x = group_var, y = climate_var) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text  = ggplot2::element_text(size = 12)
      )

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

  # -----------------------------------------------------
  # Climate-response relationship
  # -----------------------------------------------------
  if (type == "relation") {

    if (is.null(metric_var)) {
      if (obj_kind == "daily" && "amplitude" %in% names(dat)) {
        metric_var <- "amplitude"
      } else if (obj_kind == "ZG") {
        if ("ABr.value" %in% names(dat) && any(!is.na(dat$ABr.value))) {
          metric_var <- "ABr.value"
        } else if ("Magnitude" %in% names(dat)) {
          metric_var <- "Magnitude"
        } else {
          metric_var <- names(dat)[2]
        }
      } else if (obj_kind == "SC") {
        if ("Magnitude" %in% names(dat)) {
          metric_var <- "Magnitude"
        } else {
          metric_var <- names(dat)[2]
        }
      } else {
        stop("Please provide 'metric_var' for type = 'relation'.")
      }
    }

    if (!metric_var %in% names(dat)) {
      stop(sprintf("Metric variable '%s' was not found in the selected data.", metric_var))
    }

    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = .data[[climate_var]], y = .data[[metric_var]])
    )

    if (!is.null(group_var)) {
      p <- p +
        ggplot2::geom_point(
          ggplot2::aes(color = .data[[group_var]]),
          alpha = point_alpha
        )
    } else {
      p <- p +
        ggplot2::geom_point(alpha = point_alpha)
    }

    if (isTRUE(add_smooth)) {
      p <- p +
        ggplot2::geom_smooth(method = "lm", se = TRUE, color = "black")
    }

    p <- p +
      ggplot2::labs(x = climate_var, y = metric_var, color = group_var) +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text  = ggplot2::element_text(size = 12)
      )

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


# =========================================================
# Compare multiple climate variables
# =========================================================

#' Compare multiple climate variables attached to dendrometer outputs
#'
#' @description
#' Plots and compares multiple climate variables attached to climate-augmented
#' dendrometer outputs. It works with objects of class \code{daily_output_clim},
#' \code{ZG_output_clim}, and \code{SC_output_clim}.
#'
#' Supported plot types include faceted time-series, grouped boxplots,
#' violin plots, combined violin-boxplots, correlation heatmaps, and
#' regression heatmaps.
#'
#' @param x A climate-augmented dendrometer object.
#' @param climate_vars Character vector of climate variables to compare. If
#'   \code{NULL}, all detected climate variables at the selected scale are used.
#' @param numeric_vars Optional character vector of numeric variables used when
#'   \code{type = "cor_heatmap"} or \code{type = "regression_heatmap"}. If
#'   \code{NULL}, all numeric variables in the selected data are used.
#' @param scale Character. Data level to use. One of \code{"auto"},
#'   \code{"daily"}, \code{"cycle"}, or \code{"point"}.
#' @param type Character. Plot type. One of \code{"timeseries"},
#'   \code{"boxplot"}, \code{"violin"}, \code{"both"},
#'   \code{"cor_heatmap"}, or \code{"regression_heatmap"}.
#' @param group_var Optional character. Grouping variable. If \code{NULL},
#'   \code{Day_status} is used for daily outputs when available, and
#'   \code{Phases} is used for phase outputs.
#' @param Year Optional numeric year or vector of years for subsetting.
#' @param DOY Optional numeric vector of length 2 giving the day-of-year range.
#' @param box_scales Character. Facet scale behavior for grouped plots. One of
#'   \code{"free_y"} or \code{"fixed"}.
#' @param cor_method Character. Correlation method. One of \code{"pearson"},
#'   \code{"spearman"}, or \code{"kendall"}.
#' @param use_pairwise Logical. If \code{TRUE}, pairwise complete observations
#'   are used in correlations.
#' @param regression_stat Character. Statistic shown in regression heatmaps.
#'   One of \code{"r.squared"}, \code{"slope"}, or \code{"p.value"}.
#' @param point_alpha Numeric. Alpha value for points.
#' @param show_values Logical. If \code{TRUE}, numeric values are printed inside
#'   heatmap tiles.
#' @param temporal Deprecated argument kept for compatibility with older
#'   examples. It is ignored.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @details
#' For \code{type = "regression_heatmap"}, pairwise simple linear regressions
#' are fitted as \code{y ~ x} for every variable combination. The heatmap can
#' display \eqn{R^2}, slope, or p-value.
#'
#' @examples
#' \donttest{
#' data(gf_nepa17)
#' data(ktm_clim_hourly)
#'
#' zg <- phase.zg(df = gf_nepa17[1:800, ], TreeNum = 1)
#' zg_clim <- dm_add_climate(
#'   zg,
#'   ktm_clim_hourly,
#'   scale = "phase",
#'   mean_vars = c("temp", "VPD", "RH"),
#'   max_vars  = c("temp", "VPD"),
#'   sum_vars  = c("prec")
#' )
#'
#' dm_plot_climate_compare(
#'   zg_clim,
#'   climate_vars = c("temp_mean_phase", "VPD_mean_phase", "RH_mean_phase"),
#'   scale = "cycle",
#'   type = "boxplot"
#' )
#'
#' plot(
#'   zg_clim,
#'   climate_vars = c("temp_mean_phase", "VPD_mean_phase", "RH_mean_phase"),
#'   scale = "cycle",
#'   type = "boxplot"
#' )
#' }
#'
#' @importFrom ggplot2 ggplot aes geom_line geom_point geom_boxplot geom_violin
#'   geom_tile geom_text facet_wrap theme_minimal theme element_text labs
#'   position_jitter scale_fill_gradient2 scale_fill_gradient
#' @importFrom tidyr pivot_longer
#' @importFrom tibble as_tibble
#' @importFrom rlang .data
#' @importFrom stats cor lm var coef
#'
#' @export
dm_plot_climate_compare <- function(x,
                                    climate_vars = NULL,
                                    numeric_vars = NULL,
                                    scale = c("auto", "daily", "cycle", "point"),
                                    type = c("timeseries", "boxplot", "violin", "both",
                                             "cor_heatmap", "regression_heatmap"),
                                    group_var = NULL,
                                    Year = NULL,
                                    DOY = NULL,
                                    box_scales = c("free_y", "fixed"),
                                    cor_method = c("pearson", "spearman", "kendall"),
                                    use_pairwise = TRUE,
                                    regression_stat = c("r.squared", "slope", "p.value"),
                                    point_alpha = 0.5,
                                    show_values = FALSE,
                                    temporal = NULL) {

  scale <- match.arg(scale)
  type <- match.arg(type)
  box_scales <- match.arg(box_scales)
  cor_method <- match.arg(cor_method)
  regression_stat <- match.arg(regression_stat)

  dat <- NULL
  time_col <- NULL
  default_group <- NULL
  clim_candidates <- NULL

  value <- NULL
  Var1 <- NULL
  Var2 <- NULL

  # -----------------------------------------------------
  # Choose data source
  # -----------------------------------------------------
  if (inherits(x, "daily_output_clim")) {

    dat <- tibble::as_tibble(x)
    time_col <- "DATE"
    default_group <- if ("Day_status" %in% names(dat)) "Day_status" else NULL

    clim_candidates <- setdiff(
      names(dat),
      c(
        "DATE", "Min", "Time_min", "Max", "Time_max", "mean", "median",
        "amplitude", "Time_min_h", "Time_max_h", "lag_h", "Remarks",
        "Max_diff", "Day_status"
      )
    )

  } else if (inherits(x, "ZG_output_clim") || inherits(x, "SC_output_clim")) {

    kind <- if (inherits(x, "ZG_output_clim") || inherits(x, "ZG_output")) {
      "ZG"
    } else {
      "SC"
    }

    cyc_name <- if (kind == "ZG") "ZG_cycle" else "SC_cycle"
    ph_name  <- if (kind == "ZG") "ZG_phase" else "SC_phase"

    cyc_dat <- tibble::as_tibble(x[[cyc_name]])
    ph_dat  <- tibble::as_tibble(x[[ph_name]])

    cyc_clim <- .dm_detect_cycle_climate_cols(cyc_dat, kind = kind)
    ph_clim  <- .dm_detect_point_climate_cols(ph_dat, kind = kind)

    if (scale == "auto") {
      if (!is.null(climate_vars)) {
        if (all(climate_vars %in% cyc_clim)) {
          scale <- "cycle"
        } else if (all(climate_vars %in% ph_clim)) {
          scale <- "point"
        } else {
          stop(
            "Selected climate_vars were not found consistently in ",
            "cycle-level or point-level data."
          )
        }
      } else if (type %in% c("cor_heatmap", "regression_heatmap")) {
        scale <- "cycle"
      } else if (length(cyc_clim) > 0) {
        scale <- "cycle"
      } else {
        scale <- "point"
      }
    }

    if (scale == "daily") {
      stop("For ZG/SC outputs, use scale = 'cycle' or scale = 'point'.")
    }

    if (scale == "cycle") {
      dat <- cyc_dat
      time_col <- "Start"
      clim_candidates <- cyc_clim
    } else {
      dat <- ph_dat
      time_col <- "TIME"
      clim_candidates <- ph_clim
    }

    default_group <- "Phases"

    if ("Phases" %in% names(dat)) {
      dat$Phases <- .dm_label_phase(dat$Phases, obj_class = kind)
    }

  } else {
    stop(
      "'x' must be a climate-augmented object returned by dm_add_climate(), ",
      "dm_join_daily_clim(), dm_join_phase_clim(), or dm_join_subdaily_clim()."
    )
  }

  if (is.null(group_var)) {
    group_var <- default_group
  }

  if (!is.null(group_var) && !group_var %in% names(dat)) {
    stop(sprintf("Grouping variable '%s' was not found in the selected data.", group_var))
  }

  dat <- .dm_subset_window(dat, time_col = time_col, Year = Year, DOY = DOY)

  # -----------------------------------------------------
  # Select climate variables for faceted plots
  # -----------------------------------------------------
  if (type %in% c("timeseries", "boxplot", "violin", "both")) {

    if (is.null(climate_vars)) {
      climate_vars <- clim_candidates
    }

    if (length(climate_vars) == 0) {
      stop("No climate variables were detected for plotting.")
    }

    miss <- setdiff(climate_vars, names(dat))
    if (length(miss) > 0) {
      stop(sprintf(
        "The following climate_vars were not found: %s",
        paste(miss, collapse = ", ")
      ))
    }

    pdat <- tidyr::pivot_longer(
      dat,
      cols = climate_vars,
      names_to = "Climate_var",
      values_to = "value"
    )

    pdat <- pdat[!is.na(pdat$value), , drop = FALSE]
  }

  # -----------------------------------------------------
  # Time series
  # -----------------------------------------------------
  if (type == "timeseries") {

    p <- ggplot2::ggplot(
      pdat,
      ggplot2::aes(x = .data[[time_col]], y = value)
    )

    if (!is.null(group_var)) {
      p <- p +
        ggplot2::geom_line(
          ggplot2::aes(color = .data[[group_var]]),
          linewidth = 0.45
        ) +
        ggplot2::geom_point(
          ggplot2::aes(color = .data[[group_var]]),
          size = 1.2,
          alpha = point_alpha
        ) +
        ggplot2::labs(color = group_var)
    } else {
      p <- p +
        ggplot2::geom_line(linewidth = 0.45) +
        ggplot2::geom_point(size = 1.2, alpha = point_alpha)
    }

    p <- p +
      ggplot2::facet_wrap(~ Climate_var, scales = "free_y", ncol = 1) +
      ggplot2::labs(x = time_col, y = "Value") +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text  = ggplot2::element_text(size = 11)
      )

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

  # -----------------------------------------------------
  # Boxplot / violin / both
  # -----------------------------------------------------
  if (type %in% c("boxplot", "violin", "both")) {

    if (is.null(group_var)) {
      stop("A grouping variable is required for grouped comparison plots.")
    }

    p <- ggplot2::ggplot(
      pdat,
      ggplot2::aes(x = .data[[group_var]], y = value)
    )

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

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

    p <- p +
      ggplot2::geom_point(
        position = ggplot2::position_jitter(width = 0.12, height = 0),
        size = 1.1,
        alpha = point_alpha,
        color = "grey25"
      ) +
      ggplot2::facet_wrap(~ Climate_var, scales = box_scales) +
      ggplot2::labs(x = group_var, y = "Value") +
      ggplot2::theme_minimal() +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text  = ggplot2::element_text(size = 11)
      )

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

  # -----------------------------------------------------
  # Heatmap inputs
  # -----------------------------------------------------
  if (type %in% c("cor_heatmap", "regression_heatmap")) {

    if (is.null(numeric_vars)) {
      numeric_vars <- names(dat)[vapply(dat, is.numeric, logical(1))]
    }

    miss <- setdiff(numeric_vars, names(dat))
    if (length(miss) > 0) {
      stop(sprintf(
        "The following numeric_vars were not found: %s",
        paste(miss, collapse = ", ")
      ))
    }

    numeric_vars <- numeric_vars[vapply(dat[numeric_vars], is.numeric, logical(1))]

    if (length(numeric_vars) < 2) {
      stop("At least two numeric variables are required for heatmap plotting.")
    }

    mdat <- dat[, numeric_vars, drop = FALSE]
  }

  # -----------------------------------------------------
  # Correlation heatmap
  # -----------------------------------------------------
  if (type == "cor_heatmap") {

    cm <- stats::cor(
      mdat,
      use = if (isTRUE(use_pairwise)) "pairwise.complete.obs" else "complete.obs",
      method = cor_method
    )

    cdat <- as.data.frame(as.table(cm))
    names(cdat) <- c("Var1", "Var2", "value")

    p <- ggplot2::ggplot(
      cdat,
      ggplot2::aes(x = Var1, y = Var2, fill = value)
    ) +
      ggplot2::geom_tile() +
      ggplot2::scale_fill_gradient2(
        low = "firebrick",
        mid = "white",
        high = "steelblue",
        midpoint = 0,
        limits = c(-1, 1)
      ) +
      ggplot2::labs(
        x = NULL,
        y = NULL,
        fill = paste0("Correlation\n(", cor_method, ")")
      ) +
      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)
      )

    if (isTRUE(show_values)) {
      p <- p +
        ggplot2::geom_text(
          ggplot2::aes(label = sprintf("%.2f", value)),
          size = 3
        )
    }

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

  # -----------------------------------------------------
  # Regression heatmap
  # -----------------------------------------------------
  if (type == "regression_heatmap") {

    pair_reg_stat <- function(xv, yv, what = c("r.squared", "slope", "p.value")) {

      what <- match.arg(what)
      ok <- is.finite(xv) & is.finite(yv)

      if (sum(ok) < 3) {
        return(NA_real_)
      }

      if (stats::var(xv[ok]) == 0 || stats::var(yv[ok]) == 0) {
        return(NA_real_)
      }

      fit <- try(stats::lm(yv[ok] ~ xv[ok]), silent = TRUE)

      if (inherits(fit, "try-error")) {
        return(NA_real_)
      }

      sm <- summary(fit)

      if (what == "r.squared") {
        return(unname(sm$r.squared))
      }

      if (what == "slope") {
        return(unname(stats::coef(fit)[2]))
      }

      if (what == "p.value") {
        return(unname(sm$coefficients[2, 4]))
      }

      NA_real_
    }

    res <- expand.grid(
      Var1 = numeric_vars,
      Var2 = numeric_vars,
      stringsAsFactors = FALSE
    )

    res$value <- mapply(
      FUN = function(v1, v2) {
        pair_reg_stat(mdat[[v1]], mdat[[v2]], what = regression_stat)
      },
      v1 = res$Var1,
      v2 = res$Var2
    )

    p <- ggplot2::ggplot(
      res,
      ggplot2::aes(x = Var1, y = Var2, fill = value)
    ) +
      ggplot2::geom_tile()

    if (regression_stat == "r.squared") {
      p <- p +
        ggplot2::scale_fill_gradient(
          low = "white",
          high = "steelblue",
          limits = c(0, 1)
        )
    } else if (regression_stat == "p.value") {
      p <- p +
        ggplot2::scale_fill_gradient(
          low = "firebrick",
          high = "white"
        )
    } else {
      p <- p +
        ggplot2::scale_fill_gradient2(
          low = "firebrick",
          mid = "white",
          high = "steelblue",
          midpoint = 0
        )
    }

    p <- p +
      ggplot2::labs(x = NULL, y = NULL, fill = regression_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)
      )

    if (isTRUE(show_values)) {
      lab_fmt <- if (regression_stat == "p.value") "%.3g" else "%.2f"
      p <- p +
        ggplot2::geom_text(
          ggplot2::aes(label = sprintf(lab_fmt, value)),
          size = 3
        )
    }

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


# =========================================================
# S3 plot methods
# =========================================================

.plot_climate_s3 <- function(x,
                             y = NULL,
                             ...,
                             climate_var = NULL,
                             climate_vars = NULL,
                             numeric_vars = NULL,
                             compare = FALSE,
                             temporal = NULL) {

  if (!is.null(y)) {
    if (isTRUE(compare)) {
      if (is.null(climate_vars)) {
        climate_vars <- y
      }
    } else {
      if (is.null(climate_var)) {
        climate_var <- y
      }
    }
  }

  if (!isTRUE(compare) && (!is.null(climate_vars) || !is.null(numeric_vars))) {
    compare <- TRUE
  }

  if (isTRUE(compare)) {
    return(
      dm_plot_climate_compare(
        x = x,
        climate_vars = climate_vars,
        numeric_vars = numeric_vars,
        ...
      )
    )
  }

  if (is.null(climate_var)) {
    stop(
      "Please provide 'climate_var', for example: ",
      "plot(x, climate_var = 'VPD_mean_phase'). ",
      "For multiple variables, use compare = TRUE or provide 'climate_vars'."
    )
  }

  dm_plot_climate(
    x = x,
    climate_var = climate_var,
    ...
  )
}


#' Plot climate-augmented daily dendrometer output
#'
#' @description
#' S3 plotting method for objects of class \code{daily_output_clim}. This allows
#' climate-augmented daily dendrometer outputs to be visualized directly with
#' the generic \code{plot()} function.
#'
#' @param x An object of class \code{daily_output_clim}.
#' @param y Optional climate variable name passed as the second argument.
#' @param ... Additional arguments passed to \code{dm_plot_climate()} or
#'   \code{dm_plot_climate_compare()}.
#' @param climate_var Character. Name of one climate variable to plot.
#' @param climate_vars Character vector of climate variables for comparison
#'   plots.
#' @param numeric_vars Character vector of numeric variables used for correlation
#'   or regression heatmaps.
#' @param compare Logical. If \code{TRUE}, \code{dm_plot_climate_compare()} is
#'   called.
#' @param temporal Deprecated argument kept for compatibility with older
#'   examples. It is ignored.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @examples
#' \donttest{
#' # plot(daily_clim_output, climate_var = "temp_mean", type = "timeseries")
#' }
#'
#' @method plot daily_output_clim
#' @export
plot.daily_output_clim <- function(x,
                                   y = NULL,
                                   ...,
                                   climate_var = NULL,
                                   climate_vars = NULL,
                                   numeric_vars = NULL,
                                   compare = FALSE,
                                   temporal = NULL) {

  .plot_climate_s3(
    x = x,
    y = y,
    ...,
    climate_var = climate_var,
    climate_vars = climate_vars,
    numeric_vars = numeric_vars,
    compare = compare,
    temporal = temporal
  )
}


#' Plot climate-augmented zero-growth output
#'
#' @description
#' S3 plotting method for objects of class \code{ZG_output_clim}. This allows
#' climate-augmented zero-growth outputs to be visualized directly with the
#' generic \code{plot()} function.
#'
#' @param x An object of class \code{ZG_output_clim}.
#' @param y Optional climate variable name passed as the second argument.
#' @param ... Additional arguments passed to \code{dm_plot_climate()} or
#'   \code{dm_plot_climate_compare()}.
#' @param climate_var Character. Name of one climate variable to plot.
#' @param climate_vars Character vector of climate variables for comparison
#'   plots.
#' @param numeric_vars Character vector of numeric variables used for correlation
#'   or regression heatmaps.
#' @param compare Logical. If \code{TRUE}, \code{dm_plot_climate_compare()} is
#'   called.
#' @param temporal Deprecated argument kept for compatibility with older
#'   examples. It is ignored.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @examples
#' \donttest{
#' # plot(zg_clim, climate_var = "VPD_mean_phase", scale = "cycle",
#' #      type = "boxplot")
#' }
#'
#' @method plot ZG_output_clim
#' @export
plot.ZG_output_clim <- function(x,
                                y = NULL,
                                ...,
                                climate_var = NULL,
                                climate_vars = NULL,
                                numeric_vars = NULL,
                                compare = FALSE,
                                temporal = NULL) {

  .plot_climate_s3(
    x = x,
    y = y,
    ...,
    climate_var = climate_var,
    climate_vars = climate_vars,
    numeric_vars = numeric_vars,
    compare = compare,
    temporal = temporal
  )
}


#' Plot climate-augmented stem-cycle output
#'
#' @description
#' S3 plotting method for objects of class \code{SC_output_clim}. This allows
#' climate-augmented stem-cycle outputs to be visualized directly with the
#' generic \code{plot()} function.
#'
#' @param x An object of class \code{SC_output_clim}.
#' @param y Optional climate variable name passed as the second argument.
#' @param ... Additional arguments passed to \code{dm_plot_climate()} or
#'   \code{dm_plot_climate_compare()}.
#' @param climate_var Character. Name of one climate variable to plot.
#' @param climate_vars Character vector of climate variables for comparison
#'   plots.
#' @param numeric_vars Character vector of numeric variables used for correlation
#'   or regression heatmaps.
#' @param compare Logical. If \code{TRUE}, \code{dm_plot_climate_compare()} is
#'   called.
#' @param temporal Deprecated argument kept for compatibility with older
#'   examples. It is ignored.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @examples
#' \donttest{
#' # plot(sc_clim, climate_var = "temp_mean_phase", scale = "cycle",
#' #      type = "violin")
#' }
#'
#' @method plot SC_output_clim
#' @export
plot.SC_output_clim <- function(x,
                                y = NULL,
                                ...,
                                climate_var = NULL,
                                climate_vars = NULL,
                                numeric_vars = NULL,
                                compare = FALSE,
                                temporal = NULL) {

  .plot_climate_s3(
    x = x,
    y = y,
    ...,
    climate_var = climate_var,
    climate_vars = climate_vars,
    numeric_vars = numeric_vars,
    compare = compare,
    temporal = temporal
  )
}