R/phase.sc.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Application of the stem-cycle approach to classify dendrometer phases
#'
#' @description
#' Implements the \strong{stem-cycle approach} (Downes et al., 1999; Deslauriers et al., 2011)
#' to divide a dendrometer time series into three biologically meaningful phases:
#' \enumerate{
#'   \item \strong{Shrinkage} (phase = 1): the dendrometer reading decreases
#'         compared to the previous reading.
#'   \item \strong{Expansion} (phase = 2): the dendrometer reading increases
#'         compared to the previous reading, but remains below the previous maximum.
#'   \item \strong{Increment} (phase = 3): the dendrometer reading exceeds
#'         the previous maximum (irreversible stem growth).
#' }
#'
#' For each contiguous phase, the function calculates duration, magnitude, rate,
#' and assigns day-of-year information. Optionally, the dendrometer series may be
#' smoothed before phase calculation to reduce noise and spurious phase changes.
#'
#' @details
#' Classification uses the cumulative maximum of the dendrometer series:
#' \itemize{
#'   \item If the cumulative maximum increases, the phase is labeled \emph{Increment} (3).
#'   \item If the cumulative maximum is constant and the first difference is positive,
#'         the phase is \emph{Expansion} (2).
#'   \item If the cumulative maximum is constant and the first difference is negative,
#'         the phase is \emph{Shrinkage} (1).
#' }
#'
#' The function returns both phase-level summaries (\code{SC_cycle}) and point-level
#' labels (\code{SC_phase}). Optional smoothing uses \code{smooth_dm} with
#' \code{method = "median_mean"} and a window length between 1–24 hours.
#'
#' @references
#' Deslauriers A, Rossi S, Turcotte A, Morin H, Krause C (2011)
#' A three-step procedure in SAS to analyze the time series from automatic dendrometers.
#' \emph{Dendrochronologia} 29:151–161. \doi{10.1016/j.dendro.2011.01.008}
#'
#' Downes G, Beadle C, Worledge D (1999)
#' Daily stem growth patterns in irrigated \emph{Eucalyptus globulus} and \emph{E. nitens} in relation to climate.
#' \emph{Trees} 14:102–111. \doi{10.1007/PL00009752}
#'
#' @param df A data frame with the first column containing date-time in the format
#'   \code{"yyyy-mm-dd HH:MM:SS"} (or convertible to \code{POSIXct}), followed by one or more
#'   dendrometer measurement columns (mm).
#' @param TreeNum Integer. The index of the dendrometer column to analyze.
#'   For example, \code{TreeNum = 1} selects the first dendrometer series after the time column.
#' @param smoothing Numeric or \code{NULL}. Length of the smoothing window in hours
#'   (1–24). If \code{NULL} (default), no smoothing is applied. If provided, the
#'   dendrometer series is smoothed using \code{smooth_dm(method = "median_mean")}
#'   prior to phase classification.
#'
#' @return
#' A list of class \code{"SC_output"} containing:
#' \describe{
#'   \item{SC_cycle}{A tibble with one row per contiguous phase, including:
#'     \itemize{
#'       \item \code{Phases} – Phase type (1 = Shrinkage, 2 = Expansion, 3 = Increment)
#'       \item \code{Start}, \code{End} – \code{POSIXct} start and end time of the phase
#'       \item \code{Duration_h}, \code{Duration_m} – Phase duration (hours, minutes)
#'       \item \code{Magnitude} – Change in dendrometer value during the phase (measurement unit)
#'       \item \code{rate} – Rate of change expressed in (\code{Magnitude*1000/Duration_h}) (eg. \eqn{\mu}m/hour)
#'       \item \code{DOY} – Day-of-year at phase start
#'     }}
#'   \item{SC_phase}{A tibble of point-level values including:
#'     \itemize{
#'       \item \code{TIME} – timestamp
#'       \item \code{dm} – dendrometer measurement
#'       \item \code{Phases} – phase assignment for each timestamp
#'     }}
#' }
#'
#' @seealso
#' \code{\link{phase.zg}} for the zero-growth approach;
#' \code{\link{smooth_dm}} for smoothing dendrometer series.
#'
#' @examples
#' \donttest{
#' library(dendRoAnalyst)
#' data(gf_nepa17)
#'
#' # Apply stem-cycle approach without smoothing
#' sc1 <- phase.sc(df = gf_nepa17, TreeNum = 1)
#' head(sc1$SC_cycle, 5)
#' head(sc1$SC_phase, 5)
#'
#' # Apply with 12-hour smoothing to reduce noise
#' sc2 <- phase.sc(df = gf_nepa17, TreeNum = 1, smoothing = 12)
#' head(sc2$SC_cycle, 5)
#' }
#'
#' @importFrom lubridate ymd_hms yday
#' @importFrom dplyr mutate select rename
#' @importFrom tibble tibble as_tibble
#' @importFrom tidyr pivot_longer
#' @importFrom pspline smooth.Pspline
#'
#' @export

phase.sc<-function(df, TreeNum, smoothing=NULL){
  ##################### function to calculate phases################
  phase_cal <- function(y){
    all_max <- cummax(y)
    all_max_diff <-diff(all_max)
    #all_max_diff <-all_max_diff[2:length(all_max_diff)]
    phases<- vector('integer', length(y)-1)
    # phases[all_max_diff > 0] <- 3
    y_diff <- diff(y)
    # phases[all_max_diff = 0 & y_diff>=0] <- 2
    # phases[all_max_diff = 0 & y_diff<0] <- 1
    phases<-sapply(seq_along(all_max_diff), function(i){
      if(all_max_diff[i]>0){
        ph <- 3
      }else{
        if(y_diff[i]>=0){
          ph <- 2
        }else{
          ph <- 1
        }
      }
      #print(length(phases))
      return(ph)
    })
    return(phases)
  }
  #####################function to cal phase stats #################
  phase_sats<-function(dm_data, y){
    ph_diff<-c(0,diff(y))
    dy<-yday(c(dm_data$TIME[1], dm_data$TIME[ph_diff!=0]))
    ph<-c(y[1],y[ph_diff!=0])
    magn<-c(dm_data$dm[1], dm_data$dm[ph_diff!=0], dm_data$dm[nrow(dm_data)])
    magn<-diff(magn)
    strt_t<-c(dm_data$TIME[1], dm_data$TIME[ph_diff!=0])
    #strt_t<-strt_t[1:length(strt_t)-1]
    end_t <- c(strt_t[2:length(strt_t)], dm_data$TIME[nrow(dm_data)])
    dur_h<-as.numeric(difftime(end_t, strt_t), units='hours')
    dur_m<-as.numeric(difftime(end_t, strt_t), units='mins')
    out <- tibble(
      'Phases' = ph,
      'Start' = strt_t,
      'End' = end_t,
      'Duration_h' = dur_h,
      'Duration_m' = dur_m,
      'Magnitude' = magn,
      'rate' = magn*1000/dur_h,
      'DOY' = dy
    )
    return(out)
  }
  #####################function to calculate resolution ############
  reso_den<-function(input_time){
    time1<-input_time
    reference<-time1[1]
    time_min<-as.integer(difftime(time1,reference, units = "mins"))
    diff_time<-diff(time_min)
    diff2_time<-unique(diff_time)
    reso<-mean(diff2_time)
    if (length(diff2_time) > 1) {
      print(diff2_time)
      warning('Warning: The temporal resolution of dendrometer data is not consistent, For better result, please use dendrometer data with consistent resolution. There may be NA values in the dataset.')
      #cat(paste('Mean temporal resolution is :', round(reso),' minutes.'))
    }else{
      reso<-mean(diff2_time)
      #cat(paste('Temporal resolution is :', round(reso),' minutes.'))
      return(round(reso))
    }
  }
  ##################################################################
  if (!inherits(df[[1]], 'Date') && !inherits(df[[1]], 'POSIXct')) {
    df[[1]] <- ymd_hms(df[[1]])
  }
  dm_data <- as_tibble(df)%>%
    dplyr::select(c(1,TreeNum+1))
  sf<-smoothing
  dm_data <- dm_data %>%
    dplyr::rename(TIME = 1,
           dm = 2)
  r.denro<-reso_den(dm_data$TIME)
  if(is.null(sf)==T){
    y_sm <- dm_data$dm
  }else{
    if(sf<1|sf>24){
      stop('smoothing must be between 1 and 24.')
    }else{
      warning(paste('You are applying smoothing to raw DM data! The smoothing value is ',sf,' hour(s).'))
      y_sm<-smooth_dm(time = dm_data$TIME, dm = dm_data$dm, method = 'median_mean', window_hours = sf)
    }
  }
  y<-phase_cal(y_sm)
  ph_st<-phase_sats(dm_data, y)
  dm_data <- dm_data%>%
    dplyr::mutate('Phases' = c(NA,y))
  out<-list(SC_cycle = ph_st, SC_phase = tibble(dm_data))
  class(out) <- "SC_output"
  return(out)
}