R/pick_stalta.R

In coffeemuggler/eseis: Environmental Seismology Toolbox

#' Calculate stal-lta-ratio.
#' 
#' The function calculates the ratio of the short-term-average and 
#' long-term-average of the input signal.
#' 
#' @param data \code{eseis} object, \code{numeric} vector or list of 
#' objects, data set to be processed.
#' 
#' @param time \code{POSIXct} vector, time vector of the signal(s). If not 
#' provided, a synthetic time vector will be created.
#' 
#' @param dt \code{Numeric} value, sampling period. If omitted, either 
#' estimated from \code{time} or set to 0.01 s (i.e., f = 100 Hz). 
#' 
#' @param sta \code{Numeric} value, number of samples for short-term window.
#' 
#' @param lta \code{Numeric} value, number of samples for long-term window.
#' 
#' @param freeze \code{Logical} value, option to freeze lta value at start of 
#' an event. Useful to avoid self-adjustment of lta for long-duration events. 
#' 
#' @param on \code{Numeric} value, threshold value for event onset.
#' 
#' @param off \code{Numeric} value, threshold value for event end.
#' 
#' @return \code{data frame}, detected events (ID, start time, duration in 
#' seconds, STA-LTA vaue).
#' 
#' @author Michael Dietze
#' 
#' @keywords eseis
#' 
#' @examples
#' 
#' ## load example data
#' data(rockfall)
#' 
#' ## filter signal
#' rockfall_f <- signal_filter(data = rockfall_eseis,
#'                             f = c(1, 90), 
#'                             p = 0.05)
#'                     
#' ## calculate signal envelope
#' rockfall_e <- signal_envelope(data = rockfall_f)
#' 
#' ## pick earthquake and rockfall event
#' p <- pick_stalta(data = rockfall_e,
#'                  sta = 100, 
#'                  lta = 18000, 
#'                  freeze = TRUE, 
#'                  on = 5, 
#'                  off = 3)
#'                  
#' p$picks
#'                      
#' @export pick_stalta
pick_stalta <- function(
  data,
  time,
  dt,
  sta,
  lta,
  freeze = FALSE,
  on,
  off
) {
  
  ## check/store presence of time vector
  if(missing(time) == TRUE) {
    
    time_in <- NULL
  } else {
    
    time_in <- time
  }
  
  if(missing(dt) == TRUE) {
    
    dt <- NULL
  } else {
    
    dt <- dt
  }

  ## check data structure
  if(class(data)[1] == "list") {
    
    ## apply function to list
    data_out <- lapply(X = data, 
                       FUN = eseis::pick_stalta,
                       time = time_in,
                       dt = dt,
                       sta = sta,
                       lta = lta,
                       freeze = freeze,
                       on = on,
                       off = off)
    
    ## return output
    return(data_out)
    
  } else {
    
    ## get start time
    eseis_t_0 <- Sys.time()
    
    ## collect function arguments
    eseis_arguments <- list(data = "",
                            time = time_in,
                            dt = dt,
                            sta = sta,
                            lta = lta,
                            freeze = freeze,
                            on = on,
                            off = off)
    
    ## homogenise data structure
    if(class(data)[1] == "eseis") {
      
      ## set eseis flag
      eseis_class <- TRUE
      
      ## store initial object
      eseis_data <- data
      
      ## extract signal vector
      data <- eseis_data$signal
      
      ## extract sampling period
      dt <- eseis_data$meta$dt
      
      ## generate time vector
      time <- seq(from = eseis_data$meta$starttime, 
                  by = dt, 
                  length.out = eseis_data$meta$n)
      
    } else {
      
      ## set eseis flag
      eseis_class <- FALSE
      
      ## assign dt and time
      if(is.null(dt) == TRUE) {
        
        if(is.null(time) == TRUE) {
          
          dt <- 1 / 200
          
          warning("No dt provided. Set to 1 / 200 s by default!")
        } else {
          
          if(length(time) > 1000) {
            
            time_dt <- time[1:1000]
            
          } else {
            
            time_dt <- time
            
          }
          
          dt_estimate <- mean(x = diff(x = as.numeric(time)), 
                              na.rm = TRUE)
          print(paste("dt estimated from time vector (", 
                      signif(x = dt_estimate),
                      " s)",
                      sep = ""))
          
          dt <- dt_estimate
        }
        
        ## handle missing time vector
        if(is.null(time) == TRUE) {
          
          time <- seq(from = 
                        as.POSIXct(x = strptime(x = "0000-01-01 00:00:00",
                                                format = "%Y-%m-%d %H:%M:%S", 
                                                tz = "UTC"), tz = "UTC"),
                      by = dt,
                      length.out = length(data))
          
          print("No or non-POSIXct time data provided. Default data created!")
          
        }
        
      }
    }
 
    ## calculate sta vector
    data_sta <- caTools::runmean(x = data, 
                                 k = sta, 
                                 alg = "fast", 
                                 endrule = "NA",
                                 align = "right")
    
    ## calculate lta vector
    data_lta <- caTools::runmean(x = data, 
                                 k = lta, 
                                 alg = "fast", 
                                 endrule = "NA",
                                 align = "right")
    
    ## calculate stalta vector
    data_stalta <- data_sta / data_lta
    
    ## create output data set
    event <- rep(x = 0, 
                 times = length(data))
    
    if(freeze == TRUE) {
      
      event <- .stalta_event_freeze(event_length = length(event),
                                    data_sta = data_sta,
                                    data_lta = data_lta,
                                    on = on,
                                    off = off)
    } else {
      
      ## calculate ratio
      ratio <- data_sta / data_lta
      ratio[is.na(ratio)] <- 0
      
      ## assign event trigger
      event <- .stalta_event_nofreeze(event_length = length(event),
                                      ratio = ratio,
                                      on = on,
                                      off = off)
    }

    ## calculate event state switches
    event_diff <- diff(x = event)

    ## asign on and off states
    event_on <- time[event_diff == 1]
    event_off <- time[event_diff == -1]

    ## calculate event duration
    event_duration <- as.numeric(difftime(time1 = event_off, 
                                          time2 = event_on, 
                                          units = "sec"))
    
    # ## get maximum sta-lta value
    # stalta_max <- rep(NA, length = length(event_on))
    # for(i in 1:length(stalta_max)) {
    #   stalta_max[i] <- max((data_sta / data_lta)[time >= event_on[i] & 
    #                                time <= event_off[i]])
    # }

    ## get sta-lta value at event onset
    event_stalta <- data_stalta[event_diff == 1]
    
    ## create ID vector
    ID <- seq(from = 1, 
              along.with = event_duration)
    
    ## build output data set
    if(length(ID) > 0) {
      
      data_out <- data.frame(ID = ID,
                             start = event_on,
                             duration = event_duration,
                             stalta = event_stalta)
      
    } else {
      
      data_out <- data.frame(ID = NA,
                             start = NA,
                             duration = NA,
                             stalta = NA)[-1,]
    }
    
    ## optionally rebuild eseis object
    if(eseis_class == TRUE) {
      
      ## assign aggregated signal vector
      eseis_data <- list(picks = data_out,
                         history = eseis_data$history)
      
      ## calculate function call duration
      eseis_duration <- as.numeric(difftime(time1 = Sys.time(), 
                                            time2 = eseis_t_0, 
                                            units = "secs"))
      
      ## update object history
      eseis_data$history[[length(eseis_data$history) + 1]] <- 
        list(time = Sys.time(),
             call = "pick_stalta()",
             arguments = eseis_arguments,
             duration = eseis_duration)
      names(eseis_data$history)[length(eseis_data$history)] <- 
        as.character(length(eseis_data$history))
      
      ## set S3 class name
      class(eseis_data)[1] <- "eseis"
      
      ## assign eseis object to output data set
      data_out <- eseis_data
    }

    ## return output
    return(data_out)
  }
}