R/dm_smooth.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Smoothing of Dendrometer Time Series
#'
#' @description
#' Applies various smoothing techniques to dendrometer (dm) time series data
#' using a user-defined or automatically detected temporal resolution. The function
#' supports several smoothing methods: robust median+mean, penalized spline,
#' Savitzky-Golay filter, exponential moving average (EMA), and LOESS.
#'
#' @details
#' The function is designed to smooth dendrometer time series (e.g., 10–60 min resolution)
#' while preserving key features such as diurnal fluctuations or long-term growth,
#' depending on the window size. It auto-detects the resolution (in minutes) if not provided.
#'
#' @param time A POSIXct vector representing the time column.
#' @param dm A numeric vector of dendrometer values corresponding to \code{time}.
#' @param resolution_min Integer. The resolution of the time series in minutes. If \code{NULL},
#'   the resolution is auto-detected based on median time difference.
#' @param method Smoothing method. One of: \code{"median_mean"}, \code{"pspline"},
#'   \code{"sg"} (Savitzky-Golay), \code{"ema"} (exponential moving average), or \code{"loess"}.
#' @param window_hours Numeric. Smoothing window length in hours. Converted to points using resolution.
#' @param sg_order Integer. Polynomial order for Savitzky-Golay smoothing. Ignored unless \code{method = "sg"}.
#' @param ema_alpha Numeric. Smoothing factor (0–1) for exponential moving average. If \code{NULL},
#'   it is automatically calculated from the window size.
#'
#' @return A numeric vector of the same length as \code{dm} containing the smoothed values.
#'
#' @examples
#' \dontrun{
#' # Example: Create synthetic dendrometer time series (30-min resolution)
#' time_seq <- seq.POSIXt(from = as.POSIXct("2023-06-01 00:00:00"),
#'                        to   = as.POSIXct("2023-06-03 00:00:00"),
#'                        by   = "30 mins")
#' set.seed(123)
#' dm_raw <- cumsum(rnorm(length(time_seq), mean = 0.005, sd = 0.01)) +
#'           0.1 * sin(2 * pi * as.numeric(difftime(time_seq, min(time_seq), units = "hours")) / 24)
#'
#' # Median plus moving mean smoothing (default robust filter)
#' dm_medmean <- smooth_dm(time = time_seq, dm = dm_raw,
#'                         method = "median_mean", window_hours = 3)
#'
#' # Penalized spline smoothing
#' dm_pspline <- smooth_dm(time = time_seq, dm = dm_raw,
#'                         method = "pspline", window_hours = 6)
#'
#' # Savitzky-Golay filter (requires 'signal' package)
#' if (requireNamespace("signal", quietly = TRUE)) {
#'   dm_sg <- smooth_dm(time = time_seq, dm = dm_raw,
#'                      method = "sg", window_hours = 2, sg_order = 2)
#' }
#'
#' # Exponential moving average smoothing
#' dm_ema <- smooth_dm(time = time_seq, dm = dm_raw,
#'                     method = "ema", window_hours = 4)
#'
#' # LOESS smoothing
#' dm_loess <- smooth_dm(time = time_seq, dm = dm_raw,
#'                       method = "loess", window_hours = 3)
#'
#' # Plot raw and smoothed series
#' plot(time_seq, dm_raw, type = "l", col = "gray", lwd = 1, main = "Smoothed Dendrometer Series",
#'      ylab = "DM", xlab = "Time")
#' lines(time_seq, dm_medmean, col = "blue", lwd = 2)
#' lines(time_seq, dm_pspline, col = "green", lwd = 2)
#' lines(time_seq, dm_ema, col = "orange", lwd = 2)
#' lines(time_seq, dm_loess, col = "purple", lwd = 2)
#' legend("topright", legend = c("Raw", "Median+Mean", "P-spline", "EMA", "LOESS"),
#'        col = c("gray", "blue", "green", "orange", "purple"), lwd = 2)
#' }
#'
#' @importFrom zoo rollapply
#' @importFrom stats runmed loess
#' @importFrom utils tail
#' @export

smooth_dm <- function(time, dm,
                      resolution_min = NULL,
                      method = c("median_mean","pspline","sg","ema","loess"),
                      window_hours = 3,
                      sg_order = 2,
                      ema_alpha = NULL) {
  stopifnot(length(time) == length(dm))
  method <- match.arg(method)

  # auto-detect resolution if needed
  if (is.null(resolution_min)) {
    dt <- as.numeric(diff(as.POSIXct(time)), units = "mins")
    dt <- dt[is.finite(dt) & dt > 0]
    if (!length(dt)) stop("Cannot detect resolution")
    resolution_min <- as.integer(round(stats::median(dt)))
  }
  # window length in points
  k <- max(3L, as.integer(round(window_hours * 60 / resolution_min)))
  if (k %% 2 == 0) k <- k + 1L  # odd window works better for centered filters

  x <- dm

  if (method == "median_mean") {
    x_med <- stats::runmed(x, k = k, endrule = "median")
    y_sm  <- zoo::rollapply(x_med, width = k, align = "center",
                            FUN = function(z) mean(z, na.rm = TRUE),
                            partial = TRUE, fill = NA)
    na_idx <- which(is.na(y_sm))
    if (length(na_idx)) y_sm[na_idx] <- x_med[na_idx]
    return(y_sm)
  }

  if (method == "pspline") {
    spar_pts <- k
    fit <- pspline::smooth.Pspline(x = seq_along(x), y = x, norder = 2, spar = spar_pts)
    return(as.numeric(fit$ysmth))
  }

  if (method == "sg") {
    if (!requireNamespace("signal", quietly = TRUE))
      stop("Install 'signal' for Savitzky Golay")
    if (k <= sg_order) k <- sg_order + 1 + (sg_order %% 2 == 0)
    if (k %% 2 == 0) k <- k + 1
    y_sm <- signal::sgolayfilt(x, p = sg_order, n = k)
    return(as.numeric(y_sm))
  }

  if (method == "ema") {
    if (is.null(ema_alpha)) {
      ema_alpha <- 1 - exp(-log(2) / (k/2))
    }
    y_fwd <- stats::filter(x, ema_alpha, method = "recursive", init = x[1])
    y_bwd <- rev(stats::filter(rev(x), ema_alpha, method = "recursive", init = tail(x,1)))
    y_sm  <- (as.numeric(y_fwd) + as.numeric(y_bwd)) / 2
    return(y_sm)
  }

  if (method == "loess") {
    span <- min(0.99, max(0.05, k / length(x)))
    y_sm <- stats::loess(x ~ seq_along(x), span = span, degree = 1)$fitted
    return(as.numeric(y_sm))
  }
}