R/plot_signal.R

In eseis: Environmental Seismology Toolbox

#' Plot a seismic signal
#' 
#' This function plots a line graph of a seismic signal. To avoid long plot 
#' preparation times the signal is reduced to a given number of points.
#' 
#' The \code{format} argument is based on hints provided by Sebastian 
#' Kreutzer and Christoph Burow. It allows plotting time axis units in 
#' user defined formats. The time format must be provided as character string 
#' using the POSIX standard (see documentation of \code{strptime} for a list 
#' of available keywords), e.g., "%H:%M:%S" (hour:minute:second),
#' "%Y-%m-%d" (year-month-day).
#' 
#' @param data \code{eseis} object or \code{numeric} vector, data set to 
#' be plotted.
#' 
#' @param time \code{POSIXct} vector, corresponding time vector.
#' 
#' @param n \code{Numeric} value, number of values to which the dataset 
#' is reduced. Default is \code{10000}.
#' 
#' @param \dots Further arguments passed to the plot function.
#' 
#' @return A line plot of a seismic wave form.
#' 
#' @author Michael Dietze
#' 
#' @keywords eseis
#' 
#' @examples
#' 
#' ## load example data set
#' data(rockfall)
#' 
#' ## plot data set straightforward
#' plot_signal(data = rockfall_eseis)
#' 
#' ## plot data set with lower resolution
#' plot_signal(data = rockfall_eseis, n = 100)
#' 
#' ## plot data set but not as an eseis object
#' plot_signal(data = rockfall_z, time = rockfall_t)
#' 
#' ## load earthquake data set
#' data(earthquake)
#' 
#' ## plot all three components (after changing plot options)
#' pars <- par(no.readonly = TRUE)
#' par(mfcol = c(3, 1))
#' 
#' plt <- lapply(s, plot_signal, t = t)
#' 
#' par(pars)
#' 
#' @export plot_signal
#' 
plot_signal <- function(
    
  data,
  time,
  n = 10000,
  ...
) {
  
  ## check if time vector is provided
  if(missing(time) == TRUE) {
    
    time <- NA
  }
  
  ## extract additional plot arguments
  args <- list(...)
  
  ## check/set plot arguments
  if ("type" %in% names(args)) {
    type <- args$type
  }
  else {
    
    type = "l"
  }
  
  if ("xlab" %in% names(args)) {
    xlab <- args$xlab
  }
  else {
    
    xlab = "Time"
  }
  
  if ("ylab" %in% names(args)) {
    ylab <- args$ylab
  }
  else {
    
    ylab = "Amplitude"
  }
  
  if ("axes" %in% names(args)) {
    axes <- args$axes
  }
  else {
    
    axes <- TRUE
  }
  
  if ("format" %in% names(args)) {
    format <- args$format
  }
  else {
    
    format <- ""
  }
  
  ## remove keywords
  keywords <- c("format", "axes", "xlab", "ylab")
  args <- args[!names(args)%in%keywords]
  
  if(class(data)[1] == "eseis") {
    
    ## put eseis object into list
    data <- list(data = data)
  } else if(class(data)[1] == "numeric") {
    
    ## put data vectors into list
    data <- list(data)
  }
  
  ## process list elements
  for(nn in 1:length(data)) {
    
    ## check if input object is of class eseis
    if(class(data[[nn]])[1] == "eseis") {
      
      ## set eseis flag
      eseis_class <- TRUE
      
      ## store initial object
      eseis_data <- data[[nn]]
      
      ## extract signal vector
      data[[nn]] <- eseis_data$signal
      
      ## create time vector
      time <- seq(from = eseis_data$meta$starttime, 
                  by = eseis_data$meta$dt, 
                  length.out = eseis_data$meta$n)
      
    } else {
      
      ## set eseis flag
      eseis_class <- FALSE
    }
    
    ## account for data sets smaller than n
    if(length(data[[nn]]) < n) {
      
      n <- length(data[[nn]])
    }
    
    ## get NA values to pad data set
    n_pad <- ceiling(x = length(data[[nn]]) / n) * n - length(data[[nn]])
    
    ## pad data set
    s_pad <- c(data[[nn]], rep(NA, 
                         times = n_pad))
    
    t_pad <- c(time, rep(NA, 
                         times = n_pad))
    
    ## convert signal vector to matrix
    S <- matrix(data = s_pad, 
                ncol = n)
    
    ## calculate columnwise min and max
    S_min <- matrixStats::colMins(x = S, 
                                  na.rm = TRUE)
    
    S_max <- matrixStats::colMaxs(x = S, 
                                  na.rm = TRUE)
    
    ## convert min and max alternatingly to vector
    s_plot <- as.numeric(rbind(S_min, S_max))
    
    ## get time vector subset
    i <- seq(from = 1, 
             to = length(t_pad), 
             by = nrow(S))
    
    ## make pairs of time vector subsets
    t_plot <- rep(t_pad[i], each = 2)
    
    
    ## generate plot
    do.call(what = graphics::plot, 
            args = c(list(t_plot, 
                          s_plot, 
                          type = type, 
                          xlab = xlab,
                          ylab = ylab,
                          axes = FALSE), 
                     args))
    
    ## optionally add axes
    if(axes == TRUE) {
      
      graphics::axis.POSIXct(side = 1, 
                             x = t_plot, 
                             format = format)
      
      graphics::axis(side = 2)
      
      graphics::box(which = "plot")
    }
  }
}