R/dm_wavelet.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' Wavelet analysis of dendrometer series
#'
#' @description
#' Performs continuous wavelet analysis on dendrometer time series.
#'
#' The function can work with:
#' \itemize{
#'   \item raw dendrometer series from a data frame,
#'   \item detrended dendrometer series from \code{dm.detrend.fit()} output,
#'   \item first differences of either raw or detrended series.
#' }
#'
#' It automatically:
#' \itemize{
#'   \item identifies the temporal resolution of the input series,
#'   \item regularizes the series to a complete time grid,
#'   \item optionally interpolates missing values,
#'   \item performs wavelet analysis separately for each selected tree/series,
#'   \item converts wavelet periods to hours for summary and plotting.
#' }
#'
#' @param x Input data. Either:
#'   \describe{
#'     \item{data.frame}{First column must be time, remaining selected numeric
#'       columns are dendrometer series.}
#'     \item{dm_detrended object}{Uses \code{x$detrended_data} as input.}
#'   }
#' @param TreeNum Either \code{"all"} to use all available dendrometer series,
#'   a numeric vector selecting series by position, or a character vector of
#'   series names.
#' @param source Which series to analyze. One of:
#'   \describe{
#'     \item{`"auto"`}{Uses raw series for a data frame input and detrended
#'       series for a \code{dm_detrended} object.}
#'     \item{`"raw"`}{Use raw dendrometer series from a data frame input.}
#'     \item{`"detrended"`}{Use detrended series from a \code{dm_detrended} object.}
#'     \item{`"first_diff"`}{Use first differences of the available input series.}
#'   }
#' @param detrended_col For \code{dm_detrended} input, which table to use:
#'   \code{"detrended_data"}.
#' @param na_action How to handle missing values after completing the regular time
#'   grid. One of:
#'   \describe{
#'     \item{`"interpolate"`}{Linearly interpolate internal gaps and carry ends
#'       forward/backward.}
#'     \item{`"fail"`}{Stop if missing values are found.}
#'   }
#' @param loess_span Smoothing span passed to \code{WaveletComp::analyze.wavelet()}.
#' @param dj Frequency resolution parameter passed to
#'   \code{WaveletComp::analyze.wavelet()}.
#' @param lowerPeriod Optional lower period bound in native time units of the
#'   detected input resolution. If \code{NULL}, the default from
#'   \code{WaveletComp::analyze.wavelet()} is used.
#' @param upperPeriod Optional upper period bound in native time units of the
#'   detected input resolution. If \code{NULL}, the default from
#'   \code{WaveletComp::analyze.wavelet()} is used.
#' @param make_pval Logical; if \code{TRUE}, compute Monte Carlo significance.
#' @param n_sim Number of simulations used when \code{make_pval = TRUE}.
#' @param verbose Logical; if \code{TRUE}, prints a completion message.
#'
#' @return
#' An object of class \code{"dm_wavelet"} with elements:
#' \describe{
#'   \item{call}{The matched function call.}
#'   \item{input_type}{Either \code{"raw"} or \code{"dm_detrended"}.}
#'   \item{source}{Series source actually used.}
#'   \item{series}{Character vector of analyzed series names.}
#'   \item{time_unit}{Detected time unit of the input series.}
#'   \item{dt}{Detected time step in native units.}
#'   \item{dt_hours}{Detected time step expressed in hours.}
#'   \item{resolution_seconds}{Detected time step in seconds.}
#'   \item{data_used}{Regularized time series used for analysis.}
#'   \item{results}{Named list of wavelet results, one per series. Each entry
#'     contains the original \code{WaveletComp} object plus time and
#'     period-in-hours information.}
#'   \item{settings}{Analysis settings.}
#' }
#'
#' @importFrom dplyr left_join arrange %>%
#' @importFrom tibble as_tibble tibble
#' @importFrom lubridate ymd_hms parse_date_time
#' @importFrom zoo na.approx na.locf
#' @export
dm_wavelet <- function(x,
                       TreeNum = "all",
                       source = c("auto", "raw", "detrended", "first_diff"),
                       detrended_col = c("detrended_data"),
                       na_action = c("interpolate", "fail"),
                       loess_span = 0.75,
                       dj = 1 / 20,
                       lowerPeriod = NULL,
                       upperPeriod = NULL,
                       make_pval = TRUE,
                       n_sim = 10,
                       verbose = TRUE) {

  source <- match.arg(source)
  detrended_col <- match.arg(detrended_col)
  na_action <- match.arg(na_action)

  if (!requireNamespace("WaveletComp", quietly = TRUE)) {
    stop("Package 'WaveletComp' is required for dm_wavelet().")
  }

  cl <- match.call()

  dmw_parse_time <- function(tt) {
    if (inherits(tt, "POSIXct") || inherits(tt, "POSIXt")) return(as.POSIXct(tt))
    if (inherits(tt, "Date")) return(as.POSIXct(tt))
    out <- suppressWarnings(lubridate::ymd_hms(tt, quiet = TRUE))
    if (all(is.na(out))) {
      out <- suppressWarnings(
        lubridate::parse_date_time(
          tt,
          orders = c("ymd HMS", "ymd HM", "ymd",
                     "dmy HMS", "dmy HM", "dmy",
                     "mdy HMS", "mdy HM", "mdy"),
          quiet = TRUE
        )
      )
    }
    out
  }

  dmw_select_cols <- function(dat, TreeNum, meta_cols = character(0)) {
    series_cols <- setdiff(names(dat), meta_cols)
    series_cols <- series_cols[vapply(dat[series_cols], is.numeric, logical(1))]

    if (length(series_cols) == 0) {
      stop("No numeric dendrometer series available.")
    }

    if (is.character(TreeNum) && length(TreeNum) == 1 && tolower(TreeNum) == "all") {
      keep <- series_cols
    } else if (is.numeric(TreeNum)) {
      idx <- as.integer(TreeNum)
      if (any(is.na(idx)) || any(idx < 1) || any(idx > length(series_cols))) {
        stop("'TreeNum' numeric values are out of range.")
      }
      keep <- series_cols[idx]
    } else if (is.character(TreeNum)) {
      miss <- setdiff(TreeNum, series_cols)
      if (length(miss) > 0) {
        warning("These requested series were not found and were ignored: ",
                paste(miss, collapse = ", "))
      }
      keep <- intersect(TreeNum, series_cols)
      if (length(keep) == 0) {
        stop("None of the requested TreeNum names were found.")
      }
    } else {
      stop("'TreeNum' must be 'all', numeric indices, or character names.")
    }

    keep
  }

  dmw_infer_dt <- function(time_vec) {
    dsec <- as.numeric(diff(time_vec), units = "secs")
    dsec <- dsec[is.finite(dsec) & dsec > 0]

    if (length(dsec) == 0) {
      stop("Could not infer time step from TIME column.")
    }

    med <- stats::median(dsec, na.rm = TRUE)

    if (med >= 86400 && abs(med / 86400 - round(med / 86400)) < 1e-6) {
      return(list(dt = med / 86400, unit = "days", sec = med, hours = med / 3600))
    }
    if (med >= 3600 && abs(med / 3600 - round(med / 3600)) < 1e-6) {
      return(list(dt = med / 3600, unit = "hours", sec = med, hours = med / 3600))
    }
    if (med >= 60 && abs(med / 60 - round(med / 60)) < 1e-6) {
      return(list(dt = med / 60, unit = "mins", sec = med, hours = med / 3600))
    }

    list(dt = med, unit = "secs", sec = med, hours = med / 3600)
  }

  dmw_period_to_hours <- function(period, time_unit) {
    if (time_unit == "days") return(period * 24)
    if (time_unit == "hours") return(period)
    if (time_unit == "mins") return(period / 60)
    if (time_unit == "secs") return(period / 3600)
    period
  }

  dmw_complete_regular <- function(dat, dt_sec) {
    tz_use <- attr(dat$TIME, "tzone")
    if (is.null(tz_use) || length(tz_use) == 0) tz_use <- "UTC"

    full_time <- seq(
      from = min(dat$TIME, na.rm = TRUE),
      to = max(dat$TIME, na.rm = TRUE),
      by = dt_sec
    )

    tibble::tibble(TIME = as.POSIXct(full_time, origin = "1970-01-01", tz = tz_use)) %>%
      dplyr::left_join(dat, by = "TIME") %>%
      dplyr::arrange(.data$TIME)
  }

  dmw_fill_na <- function(z) {
    z2 <- zoo::na.approx(z, na.rm = FALSE)
    z2 <- zoo::na.locf(z2, na.rm = FALSE)
    z2 <- zoo::na.locf(z2, fromLast = TRUE, na.rm = FALSE)
    z2
  }

  if (inherits(x, "dm_detrended")) {
    input_type <- "dm_detrended"
    if (source == "raw") {
      stop("source = 'raw' is not valid for a dm_detrended object.")
    }
    dat0 <- tibble::as_tibble(x[[detrended_col]])
    meta_cols <- intersect(
      c("TIME", "doy", "season_label", "season_start", "season_end", "season_day"),
      names(dat0)
    )
    keep_cols <- dmw_select_cols(dat0, TreeNum = TreeNum, meta_cols = meta_cols)
    dat <- dat0[, c("TIME", keep_cols), drop = FALSE]
    source_used <- if (source == "auto") "detrended" else source
  } else if (is.data.frame(x)) {
    input_type <- "raw"
    if (source == "detrended") {
      stop("source = 'detrended' requires a dm_detrended object.")
    }
    dat0 <- tibble::as_tibble(x)
    names(dat0)[1] <- "TIME"
    keep_cols <- dmw_select_cols(dat0, TreeNum = TreeNum, meta_cols = "TIME")
    dat <- dat0[, c("TIME", keep_cols), drop = FALSE]
    source_used <- if (source == "auto") "raw" else source
  } else {
    stop("'x' must be either a data.frame or an object of class 'dm_detrended'.")
  }

  dat$TIME <- dmw_parse_time(dat$TIME)
  if (any(is.na(dat$TIME))) {
    stop("Some timestamps in the first column could not be parsed.")
  }

  dat <- dat[order(dat$TIME), , drop = FALSE]

  dt_info <- dmw_infer_dt(dat$TIME)
  dat_reg <- dmw_complete_regular(dat, dt_sec = dt_info$sec)

  series_names <- setdiff(names(dat_reg), "TIME")

  for (ss in series_names) {
    if (na_action == "interpolate") {
      dat_reg[[ss]] <- dmw_fill_na(dat_reg[[ss]])
    } else {
      if (any(is.na(dat_reg[[ss]]))) {
        stop("Missing values found after regularization. Use na_action = 'interpolate' or pre-clean the series.")
      }
    }
  }

  if (source_used == "first_diff") {
    for (ss in series_names) {
      dat_reg[[ss]] <- c(NA_real_, diff(dat_reg[[ss]]))
    }
    dat_reg <- dat_reg[-1, , drop = FALSE]
  }

  for (ss in series_names) {
    if (any(is.na(dat_reg[[ss]]))) {
      stop("Series '", ss, "' still contains missing values after preprocessing.")
    }
    if (sum(is.finite(dat_reg[[ss]])) < 8) {
      stop("Series '", ss, "' is too short for wavelet analysis after preprocessing.")
    }
  }

  res_list <- vector("list", length(series_names))
  names(res_list) <- series_names

  for (ss in series_names) {
    tmp <- data.frame(series = as.numeric(dat_reg[[ss]]))

    wave_args <- list(
      my.data = tmp,
      my.series = "series",
      loess.span = loess_span,
      dt = dt_info$dt,
      dj = dj,
      make.pval = make_pval,
      n.sim = n_sim
    )

    if (!is.null(lowerPeriod)) wave_args$lowerPeriod <- lowerPeriod
    if (!is.null(upperPeriod)) wave_args$upperPeriod <- upperPeriod

    wt <- do.call(WaveletComp::analyze.wavelet, wave_args)

    period_hours <- dmw_period_to_hours(wt$Period, dt_info$unit)

    res_list[[ss]] <- list(
      wavelet = wt,
      time = dat_reg$TIME,
      series = dat_reg[[ss]],
      period_native = wt$Period,
      period_hours = period_hours,
      avg_power = if (!is.null(wt$Power.avg)) {
        as.numeric(wt$Power.avg)
      } else {
        rowMeans(wt$Power, na.rm = TRUE)
      }
    )
  }

  out <- list(
    call = cl,
    input_type = input_type,
    source = source_used,
    series = series_names,
    time_unit = dt_info$unit,
    dt = dt_info$dt,
    dt_hours = dt_info$hours,
    resolution_seconds = dt_info$sec,
    data_used = dat_reg,
    results = res_list,
    settings = list(
      TreeNum = TreeNum,
      loess_span = loess_span,
      dj = dj,
      lowerPeriod = lowerPeriod,
      upperPeriod = upperPeriod,
      make_pval = make_pval,
      n_sim = n_sim,
      na_action = na_action
    )
  )

  class(out) <- "dm_wavelet"

  if (isTRUE(verbose)) {
    message(
      "dm_wavelet completed for ", length(series_names),
      " series. Detected resolution: ", dt_info$dt, " ", dt_info$unit,
      " (", round(dt_info$hours, 4), " hours)."
    )
  }

  out
}


#' Summarize a dm_wavelet object
#'
#' @description
#' Summarizes wavelet-analysis results produced by \code{dm_wavelet()}.
#'
#' For each analyzed series, the summary reports:
#' \itemize{
#'   \item the dominant period in hours,
#'   \item the corresponding maximum average wavelet power,
#'   \item the analyzed period range in hours,
#'   \item the detected temporal resolution.
#' }
#'
#' @param object An object of class \code{"dm_wavelet"}.
#' @param top_n Integer. Number of strongest periods to report per series.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return An object of class \code{"summary.dm_wavelet"}.
#' @method summary dm_wavelet
#' @export
summary.dm_wavelet <- function(object, top_n = 3, ...) {

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

  if (!is.numeric(top_n) || length(top_n) != 1 || is.na(top_n) || top_n < 1) {
    stop("'top_n' must be a positive integer.")
  }
  top_n <- as.integer(top_n)

  series_rows <- vector("list", length(object$results))
  top_rows <- vector("list", length(object$results))
  nm <- names(object$results)

  for (i in seq_along(object$results)) {
    ss <- nm[i]
    res <- object$results[[i]]

    per_h <- as.numeric(res$period_hours)
    avg_pow <- as.numeric(res$avg_power)

    ok <- is.finite(per_h) & is.finite(avg_pow)

    if (sum(ok) == 0) {
      series_rows[[i]] <- tibble::tibble(
        series = ss,
        n_obs = length(res$series),
        dominant_period_hours = NA_real_,
        max_average_power = NA_real_,
        min_period_hours = NA_real_,
        max_period_hours = NA_real_,
        dt = object$dt,
        dt_hours = object$dt_hours,
        time_unit = object$time_unit
      )

      top_rows[[i]] <- tibble::tibble(
        series = ss,
        rank = integer(0),
        period_hours = numeric(0),
        average_power = numeric(0)
      )
      next
    }

    per_h <- per_h[ok]
    avg_pow <- avg_pow[ok]
    ord <- order(avg_pow, decreasing = TRUE)
    ord_top <- head(ord, top_n)

    series_rows[[i]] <- tibble::tibble(
      series = ss,
      n_obs = length(res$series),
      dominant_period_hours = per_h[ord[1]],
      max_average_power = avg_pow[ord[1]],
      min_period_hours = min(per_h, na.rm = TRUE),
      max_period_hours = max(per_h, na.rm = TRUE),
      dt = object$dt,
      dt_hours = object$dt_hours,
      time_unit = object$time_unit
    )

    top_rows[[i]] <- tibble::tibble(
      series = ss,
      rank = seq_along(ord_top),
      period_hours = per_h[ord_top],
      average_power = avg_pow[ord_top]
    )
  }

  out <- list(
    overview = tibble::tibble(
      input_type = object$input_type,
      source = object$source,
      n_series = length(object$series),
      dt = object$dt,
      dt_hours = object$dt_hours,
      time_unit = object$time_unit
    ),
    series_summary = dplyr::bind_rows(series_rows),
    top_periods = dplyr::bind_rows(top_rows),
    settings = object$settings
  )

  class(out) <- "summary.dm_wavelet"
  out
}


#' @method print summary.dm_wavelet
#' @export
print.summary.dm_wavelet <- function(x, ...) {
  cat("dm_wavelet summary\n")
  cat("------------------\n")
  print(x$overview)

  cat("\nSeries summary:\n")
  print(x$series_summary)

  cat("\nTop periods (hours):\n")
  print(x$top_periods)

  invisible(x)
}


#' @title Plot method for wavelet analysis output
#'
#' @description
#' S3 plotting method for objects returned by \code{dm_wavelet()}.
#'
#' The function provides three main visualization types:
#' \itemize{
#'   \item \code{"series"}: plots the analyzed input dendrometer series over time,
#'   \item \code{"average"}: plots the average wavelet power spectrum against period,
#'   \item \code{"power"}: plots the full wavelet power spectrum as a time-period image.
#' }
#'
#' Periods are shown in \strong{hours}, regardless of the original temporal
#' resolution of the input series. For power plots, the x-axis shows the actual
#' time series timestamps and the y-axis shows the wavelet period in hours.
#'
#' @param x An object of class \code{"dm_wavelet"} returned by
#'   \code{dm_wavelet()}.
#' @param y Unused.
#' @param series Optional character vector giving one or more series names to
#'   plot. If \code{NULL}, the first available series is used.
#' @param type Character string specifying the plot type. One of:
#'   \describe{
#'     \item{`"power"`}{Wavelet power spectrum as a raster image.}
#'     \item{`"average"`}{Average wavelet power spectrum across the full time series.}
#'     \item{`"series"`}{Original analyzed series over time.}
#'   }
#' @param facet Logical. If \code{TRUE} and more than one series is selected,
#'   panels are faceted by series.
#' @param log_period Logical. If \code{TRUE}, the period axis is shown on a
#'   log10 scale where applicable.
#' @param log_power Logical. If \code{TRUE} and \code{type = "power"}, the
#'   plotted wavelet power is transformed using \code{log10(power + eps)} to
#'   improve contrast.
#' @param clip_quantile Optional numeric vector of length 2 giving lower and
#'   upper quantiles used to clip wavelet power before plotting, for example
#'   \code{c(0.01, 0.99)}. This is useful when a few extreme values dominate the
#'   colour scale. Use \code{NULL} to disable clipping.
#' @param show_sig Logical. If \code{TRUE}, overlays significance information
#'   when available. For \code{type = "average"}, significant periods are marked
#'   as points. For \code{type = "power"}, significant regions are shown as
#'   contour lines.
#' @param show_coi Logical. If \code{TRUE} and \code{type = "power"}, the cone of
#'   influence (COI) is added as a shaded overlay.
#' @param coi_fill Fill colour used for the cone of influence shading.
#' @param coi_alpha Numeric transparency of the cone of influence shading.
#' @param siglvl Numeric significance threshold between 0 and 1. Default is
#'   \code{0.05}.
#' @param sig_color Colour used for significance overlays.
#' @param sig_size Line width or point size used for significance overlays.
#' @param main Optional plot title. If \code{NULL}, a default title is used.
#' @param ... Further arguments passed to or from other methods. Currently unused.
#'
#' @details
#' For \code{type = "power"}, the function uses the wavelet power matrix stored
#' in the \code{"dm_wavelet"} object and converts the Fourier periods to hours.
#' The power plot may also show:
#' \itemize{
#'   \item significance contours, when p-values are available,
#'   \item the cone of influence (COI), when COI information is available.
#' }
#'
#' For \code{type = "average"}, the function plots the average wavelet power
#' spectrum and may optionally indicate significant periods when average-spectrum
#' p-values are available.
#'
#' For \code{type = "series"}, the original analyzed series are plotted without
#' any wavelet transformation.
#'
#' @return
#' A \code{ggplot2} object.
#'
#' @seealso
#' \code{\link{dm_wavelet}}, \code{\link[ggplot2]{ggplot}}
#'
#' @examples
#' \donttest{
#' wv <- dm_wavelet(
#'   x = gf_nepa17,
#'   TreeNum = 1:2,
#'   source = "raw",
#'   make_pval = TRUE,
#'   verbose = FALSE
#' )
#'
#' # original series
#' plot(wv, type = "series")
#'
#' # average wavelet power
#' plot(wv, type = "average")
#'
#' # full power spectrum
#' plot(wv, type = "power")
#'
#' # one selected series
#' plot(wv, series = names(wv$results)[1], type = "power")
#'
#' # stronger contrast in the power plot
#' plot(
#'   wv,
#'   type = "power",
#'   log_power = TRUE,
#'   clip_quantile = c(0.05, 0.95)
#' )
#' }
#'
#' @method plot dm_wavelet
#' @importFrom dplyr bind_rows
#' @importFrom tibble tibble
#' @importFrom ggplot2 ggplot aes geom_raster geom_line geom_point geom_contour
#'   geom_polygon facet_wrap theme_bw theme element_text labs scale_y_log10
#'   scale_x_log10 scale_fill_viridis_c coord_cartesian
#' @export
plot.dm_wavelet <- function(x,
                            y = NULL,
                            series = NULL,
                            type = c("power", "average", "series"),
                            facet = TRUE,
                            log_period = TRUE,
                            log_power = TRUE,
                            clip_quantile = c(0.01, 0.99),
                            show_sig = TRUE,
                            show_coi = TRUE,
                            coi_fill = "white",
                            coi_alpha = 0.45,
                            siglvl = 0.05,
                            sig_color = "black",
                            sig_size = 0.4,
                            main = NULL,
                            ...) {
  TIME <- value <- period_hours <- average_power <- power_plot <- pval <- NULL

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

  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("Package 'ggplot2' is required for plot.dm_wavelet().")
  }

  type <- match.arg(type)

  if (is.null(series)) {
    series <- x$series[1]
  }

  miss <- setdiff(series, x$series)
  if (length(miss) > 0) {
    stop("These requested series were not found: ", paste(miss, collapse = ", "))
  }

  if (!is.numeric(siglvl) || length(siglvl) != 1 || is.na(siglvl) ||
      siglvl <= 0 || siglvl >= 1) {
    stop("'siglvl' must be a single number between 0 and 1.")
  }

  dmw_period_to_hours <- function(period, time_unit) {
    if (time_unit == "days") return(period * 24)
    if (time_unit == "hours") return(period)
    if (time_unit == "mins") return(period / 60)
    if (time_unit == "secs") return(period / 3600)
    period
  }

  if (type == "series") {
    df_list <- lapply(series, function(ss) {
      res <- x$results[[ss]]
      tibble::tibble(
        TIME = res$time,
        value = as.numeric(res$series),
        series = ss
      )
    })

    dat <- dplyr::bind_rows(df_list)

    p <- ggplot2::ggplot(dat, ggplot2::aes(x = TIME, y = value, colour = series)) +
      ggplot2::geom_line() +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = if (length(series) > 1) "right" else "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = "Time",
        y = "Series value",
        colour = "Series",
        title = if (is.null(main)) "Input series used for wavelet analysis" else main
      )

    if (length(series) > 1 && isTRUE(facet)) {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ series"), scales = "free_y", ncol = 1)
    }

    return(p)
  }

  if (type == "average") {
    df_list <- lapply(series, function(ss) {
      res <- x$results[[ss]]
      wt <- res$wavelet

      out <- tibble::tibble(
        period_hours = as.numeric(res$period_hours),
        average_power = as.numeric(res$avg_power),
        series = ss
      )

      if (!is.null(wt$Power.avg.pval)) {
        pavg <- as.numeric(wt$Power.avg.pval)
        n <- min(length(pavg), nrow(out))
        out$avg_pval <- NA_real_
        out$avg_pval[seq_len(n)] <- pavg[seq_len(n)]
      } else {
        out$avg_pval <- NA_real_
      }

      out
    })

    dat <- dplyr::bind_rows(df_list)

    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = period_hours, y = average_power, colour = series)
    ) +
      ggplot2::geom_line() +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = if (length(series) > 1) "right" else "none",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = "Period (hours)",
        y = "Average wavelet power",
        colour = "Series",
        title = if (is.null(main)) "Average wavelet power spectrum" else main
      )

    if (isTRUE(log_period) && all(dat$period_hours > 0, na.rm = TRUE)) {
      p <- p + ggplot2::scale_x_log10()
    }

    if (isTRUE(show_sig) && "avg_pval" %in% names(dat) && any(is.finite(dat$avg_pval))) {
      sig_dat <- dat[is.finite(dat$avg_pval) & dat$avg_pval <= siglvl, , drop = FALSE]
      if (nrow(sig_dat) > 0) {
        p <- p + ggplot2::geom_point(
          data = sig_dat,
          mapping = ggplot2::aes(x = period_hours, y = average_power),
          inherit.aes = FALSE,
          colour = sig_color,
          size = max(1, sig_size * 2)
        )
      }
    }

    if (length(series) > 1 && isTRUE(facet)) {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ series"), scales = "free_y", ncol = 1)
    }

    return(p)
  }

  if (type == "power") {
    df_list <- lapply(series, function(ss) {
      res <- x$results[[ss]]
      wt <- res$wavelet
      pow <- as.matrix(wt$Power)

      n_p <- min(nrow(pow), length(res$period_hours))
      n_t <- min(ncol(pow), length(res$time))

      pow <- pow[seq_len(n_p), seq_len(n_t), drop = FALSE]
      per <- as.numeric(res$period_hours[seq_len(n_p)])
      tt <- res$time[seq_len(n_t)]

      df <- expand.grid(
        period_hours = per,
        TIME = tt,
        KEEP.OUT.ATTRS = FALSE,
        stringsAsFactors = FALSE
      )
      df$power <- c(pow)
      df$series <- ss

      if (!is.null(wt$Power.pval)) {
        pmat <- as.matrix(wt$Power.pval)
        pmat <- pmat[seq_len(n_p), seq_len(n_t), drop = FALSE]
        df$pval <- c(pmat)
      } else {
        df$pval <- NA_real_
      }

      df
    })

    dat <- dplyr::bind_rows(df_list)

    # keep only valid positive periods
    dat <- dat[is.finite(dat$period_hours) & dat$period_hours > 0, , drop = FALSE]

    # fixed y range from actual wavelet periods only
    y_rng <- range(dat$period_hours, na.rm = TRUE)

    if (!is.null(clip_quantile)) {
      if (!is.numeric(clip_quantile) || length(clip_quantile) != 2 ||
          any(is.na(clip_quantile)) ||
          clip_quantile[1] < 0 || clip_quantile[2] > 1 ||
          clip_quantile[1] >= clip_quantile[2]) {
        stop("'clip_quantile' must be NULL or a numeric vector like c(0.01, 0.99).")
      }

      qv <- stats::quantile(dat$power, probs = clip_quantile, na.rm = TRUE, names = FALSE)
      dat$power_plot <- pmin(pmax(dat$power, qv[1]), qv[2])
    } else {
      dat$power_plot <- dat$power
    }

    if (isTRUE(log_power)) {
      eps <- max(min(dat$power_plot[dat$power_plot > 0], na.rm = TRUE) / 10, 1e-12)
      if (!is.finite(eps)) eps <- 1e-12
      dat$power_plot <- log10(dat$power_plot + eps)
      fill_lab <- "log10(Power)"
    } else {
      fill_lab <- "Power"
    }

    p <- ggplot2::ggplot(
      dat,
      ggplot2::aes(x = TIME, y = period_hours, fill = power_plot)
    ) +
      ggplot2::geom_raster() +
      ggplot2::scale_fill_viridis_c(option = "D", na.value = "grey85") +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = "right",
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      ) +
      ggplot2::labs(
        x = "Time",
        y = "Period (hours)",
        fill = fill_lab,
        title = if (is.null(main)) "Wavelet power spectrum" else main
      )

    if (isTRUE(show_coi)) {
      coi_list <- lapply(series, function(ss) {
        res <- x$results[[ss]]
        wt <- res$wavelet

        if (is.null(wt$coi.1) || is.null(wt$coi.2) ||
            is.null(wt$axis.1) || is.null(wt$axis.2)) {
          return(NULL)
        }

        tz_use <- attr(res$time, "tzone")
        if (is.null(tz_use) || length(tz_use) == 0) tz_use <- "UTC"

        x_num <- stats::approx(
          x = as.numeric(wt$axis.1),
          y = as.numeric(res$time),
          xout = as.numeric(wt$coi.1),
          rule = 2
        )$y

        period_native <- 2^as.numeric(wt$coi.2)
        period_hours <- dmw_period_to_hours(period_native, x$time_unit)

        # clip COI to actual plotted y range
        period_hours <- pmin(pmax(period_hours, y_rng[1]), y_rng[2])

        tibble::tibble(
          TIME = as.POSIXct(x_num, origin = "1970-01-01", tz = tz_use),
          period_hours = period_hours,
          series = ss
        )
      })

      coi_dat <- dplyr::bind_rows(coi_list)

      if (nrow(coi_dat) > 0) {
        p <- p + ggplot2::geom_polygon(
          data = coi_dat,
          mapping = ggplot2::aes(x = TIME, y = period_hours, group = series),
          inherit.aes = FALSE,
          fill = coi_fill,
          alpha = coi_alpha,
          colour = NA
        )
      }
    }

    if (isTRUE(show_sig) && "pval" %in% names(dat) && any(is.finite(dat$pval))) {
      sig_dat <- dat[is.finite(dat$pval), , drop = FALSE]
      if (nrow(sig_dat) > 0) {
        p <- p + ggplot2::geom_contour(
          data = sig_dat,
          mapping = ggplot2::aes(
            x = TIME,
            y = period_hours,
            z = pval
          ),
          breaks = siglvl,
          colour = sig_color,
          linewidth = sig_size,
          inherit.aes = FALSE
        )
      }
    }

    # force the y axis to the real period range
    if (isTRUE(log_period) && all(y_rng > 0)) {
      p <- p + ggplot2::scale_y_log10(limits = y_rng)
    } else {
      p <- p + ggplot2::coord_cartesian(ylim = y_rng, expand = FALSE)
    }

    if (length(series) > 1 && isTRUE(facet)) {
      p <- p + ggplot2::facet_wrap(stats::as.formula("~ series"), scales = "free_x", ncol = 1)
    }

    return(p)
  }

  stop("Unknown plot type.")
}