Nothing
R/plot_ZG_output.R
In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data
Defines functions
plot.ZG_output
Documented in
plot.ZG_output
#' @title Plot method for zero-growth output
#'
#' @description
#' Unified S3 plotting method for objects returned by \code{phase.zg()}.
#' It supports time-series views of GRO and TWD, phase ribbons and transitions,
#' TWD-only plots, ABr summaries, phase-balance views, and boxplots of
#' phase-level statistics. For boxplots, phase summaries are first assembled to
#' the requested daily or monthly time scale, and rows with missing x-axis
#' grouping labels are removed before plotting.
#'
#' The plotting method operates on the two tables returned by
#' \code{phase.zg()}:
#' \itemize{
#' \item \code{ZG_phase}: point-level time series with \code{TIME},
#' \code{dm}, \code{Phases}, \code{TWD}, and \code{GRO};
#' \item \code{ZG_cycle}: phase-level summaries including
#' \code{max.twd}, \code{AUC.load}, \code{AUC.total}, and
#' \code{ABr.value}.
#' }
#'
#' @details
#' \strong{ABr formula shown by this plot method.}
#'
#' When \code{type = "abr"} or when \code{stat = "ABr.value"} is used in
#' boxplots, the displayed values come directly from \code{phase.zg()}:
#'
#' \deqn{
#' ABr = \max(TWD) \times \left(\frac{AUC_{load}}{\max(TWD)}\right)^{\beta}
#' }
#'
#' where
#'
#' \deqn{AUC_{load} = \int_{t_{start}}^{t_{max}} TWD(t)\,dt}
#'
#' and
#'
#' \deqn{AUC_{total} = \int_{t_{start}}^{t_{end}} TWD(t)\,dt.}
#'
#' This function does not recompute ABr or AUC values; it only visualizes
#' the values already stored in the \code{ZG_output} object.
#'
#' \strong{Aggregation behavior.}
#'
#' For \code{temporal = "daily"} or \code{"monthly"}, point-level series
#' are aggregated before plotting. TWD and GRO can be summarized using
#' \code{"mean"} or \code{"max"}, depending on \code{twd_fun} and
#' \code{gro_fun}.
#'
#' For phase-level summaries in \code{ZG_cycle}, the plot method aggregates
#' statistics by daily or monthly period. In those aggregated summaries:
#' \itemize{
#' \item \code{Duration_h}, \code{Magnitude}, \code{AUC.load},
#' \code{AUC.total}, and \code{ABr.value} are summed;
#' \item \code{max.twd} is taken as the maximum;
#' \item \code{Avg.twd} and \code{STD.twd} are averaged;
#' \item \code{rate} is recalculated from aggregated magnitude and duration.
#' }
#'
#' \strong{Phase balance.}
#'
#' For \code{type = "balance"}, two balance styles are available through
#' \code{balance_mode}. With \code{balance_mode = "time_of_day"}, the plot
#' draws TWD and GRO intervals on a 0--24 hour y-axis, preserving the order in
#' which phases occurred within each day, similar to the balance option in
#' \code{plot.SC_output()}. With \code{balance_mode = "duration"}, the plot
#' draws stacked bars showing the total number of hours spent in each phase.
#' Missing or irregular intervals can be handled with \code{balance_gap}.
#'
#' \strong{Boxplot statistics.}
#'
#' The following phase-level variables can be visualized in
#' \code{type = "boxplot"}. For \code{box_group = "period"}, boxplots are
#' drawn on a continuous date axis spanning the selected multi-annual data
#' window, so empty days or months remain visible as gaps rather than being
#' removed from the x-axis:
#' \itemize{
#' \item \code{"Duration_h"}
#' \item \code{"Magnitude"}
#' \item \code{"rate"}
#' \item \code{"max.twd"}
#' \item \code{"Avg.twd"}
#' \item \code{"STD.twd"}
#' \item \code{"AUC.load"}
#' \item \code{"AUC.total"}
#' \item \code{"ABr.value"}
#' }
#'
#' \strong{Phase handling.}
#'
#' For \code{phase = "auto"}, TWD-specific statistics
#' (\code{max.twd}, \code{Avg.twd}, \code{STD.twd}, \code{AUC.load},
#' \code{AUC.total}, \code{ABr.value}) automatically select TWD phases,
#' whereas GRO-specific variables (\code{Magnitude}, \code{rate}) select GRO
#' phases.
#'
#' @param x Object of class \code{"ZG_output"} returned by \code{phase.zg()}.
#' @param y Unused.
#' @param DOY Optional numeric vector of length 2 giving start and end
#' day-of-year for truncating the plotting window.
#' @param Year Optional numeric year used together with \code{DOY}.
#' @param view Optional backward-compatible argument. Use \code{"boxplot"}
#' to force boxplot mode.
#' @param type Plot type. One of \code{"gro_twd"}, \code{"phase_ribbon"},
#' \code{"transition"}, \code{"twd"}, \code{"abr"},
#' \code{"phase_summary"}, \code{"balance"}, or \code{"boxplot"}.
#' @param temporal Temporal scale. One of \code{"raw"}, \code{"daily"},
#' or \code{"monthly"}.
#' @param x_axis X-axis style for time-based plots. One of \code{"time"} or
#' \code{"doy"}. When \code{"doy"} is used, timestamps are converted to
#' day-of-year. For \code{temporal = "monthly"}, \code{"doy"} is not
#' meaningful and time is used instead.
#' @param stat For \code{type = "boxplot"}, one of
#' \code{"Duration_h"}, \code{"Magnitude"}, \code{"rate"},
#' \code{"max.twd"}, \code{"Avg.twd"}, \code{"STD.twd"},
#' \code{"AUC.load"}, \code{"AUC.total"}, or \code{"ABr.value"}.
#' @param phase For \code{type = "boxplot"}, one of \code{"auto"},
#' \code{"TWD"}, \code{"GRO"}, or \code{"all"}.
#' @param box_group For \code{type = "boxplot"}, grouping mode:
#' \code{"period"}, \code{"month_of_year"}, or \code{"doy"}.
#' \code{"period"} uses exact dates/months from the already aggregated
#' daily or monthly phase summaries; the other two pool values across
#' repeated seasonal positions. Rows with missing grouping labels are
#' removed before plotting.
#' @param twd_fun Aggregation function for point-level TWD when
#' \code{temporal != "raw"}. One of \code{"mean"} or \code{"max"}.
#' @param gro_fun Aggregation function for point-level GRO when
#' \code{temporal != "raw"}. One of \code{"max"} or \code{"mean"}.
#' @param transition_linetype Line type used for transition markers.
#' @param transition_alpha Transparency of transition markers.
#' @param cols Vector of two colours for TWD and GRO. It may be named
#' (\code{c(TWD = "red", GRO = "blue")}) or unnamed length 2.
#' @param singleton_as_points Logical. If \code{TRUE}, the boxplot panel
#' falls back to a point plot when every group contains only one value.
#' @param balance_mode For \code{type = "balance"}, one of
#' \code{"time_of_day"} or \code{"duration"}. \code{"time_of_day"} draws
#' continuous TWD/GRO intervals along the y-axis from 0 to 24 h, preserving
#' when phases occurred within each day. \code{"duration"} draws stacked bars
#' of total hours spent in each phase per day or month.
#' @param balance_gap For \code{type = "balance"}, one of
#' \code{"carry_forward"} or \code{"observed_only"}.
#' \code{"carry_forward"} assigns gaps to the previous observed phase;
#' \code{"observed_only"} counts only observed intervals.
#' @param ... Unused.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @examples
#' \donttest{
#' # data(gf_nepa17)
#' # zg <- phase.zg(df = gf_nepa17[1:500, ], TreeNum = 1, beta = 0.1)
#'
#' # Raw GRO and TWD time series
#' # plot(zg)
#'
#' # Daily aggregated series
#' # plot(zg, temporal = "daily")
#'
#' # Monthly TWD plot
#' # plot(zg, temporal = "monthly", type = "twd")
#'
#' # ABr bar plot
#' # plot(zg, temporal = "monthly", type = "abr")
#'
#' # Daily phase timing, similar to plot.SC_output balance
#' # plot(zg, type = "balance", temporal = "daily",
#' # balance_mode = "time_of_day")
#'
#' # Stacked phase-duration balance
#' # plot(zg, type = "balance", temporal = "daily",
#' # balance_mode = "duration")
#'
#' # Transition plot in day-of-year coordinates
#' # plot(zg, type = "transition", x_axis = "doy",
#' # DOY = c(50, 100), Year = 2017)
#'
#' # Boxplots of peak TWD
#' # plot(zg, type = "boxplot", stat = "max.twd",
#' # temporal = "daily", box_group = "doy")
#'
#' # Boxplots of AUC.load
#' # plot(zg, type = "boxplot", stat = "AUC.load",
#' # temporal = "monthly", box_group = "month_of_year")
#'
#' # Boxplots of ABr.value
#' # plot(zg, type = "boxplot", stat = "ABr.value",
#' # temporal = "monthly", box_group = "month_of_year")
#' # plot(zg, view = "boxplot", stat = "ABr.value", temporal = "monthly")
#' }
#'
#' @method plot ZG_output
#' @importFrom dplyr %>% mutate rename group_by summarise arrange if_else filter bind_rows
#' @importFrom tidyr pivot_longer
#' @importFrom tibble tibble
#' @importFrom lubridate floor_date yday month
#' @importFrom ggplot2 ggplot aes geom_line geom_area geom_rect geom_point
#' geom_col geom_boxplot geom_vline facet_wrap theme_minimal theme theme_bw
#' element_text labs scale_color_manual scale_fill_manual position_jitter
#' scale_x_date scale_x_discrete scale_y_continuous coord_cartesian
#' @export
plot.ZG_output <- function(x,
y = NULL,
DOY = NULL,
Year = NULL,
view = NULL,
type = c("gro_twd", "phase_ribbon", "transition",
"twd", "abr", "phase_summary",
"balance", "boxplot"),
temporal = c("raw", "daily", "monthly"),
x_axis = c("time", "doy"),
stat = c("Duration_h", "Magnitude", "rate",
"max.twd", "Avg.twd", "STD.twd",
"AUC.load", "AUC.total", "ABr.value"),
phase = c("auto", "TWD", "GRO", "all"),
box_group = c("period", "month_of_year", "doy"),
twd_fun = c("mean", "max"),
gro_fun = c("max", "mean"),
transition_linetype = "dashed",
transition_alpha = 0.55,
cols = c(TWD = "red", GRO = "blue"),
singleton_as_points = TRUE,
balance_mode = c("time_of_day", "duration"),
balance_gap = c("carry_forward", "observed_only"),
...) {
TIME <- dm <- GRO <- TWD <- Variable <- Value <- NULL
Start <- End <- Duration_h <- Magnitude <- rate <- max.twd <- Max.twd.time <- NULL
ABr.value <- Avg.twd <- STD.twd <- AUC.load <- AUC.total <- period <- Phase <- NULL
group_label <- value <- hours <- Phases <- Metric <- xmin <- xmax <- X <- NULL
group_date <- seg_start <- seg_end <- day <- ymin <- ymax <- unit <- total_hours <- NULL
if (!inherits(x, "ZG_output")) {
stop("Input must be an object of class 'ZG_output'.", call. = FALSE)
}
if (is.null(names(cols))) {
if (length(cols) != 2) {
stop(
"'cols' must be a named vector with TWD and GRO, ",
"or an unnamed vector of length 2.",
call. = FALSE
)
}
names(cols) <- c("TWD", "GRO")
}
if (!all(c("TWD", "GRO") %in% names(cols))) {
stop("'cols' must contain names 'TWD' and 'GRO'.", call. = FALSE)
}
if (!is.null(view)) {
if (!view %in% c("series", "boxplot")) {
stop("'view' must be either 'series' or 'boxplot'.", call. = FALSE)
}
if (view == "boxplot") {
type <- "boxplot"
}
}
type <- match.arg(type)
temporal <- match.arg(temporal)
x_axis <- match.arg(x_axis)
stat <- match.arg(stat)
phase <- match.arg(phase)
box_group <- match.arg(box_group)
twd_fun <- match.arg(twd_fun)
gro_fun <- match.arg(gro_fun)
balance_mode <- match.arg(balance_mode)
balance_gap <- match.arg(balance_gap)
if (type == "boxplot" && temporal == "raw") {
stop("For type = 'boxplot', temporal must be 'daily' or 'monthly'.", call. = FALSE)
}
.time_unit <- function(x) {
switch(
x,
raw = NA_character_,
daily = "day",
monthly = "month"
)
}
.next_period_end <- function(period_start, unit) {
if (unit == "month") {
seq(period_start, length.out = 2, by = "1 month")[2]
} else {
seq(period_start, length.out = 2, by = "1 day")[2]
}
}
.x_values <- function(tt, temporal, x_axis) {
if (x_axis == "time") {
return(list(x = tt, xlab = "Time"))
}
if (temporal == "monthly") {
warning(
"x_axis = 'doy' is not meaningful for temporal = 'monthly'; using time instead.",
call. = FALSE
)
return(list(x = tt, xlab = "Time"))
}
list(
x = lubridate::yday(tt),
xlab = "Day of year"
)
}
.phase_factor <- function(x) {
factor(x, levels = c(1, 2), labels = c("TWD", "GRO"))
}
.safe_sum <- function(z) {
if (all(is.na(z))) NA_real_ else sum(z, na.rm = TRUE)
}
.safe_mean <- function(z) {
if (all(is.na(z))) NA_real_ else mean(z, na.rm = TRUE)
}
.safe_max <- function(z) {
if (all(is.na(z))) NA_real_ else max(z, na.rm = TRUE)
}
.safe_time_of_max <- function(values, times) {
if (all(is.na(values))) {
return(as.POSIXct(NA, tz = attr(times, "tzone")))
}
idx <- which.max(replace(values, is.na(values), -Inf))
times[idx][1]
}
.dominant_phase <- function(z) {
if (sum(z == 2L, na.rm = TRUE) >= sum(z == 1L, na.rm = TRUE)) {
2L
} else {
1L
}
}
.hour_decimal <- function(tt) {
day_start <- lubridate::floor_date(tt, unit = "day")
as.numeric(difftime(tt, day_start, units = "hours"))
}
.truncate_zg <- function(obj, DOY = NULL, Year = NULL) {
if (is.null(DOY) && is.null(Year)) {
td <- obj$ZG_phase %>%
dplyr::rename(TIME = 1, dm = 2, Phases = 3, TWD = 4, GRO = 5)
} else {
if (is.null(DOY) || is.null(Year) || length(DOY) != 2) {
stop("Provide both DOY = c(start, end) and Year, or neither.", call. = FALSE)
}
td <- dendro.truncate(obj$ZG_phase, CalYear = Year, DOY = DOY) %>%
dplyr::rename(TIME = 1, dm = 2, Phases = 3, TWD = 4, GRO = 5)
}
if (nrow(td) == 0) {
stop("No data available for the requested plotting period.", call. = FALSE)
}
cyc <- obj$ZG_cycle %>%
dplyr::filter(
End >= min(td$TIME, na.rm = TRUE),
Start <= max(td$TIME, na.rm = TRUE)
)
list(td = td, cyc = cyc)
}
.aggregate_points <- function(td, temporal, twd_fun, gro_fun) {
if (temporal == "raw") {
td$Phase <- .phase_factor(td$Phases)
return(td)
}
unit <- .time_unit(temporal)
agg_twd <- function(z) {
if (all(is.na(z))) {
return(NA_real_)
}
if (twd_fun == "mean") {
mean(z, na.rm = TRUE)
} else {
max(z, na.rm = TRUE)
}
}
agg_gro <- function(z) {
if (all(is.na(z))) {
return(NA_real_)
}
if (gro_fun == "max") {
max(z, na.rm = TRUE)
} else {
mean(z, na.rm = TRUE)
}
}
td %>%
dplyr::mutate(period = lubridate::floor_date(TIME, unit = unit)) %>%
dplyr::group_by(period) %>%
dplyr::summarise(
dm = mean(dm, na.rm = TRUE),
GRO = agg_gro(GRO),
TWD = agg_twd(TWD),
Phases = .dominant_phase(Phases),
.groups = "drop"
) %>%
dplyr::rename(TIME = period) %>%
dplyr::mutate(Phase = .phase_factor(Phases))
}
.aggregate_cycles <- function(cyc, temporal) {
if (temporal == "raw") {
return(cyc)
}
unit <- .time_unit(temporal)
cyc %>%
dplyr::mutate(period = lubridate::floor_date(Start, unit = unit)) %>%
dplyr::group_by(period, Phases) %>%
dplyr::summarise(
Start = min(Start, na.rm = TRUE),
End = max(End, na.rm = TRUE),
Duration_h = sum(Duration_h, na.rm = TRUE),
Magnitude = .safe_sum(Magnitude),
rate = {
mag <- .safe_sum(Magnitude)
dur <- sum(Duration_h, na.rm = TRUE)
if (is.na(mag) || dur == 0) {
NA_real_
} else {
mag * 1000 / dur
}
},
max.twd = .safe_max(max.twd),
Max.twd.time = .safe_time_of_max(max.twd, Max.twd.time),
AUC.load = .safe_sum(AUC.load),
AUC.total = .safe_sum(AUC.total),
ABr.value = .safe_sum(ABr.value),
Avg.twd = .safe_mean(Avg.twd),
STD.twd = .safe_mean(STD.twd),
DOY = lubridate::yday(min(Start, na.rm = TRUE)),
.groups = "drop"
)
}
.balance_segments <- function(td_raw, temporal, gap_mode = "carry_forward") {
unit <- if (temporal == "monthly") "month" else "day"
td_i <- td_raw %>%
dplyr::filter(!is.na(TIME), !is.na(Phases)) %>%
dplyr::arrange(TIME)
if (nrow(td_i) < 2) {
stop("At least two time steps are required to calculate phase balance.", call. = FALSE)
}
diffs_sec <- as.numeric(diff(td_i$TIME), units = "secs")
diffs_sec <- diffs_sec[is.finite(diffs_sec) & diffs_sec > 0]
if (length(diffs_sec) == 0) {
stop("Could not determine the temporal resolution of the phase series.", call. = FALSE)
}
res_sec <- stats::median(diffs_sec, na.rm = TRUE)
td_i$interval_start <- td_i$TIME
td_i$interval_end <- dplyr::lead(td_i$TIME)
td_i$interval_end[nrow(td_i)] <- td_i$TIME[nrow(td_i)] + res_sec
if (gap_mode == "observed_only") {
interval_sec <- as.numeric(
difftime(td_i$interval_end, td_i$interval_start, units = "secs")
)
long_gap <- is.finite(interval_sec) & interval_sec > res_sec * 1.5
td_i$interval_end[long_gap] <- td_i$interval_start[long_gap] + res_sec
}
rows <- list()
k <- 0L
for (i in seq_len(nrow(td_i))) {
st <- td_i$interval_start[i]
en <- td_i$interval_end[i]
ph <- td_i$Phases[i]
if (is.na(st) || is.na(en) || is.na(ph) || en <= st) {
next
}
while (st < en) {
period_start <- lubridate::floor_date(st, unit = unit)
period_end <- .next_period_end(period_start, unit = unit)
seg_end <- if (en < period_end) {
en
} else {
period_end
}
if (seg_end <= st) {
break
}
k <- k + 1L
rows[[k]] <- tibble::tibble(
period = period_start,
seg_start = st,
seg_end = seg_end,
Phases = ph,
hours = as.numeric(difftime(seg_end, st, units = "hours"))
)
st <- seg_end
}
}
if (length(rows) == 0) {
stop("No valid phase-duration segments could be calculated.", call. = FALSE)
}
dplyr::bind_rows(rows) %>%
dplyr::mutate(
Phase = .phase_factor(Phases),
unit = unit
)
}
out <- .truncate_zg(x, DOY = DOY, Year = Year)
td_raw <- out$td
cyc_raw <- out$cyc
td <- .aggregate_points(
td_raw,
temporal = temporal,
twd_fun = twd_fun,
gro_fun = gro_fun
)
cyc <- .aggregate_cycles(cyc_raw, temporal = temporal)
if (type %in% c("gro_twd", "phase_ribbon", "transition", "twd")) {
td_plot <- td
}
if (type == "gro_twd") {
data_long <- tidyr::pivot_longer(
td_plot,
cols = c(GRO, TWD),
names_to = "Variable",
values_to = "Value"
) %>%
dplyr::mutate(
Value = dplyr::if_else(Variable == "TWD", -Value, Value),
Variable = factor(Variable, levels = c("GRO", "TWD"))
) %>%
dplyr::arrange(TIME)
xv <- .x_values(data_long$TIME, temporal, x_axis)
data_long$X <- xv$x
data_long <- dplyr::filter(data_long, !is.na(X))
p <- ggplot2::ggplot(data_long, ggplot2::aes(x = X, y = Value)) +
ggplot2::geom_line(
data = dplyr::filter(data_long, Variable == "GRO"),
color = unname(cols["GRO"]),
linewidth = 0.5
) +
ggplot2::geom_area(
data = dplyr::filter(data_long, Variable == "TWD"),
fill = unname(cols["TWD"]),
alpha = 0.5
) +
ggplot2::facet_wrap(~ Variable, scales = "free_y", ncol = 1) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Value"
} else {
paste0("Value (", temporal, ")")
}
) +
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 == "phase_ribbon") {
td_r <- td_plot %>% dplyr::filter(!is.na(Phases))
xv <- .x_values(td_r$TIME, temporal, x_axis)
td_r$X <- xv$x
td_r <- dplyr::filter(td_r, !is.na(X))
rr <- rle(td_r$Phases)
idx_end <- cumsum(rr$lengths)
idx_start <- c(1, head(idx_end + 1, -1))
ribbons <- tibble::tibble(
xmin = td_r$X[idx_start],
xmax = td_r$X[idx_end],
Phase = factor(rr$values, levels = c(1, 2), labels = c("TWD", "GRO"))
)
p <- ggplot2::ggplot(td_r, ggplot2::aes(x = X)) +
ggplot2::geom_rect(
data = ribbons,
ggplot2::aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf, fill = Phase),
inherit.aes = FALSE,
alpha = 0.15
) +
ggplot2::geom_line(ggplot2::aes(y = dm), linewidth = 0.4) +
ggplot2::geom_line(
ggplot2::aes(y = GRO),
linetype = 2,
linewidth = 0.5,
color = unname(cols["GRO"])
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = "forestgreen")
) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Stem radius / growth line (mm)"
} else {
paste0("Stem radius / growth line (", temporal, ")")
},
fill = "Phase"
) +
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 == "transition") {
td_tr <- td_plot %>%
dplyr::filter(!is.na(Phases)) %>%
dplyr::arrange(TIME)
xv <- .x_values(td_tr$TIME, temporal, x_axis)
td_tr$X <- xv$x
td_tr <- dplyr::filter(td_tr, !is.na(X))
trans_idx <- which(diff(td_tr$Phases) != 0) + 1L
transitions <- if (length(trans_idx) > 0) {
td_tr[trans_idx, , drop = FALSE]
} else {
td_tr[0, , drop = FALSE]
}
p <- ggplot2::ggplot(td_tr, ggplot2::aes(x = X, y = dm)) +
ggplot2::geom_line(color = "grey35", linewidth = 0.5) +
ggplot2::geom_point(ggplot2::aes(color = Phase), size = 1.8) +
ggplot2::geom_vline(
data = transitions,
ggplot2::aes(xintercept = X, color = Phase),
linetype = transition_linetype,
alpha = transition_alpha,
show.legend = FALSE
) +
ggplot2::geom_point(
data = transitions,
ggplot2::aes(x = X, y = dm, color = Phase),
size = 3,
inherit.aes = FALSE
) +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"])),
name = "New phase"
) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Stem radius (mm)"
} else {
paste0("Stem radius (", temporal, ")")
},
title = "Zero-growth phase transitions"
) +
ggplot2::theme_minimal() +
ggplot2::theme_bw() +
ggplot2::theme(
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 12)
)
print(p)
return(invisible(p))
}
if (type == "twd") {
twd_peaks <- cyc %>%
dplyr::filter(Phases == 1L, !is.na(max.twd), !is.na(Max.twd.time))
xv <- .x_values(td_plot$TIME, temporal, x_axis)
td_plot$X <- xv$x
td_plot <- dplyr::filter(td_plot, !is.na(X))
if (nrow(twd_peaks) > 0) {
twd_peaks$X <- .x_values(twd_peaks$Max.twd.time, temporal, x_axis)$x
twd_peaks <- dplyr::filter(twd_peaks, !is.na(X))
}
p <- ggplot2::ggplot(td_plot, ggplot2::aes(x = X, y = TWD)) +
ggplot2::geom_area(fill = unname(cols["TWD"]), alpha = 0.35) +
ggplot2::geom_line(color = unname(cols["TWD"]), linewidth = 0.4) +
ggplot2::geom_point(
data = twd_peaks,
ggplot2::aes(x = X, y = max.twd),
inherit.aes = FALSE,
color = "black",
size = 1.8
) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Tree water deficit (mm)"
} else {
paste0("Tree water deficit (", temporal, ")")
}
) +
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 == "abr") {
abr_dat <- cyc %>%
dplyr::filter(Phases == 1L, !is.na(ABr.value))
abr_time <- if (temporal != "raw" && "period" %in% names(abr_dat)) {
abr_dat$period
} else {
abr_dat$Start
}
xv <- .x_values(abr_time, temporal, x_axis)
abr_dat$X <- xv$x
abr_dat <- dplyr::filter(abr_dat, !is.na(X))
p <- ggplot2::ggplot(abr_dat, ggplot2::aes(x = X, y = ABr.value)) +
ggplot2::geom_col(fill = "firebrick", alpha = 0.8) +
ggplot2::labs(
x = if (x_axis == "doy" && temporal != "monthly") {
"Day of year"
} else if (temporal == "raw") {
"Start of TWD phase"
} else {
"Aggregation period"
},
y = if (temporal == "raw") {
"ABr value"
} else {
paste0("ABr value (", temporal, ")")
}
) +
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 == "phase_summary") {
sum_dat <- cyc %>%
dplyr::mutate(
Phase = .phase_factor(Phases),
Metric = ifelse(Phases == 1L, max.twd, Magnitude)
)
p <- ggplot2::ggplot(sum_dat, ggplot2::aes(x = Duration_h, y = Metric, color = Phase)) +
ggplot2::geom_point(size = 2, alpha = 0.8) +
ggplot2::facet_wrap(~ Phase, scales = "free_y") +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::labs(
x = "Phase duration (h)",
y = "Phase metric (max.twd for TWD, Magnitude for GRO)",
color = "Phase"
) +
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 == "balance") {
unit <- if (temporal == "monthly") "month" else "day"
if (temporal == "raw") {
warning("For type = 'balance', raw data are summarised to daily balance.", call. = FALSE)
unit <- "day"
}
if (balance_mode == "duration") {
bal_segments <- .balance_segments(
td_raw = td_raw,
temporal = if (unit == "month") "monthly" else "daily",
gap_mode = balance_gap
)
bal <- bal_segments %>%
dplyr::group_by(period, Phase) %>%
dplyr::summarise(hours = sum(hours, na.rm = TRUE), .groups = "drop")
if (unit == "day") {
totals <- bal %>%
dplyr::group_by(period) %>%
dplyr::summarise(total_hours = sum(hours, na.rm = TRUE), .groups = "drop")
incomplete_days <- totals %>%
dplyr::filter(abs(total_hours - 24) > 1e-6)
if (nrow(incomplete_days) > 0) {
warning(
"Some daily balance totals do not equal 24 h. ",
"This usually indicates an incomplete first/last day, missing timestamps, ",
"or truncated data. Use balance_gap = 'carry_forward' to allocate gaps ",
"to the previous phase, or balance_gap = 'observed_only' to count only observed intervals.",
call. = FALSE
)
}
}
if (x_axis == "doy" && unit == "day") {
bal$X <- lubridate::yday(bal$period)
x_lab <- "Day of year"
use_date_scale <- FALSE
} else {
if (x_axis == "doy" && unit == "month") {
warning(
"x_axis = 'doy' is not meaningful for monthly balance; using time instead.",
call. = FALSE
)
}
bal$X <- as.Date(bal$period)
x_lab <- if (unit == "month") "Month" else "Date"
use_date_scale <- TRUE
}
bal <- dplyr::filter(bal, !is.na(X))
p <- ggplot2::ggplot(bal, ggplot2::aes(x = X, y = hours, fill = Phase)) +
ggplot2::geom_col() +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"])),
name = "Phase"
) +
ggplot2::labs(
title = "Phase balance",
x = x_lab,
y = if (unit == "month") "Hours per month" else "Hours per day"
) +
ggplot2::theme_minimal() +
ggplot2::theme(
legend.position = "top",
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 12)
)
if (use_date_scale) {
p <- p + ggplot2::scale_x_date()
}
if (unit == "day") {
p <- p + ggplot2::coord_cartesian(ylim = c(0, 24))
}
print(p)
return(invisible(p))
}
if (balance_mode == "time_of_day") {
if (unit == "month") {
warning(
"balance_mode = 'time_of_day' is designed for daily phase timing. ",
"Using daily phase timing even though temporal = 'monthly'.",
call. = FALSE
)
}
seg <- .balance_segments(
td_raw = td_raw,
temporal = "daily",
gap_mode = balance_gap
) %>%
dplyr::mutate(
day = as.Date(period),
ymin = .hour_decimal(seg_start),
ymax = dplyr::if_else(
as.Date(seg_end) > day & abs(.hour_decimal(seg_end)) < 1e-8,
24,
.hour_decimal(seg_end)
)
) %>%
dplyr::filter(!is.na(ymin), !is.na(ymax), ymax > ymin)
totals <- seg %>%
dplyr::group_by(day) %>%
dplyr::summarise(total_hours = sum(ymax - ymin, na.rm = TRUE), .groups = "drop")
incomplete_days <- totals %>%
dplyr::filter(abs(total_hours - 24) > 1e-6)
if (nrow(incomplete_days) > 0) {
warning(
"Some daily balance totals do not equal 24 h. ",
"This usually indicates an incomplete first/last day, missing timestamps, ",
"or truncated data. The plot still preserves the actual within-day timing ",
"of observed or gap-filled phase intervals.",
call. = FALSE
)
}
if (x_axis == "doy") {
seg$X <- lubridate::yday(seg$day)
seg$xmin <- seg$X - 0.45
seg$xmax <- seg$X + 0.45
x_lab <- "Day of year"
use_date_scale <- FALSE
} else {
seg$X <- seg$day
seg$xmin <- seg$day - 0.45
seg$xmax <- seg$day + 0.45
x_lab <- "Date"
use_date_scale <- TRUE
}
p <- ggplot2::ggplot(seg) +
ggplot2::geom_rect(
ggplot2::aes(
xmin = xmin,
xmax = xmax,
ymin = ymin,
ymax = ymax,
fill = Phase
),
color = NA
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"])),
name = "Phase"
) +
ggplot2::scale_y_continuous(
limits = c(0, 24),
breaks = seq(0, 24, by = 6)
) +
ggplot2::labs(
title = "Daily phase timing",
x = x_lab,
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)
)
if (use_date_scale) {
p <- p + ggplot2::scale_x_date()
}
print(p)
return(invisible(p))
}
}
if (type == "boxplot") {
cyc_box <- cyc
if (phase == "auto") {
if (stat %in% c("max.twd", "Avg.twd", "STD.twd", "AUC.load", "AUC.total", "ABr.value")) {
phase <- "TWD"
} else if (stat %in% c("Magnitude", "rate")) {
phase <- "GRO"
} else {
phase <- "all"
}
}
if (phase == "TWD") {
cyc_box <- dplyr::filter(cyc_box, Phases == 1L)
}
if (phase == "GRO") {
cyc_box <- dplyr::filter(cyc_box, Phases == 2L)
}
cyc_box <- cyc_box %>%
dplyr::mutate(Phase = .phase_factor(Phases))
box_time <- if ("period" %in% names(cyc_box)) {
cyc_box$period
} else {
cyc_box$Start
}
use_continuous_time <- FALSE
if (box_group == "period") {
use_continuous_time <- TRUE
if (temporal == "daily") {
full_dates <- seq(
as.Date(min(td_raw$TIME, na.rm = TRUE)),
as.Date(max(td_raw$TIME, na.rm = TRUE)),
by = "day"
)
cyc_box$group_date <- as.Date(box_time)
x_lab <- "Date"
x_limits <- range(full_dates)
box_width <- 0.8
point_jitter_width <- 0.2
date_breaks <- if (length(full_dates) > 120) "3 months" else "1 month"
date_labels <- if (length(unique(format(full_dates, "%Y"))) > 1) "%Y-%m" else "%m-%d"
} else {
full_dates <- seq(
as.Date(lubridate::floor_date(min(td_raw$TIME, na.rm = TRUE), "month")),
as.Date(lubridate::floor_date(max(td_raw$TIME, na.rm = TRUE), "month")),
by = "month"
)
cyc_box$group_date <- as.Date(lubridate::floor_date(box_time, "month"))
x_lab <- "Date"
x_limits <- range(full_dates)
box_width <- 20
point_jitter_width <- 5
date_breaks <- if (length(full_dates) > 24) "6 months" else "1 month"
date_labels <- if (length(unique(format(full_dates, "%Y"))) > 1) "%Y-%m" else "%b"
}
} else if (box_group == "month_of_year") {
cyc_box$group_label <- factor(
as.character(lubridate::month(box_time, label = TRUE, abbr = TRUE)),
levels = month.abb
)
x_lab <- "Month of year"
} else {
full_levels <- as.character(1:366)
cyc_box$group_label <- factor(
as.character(lubridate::yday(box_time)),
levels = full_levels
)
x_lab <- "Day of year"
}
cyc_box$value <- cyc_box[[stat]]
if (use_continuous_time) {
cyc_box <- cyc_box %>%
dplyr::filter(!is.na(value), !is.na(group_date))
} else {
cyc_box <- cyc_box %>%
dplyr::filter(!is.na(value), !is.na(group_label))
}
if (nrow(cyc_box) == 0) {
stop("No non-missing values available for the selected statistic/phase.", call. = FALSE)
}
if (use_continuous_time) {
count_key <- format(cyc_box$group_date, "%Y-%m-%d")
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
count_key <- paste(count_key, cyc_box$Phase, sep = "_")
}
} else {
count_key <- as.character(cyc_box$group_label)
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
count_key <- paste(count_key, cyc_box$Phase, sep = "_")
}
}
counts <- table(count_key)
singleton_groups <- all(counts <= 1)
if (singleton_groups) {
if (box_group == "period") {
warning(
"Each group has only one value, so boxplots collapse to lines. ",
"This is expected for box_group = 'period'. ",
"Use box_group = 'month_of_year' for pooled seasonal boxplots.",
call. = FALSE
)
} else if (box_group == "doy" && length(unique(format(cyc_box$Start, "%Y"))) <= 1) {
warning(
"box_group = 'doy' needs repeated day-of-year values across multiple years. ",
"Your selected data appear to cover only one year, so each DOY has one value.",
call. = FALSE
)
} else {
warning("Each boxplot group has only one value; boxes collapse to lines.", call. = FALSE)
}
}
if (singleton_groups && isTRUE(singleton_as_points)) {
if (use_continuous_time) {
p <- ggplot2::ggplot(
cyc_box,
ggplot2::aes(x = group_date, y = value, color = Phase)
) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = point_jitter_width, height = 0),
size = 2,
alpha = 0.8
) +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::scale_x_date(
limits = x_limits,
date_breaks = date_breaks,
date_labels = date_labels
) +
ggplot2::labs(
x = x_lab,
y = stat,
color = "Phase"
) +
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(
cyc_box,
ggplot2::aes(x = group_label, y = value, color = Phase)
) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = 0.12, height = 0),
size = 2,
alpha = 0.8
) +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::labs(
x = x_lab,
y = stat,
color = "Phase"
) +
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)
) +
ggplot2::scale_x_discrete(drop = FALSE)
}
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
p <- p + ggplot2::facet_wrap(~ Phase, scales = "free_y")
}
print(p)
return(invisible(p))
}
if (use_continuous_time) {
p <- ggplot2::ggplot(
cyc_box,
ggplot2::aes(
x = group_date,
y = value,
fill = Phase,
group = interaction(group_date, Phase, drop = TRUE)
)
) +
ggplot2::geom_boxplot(width = box_width, outlier.shape = NA) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = point_jitter_width, height = 0),
size = 1.2,
alpha = 0.6
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::scale_x_date(
limits = x_limits,
date_breaks = date_breaks,
date_labels = date_labels
) +
ggplot2::labs(
x = x_lab,
y = stat,
fill = "Phase"
) +
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(
cyc_box,
ggplot2::aes(x = group_label, y = value, fill = Phase)
) +
ggplot2::geom_boxplot(outlier.shape = NA) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = 0.15, height = 0),
size = 1.2,
alpha = 0.6
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::labs(
x = x_lab,
y = stat,
fill = "Phase"
) +
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)
) +
ggplot2::scale_x_discrete(drop = FALSE)
}
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
p <- p + ggplot2::facet_wrap(~ Phase, scales = "free_y")
}
print(p)
return(invisible(p))
}
}
Try the dendRoAnalyst package in your browser
Any scripts or data that you put into this service are public.
dendRoAnalyst documentation built on May 20, 2026, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.ZG_output
Documented in plot.ZG_output
#' @title Plot method for zero-growth output
#'
#' @description
#' Unified S3 plotting method for objects returned by \code{phase.zg()}.
#' It supports time-series views of GRO and TWD, phase ribbons and transitions,
#' TWD-only plots, ABr summaries, phase-balance views, and boxplots of
#' phase-level statistics. For boxplots, phase summaries are first assembled to
#' the requested daily or monthly time scale, and rows with missing x-axis
#' grouping labels are removed before plotting.
#'
#' The plotting method operates on the two tables returned by
#' \code{phase.zg()}:
#' \itemize{
#' \item \code{ZG_phase}: point-level time series with \code{TIME},
#' \code{dm}, \code{Phases}, \code{TWD}, and \code{GRO};
#' \item \code{ZG_cycle}: phase-level summaries including
#' \code{max.twd}, \code{AUC.load}, \code{AUC.total}, and
#' \code{ABr.value}.
#' }
#'
#' @details
#' \strong{ABr formula shown by this plot method.}
#'
#' When \code{type = "abr"} or when \code{stat = "ABr.value"} is used in
#' boxplots, the displayed values come directly from \code{phase.zg()}:
#'
#' \deqn{
#' ABr = \max(TWD) \times \left(\frac{AUC_{load}}{\max(TWD)}\right)^{\beta}
#' }
#'
#' where
#'
#' \deqn{AUC_{load} = \int_{t_{start}}^{t_{max}} TWD(t)\,dt}
#'
#' and
#'
#' \deqn{AUC_{total} = \int_{t_{start}}^{t_{end}} TWD(t)\,dt.}
#'
#' This function does not recompute ABr or AUC values; it only visualizes
#' the values already stored in the \code{ZG_output} object.
#'
#' \strong{Aggregation behavior.}
#'
#' For \code{temporal = "daily"} or \code{"monthly"}, point-level series
#' are aggregated before plotting. TWD and GRO can be summarized using
#' \code{"mean"} or \code{"max"}, depending on \code{twd_fun} and
#' \code{gro_fun}.
#'
#' For phase-level summaries in \code{ZG_cycle}, the plot method aggregates
#' statistics by daily or monthly period. In those aggregated summaries:
#' \itemize{
#' \item \code{Duration_h}, \code{Magnitude}, \code{AUC.load},
#' \code{AUC.total}, and \code{ABr.value} are summed;
#' \item \code{max.twd} is taken as the maximum;
#' \item \code{Avg.twd} and \code{STD.twd} are averaged;
#' \item \code{rate} is recalculated from aggregated magnitude and duration.
#' }
#'
#' \strong{Phase balance.}
#'
#' For \code{type = "balance"}, two balance styles are available through
#' \code{balance_mode}. With \code{balance_mode = "time_of_day"}, the plot
#' draws TWD and GRO intervals on a 0--24 hour y-axis, preserving the order in
#' which phases occurred within each day, similar to the balance option in
#' \code{plot.SC_output()}. With \code{balance_mode = "duration"}, the plot
#' draws stacked bars showing the total number of hours spent in each phase.
#' Missing or irregular intervals can be handled with \code{balance_gap}.
#'
#' \strong{Boxplot statistics.}
#'
#' The following phase-level variables can be visualized in
#' \code{type = "boxplot"}. For \code{box_group = "period"}, boxplots are
#' drawn on a continuous date axis spanning the selected multi-annual data
#' window, so empty days or months remain visible as gaps rather than being
#' removed from the x-axis:
#' \itemize{
#' \item \code{"Duration_h"}
#' \item \code{"Magnitude"}
#' \item \code{"rate"}
#' \item \code{"max.twd"}
#' \item \code{"Avg.twd"}
#' \item \code{"STD.twd"}
#' \item \code{"AUC.load"}
#' \item \code{"AUC.total"}
#' \item \code{"ABr.value"}
#' }
#'
#' \strong{Phase handling.}
#'
#' For \code{phase = "auto"}, TWD-specific statistics
#' (\code{max.twd}, \code{Avg.twd}, \code{STD.twd}, \code{AUC.load},
#' \code{AUC.total}, \code{ABr.value}) automatically select TWD phases,
#' whereas GRO-specific variables (\code{Magnitude}, \code{rate}) select GRO
#' phases.
#'
#' @param x Object of class \code{"ZG_output"} returned by \code{phase.zg()}.
#' @param y Unused.
#' @param DOY Optional numeric vector of length 2 giving start and end
#' day-of-year for truncating the plotting window.
#' @param Year Optional numeric year used together with \code{DOY}.
#' @param view Optional backward-compatible argument. Use \code{"boxplot"}
#' to force boxplot mode.
#' @param type Plot type. One of \code{"gro_twd"}, \code{"phase_ribbon"},
#' \code{"transition"}, \code{"twd"}, \code{"abr"},
#' \code{"phase_summary"}, \code{"balance"}, or \code{"boxplot"}.
#' @param temporal Temporal scale. One of \code{"raw"}, \code{"daily"},
#' or \code{"monthly"}.
#' @param x_axis X-axis style for time-based plots. One of \code{"time"} or
#' \code{"doy"}. When \code{"doy"} is used, timestamps are converted to
#' day-of-year. For \code{temporal = "monthly"}, \code{"doy"} is not
#' meaningful and time is used instead.
#' @param stat For \code{type = "boxplot"}, one of
#' \code{"Duration_h"}, \code{"Magnitude"}, \code{"rate"},
#' \code{"max.twd"}, \code{"Avg.twd"}, \code{"STD.twd"},
#' \code{"AUC.load"}, \code{"AUC.total"}, or \code{"ABr.value"}.
#' @param phase For \code{type = "boxplot"}, one of \code{"auto"},
#' \code{"TWD"}, \code{"GRO"}, or \code{"all"}.
#' @param box_group For \code{type = "boxplot"}, grouping mode:
#' \code{"period"}, \code{"month_of_year"}, or \code{"doy"}.
#' \code{"period"} uses exact dates/months from the already aggregated
#' daily or monthly phase summaries; the other two pool values across
#' repeated seasonal positions. Rows with missing grouping labels are
#' removed before plotting.
#' @param twd_fun Aggregation function for point-level TWD when
#' \code{temporal != "raw"}. One of \code{"mean"} or \code{"max"}.
#' @param gro_fun Aggregation function for point-level GRO when
#' \code{temporal != "raw"}. One of \code{"max"} or \code{"mean"}.
#' @param transition_linetype Line type used for transition markers.
#' @param transition_alpha Transparency of transition markers.
#' @param cols Vector of two colours for TWD and GRO. It may be named
#' (\code{c(TWD = "red", GRO = "blue")}) or unnamed length 2.
#' @param singleton_as_points Logical. If \code{TRUE}, the boxplot panel
#' falls back to a point plot when every group contains only one value.
#' @param balance_mode For \code{type = "balance"}, one of
#' \code{"time_of_day"} or \code{"duration"}. \code{"time_of_day"} draws
#' continuous TWD/GRO intervals along the y-axis from 0 to 24 h, preserving
#' when phases occurred within each day. \code{"duration"} draws stacked bars
#' of total hours spent in each phase per day or month.
#' @param balance_gap For \code{type = "balance"}, one of
#' \code{"carry_forward"} or \code{"observed_only"}.
#' \code{"carry_forward"} assigns gaps to the previous observed phase;
#' \code{"observed_only"} counts only observed intervals.
#' @param ... Unused.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @examples
#' \donttest{
#' # data(gf_nepa17)
#' # zg <- phase.zg(df = gf_nepa17[1:500, ], TreeNum = 1, beta = 0.1)
#'
#' # Raw GRO and TWD time series
#' # plot(zg)
#'
#' # Daily aggregated series
#' # plot(zg, temporal = "daily")
#'
#' # Monthly TWD plot
#' # plot(zg, temporal = "monthly", type = "twd")
#'
#' # ABr bar plot
#' # plot(zg, temporal = "monthly", type = "abr")
#'
#' # Daily phase timing, similar to plot.SC_output balance
#' # plot(zg, type = "balance", temporal = "daily",
#' # balance_mode = "time_of_day")
#'
#' # Stacked phase-duration balance
#' # plot(zg, type = "balance", temporal = "daily",
#' # balance_mode = "duration")
#'
#' # Transition plot in day-of-year coordinates
#' # plot(zg, type = "transition", x_axis = "doy",
#' # DOY = c(50, 100), Year = 2017)
#'
#' # Boxplots of peak TWD
#' # plot(zg, type = "boxplot", stat = "max.twd",
#' # temporal = "daily", box_group = "doy")
#'
#' # Boxplots of AUC.load
#' # plot(zg, type = "boxplot", stat = "AUC.load",
#' # temporal = "monthly", box_group = "month_of_year")
#'
#' # Boxplots of ABr.value
#' # plot(zg, type = "boxplot", stat = "ABr.value",
#' # temporal = "monthly", box_group = "month_of_year")
#' # plot(zg, view = "boxplot", stat = "ABr.value", temporal = "monthly")
#' }
#'
#' @method plot ZG_output
#' @importFrom dplyr %>% mutate rename group_by summarise arrange if_else filter bind_rows
#' @importFrom tidyr pivot_longer
#' @importFrom tibble tibble
#' @importFrom lubridate floor_date yday month
#' @importFrom ggplot2 ggplot aes geom_line geom_area geom_rect geom_point
#' geom_col geom_boxplot geom_vline facet_wrap theme_minimal theme theme_bw
#' element_text labs scale_color_manual scale_fill_manual position_jitter
#' scale_x_date scale_x_discrete scale_y_continuous coord_cartesian
#' @export
plot.ZG_output <- function(x,
y = NULL,
DOY = NULL,
Year = NULL,
view = NULL,
type = c("gro_twd", "phase_ribbon", "transition",
"twd", "abr", "phase_summary",
"balance", "boxplot"),
temporal = c("raw", "daily", "monthly"),
x_axis = c("time", "doy"),
stat = c("Duration_h", "Magnitude", "rate",
"max.twd", "Avg.twd", "STD.twd",
"AUC.load", "AUC.total", "ABr.value"),
phase = c("auto", "TWD", "GRO", "all"),
box_group = c("period", "month_of_year", "doy"),
twd_fun = c("mean", "max"),
gro_fun = c("max", "mean"),
transition_linetype = "dashed",
transition_alpha = 0.55,
cols = c(TWD = "red", GRO = "blue"),
singleton_as_points = TRUE,
balance_mode = c("time_of_day", "duration"),
balance_gap = c("carry_forward", "observed_only"),
...) {
TIME <- dm <- GRO <- TWD <- Variable <- Value <- NULL
Start <- End <- Duration_h <- Magnitude <- rate <- max.twd <- Max.twd.time <- NULL
ABr.value <- Avg.twd <- STD.twd <- AUC.load <- AUC.total <- period <- Phase <- NULL
group_label <- value <- hours <- Phases <- Metric <- xmin <- xmax <- X <- NULL
group_date <- seg_start <- seg_end <- day <- ymin <- ymax <- unit <- total_hours <- NULL
if (!inherits(x, "ZG_output")) {
stop("Input must be an object of class 'ZG_output'.", call. = FALSE)
}
if (is.null(names(cols))) {
if (length(cols) != 2) {
stop(
"'cols' must be a named vector with TWD and GRO, ",
"or an unnamed vector of length 2.",
call. = FALSE
)
}
names(cols) <- c("TWD", "GRO")
}
if (!all(c("TWD", "GRO") %in% names(cols))) {
stop("'cols' must contain names 'TWD' and 'GRO'.", call. = FALSE)
}
if (!is.null(view)) {
if (!view %in% c("series", "boxplot")) {
stop("'view' must be either 'series' or 'boxplot'.", call. = FALSE)
}
if (view == "boxplot") {
type <- "boxplot"
}
}
type <- match.arg(type)
temporal <- match.arg(temporal)
x_axis <- match.arg(x_axis)
stat <- match.arg(stat)
phase <- match.arg(phase)
box_group <- match.arg(box_group)
twd_fun <- match.arg(twd_fun)
gro_fun <- match.arg(gro_fun)
balance_mode <- match.arg(balance_mode)
balance_gap <- match.arg(balance_gap)
if (type == "boxplot" && temporal == "raw") {
stop("For type = 'boxplot', temporal must be 'daily' or 'monthly'.", call. = FALSE)
}
.time_unit <- function(x) {
switch(
x,
raw = NA_character_,
daily = "day",
monthly = "month"
)
}
.next_period_end <- function(period_start, unit) {
if (unit == "month") {
seq(period_start, length.out = 2, by = "1 month")[2]
} else {
seq(period_start, length.out = 2, by = "1 day")[2]
}
}
.x_values <- function(tt, temporal, x_axis) {
if (x_axis == "time") {
return(list(x = tt, xlab = "Time"))
}
if (temporal == "monthly") {
warning(
"x_axis = 'doy' is not meaningful for temporal = 'monthly'; using time instead.",
call. = FALSE
)
return(list(x = tt, xlab = "Time"))
}
list(
x = lubridate::yday(tt),
xlab = "Day of year"
)
}
.phase_factor <- function(x) {
factor(x, levels = c(1, 2), labels = c("TWD", "GRO"))
}
.safe_sum <- function(z) {
if (all(is.na(z))) NA_real_ else sum(z, na.rm = TRUE)
}
.safe_mean <- function(z) {
if (all(is.na(z))) NA_real_ else mean(z, na.rm = TRUE)
}
.safe_max <- function(z) {
if (all(is.na(z))) NA_real_ else max(z, na.rm = TRUE)
}
.safe_time_of_max <- function(values, times) {
if (all(is.na(values))) {
return(as.POSIXct(NA, tz = attr(times, "tzone")))
}
idx <- which.max(replace(values, is.na(values), -Inf))
times[idx][1]
}
.dominant_phase <- function(z) {
if (sum(z == 2L, na.rm = TRUE) >= sum(z == 1L, na.rm = TRUE)) {
2L
} else {
1L
}
}
.hour_decimal <- function(tt) {
day_start <- lubridate::floor_date(tt, unit = "day")
as.numeric(difftime(tt, day_start, units = "hours"))
}
.truncate_zg <- function(obj, DOY = NULL, Year = NULL) {
if (is.null(DOY) && is.null(Year)) {
td <- obj$ZG_phase %>%
dplyr::rename(TIME = 1, dm = 2, Phases = 3, TWD = 4, GRO = 5)
} else {
if (is.null(DOY) || is.null(Year) || length(DOY) != 2) {
stop("Provide both DOY = c(start, end) and Year, or neither.", call. = FALSE)
}
td <- dendro.truncate(obj$ZG_phase, CalYear = Year, DOY = DOY) %>%
dplyr::rename(TIME = 1, dm = 2, Phases = 3, TWD = 4, GRO = 5)
}
if (nrow(td) == 0) {
stop("No data available for the requested plotting period.", call. = FALSE)
}
cyc <- obj$ZG_cycle %>%
dplyr::filter(
End >= min(td$TIME, na.rm = TRUE),
Start <= max(td$TIME, na.rm = TRUE)
)
list(td = td, cyc = cyc)
}
.aggregate_points <- function(td, temporal, twd_fun, gro_fun) {
if (temporal == "raw") {
td$Phase <- .phase_factor(td$Phases)
return(td)
}
unit <- .time_unit(temporal)
agg_twd <- function(z) {
if (all(is.na(z))) {
return(NA_real_)
}
if (twd_fun == "mean") {
mean(z, na.rm = TRUE)
} else {
max(z, na.rm = TRUE)
}
}
agg_gro <- function(z) {
if (all(is.na(z))) {
return(NA_real_)
}
if (gro_fun == "max") {
max(z, na.rm = TRUE)
} else {
mean(z, na.rm = TRUE)
}
}
td %>%
dplyr::mutate(period = lubridate::floor_date(TIME, unit = unit)) %>%
dplyr::group_by(period) %>%
dplyr::summarise(
dm = mean(dm, na.rm = TRUE),
GRO = agg_gro(GRO),
TWD = agg_twd(TWD),
Phases = .dominant_phase(Phases),
.groups = "drop"
) %>%
dplyr::rename(TIME = period) %>%
dplyr::mutate(Phase = .phase_factor(Phases))
}
.aggregate_cycles <- function(cyc, temporal) {
if (temporal == "raw") {
return(cyc)
}
unit <- .time_unit(temporal)
cyc %>%
dplyr::mutate(period = lubridate::floor_date(Start, unit = unit)) %>%
dplyr::group_by(period, Phases) %>%
dplyr::summarise(
Start = min(Start, na.rm = TRUE),
End = max(End, na.rm = TRUE),
Duration_h = sum(Duration_h, na.rm = TRUE),
Magnitude = .safe_sum(Magnitude),
rate = {
mag <- .safe_sum(Magnitude)
dur <- sum(Duration_h, na.rm = TRUE)
if (is.na(mag) || dur == 0) {
NA_real_
} else {
mag * 1000 / dur
}
},
max.twd = .safe_max(max.twd),
Max.twd.time = .safe_time_of_max(max.twd, Max.twd.time),
AUC.load = .safe_sum(AUC.load),
AUC.total = .safe_sum(AUC.total),
ABr.value = .safe_sum(ABr.value),
Avg.twd = .safe_mean(Avg.twd),
STD.twd = .safe_mean(STD.twd),
DOY = lubridate::yday(min(Start, na.rm = TRUE)),
.groups = "drop"
)
}
.balance_segments <- function(td_raw, temporal, gap_mode = "carry_forward") {
unit <- if (temporal == "monthly") "month" else "day"
td_i <- td_raw %>%
dplyr::filter(!is.na(TIME), !is.na(Phases)) %>%
dplyr::arrange(TIME)
if (nrow(td_i) < 2) {
stop("At least two time steps are required to calculate phase balance.", call. = FALSE)
}
diffs_sec <- as.numeric(diff(td_i$TIME), units = "secs")
diffs_sec <- diffs_sec[is.finite(diffs_sec) & diffs_sec > 0]
if (length(diffs_sec) == 0) {
stop("Could not determine the temporal resolution of the phase series.", call. = FALSE)
}
res_sec <- stats::median(diffs_sec, na.rm = TRUE)
td_i$interval_start <- td_i$TIME
td_i$interval_end <- dplyr::lead(td_i$TIME)
td_i$interval_end[nrow(td_i)] <- td_i$TIME[nrow(td_i)] + res_sec
if (gap_mode == "observed_only") {
interval_sec <- as.numeric(
difftime(td_i$interval_end, td_i$interval_start, units = "secs")
)
long_gap <- is.finite(interval_sec) & interval_sec > res_sec * 1.5
td_i$interval_end[long_gap] <- td_i$interval_start[long_gap] + res_sec
}
rows <- list()
k <- 0L
for (i in seq_len(nrow(td_i))) {
st <- td_i$interval_start[i]
en <- td_i$interval_end[i]
ph <- td_i$Phases[i]
if (is.na(st) || is.na(en) || is.na(ph) || en <= st) {
next
}
while (st < en) {
period_start <- lubridate::floor_date(st, unit = unit)
period_end <- .next_period_end(period_start, unit = unit)
seg_end <- if (en < period_end) {
en
} else {
period_end
}
if (seg_end <= st) {
break
}
k <- k + 1L
rows[[k]] <- tibble::tibble(
period = period_start,
seg_start = st,
seg_end = seg_end,
Phases = ph,
hours = as.numeric(difftime(seg_end, st, units = "hours"))
)
st <- seg_end
}
}
if (length(rows) == 0) {
stop("No valid phase-duration segments could be calculated.", call. = FALSE)
}
dplyr::bind_rows(rows) %>%
dplyr::mutate(
Phase = .phase_factor(Phases),
unit = unit
)
}
out <- .truncate_zg(x, DOY = DOY, Year = Year)
td_raw <- out$td
cyc_raw <- out$cyc
td <- .aggregate_points(
td_raw,
temporal = temporal,
twd_fun = twd_fun,
gro_fun = gro_fun
)
cyc <- .aggregate_cycles(cyc_raw, temporal = temporal)
if (type %in% c("gro_twd", "phase_ribbon", "transition", "twd")) {
td_plot <- td
}
if (type == "gro_twd") {
data_long <- tidyr::pivot_longer(
td_plot,
cols = c(GRO, TWD),
names_to = "Variable",
values_to = "Value"
) %>%
dplyr::mutate(
Value = dplyr::if_else(Variable == "TWD", -Value, Value),
Variable = factor(Variable, levels = c("GRO", "TWD"))
) %>%
dplyr::arrange(TIME)
xv <- .x_values(data_long$TIME, temporal, x_axis)
data_long$X <- xv$x
data_long <- dplyr::filter(data_long, !is.na(X))
p <- ggplot2::ggplot(data_long, ggplot2::aes(x = X, y = Value)) +
ggplot2::geom_line(
data = dplyr::filter(data_long, Variable == "GRO"),
color = unname(cols["GRO"]),
linewidth = 0.5
) +
ggplot2::geom_area(
data = dplyr::filter(data_long, Variable == "TWD"),
fill = unname(cols["TWD"]),
alpha = 0.5
) +
ggplot2::facet_wrap(~ Variable, scales = "free_y", ncol = 1) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Value"
} else {
paste0("Value (", temporal, ")")
}
) +
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 == "phase_ribbon") {
td_r <- td_plot %>% dplyr::filter(!is.na(Phases))
xv <- .x_values(td_r$TIME, temporal, x_axis)
td_r$X <- xv$x
td_r <- dplyr::filter(td_r, !is.na(X))
rr <- rle(td_r$Phases)
idx_end <- cumsum(rr$lengths)
idx_start <- c(1, head(idx_end + 1, -1))
ribbons <- tibble::tibble(
xmin = td_r$X[idx_start],
xmax = td_r$X[idx_end],
Phase = factor(rr$values, levels = c(1, 2), labels = c("TWD", "GRO"))
)
p <- ggplot2::ggplot(td_r, ggplot2::aes(x = X)) +
ggplot2::geom_rect(
data = ribbons,
ggplot2::aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf, fill = Phase),
inherit.aes = FALSE,
alpha = 0.15
) +
ggplot2::geom_line(ggplot2::aes(y = dm), linewidth = 0.4) +
ggplot2::geom_line(
ggplot2::aes(y = GRO),
linetype = 2,
linewidth = 0.5,
color = unname(cols["GRO"])
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = "forestgreen")
) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Stem radius / growth line (mm)"
} else {
paste0("Stem radius / growth line (", temporal, ")")
},
fill = "Phase"
) +
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 == "transition") {
td_tr <- td_plot %>%
dplyr::filter(!is.na(Phases)) %>%
dplyr::arrange(TIME)
xv <- .x_values(td_tr$TIME, temporal, x_axis)
td_tr$X <- xv$x
td_tr <- dplyr::filter(td_tr, !is.na(X))
trans_idx <- which(diff(td_tr$Phases) != 0) + 1L
transitions <- if (length(trans_idx) > 0) {
td_tr[trans_idx, , drop = FALSE]
} else {
td_tr[0, , drop = FALSE]
}
p <- ggplot2::ggplot(td_tr, ggplot2::aes(x = X, y = dm)) +
ggplot2::geom_line(color = "grey35", linewidth = 0.5) +
ggplot2::geom_point(ggplot2::aes(color = Phase), size = 1.8) +
ggplot2::geom_vline(
data = transitions,
ggplot2::aes(xintercept = X, color = Phase),
linetype = transition_linetype,
alpha = transition_alpha,
show.legend = FALSE
) +
ggplot2::geom_point(
data = transitions,
ggplot2::aes(x = X, y = dm, color = Phase),
size = 3,
inherit.aes = FALSE
) +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"])),
name = "New phase"
) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Stem radius (mm)"
} else {
paste0("Stem radius (", temporal, ")")
},
title = "Zero-growth phase transitions"
) +
ggplot2::theme_minimal() +
ggplot2::theme_bw() +
ggplot2::theme(
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 12)
)
print(p)
return(invisible(p))
}
if (type == "twd") {
twd_peaks <- cyc %>%
dplyr::filter(Phases == 1L, !is.na(max.twd), !is.na(Max.twd.time))
xv <- .x_values(td_plot$TIME, temporal, x_axis)
td_plot$X <- xv$x
td_plot <- dplyr::filter(td_plot, !is.na(X))
if (nrow(twd_peaks) > 0) {
twd_peaks$X <- .x_values(twd_peaks$Max.twd.time, temporal, x_axis)$x
twd_peaks <- dplyr::filter(twd_peaks, !is.na(X))
}
p <- ggplot2::ggplot(td_plot, ggplot2::aes(x = X, y = TWD)) +
ggplot2::geom_area(fill = unname(cols["TWD"]), alpha = 0.35) +
ggplot2::geom_line(color = unname(cols["TWD"]), linewidth = 0.4) +
ggplot2::geom_point(
data = twd_peaks,
ggplot2::aes(x = X, y = max.twd),
inherit.aes = FALSE,
color = "black",
size = 1.8
) +
ggplot2::labs(
x = xv$xlab,
y = if (temporal == "raw") {
"Tree water deficit (mm)"
} else {
paste0("Tree water deficit (", temporal, ")")
}
) +
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 == "abr") {
abr_dat <- cyc %>%
dplyr::filter(Phases == 1L, !is.na(ABr.value))
abr_time <- if (temporal != "raw" && "period" %in% names(abr_dat)) {
abr_dat$period
} else {
abr_dat$Start
}
xv <- .x_values(abr_time, temporal, x_axis)
abr_dat$X <- xv$x
abr_dat <- dplyr::filter(abr_dat, !is.na(X))
p <- ggplot2::ggplot(abr_dat, ggplot2::aes(x = X, y = ABr.value)) +
ggplot2::geom_col(fill = "firebrick", alpha = 0.8) +
ggplot2::labs(
x = if (x_axis == "doy" && temporal != "monthly") {
"Day of year"
} else if (temporal == "raw") {
"Start of TWD phase"
} else {
"Aggregation period"
},
y = if (temporal == "raw") {
"ABr value"
} else {
paste0("ABr value (", temporal, ")")
}
) +
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 == "phase_summary") {
sum_dat <- cyc %>%
dplyr::mutate(
Phase = .phase_factor(Phases),
Metric = ifelse(Phases == 1L, max.twd, Magnitude)
)
p <- ggplot2::ggplot(sum_dat, ggplot2::aes(x = Duration_h, y = Metric, color = Phase)) +
ggplot2::geom_point(size = 2, alpha = 0.8) +
ggplot2::facet_wrap(~ Phase, scales = "free_y") +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::labs(
x = "Phase duration (h)",
y = "Phase metric (max.twd for TWD, Magnitude for GRO)",
color = "Phase"
) +
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 == "balance") {
unit <- if (temporal == "monthly") "month" else "day"
if (temporal == "raw") {
warning("For type = 'balance', raw data are summarised to daily balance.", call. = FALSE)
unit <- "day"
}
if (balance_mode == "duration") {
bal_segments <- .balance_segments(
td_raw = td_raw,
temporal = if (unit == "month") "monthly" else "daily",
gap_mode = balance_gap
)
bal <- bal_segments %>%
dplyr::group_by(period, Phase) %>%
dplyr::summarise(hours = sum(hours, na.rm = TRUE), .groups = "drop")
if (unit == "day") {
totals <- bal %>%
dplyr::group_by(period) %>%
dplyr::summarise(total_hours = sum(hours, na.rm = TRUE), .groups = "drop")
incomplete_days <- totals %>%
dplyr::filter(abs(total_hours - 24) > 1e-6)
if (nrow(incomplete_days) > 0) {
warning(
"Some daily balance totals do not equal 24 h. ",
"This usually indicates an incomplete first/last day, missing timestamps, ",
"or truncated data. Use balance_gap = 'carry_forward' to allocate gaps ",
"to the previous phase, or balance_gap = 'observed_only' to count only observed intervals.",
call. = FALSE
)
}
}
if (x_axis == "doy" && unit == "day") {
bal$X <- lubridate::yday(bal$period)
x_lab <- "Day of year"
use_date_scale <- FALSE
} else {
if (x_axis == "doy" && unit == "month") {
warning(
"x_axis = 'doy' is not meaningful for monthly balance; using time instead.",
call. = FALSE
)
}
bal$X <- as.Date(bal$period)
x_lab <- if (unit == "month") "Month" else "Date"
use_date_scale <- TRUE
}
bal <- dplyr::filter(bal, !is.na(X))
p <- ggplot2::ggplot(bal, ggplot2::aes(x = X, y = hours, fill = Phase)) +
ggplot2::geom_col() +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"])),
name = "Phase"
) +
ggplot2::labs(
title = "Phase balance",
x = x_lab,
y = if (unit == "month") "Hours per month" else "Hours per day"
) +
ggplot2::theme_minimal() +
ggplot2::theme(
legend.position = "top",
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 12)
)
if (use_date_scale) {
p <- p + ggplot2::scale_x_date()
}
if (unit == "day") {
p <- p + ggplot2::coord_cartesian(ylim = c(0, 24))
}
print(p)
return(invisible(p))
}
if (balance_mode == "time_of_day") {
if (unit == "month") {
warning(
"balance_mode = 'time_of_day' is designed for daily phase timing. ",
"Using daily phase timing even though temporal = 'monthly'.",
call. = FALSE
)
}
seg <- .balance_segments(
td_raw = td_raw,
temporal = "daily",
gap_mode = balance_gap
) %>%
dplyr::mutate(
day = as.Date(period),
ymin = .hour_decimal(seg_start),
ymax = dplyr::if_else(
as.Date(seg_end) > day & abs(.hour_decimal(seg_end)) < 1e-8,
24,
.hour_decimal(seg_end)
)
) %>%
dplyr::filter(!is.na(ymin), !is.na(ymax), ymax > ymin)
totals <- seg %>%
dplyr::group_by(day) %>%
dplyr::summarise(total_hours = sum(ymax - ymin, na.rm = TRUE), .groups = "drop")
incomplete_days <- totals %>%
dplyr::filter(abs(total_hours - 24) > 1e-6)
if (nrow(incomplete_days) > 0) {
warning(
"Some daily balance totals do not equal 24 h. ",
"This usually indicates an incomplete first/last day, missing timestamps, ",
"or truncated data. The plot still preserves the actual within-day timing ",
"of observed or gap-filled phase intervals.",
call. = FALSE
)
}
if (x_axis == "doy") {
seg$X <- lubridate::yday(seg$day)
seg$xmin <- seg$X - 0.45
seg$xmax <- seg$X + 0.45
x_lab <- "Day of year"
use_date_scale <- FALSE
} else {
seg$X <- seg$day
seg$xmin <- seg$day - 0.45
seg$xmax <- seg$day + 0.45
x_lab <- "Date"
use_date_scale <- TRUE
}
p <- ggplot2::ggplot(seg) +
ggplot2::geom_rect(
ggplot2::aes(
xmin = xmin,
xmax = xmax,
ymin = ymin,
ymax = ymax,
fill = Phase
),
color = NA
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"])),
name = "Phase"
) +
ggplot2::scale_y_continuous(
limits = c(0, 24),
breaks = seq(0, 24, by = 6)
) +
ggplot2::labs(
title = "Daily phase timing",
x = x_lab,
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)
)
if (use_date_scale) {
p <- p + ggplot2::scale_x_date()
}
print(p)
return(invisible(p))
}
}
if (type == "boxplot") {
cyc_box <- cyc
if (phase == "auto") {
if (stat %in% c("max.twd", "Avg.twd", "STD.twd", "AUC.load", "AUC.total", "ABr.value")) {
phase <- "TWD"
} else if (stat %in% c("Magnitude", "rate")) {
phase <- "GRO"
} else {
phase <- "all"
}
}
if (phase == "TWD") {
cyc_box <- dplyr::filter(cyc_box, Phases == 1L)
}
if (phase == "GRO") {
cyc_box <- dplyr::filter(cyc_box, Phases == 2L)
}
cyc_box <- cyc_box %>%
dplyr::mutate(Phase = .phase_factor(Phases))
box_time <- if ("period" %in% names(cyc_box)) {
cyc_box$period
} else {
cyc_box$Start
}
use_continuous_time <- FALSE
if (box_group == "period") {
use_continuous_time <- TRUE
if (temporal == "daily") {
full_dates <- seq(
as.Date(min(td_raw$TIME, na.rm = TRUE)),
as.Date(max(td_raw$TIME, na.rm = TRUE)),
by = "day"
)
cyc_box$group_date <- as.Date(box_time)
x_lab <- "Date"
x_limits <- range(full_dates)
box_width <- 0.8
point_jitter_width <- 0.2
date_breaks <- if (length(full_dates) > 120) "3 months" else "1 month"
date_labels <- if (length(unique(format(full_dates, "%Y"))) > 1) "%Y-%m" else "%m-%d"
} else {
full_dates <- seq(
as.Date(lubridate::floor_date(min(td_raw$TIME, na.rm = TRUE), "month")),
as.Date(lubridate::floor_date(max(td_raw$TIME, na.rm = TRUE), "month")),
by = "month"
)
cyc_box$group_date <- as.Date(lubridate::floor_date(box_time, "month"))
x_lab <- "Date"
x_limits <- range(full_dates)
box_width <- 20
point_jitter_width <- 5
date_breaks <- if (length(full_dates) > 24) "6 months" else "1 month"
date_labels <- if (length(unique(format(full_dates, "%Y"))) > 1) "%Y-%m" else "%b"
}
} else if (box_group == "month_of_year") {
cyc_box$group_label <- factor(
as.character(lubridate::month(box_time, label = TRUE, abbr = TRUE)),
levels = month.abb
)
x_lab <- "Month of year"
} else {
full_levels <- as.character(1:366)
cyc_box$group_label <- factor(
as.character(lubridate::yday(box_time)),
levels = full_levels
)
x_lab <- "Day of year"
}
cyc_box$value <- cyc_box[[stat]]
if (use_continuous_time) {
cyc_box <- cyc_box %>%
dplyr::filter(!is.na(value), !is.na(group_date))
} else {
cyc_box <- cyc_box %>%
dplyr::filter(!is.na(value), !is.na(group_label))
}
if (nrow(cyc_box) == 0) {
stop("No non-missing values available for the selected statistic/phase.", call. = FALSE)
}
if (use_continuous_time) {
count_key <- format(cyc_box$group_date, "%Y-%m-%d")
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
count_key <- paste(count_key, cyc_box$Phase, sep = "_")
}
} else {
count_key <- as.character(cyc_box$group_label)
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
count_key <- paste(count_key, cyc_box$Phase, sep = "_")
}
}
counts <- table(count_key)
singleton_groups <- all(counts <= 1)
if (singleton_groups) {
if (box_group == "period") {
warning(
"Each group has only one value, so boxplots collapse to lines. ",
"This is expected for box_group = 'period'. ",
"Use box_group = 'month_of_year' for pooled seasonal boxplots.",
call. = FALSE
)
} else if (box_group == "doy" && length(unique(format(cyc_box$Start, "%Y"))) <= 1) {
warning(
"box_group = 'doy' needs repeated day-of-year values across multiple years. ",
"Your selected data appear to cover only one year, so each DOY has one value.",
call. = FALSE
)
} else {
warning("Each boxplot group has only one value; boxes collapse to lines.", call. = FALSE)
}
}
if (singleton_groups && isTRUE(singleton_as_points)) {
if (use_continuous_time) {
p <- ggplot2::ggplot(
cyc_box,
ggplot2::aes(x = group_date, y = value, color = Phase)
) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = point_jitter_width, height = 0),
size = 2,
alpha = 0.8
) +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::scale_x_date(
limits = x_limits,
date_breaks = date_breaks,
date_labels = date_labels
) +
ggplot2::labs(
x = x_lab,
y = stat,
color = "Phase"
) +
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(
cyc_box,
ggplot2::aes(x = group_label, y = value, color = Phase)
) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = 0.12, height = 0),
size = 2,
alpha = 0.8
) +
ggplot2::scale_color_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::labs(
x = x_lab,
y = stat,
color = "Phase"
) +
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)
) +
ggplot2::scale_x_discrete(drop = FALSE)
}
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
p <- p + ggplot2::facet_wrap(~ Phase, scales = "free_y")
}
print(p)
return(invisible(p))
}
if (use_continuous_time) {
p <- ggplot2::ggplot(
cyc_box,
ggplot2::aes(
x = group_date,
y = value,
fill = Phase,
group = interaction(group_date, Phase, drop = TRUE)
)
) +
ggplot2::geom_boxplot(width = box_width, outlier.shape = NA) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = point_jitter_width, height = 0),
size = 1.2,
alpha = 0.6
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::scale_x_date(
limits = x_limits,
date_breaks = date_breaks,
date_labels = date_labels
) +
ggplot2::labs(
x = x_lab,
y = stat,
fill = "Phase"
) +
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(
cyc_box,
ggplot2::aes(x = group_label, y = value, fill = Phase)
) +
ggplot2::geom_boxplot(outlier.shape = NA) +
ggplot2::geom_point(
position = ggplot2::position_jitter(width = 0.15, height = 0),
size = 1.2,
alpha = 0.6
) +
ggplot2::scale_fill_manual(
values = c(TWD = unname(cols["TWD"]), GRO = unname(cols["GRO"]))
) +
ggplot2::labs(
x = x_lab,
y = stat,
fill = "Phase"
) +
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)
) +
ggplot2::scale_x_discrete(drop = FALSE)
}
if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
p <- p + ggplot2::facet_wrap(~ Phase, scales = "free_y")
}
print(p)
return(invisible(p))
}
}
Try the dendRoAnalyst package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.