R/Class-Trace.R

In IRISSeismic: Classes and Methods for Seismic Data Analysis

##
##    S4 classes for handling seismic traces
##
##    Copyright (C) 2012  Mazama Science, Inc.
##    by Jonathan Callahan, jonathan@mazamascience.com
##
##    This program is free software; you can redistribute it and/or modify
##    it under the terms of the GNU General Public License as published by
##    the Free Software Foundation; either version 2 of the License, or
##    (at your option) any later version.
##
##    This program is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.
##
##    You should have received a copy of the GNU General Public License
##    along with this program; if not, write to the Free Software
##    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA


################################################################################
# R classes for a Trace obj.
#
# This is a port of some of the functionality found in obspy.core.trace.py
#
#   https://docs.obspy.org/packages/autogen/obspy.core.trace.Trace.html
################################################################################

################################################################################
# Class TraceHeader
#
# A container for additional header information for each Trace object.
#
# A TraceHeader object contains all header information (aka metadata)
# associated with a Trace object.
#
# Documentation for attributes copied verbatim from obspy.core.trace.py:
#
#    ``sampling_rate`` : float, optional
#        Sampling rate in hertz (default value is 1.0).
#    ``delta`` : float, optional
#        Sample distance in seconds (default value is 1.0).
#    ``calib`` : float, optional
#        Calibration factor (default value is 1.0).
#    ``npts`` : int, optional
#        Number of sample points (default value is 0, which implies that no data
#        is present).
#    ``network`` : string, optional
#        Network code (default is an empty string).
#    ``location`` : string, optional
#        Location code (default is an empty string).
#    ``station`` : string, optional
#        Station code (default is an empty string).
#    ``channel`` : string, optional
#        Channel code (default is an empty string).
#    ``starttime`` : :class:`~obspy.core.utcdatetime.UTCDateTime`, optional
#        Date and time of the first data sample given in UTC (default value is
#        "1970-01-01T00:00:00.0Z").
#    ``endtime`` : :class:`~obspy.core.utcdatetime.UTCDateTime`, optional
#        Date and time of the last data sample given in UTC
#        (default value is "1970-01-01T00:00:00.0Z").
#
# It can also contain some optional metadata, not included in miniseed headers:
#
#     ``latitude``  : numeric, optional
#     ``longitude`` : numeric, optional
#     ``elevation`` : numeric, optional
#     ``depth``     : numeric, optional
#     ``azimuth``   : numeric, optional
#     ``dip``       : numeric, optional
#
################################################################################

setClass(
  "TraceHeader",
  # typed slots (aka attributes) for class TraceHeader
  representation(
    sampling_rate = "numeric",
    delta = "numeric",
    calib = "numeric",
    npts = "integer",
    network = "character",
    location = "character",
    station = "character",
    channel = "character",
    quality = "character",
    starttime = "POSIXct",
    endtime = "POSIXct",
    latitude = "numeric",
    longitude = "numeric",
    elevation = "numeric",
    depth = "numeric",
    azimuth = "numeric",
    dip = "numeric",
    processing = "character"
  ),
  # default values for slots
  prototype(
    sampling_rate = 1.0,
    delta = 1.0,
    calib = 1.0,
    npts = as.integer(0),
    network = "",
    location = "",
    station = "",
    channel = "",
    quality = "X",
    starttime = as.POSIXct("1900-01-01T00:00:00", format = "%Y-%m-%dT%H:%M:%OS6", tz = "GMT"),
    endtime = as.POSIXct("1900-01-01T00:00:00", format = "%Y-%m-%dT%H:%M:%OS6", tz = "GMT"),
    latitude = as.numeric(NA),
    longitude = as.numeric(NA),
    elevation = as.numeric(NA),
    depth = as.numeric(NA),
    azimuth = as.numeric(NA),
    dip = as.numeric(NA)
  )
)

################################################################################
#
# A TraceHeader can be initialized with a headerLine as returned by the IRIS DMC timeseries webservice:
#
# TIMESERIES LD_POTS__HHZ_M, 351 samples, 100.503 sps, 2012-01-29T00:00:00.006000, SLIST, INTEGER, COUNTS
#
# or with named elements from a headerList.  Items in the headerList will override items determined
# by parsing the headerLine.
################################################################################

# TODO:  I am currently using "Trace" to mean what libmseed calls a "Trace Segment".
# TODO:  Each libmseed "Trace" has both "numsamples" and "samplecnt".
# TODO:  Need to figure out what the difference is and when to bail with zero samples in the SNCL.
# TODO:  Zero samples in a "Trace Segment" is not necessarily a deal breaker.

setMethod(
  "initialize", "TraceHeader",
  function(.Object,
           headerList = list(),
           headerLine = NULL,
           ...) {
    # Begin by parsing any headerLine that was passed in.
    if (!is.null(headerLine)) {
      headerValues <- strsplit(headerLine, ", ")[[1]]
      sncl <- strsplit(headerValues[1], " ")[[1]][2]
      snclValues <- strsplit(sncl, "_")[[1]]
      .Object@network <- snclValues[1]
      .Object@station <- snclValues[2]
      .Object@location <- snclValues[3]
      .Object@channel <- snclValues[4]
      .Object@npts <- as.integer(strsplit(headerValues[2], " ")[[1]][1])
      .Object@sampling_rate <- as.numeric(strsplit(headerValues[3], " ")[[1]][1])
      .Object@starttime <- as.POSIXct(headerValues[4], format = "%Y-%m-%dT%H:%M:%OS6", tz = "GMT") # header uses 'T' format
      .Object@latitude <- as.numeric(NA)
      .Object@longitude <- as.numeric(NA)
      .Object@elevation <- as.numeric(NA)
      .Object@depth <- as.numeric(NA)
      .Object@azimuth <- as.numeric(NA)
      .Object@dip <- as.numeric(NA)
    }

    # override with anything found in the headerList
    for (key in names(headerList)) {
      slot(.Object, key) <- headerList[[key]]
    }

    # guarantee npts and sampling_rate are non-zero numbers
    if (is.na(.Object@npts)) {
      err_msg <- paste("initialize.TraceHeader: npts=", .Object@npts, sep = "")
      if (!is.null(headerLine)) {
        err_msg <- paste(err_msg, "\nheaderLine =\"", headerLine, "\"", sep = "")
      }
      stop(err_msg)
    }
    if (is.na(.Object@sampling_rate) || (.Object@sampling_rate == 0)) {
      err_msg <- paste("initialize.TraceHeader: sampling_rate=", .Object@sampling_rate, sep = "")
      if (!is.null(headerLine)) {
        err_msg <- paste(err_msg, "\nheaderLine =\"", headerLine, "\"", sep = "")
      }
      stop(err_msg)
    }

    # calculate derived values
    delta <- 1.0 / .Object@sampling_rate
    if (!is.finite(delta)) {
      delta <- 0
    }
    .Object@delta <- delta

    if (.Object@npts == 0) {
      delta <- 0
    } else {
      delta <- (.Object@npts - 1) / .Object@sampling_rate
      if (!is.finite(delta)) {
        delta <- 0
      }
    }
    if (is.na(.Object@endtime)) {
      .Object@endtime <- .Object@starttime + delta
    }

    # return the newly created .Object
    return(.Object)
  }
)


# show method is called by show() and print() ----------------------------------

show.TraceHeader <- function(object) {
  cat("Seismic Trace TraceHeader \n")
  cat(" Network:       ", object@network, "\n")
  cat(" Station:       ", object@station, "\n")
  cat(" Location:      ", object@location, "\n")
  cat(" Channel:       ", object@channel, "\n")
  cat(" Quality:       ", object@quality, "\n")
  cat(" calib:         ", object@calib, "\n")
  cat(" npts:          ", object@npts, "\n")
  cat(" sampling rate: ", object@sampling_rate, "\n")
  cat(" delta:         ", object@delta, "\n")
  cat(" starttime:     ", format(object@starttime), "\n")
  cat(" endtime:       ", format(object@endtime), "\n")
  cat(" latitude:      ", format(object@latitude), "\n")
  cat(" longitude:     ", format(object@longitude), "\n")
  cat(" elevation:     ", format(object@elevation), "\n")
  cat(" depth:         ", format(object@depth), "\n")
  cat(" azimuth:       ", format(object@azimuth), "\n")
  cat(" dip:           ", format(object@dip), "\n")
  cat(" processing:    ", paste(object@processing, collapse = "; "), "\n")
}
# NOTE:  method signature must match generic signature for 'show' with argument: 'object'
setMethod("show", signature(object = "TraceHeader"), function(object) show.TraceHeader(object))


# as.headerLine prints TraceHeader out as a single line ------------------------
#
# TIMESERIES LD_POTS__HHZ_M, 351 samples, 100.503 sps, 2012-01-29T00:00:00.006000, SLIST, INTEGER, COUNTS

if (!isGeneric("as.headerLine")) {
  setGeneric("as.headerLine", function(obj) {
    standardGeneric("as.headerLine")
  })
}
as.headerLine.TraceHeader <- function(obj) {
  cat("TIMESERIES ", obj@network, "_", obj@station, "_", obj@location, "_", obj@channel, ", ", sep = "")
  cat(obj@npts, " samples, ", obj@sampling_rate, " sps, ", obj@starttime, ", SLIST, INTEGER, COUNTS", sep = "")
}
setMethod("as.headerLine", signature(obj = "TraceHeader"), function(obj) as.headerLine.TraceHeader(obj))

# TODO:  Setter and getter methods could be written as described in:
# TODO:    www.stat.auckland.ac.nz/S-Workshop/Gentleman/S4Objects.pdf


################################################################################
# Class Trace
#
# An object containing data of a continuous series, such as a seismic trace.
#
# Documentation for attributes copied verbatim from obspy.core.trace.py
#
#    :type data: :class:`~numpy.ndarray` or :class:`~numpy.ma.MaskedArray`
#    :param data: Array of data samples
#    :type header: dict or :class:`~obspy.core.trace.Stats`
#    :param header: Dictionary containing header fields
#
#    :var id: A SEED compatible identifier of the trace.
#    :var stats: A container :class:`~obspy.core.trace.Stats` for additional
#        header information of the trace.
#    :var Sensor:
#    :var InstrumentSensitivity:
#    :var SensitivityFrequency:
#    :var InputUnits:
#    :var data: Data samples in a :class:`~numpy.ndarray` or
#        :class:`~numpy.ma.MaskedArray`
#
################################################################################

setClass(
  "Trace",
  # typed slots (aka attributes) for class Trace
  representation(
    id = "character",
    stats = "TraceHeader",
    Sensor = "character",
    InstrumentSensitivity = "numeric",
    SensitivityFrequency = "numeric",
    InputUnits = "character",
    data = "numeric"
  ),
  # default values for slots
  prototype(
    id = "",
    stats = new("TraceHeader"),
    Sensor = "",
    # InstrumentSensitivity = 1.0,
    InstrumentSensitivity = as.numeric(NA), # if no information is available, the default should not be 1.0
    # SensitivityFrequency = 1.0,
    SensitivityFrequency = as.numeric(NA), # if no information is available, the default should not be 1.0
    InputUnits = "",
    data = c(0)
  )
)

# initialize method
setMethod(
  "initialize", "Trace",
  function(.Object,
           id = "",
           stats = new("TraceHeader"),
           Sensor = "",
           # InstrumentSensitivity = 1.0,
           InstrumentSensitivity = as.numeric(NA),
           # SensitivityFrequency = 1.0,
           SensitivityFrequency = as.numeric(NA),
           InputUnits = "",
           data = c(0),
           ...) {
    .Object@id <- id
    .Object@stats <- stats
    .Object@Sensor <- Sensor
    .Object@InstrumentSensitivity <- InstrumentSensitivity
    .Object@SensitivityFrequency <- SensitivityFrequency
    .Object@InputUnits <- InputUnits
    .Object@data <- data

    if (.Object@id == "") {
      stats <- .Object@stats
      .Object@id <- paste(stats@network, stats@station, stats@location, stats@channel, stats@quality, sep = ".")
    }

    return(.Object)
  }
)


################################################################################
# Simple methods working on the data from a Trace object
################################################################################

# NOTE:  The method signatures must match those already defined by R core functions
# NOTE:  and R standard generics.  Use "?function" to see the documentation or just
# NOTE:  "function" to see the implementation and signature.
# NOTE:
# NOTE:  When creating a new method to extend a method that has "...", always
# NOTE:  provide the S4 object as the first argument, then the "..." and then
# NOTE:  any final arguments.


# as.vector --------------------------------------------------------------------

as.vector.Trace <- function(x, mode = "any") {
  return(x@data)
}
setMethod("as.vector", signature(x = "Trace"), function(x) as.vector.Trace(x))


# isDC -------------------------------------------------------------------------

if (!isGeneric("isDC")) {
  setGeneric("isDC", function(x) {
    standardGeneric("isDC")
  })
}
isDC.Trace <- function(x) {
  x <- DDT(x, demean = TRUE, detrend = FALSE, taper = 0)
  return(all(x@data == 0, na.rm = TRUE))
}
setMethod("isDC", signature(x = "Trace"), function(x) isDC.Trace(x))


# Length -----------------------------------------------------------------------

length.Trace <- function(x) {
  return(length(x@data))
}
setMethod("length", signature(x = "Trace"), function(x) length.Trace(x))


# Maximum ----------------------------------------------------------------------

max.Trace <- function(x, ..., na.rm = TRUE) {
  return(max(x@data, na.rm = na.rm))
}
setMethod("max", signature(x = "Trace"), function(x, ...) max.Trace(x, ...))


# Mean -------------------------------------------------------------------------

mean.Trace <- function(x, ...) {
  return(mean(x@data, ...))
}
setMethod("mean", signature(x = "Trace"), function(x, ...) mean.Trace(x, ...))
# TODO:  How to get mean to use na.rm=TRUE as the default?


# Median -----------------------------------------------------------------------

median.Trace <- function(x, na.rm) {
  return(median(x@data, na.rm = na.rm))
}
# NOTE:  method signature must match generic signature for function 'median' with arguments: 'x', 'na.rm'
setMethod("median", signature(x = "Trace", na.rm = "logical"), function(x, na.rm) median.Trace(x, na.rm = na.rm))
setMethod("median", signature(x = "Trace", na.rm = "missing"), function(x, na.rm) median.Trace(x, na.rm = TRUE))


# Minimum ----------------------------------------------------------------------

min.Trace <- function(x, ..., na.rm = TRUE) {
  return(min(x@data, na.rm = na.rm))
}
setMethod("min", signature(x = "Trace"), function(x, ...) min.Trace(x, ...))


# Standard Deviation -----------------------------------------------------------

sd.Trace <- function(x, na.rm) {
  return(sd(x@data, na.rm = na.rm))
}
# NOTE:  method signature must match generic signature for for function 'sd' with arguments: 'x', 'na.rm'
setMethod("sd", signature(x = "Trace", na.rm = "logical"), function(x, na.rm) sd.Trace(x, na.rm = na.rm))
setMethod("sd", signature(x = "Trace", na.rm = "missing"), function(x, na.rm) sd.Trace(x, na.rm = TRUE))


# Root Mean Square ----------------------------------------------------

if (!isGeneric("rms")) {
  setGeneric("rms", function(x, na.rm) {
    standardGeneric("rms")
  })
}
rms.Trace <- function(x, na.rm) {
  return(sqrt(mean((x@data)^2, na.rm = na.rm)))
}
setMethod("rms", signature(x = "Trace", na.rm = "logical"), function(x, na.rm) rms.Trace(x, na.rm = na.rm))
setMethod("rms", signature(x = "Trace", na.rm = "missing"), function(x, na.rm) rms.Trace(x, na.rm = TRUE))


# Root Mean Square Variance ----------------------------------------------------

if (!isGeneric("rmsVariance")) {
  setGeneric("rmsVariance", function(x, na.rm) {
    standardGeneric("rmsVariance")
  })
}
rmsVariance.Trace <- function(x, na.rm) {
  mean <- mean(x, na.rm = na.rm)
  if (na.rm) {
    n <- length(x@data) - length(which(is.na(x@data)))
  } else {
    n <- length(x)
  }
  return(sqrt(sum((x@data - mean)^2, na.rm = na.rm) / (n)))
}
setMethod("rmsVariance", signature(x = "Trace", na.rm = "logical"), function(x, na.rm) rmsVariance.Trace(x, na.rm = na.rm))
setMethod("rmsVariance", signature(x = "Trace", na.rm = "missing"), function(x, na.rm) rmsVariance.Trace(x, na.rm = TRUE))


################################################################################
# Methods that return a new, modified Trace
################################################################################

# Multiply by a constant -------------------------------------------------

if (!isGeneric("multiplyBy")) {
  setGeneric("multiplyBy", function(x, y) {
    standardGeneric("multiplyBy")
  })
}

multiplyBy.Trace <- function(x, y) {
  id <- x@id
  stats <- x@stats
  Sensor <- x@Sensor
  InstrumentSensitivity <- x@InstrumentSensitivity
  SensitivityFrequency <- x@SensitivityFrequency
  InputUnits <- x@InputUnits
  data <- x@data * y
  stats@processing <- append(stats@processing, paste("multiply by", y))

  return(new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data))
}

setMethod("multiplyBy", signature("Trace", y = "numeric"), function(x, y) multiplyBy.Trace(x, y = y))


# Apply typical seismic DDT ----------------------------

if (!isGeneric("DDT")) {
  setGeneric("DDT", function(x, demean, detrend, taper) {
    standardGeneric("DDT")
  })
}

DDT.Trace <- function(x, demean, detrend, taper) {
  # TODO:  When traces have processingInfo we can check the processingInfo
  # TODO:  and skip any steps that have already been done.

  id <- x@id
  stats <- x@stats
  Sensor <- x@Sensor
  InstrumentSensitivity <- x@InstrumentSensitivity
  SensitivityFrequency <- x@SensitivityFrequency
  InputUnits <- x@InputUnits

  # NOTE:  Use the pracma::detrend() function.
  # NOTE:  Linear detrending also removes the mean.

  # demean only
  if (demean && !detrend) {
    data <- x@data - mean(x@data, na.rm = TRUE)
    stats@processing <- append(stats@processing, "demean")
  } else {
    data <- x@data
  }

  # detrend
  if (detrend) {
    data <- as.numeric(pracma::detrend(x@data, tt = "linear"))
    stats@processing <- append(stats@processing, "detrend")
  }

  # cosine taper
  if (taper > 0) {
    data <- stats::spec.taper(data, p = taper)
    stats@processing <- append(stats@processing, "cosine taper")
  }


  return(new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data))
}

# All parameters specified
setMethod(
  "DDT", signature("Trace", demean = "logical", detrend = "logical", taper = "numeric"),
  function(x, demean, detrend, taper) DDT.Trace(x, demean, detrend, taper)
)
# No parameters
setMethod(
  "DDT", signature("Trace", demean = "missing", detrend = "missing", taper = "missing"),
  function(x, demean, detrend, taper) DDT.Trace(x, TRUE, TRUE, 0.1)
)


# Apply Butterworth filters ------------

if (!isGeneric("butterworth")) {
  setGeneric("butterworth", function(x, n, low, high, type) {
    standardGeneric("butterworth")
  })
}

butterworth.Trace <- function(x, n, low, high, type) {
  # TODO:  When traces have processingInfo we can check the processingInfo
  # TODO:  and skip any steps that have already been done.

  data <- x@data
  id <- x@id
  stats <- x@stats
  Sensor <- x@Sensor
  InstrumentSensitivity <- x@InstrumentSensitivity
  SensitivityFrequency <- x@SensitivityFrequency
  InputUnits <- x@InputUnits

  # If neither demean nor detrend has ever been applied, demean the data
  if (!any(stringr::str_detect(x@stats@processing, "demean")) &&
    !any(stringr::str_detect(x@stats@processing, "detrend"))) {
    data <- data - mean(data, na.rm = TRUE)
    stats@processing <- append(stats@processing, "demean")
  }


  #   filt <- signal::butter(2,c(0.02/(trZ@stats@sampling_rate/2),0.04/(trZ@stats@sampling_rate/2)),type="pass")
  #   tsZ <- stats::ts(trZ@data,frequency=trZ@stats@sampling_rate)
  #   ts2 <- stats::ts(tr2@data,frequency=tr2@stats@sampling_rate)
  #   ts1 <- stats::ts(tr1@data,frequency=tr1@stats@sampling_rate)
  #
  #   trZ_f <- trZ; trZ_f@data <- as.numeric(signal::filter(filt,tsZ))
  #

  # Create normalization factor
  norm <- stats@sampling_rate / 2

  # Create Butterworth filter
  if (type == "low") {
    bf <- signal::butter(n, low / norm, type = "low")
  } else if (type == "high") {
    bf <- signal::butter(n, high / norm, type = "high")
  } else {
    bf <- signal::butter(n, c(low / norm, high / norm), type = type)
  }

  # Apply filter
  ts <- stats::ts(data, frequency = stats@sampling_rate)
  ts_filtered <- signal::filter(bf, ts)
  data <- as.numeric(ts_filtered)
  stats@processing <- append(stats@processing, paste0("Butterworth filtered with (", n, ",", low, ",", high, ",", type, ")"))

  return(new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data))
}

# All parameters specified
setMethod(
  "butterworth", signature("Trace", n = "numeric", low = "numeric", high = "numeric", type = "character"),
  function(x, n, low, high, type) butterworth.Trace(x, n, low, high, type)
)
# type missing -- band pass
setMethod(
  "butterworth", signature("Trace", n = "numeric", low = "numeric", high = "numeric", type = "missing"),
  function(x, n, low, high, type) butterworth.Trace(x, n, low, high, "pass")
)
# hi and type missing -- low pass
setMethod(
  "butterworth", signature("Trace", n = "numeric", low = "numeric", high = "missing", type = "missing"),
  function(x, n, low, high, type) butterworth.Trace(x, n, low, NULL, "high")
)
# lo and type missing -- high pass
setMethod(
  "butterworth", signature("Trace", n = "numeric", low = "missing", high = "numeric", type = "missing"),
  function(x, n, low, high, type) butterworth.Trace(x, n, NULL, high, "low")
)


# Slice out a portion of a Trace ---------------------------

# NOTE:  Unlike the ObsPy method, this method does no padding.  The returned
# NOTE:  Trace will always be a subset of the original Trace.
# NOTE:  Slicing is to the nearest sample as of version 1.4.6

if (!isGeneric("slice")) {
  setGeneric("slice", function(x, starttime, endtime) {
    standardGeneric("slice")
  })
}
slice.Trace <- function(x, starttime, endtime) {
  # id and stats will be part of the returned Trace
  id <- x@id
  stats <- x@stats
  Sensor <- x@Sensor
  InstrumentSensitivity <- x@InstrumentSensitivity
  SensitivityFrequency <- x@SensitivityFrequency
  InputUnits <- x@InputUnits
  data <- x@data

  sliced <- FALSE

  start_index <- 1
  end_index <- length(data)

  # Sanity check
  if (starttime >= endtime) {
    stop(paste("slice.Trace: requested starttime \"", starttime, "\" >= requested endtime \"", endtime, "\""))
  }
  if (starttime >= stats@endtime) {
    stop(paste("slice.Trace: requested starttime \"", starttime, "\" >= Trace endtime \"", stats@endtime, "\""))
  }
  if (endtime <= stats@starttime) {
    stop(paste("slice.Trace: requested endtime \"", endtime, "\" <= Trace starttime \"", stats@starttime, "\""))
  }

  # Modify the startime and startIndex if needed
  if (starttime > stats@starttime) {
    # NOTE:  Use difftime() instead of just subtracting to guarantee that units are "secs"
    offset_secs <- as.numeric(difftime(starttime, stats@starttime, units = "secs"))
    offset_secs <- round(offset_secs, digits = 6)
    start_index <- start_index + as.integer(floor(offset_secs * stats@sampling_rate))
    stats@starttime <- stats@starttime + offset_secs
    sliced <- TRUE
  } # else leave the start as is

  # Modify the endtime and end_index if needed
  if (endtime < stats@endtime) {
    # NOTE:  Use difftime() instead of just subtracting to guarantee that units are "secs"
    offset_secs <- as.numeric(difftime(stats@endtime, endtime, units = "secs"))
    offset_secs <- round(offset_secs, digits = 6)
    end_index <- end_index - as.integer(floor(offset_secs * stats@sampling_rate))
    stats@endtime <- stats@endtime - offset_secs
    sliced <- TRUE
  } # else leave the end as is

  if (sliced) {
    stats@npts <- as.integer(end_index - start_index + 1)
    stats@processing <- append(stats@processing, "slice")
    return(new("Trace", id, stats, Sensor, InstrumentSensitivity, SensitivityFrequency, InputUnits, data = data[start_index:end_index]))
  } else {
    return(x)
  }
}
setMethod(
  "slice", signature(x = "Trace", starttime = "POSIXct", endtime = "POSIXct"),
  function(x, starttime, endtime) slice.Trace(x, starttime = starttime, endtime = endtime)
)

################################################################################
# Method to return a new Trace assigned a specified quality code.
# This method modifies the @id slot by replacing the last character of the
# NET.STA.LOC.CHA.QUAL id string with the requested quality code value. It also
# replaces the @stats@quality slot with the requested quality code.
# The first input parameter is the Trace object. The second input
# parameter is the requested quality code.
# Quality codes must be one of the characters "R","D","M","Q".
################################################################################

if (!isGeneric("assignQuality")) {
  setGeneric(
    "assignQuality",
    function(x, quality) standardGeneric("assignQuality")
  )
}

assignQuality.Trace <- function(x, quality) {
  if (!quality %in% c("R", "D", "Q", "M")) {
    stop("Requested quality code must be one of  \"R\", \"D\", \"Q\", \"M\"")
  }
  x@id <- sub(".$", quality, x@id)
  x@stats@quality <- quality
  return(x)
}

setMethod(
  "assignQuality",
  signature(x = "Trace", quality = "character"),
  function(x, quality) assignQuality.Trace(x, quality = quality)
)

################################################################################
# Spectral Methods
################################################################################

# Hilbert transform -------------------------

if (!isGeneric("hilbert")) {
  setGeneric("hilbert", function(x) {
    standardGeneric("hilbert")
  })
}

# Code copied from seewave package hilbert() method
#   https://cran.r-project.org/web/packages/seewave/index.html
hilbert.Trace <- function(x) {
  # Always detrend, demean, cosine taper the data
  x <- DDT(x)
  hfft <- hilbertFFT(x@data)
  # The Hilbert tansform is the Imaginary part of the transformed data
  x@data <- Im(hfft)
  x@stats@processing <- append(x@stats@processing, "Hilbert transform")

  return(x)
}

setMethod(
  "hilbert", signature(x = "Trace"),
  function(x) hilbert.Trace(x)
)



# Envelope using Hilbert transform -------------------------

if (!isGeneric("envelope")) {
  setGeneric("envelope", function(x) {
    standardGeneric("envelope")
  })
}

# Code copied from seewave package env() and hilbert() methods
#   https://cran.r-project.org/web/packages/seewave/index.html
envelope.Trace <- function(x) {
  # Always detrend, demean, cosine taper the data
  x <- DDT(x)
  hfft <- hilbertFFT(x@data)
  # The envelope is the Modulus of the transformed data
  x@data <- Mod(hfft)
  x@stats@processing <- append(x@stats@processing, "envelope")

  return(x)
}

setMethod(
  "envelope", signature(x = "Trace"),
  function(x) envelope.Trace(x)
)


################################################################################
# Methods for event picking
################################################################################

# Multi-purpose SLA/TLA trigger algorithm --------------------------------------

# NOTE:  https://en.wikipedia.org/wiki/First_break_picking
# NOTE:
# NOTE:  "A Comparison of Select Trigger Algorithms for Automated Global Seismic Phase and Event Detection"
# NOTE:  https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.116.245&rep=rep1&type=pdf
# NOTE:
# NOTE:  Wong:
# NOTE:  "Automatic time-picking of first arrivals on noisy microseismic data"
# NOTE:  https://www.crewes.org/ForOurSponsors/ConferenceAbstracts/2009/CSEG/Wong_CSEG_2009.pdf
# NOTE:
# NOTE:  GEOPHYSICS, VOL. 75, NO. 4 (JULY-AUGUST 2010); P. V67-V76
# NOTE:  Automatic first-breaks picking: New strategies and algorithms"
# NOTE:  https://www.fcaglp.unlp.edu.ar/~velis/papers/PickingGeop10.pdf
# NOTE:
# NOTE:  "Adaptive microseismic event detection and automatic time picking"
# NOTE:  https://www.cspg.org/documents/Conventions/Archives/Annual/2012/279_GC2012_Adaptive_Microseismic_Event_Detection.pdf
# NOTE:
# NOTE:  Earle and Shearer:
# NOTE:  "Characterization of Global Seismograms Using an Automatic-Picking Algorithm"
# NOTE:  Bulletin of the Seismological Society of America, Vol. 84, No. 2, pp. 366-376, April 1994

if (!isGeneric("STALTA")) {
  setGeneric("STALTA", function(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment) {
    standardGeneric("STALTA")
  })
}

STALTA.Trace <- function(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment) {
  # Calculate constants
  n_sta <- staSecs * x@stats@sampling_rate
  n_lta <- ltaSecs * x@stats@sampling_rate

  if (n_sta < 1) {
    stop(paste("STALTA.Trace: STA window of", staSecs, "secs is too short for trace sampling rate of", x@stats@sampling_rate, "Hz."))
  }

  if (demean | detrend) {
    x <- DDT(x, demean, detrend, taper)
  }
  len <- length(x)

  # Stop if there is insufficient data
  if (len < n_lta) {
    stop(paste("STALTA.Trace: insufficient data."))
  }

  # DC signals will have all zeroes after demeaning.
  # Return all zeroes in this case.
  if (all(x@data == 0, na.rm = TRUE)) {
    return(x@data)
  }

  # Apply the algorithm over the entire data series

  if (algorithm == "classic_RR") {
    # "Automatic time-picking of first arrivals on noisy microseismic data"
    # NOTE:  In this reference, the LTA and STA windows use the same right edge (index is in bth windows)
    if (increment != 1) {
      stop(paste("STALTA.Trace: algorithm EarleAndShearer_envelope requires increment=1.", sep = ""))
    }
    sta <- seismicRoll::roll_mean(x@data^2, n_sta, increment, align = "right")
    lta <- seismicRoll::roll_mean(x@data^2, n_lta, increment, align = "right")
    picker <- sta / lta
  } else if (algorithm == "classic_LR") {
    # NOTE:  This is the same algorithm but with left-right alignment (index is in both windows)
    picker <- seismicRoll::roll_stalta(x@data^2, n_sta, n_lta, increment)
  } else if (algorithm == "EarleAndShearer_envelope") {
    if (increment != 1) {
      stop(paste("STALTA.Trace: algorithm EarleAndShearer_envelope requires increment=1.", sep = ""))
    }
    data <- envelope(x)@data
    sta <- seismicRoll::roll_mean(data, n_sta, increment, align = "left")
    lta <- seismicRoll::roll_mean(data, n_lta, increment, align = "right")
    picker <- sta / lta
  } else {
    stop(paste("STALTA.Trace: algorithm=\"", algorithm, "\" not recognized.", sep = ""))
  }

  return(picker)
}

# All parameters specified
setMethod(
  "STALTA", signature(
    x = "Trace", staSecs = "numeric", ltaSecs = "numeric",
    algorithm = "character", demean = "logical", detrend = "logical", taper = "numeric", increment = "numeric"
  ),
  function(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment) {
    STALTA.Trace(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment)
  }
)
# demean, detrend, taper and increment missing
setMethod(
  "STALTA", signature(
    x = "Trace", staSecs = "numeric", ltaSecs = "numeric",
    algorithm = "character", demean = "missing", detrend = "missing", taper = "missing", increment = "missing"
  ),
  function(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment) {
    STALTA.Trace(x, staSecs, ltaSecs, algorithm, TRUE, TRUE, 0.0, 1)
  }
)
# Only lstaSecs and ltaSecs
setMethod(
  "STALTA", signature(
    x = "Trace", staSecs = "numeric", ltaSecs = "numeric",
    algorithm = "missing", demean = "missing", detrend = "missing", taper = "missing", increment = "missing"
  ),
  function(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment) {
    STALTA.Trace(x, staSecs, ltaSecs, "classic_LR", TRUE, TRUE, 0.0, 1)
  }
)
# All parameters missing
setMethod(
  "STALTA", signature(
    x = "Trace", staSecs = "missing", ltaSecs = "missing",
    algorithm = "missing", demean = "missing", detrend = "missing", taper = "missing", increment = "missing"
  ),
  function(x, staSecs, ltaSecs, algorithm, demean, detrend, taper, increment) {
    STALTA.Trace(x, 3, 30, "classic_LR", TRUE, TRUE, 0.0, 1)
  }
)


# Simple triggering of p-wave onset --------------------------------------------

# NOTE:  Method to return the first timepoint at which STALTA is above a threshold
#
# NOTE:  This method only returns the first timepoint or NA.  It DOES NOT implement
# NOTE:  the full functionality found in the ObsPy method of the same name.
#
# NOTE:  https://docs.obspy.org/packages/autogen/obspy.signal.trigger.triggerOnset.html

if (!isGeneric("triggerOnset")) {
  setGeneric("triggerOnset", function(x, picker, threshold, index) {
    standardGeneric("triggerOnset")
  })
}

triggerOnset.Trace <- function(x, picker, threshold, index) {
  if (missing(threshold)) {
    threshold <- stats::quantile(picker, 0.99999, na.rm = TRUE)
  }

  if (missing(index)) {
    index <- FALSE
  }

  # At which index is the picker above the threshold?
  eventIndex <- which(picker > threshold)[1]

  # Return just the index if desired
  if (index) {
    return(eventIndex)
  }

  # Otherwise calculate the time at which that index occurs
  if (is.na(eventIndex)) {
    eventTime <- NA
  } else {
    eventTime <- x@stats@starttime + eventIndex / x@stats@sampling_rate
  }

  return(eventTime)
}

setMethod(
  "triggerOnset", signature(x = "Trace", picker = "numeric", threshold = "ANY", index = "ANY"),
  function(x, picker, threshold, index) {
    triggerOnset.Trace(x, picker, threshold, index)
  }
)


# Select the event window around p-wave onset ----------------------------------

if (!isGeneric("eventWindow")) {
  setGeneric("eventWindow", function(x, picker, threshold, windowSecs) {
    standardGeneric("eventWindow")
  })
}

eventWindow.Trace <- function(x, picker, threshold, windowSecs) {
  eventOnset <- triggerOnset(x, picker, threshold)
  windowStarttime <- eventOnset - windowSecs / 2
  windowEndtime <- eventOnset + windowSecs / 2

  windowTrace <- slice(x, windowStarttime, windowEndtime)

  return(windowTrace)
}

setMethod(
  "eventWindow", signature(x = "Trace", picker = "numeric", threshold = "numeric", windowSecs = "numeric"),
  function(x, picker, threshold, windowSecs) {
    eventWindow.Trace(x, picker, threshold, windowSecs)
  }
)
setMethod(
  "eventWindow", signature(x = "Trace", picker = "numeric", threshold = "numeric", windowSecs = "missing"),
  function(x, picker, threshold, windowSecs) {
    eventWindow.Trace(x, picker, threshold, windowSecs = 3600)
  }
)
setMethod(
  "eventWindow", signature(x = "Trace", picker = "numeric", threshold = "missing", windowSecs = "numeric"),
  function(x, picker, threshold, windowSecs) {
    eventWindow.Trace(x, picker, threshold = stats::quantile(picker, 0.999, na.rm = TRUE), windowSecs)
  }
)
setMethod(
  "eventWindow", signature(x = "Trace", picker = "numeric", threshold = "missing", windowSecs = "missing"),
  function(x, picker, threshold, windowSecs) {
    eventWindow.Trace(x, picker, threshold = stats::quantile(picker, 0.999, na.rm = TRUE), windowSecs = 3600)
  }
)


################################################################################
# Basic plotting
################################################################################

# NOTE:  I don't seem to be able to use optional arguments of type "ANY" in the
# NOTE:  function declaration.  Probably has to do with plot.~() being older
# NOTE:  style functions rather than S4 methods.

plot.Trace <- function(x, starttime = x@stats@starttime, endtime = x@stats@endtime,
                       subsampling = NULL, add = FALSE, ylab = "raw", xlab = "", ...) {
  # Create subsampled indices if necessary
  max_size <- 10000
  if (is.null(subsampling)) {
    if (length(x) <= max_size) {
      subsampling <- 1
    } else {
      subsampling <- as.integer(length(x) / (max_size / 2))
    }
  }
  indices <- seq(1, length(x), subsampling)

  # Create labels

  # remove ".M" or other quality factor designation as people don't expect it
  id <- stringr::str_sub(x@id, 1, stringr::str_length(x@id) - 2)
  main <- paste("Seismic Trace for ", id, sep = "")
  sensorText <- paste("(", x@Sensor, ")")
  # Create array of times
  times <- seq(from = x@stats@starttime, to = x@stats@endtime, length.out = length(x@data))

  if (length(x) == length(indices)) {
    xlab <- paste(xlab, "\n", " (", length(x), "points )")
  } else {
    xlab <- paste(xlab, "\n", " ( ", length(x), " points, subsampling=", subsampling, " for this plot. )", sep = "")
  }

  # Plot
  if (difftime(x@stats@endtime, x@stats@starttime, units = "days") > "1 day" & difftime(x@stats@endtime, x@stats@starttime, units = "days") < "7 days") {
    plot(x@data[indices] ~ times[indices], type = "l", main = main, xaxt = "n", xlim = c(starttime, endtime), xlab = xlab, ylab = ylab, ...)
    graphics::axis.POSIXct(1, at = seq(from = x@stats@starttime, to = x@stats@endtime, by = "day"), format = "%b %d")
  } else {
    plot(x@data[indices] ~ times[indices], type = "l", xlim = c(starttime, endtime), main = main, xlab = xlab, ylab = ylab, ...)
  }
  # title
  graphics::mtext(sensorText, line = 0.2, adj = 0.05)
}

# NOTE:  I don't seem to be able to use optional arguments of type "ANY" in the
# NOTE:  function declaration.  Probably has to do with plot.~() being older
# NOTE:  style functions rather than S4 methods.

setMethod("plot", signature(x = "Trace"), function(x, ...) plot.Trace(x, ...))


################################################################################
# END
################################################################################