R/Class-Stream.R

In IRISSeismic: Classes and Methods for Seismic Data Analysis

##
##    S4 classes for handling lists of 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 class for a Stream obj.
#
# This is a port of some of the functionality found in obspy.core.stream.py
#
#   https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.html
#
################################################################################

################################################################################
# Class Stream
#
# A Stream contains a list of Trace objects.
# The 'url' slot contains the URL of the webservice request that returned the
# stream -- useful for debugging.
#
################################################################################

setClass(
  "Stream",
  # typed slots (aka attributes) for class Stream
  representation(
    url = "character",
    requestedStarttime = "POSIXct",
    requestedEndtime = "POSIXct",
    act_flags = "integer",
    io_flags = "integer",
    dq_flags = "integer",
    timing_qual = "numeric",
    traces = "list"
  ),
  # default values for slots
  prototype(
    url = "",
    requestedStarttime = as.POSIXct("1900-01-01T00:00:00",
      format = "%Y-%m-%dT%H:%M:%OS6",
      tz = "GMT"
    ),
    requestedEndtime = as.POSIXct("1900-01-01T00:00:00",
      format = "%Y-%m-%dT%H:%M:%OS6",
      tz = "GMT"
    ),
    act_flags = rep(as.integer(0), 8),
    io_flags = rep(as.integer(0), 8),
    dq_flags = rep(as.integer(0), 8),
    timing_qual = as.numeric(NA),
    traces = list(new("Trace"))
  )
)

################################################################################
# Basic methods to operate directly on a Stream object's Trace data.
#
# The 'sd' method returns a value for each Trace in the Stream.
#
# In the obspy version of the code, the standard deviation method is called
# 'std'. In this package, the name is changed to 'sd' to match the name for the
# existing R function for standard deviation.
#
# The 'pmax', 'pmin', etc. methods return a value fore ach Trace in the Stream.
# The 'max', 'mean', etc. methods return a single value for all Traces in the
# Stream.
################################################################################

# NOTE:  The 'lapply' function returns the list that results from applying the
# NOTE:  in-line function to each element of x@traces.
# NOTE:  The 'unlist' function converts the list returned by 'lapply' into
# NOTE:  a numeric vector.

# Return unique Ids ------------------------------------------------------------

if (!isGeneric("uniqueIds")) {
  setGeneric(
    "uniqueIds",
    function(x) standardGeneric("uniqueIds")
  )
}
uniqueIds.Stream <- function(x, na.rm = FALSE) {
  ids <- unlist(lapply(x@traces, slot, "id"))
  return(unique(ids))
}
setMethod("uniqueIds", signature(x = "Stream"), function(x) uniqueIds.Stream(x))


# Parallel Length --------------------------------------------------------------

if (!isGeneric("parallelLength")) {
  setGeneric(
    "parallelLength",
    function(x) standardGeneric("parallelLength")
  )
}
parallelLength.Stream <- function(x, na.rm = FALSE) {
  return(unlist(lapply(x@traces, function(element) length(element))))
}
setMethod(
  "parallelLength",
  signature(x = "Stream"),
  function(x) parallelLength.Stream(x)
)


# Parallel Maximum -------------------------------------------------------------

if (!isGeneric("parallelMax")) {
  setGeneric(
    "parallelMax",
    function(x, na.rm) standardGeneric("parallelMax")
  )
}
parallelMax.Stream <- function(x, na.rm) {
  return(unlist(lapply(
    x@traces,
    function(element) max(element, na.rm = na.rm)
  )))
}
setMethod(
  "parallelMax",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelMax.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelMax",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelMax.Stream(x, na.rm = FALSE)
)


# Parallel Mean ----------------------------------------------------------------

if (!isGeneric("parallelMean")) {
  setGeneric(
    "parallelMean",
    function(x, na.rm) standardGeneric("parallelMean")
  )
}
parallelMean.Stream <- function(x, na.rm) {
  return(unlist(lapply(
    x@traces,
    function(element) mean(element, na.rm = na.rm)
  )))
}
setMethod(
  "parallelMean",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelMean.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelMean",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelMean.Stream(x, na.rm = FALSE)
)


# Parallel Minimum--------------------------------------------------------------

if (!isGeneric("parallelMin")) {
  setGeneric(
    "parallelMin",
    function(x, na.rm) standardGeneric("parallelMin")
  )
}
parallelMin.Stream <- function(x, na.rm) {
  return(unlist(lapply(x@traces, function(element) min(element, na.rm = na.rm))))
}
setMethod(
  "parallelMin",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelMin.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelMin",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelMin.Stream(x, na.rm = FALSE)
)


# Parallel Median --------------------------------------------------------------

if (!isGeneric("parallelMedian")) {
  setGeneric(
    "parallelMedian",
    function(x, na.rm) standardGeneric("parallelMedian")
  )
}
parallelMedian.Stream <- function(x, na.rm) {
  return(unlist(lapply(x@traces, function(element) median(element, na.rm = na.rm))))
}
setMethod(
  "parallelMedian",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelMedian.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelMedian",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelMedian.Stream(x, na.rm = FALSE)
)


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

if (!isGeneric("parallelSd")) {
  setGeneric(
    "parallelSd",
    function(x, na.rm) standardGeneric("parallelSd")
  )
}
parallelSd.Stream <- function(x, na.rm) {
  return(unlist(lapply(x@traces, function(element) sd(element, na.rm = na.rm))))
}
setMethod(
  "parallelSd",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelSd.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelSd",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelSd.Stream(x, na.rm = FALSE)
)


# Parallel RMS -----------------------------------------------------------------

if (!isGeneric("parallelRms")) {
  setGeneric(
    "parallelRms",
    function(x, na.rm) standardGeneric("parallelRms")
  )
}
parallelRms.Stream <- function(x, na.rm) {
  return(unlist(lapply(
    x@traces,
    function(element) rms(element, na.rm = na.rm)
  )))
}
setMethod(
  "parallelRms",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelRms.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelRms",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelRms.Stream(x, na.rm = FALSE)
)


# Parallel RMS variance --------------------------------------------------------

if (!isGeneric("parallelRmsVariance")) {
  setGeneric(
    "parallelRmsVariance",
    function(x, na.rm) standardGeneric("parallelRmsVariance")
  )
}
parallelRmsVariance.Stream <- function(x, na.rm) {
  return(unlist(lapply(
    x@traces,
    function(element) rmsVariance(element, na.rm = na.rm)
  )))
}
setMethod(
  "parallelRmsVariance",
  signature(x = "Stream", na.rm = "logical"),
  function(x, na.rm) parallelRmsVariance.Stream(x, na.rm = na.rm)
)
setMethod(
  "parallelRmsVariance",
  signature(x = "Stream", na.rm = "missing"),
  function(x, na.rm) parallelRmsVariance.Stream(x, na.rm = FALSE)
)


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

length.Stream <- function(x) {
  return(sum(parallelLength(x)))
}
setMethod("length", signature(x = "Stream"), function(x) length.Stream(x))


# Global Maximum ---------------------------------------------------------------

max.Stream <- function(x, ..., na.rm = FALSE) {
  return(max(parallelMax(x, na.rm = na.rm)))
}
setMethod("max", signature(x = "Stream"), function(x, ...) max.Stream(x, ...))


# Global Mean ------------------------------------------------------------------

mean.Stream <- function(x, ...) {
  data <- unlist(lapply(x@traces, slot, "data"))
  return(mean(data, ...))
}
setMethod("mean", signature(x = "Stream"), function(x, ...) mean.Stream(x, ...))


# Global Minimum ---------------------------------------------------------------

min.Stream <- function(x, ..., na.rm = FALSE) {
  return(min(parallelMin(x, na.rm = na.rm)))
}
setMethod("min", signature(x = "Stream"), function(x, ...) min.Stream(x, ...))


# Global Median ----------------------------------------------------------------

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


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

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


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

if (!isGeneric("rms")) {
  setGeneric(
    "rms",
    function(x, na.rm) standardGeneric("rms")
  )
}
rms.Stream <- function(x, na.rm) {
  data <- unlist(lapply(x@traces, slot, "data"))
  return(sqrt(mean((data)^2, na.rm = na.rm)))
}
setMethod(
  "rms",
  signature("Stream", na.rm = "logical"),
  function(x, na.rm) rms.Stream(x, na.rm = na.rm)
)
setMethod(
  "rms",
  signature("Stream", na.rm = "missing"),
  function(x, na.rm) rms.Stream(x, na.rm = FALSE)
)


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

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


################################################################################
# Method to return a new Stream object where every Trace has been multiplied
# by a constant.
################################################################################

if (!isGeneric("multiplyBy")) {
  setGeneric(
    "multiplyBy",
    function(x, y) standardGeneric("multiplyBy")
  )
}
multiplyBy.Stream <- function(x, y) {
  traces <- lapply(x@traces, function(element) multiplyBy(element, y = y))
  return(new("Stream",
    url = x@url,
    requestedStarttime = x@requestedStarttime,
    requestedEndtime = x@requestedEndtime,
    act_flags = x@act_flags,
    io_flags = x@io_flags,
    dq_flags = x@dq_flags,
    timing_qual = x@timing_qual,
    traces = traces
  ))
}
setMethod(
  "multiplyBy",
  signature("Stream", y = "numeric"),
  function(x, y) multiplyBy.Stream(x, y = y)
)

################################################################################
# Method to return a new Stream where every Trace has been assigned a specified
# quality code. This method modifies each trace @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 each @stats@quality slot with the requested
# quality code. The first input parameter is the Stream 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.Stream <- function(x, quality) {
  if (!quality %in% c("R", "D", "Q", "M")) {
    stop("Requested quality code must be one of \"R\", \"D\", \"Q\", \"M\"")
  }
  for (i in seq_along(x@traces)) {
    x@traces[[i]]@id <- sub(".$", quality, x@traces[[i]]@id)
    x@traces[[i]]@stats@quality <- quality
  }
  return(x)
}

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

################################################################################
# Method to return a list of all gaps/overlaps of Traces in the Stream object
#
# getGaps(x, min_gap)
#  min_gap -- minimum gap (sec) below which gaps will be ignored
#             (default = NULL)
#
################################################################################

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

getGaps.Stream <- function(x, min_gap) {
  # Sanity check -- single SNCL
  num_ids <- length(uniqueIds(x))
  if (num_ids > 1) {
    stop(paste("getGaps.Stream:", num_ids, "unique ids encountered in Stream."))
  }

  # NOTE:  The number of potential gaps includes an initial gap, a final gap
  # NOTE:  and a gap between each trace. Thus num_gaps = 2 + num_headers - 1.
  # NOTE:  We will set up an array for num_header+1 gaps but insert 0 if the
  # NOTE:  gap is too small to count.

  # NOTE:  Use difftime() instead of just subtracting to guarantee that units
  # are "secs"

  # Extract a list of TraceHeaders from the list of Traces
  headers <- lapply(x@traces, slot, "stats")
  num_headers <- length(headers)

  sampling_rates <- sapply(headers, slot, "sampling_rate")

  if (any(sampling_rates < 0)) {
    stop(paste("getGaps.Stream: encountered sampling rate < 0"))
  }

  # Set up arrays for information about gaps/overlaps
  gaps <- numeric(num_headers + 1)
  nsamples <- integer(num_headers + 1)

  for (i in seq(from = 1, to = num_headers)) {
    # NOTE:  The delta we calculate here has a valid datapoint at it's start and
    # end. For the purposes of calculating the correct number of missing points
    # we will calculate delta as the total time between points minus 1/sampling
    # rate. That is how many extra points could be shoved into this gap.

    if (i == 1) {
      # Initial gap (no overlap possible)
      sampling_rate <- sampling_rates[1]
      # Set min_gap and make sure it is at least 1/sampling_rate
      min_gap_new <- ifelse(is.null(min_gap), 1 / sampling_rate, min_gap)
      min_gap_new <- max(min_gap_new, 1 / sampling_rate)

      delta <- as.numeric(difftime(headers[[1]]@starttime,
        x@requestedStarttime,
        units = "secs"
      )) - 1 / sampling_rate
      if (delta > min_gap_new - 0.5 / sampling_rate) {
        gaps[1] <- delta + 1 / sampling_rate
        nsamples[1] <- as.integer(round(gaps[1] * sampling_rate))
      } else {
        gaps[1] <- 0
        nsamples[1] <- 0
      }
    } else {
      # Inter-trace gaps and overlaps
      sampling_rate <- sampling_rates[[i - 1]]
      min_gap_new <- ifelse(is.null(min_gap), 1 / sampling_rate, min_gap)
      min_gap_new <- max(min_gap_new, 1 / sampling_rate)
      h1 <- headers[[i - 1]]
      h2 <- headers[[i]]
      delta <- difftime(h2@starttime,
        h1@endtime,
        units = "secs"
      ) - 1 / sampling_rate
      if (abs(delta) > min_gap_new - 0.5 / sampling_rate) {
        gaps[i] <- delta
        nsamples[i] <- as.integer(round(abs(delta) * sampling_rate))
      } else {
        gaps[i] <- 0
        nsamples[i] <- 0
      }
    }
    if (i == num_headers) {
      # Final gap (no overlap possible)
      sampling_rate <- sampling_rates[[i]]
      min_gap_new <- ifelse(is.null(min_gap), 1 / sampling_rate, min_gap)
      min_gap_new <- max(min_gap_new, 1 / sampling_rate)
      delta <- as.numeric(difftime(x@requestedEndtime,
        headers[[num_headers]]@endtime,
        units = "secs"
      )) - 1 / sampling_rate
      if (delta > min_gap_new - 0.5 / sampling_rate) {
        gaps[num_headers + 1] <- delta
        nsamples[num_headers + 1] <- as.integer(round(delta * sampling_rate))
      } else {
        gaps[num_headers + 1] <- 0
        nsamples[num_headers + 1] <- 0
      }
    }
  }

  gap_list <- list(gaps = gaps, nsamples = nsamples)
  return(gap_list)
}

# All parameters specified
setMethod(
  "getGaps", signature(x = "Stream", min_gap = "numeric"),
  function(x, min_gap) getGaps.Stream(x, min_gap)
)
# min_gap missing
setMethod(
  "getGaps", signature(x = "Stream", min_gap = "missing"),
  function(x, min_gap) getGaps.Stream(x, NULL)
)


################################################################################
# Method to return a vector of up/down times identifying the start- and
# end-times of Traces in the Stream object.  Care is taken to identify
# overlapping Traces gaps below a minimum threshold and remove those timepoints.
#
# getUpDownTimes(Stream, min_signal, min_gap)
#   min_signal -- minimum signal (sec) below which Traces will be ignored
#                 (default = NULL)
#   min_gap -- minimum gap (sec) below which gaps will be ignored
#                 (default = NULL)
#
################################################################################

if (!isGeneric("getUpDownTimes")) {
  setGeneric(
    "getUpDownTimes",
    function(x, min_signal, min_gap) standardGeneric("getUpDownTimes")
  )
}
getUpDownTimes.Stream <- function(x, min_signal, min_gap) {
  # Extract a list of ids from the list of Traces
  num_ids <- length(uniqueIds(x))
  if (num_ids > 1) {
    stop(paste(
      "getUpDownTimes.Stream:",
      num_ids,
      "unique ids encountered in Stream."
    ))
  }

  # Extract a list of TraceHeaders
  headerList <- lapply(x@traces, slot, "stats")

  # Extract a list of start- and end-times from the list of Traces
  starttimeList <- lapply(headerList, slot, "starttime")
  endtimeList <- lapply(headerList, slot, "endtime")

  # NOTE:  unlist() will convert 'POSIXct' to 'numeric' so we need to go through
  # 'character'. This seems inefficient but we're only dealing with a small
  # number of items.
  starttimes <- as.POSIXct(
    unlist(lapply(starttimeList,
      strftime,
      format = "%Y-%m-%dT%H:%M:%OS",
      tz = "GMT"
    )),
    format = "%Y-%m-%dT%H:%M:%OS", tz = "GMT"
  )
  endtimes <- as.POSIXct(
    unlist(lapply(endtimeList,
      strftime,
      format = "%Y-%m-%dT%H:%M:%OS",
      tz = "GMT"
    )),
    format = "%Y-%m-%dT%H:%M:%OS", tz = "GMT"
  )

  # Determine signal duration for each trace
  signal_durations <- difftime(endtimes, starttimes, units = "sec")

  good_traces_flag <- signal_durations >= min_signal

  # Cull to remove those Traces whose duration is < min_signal
  headerList <- headerList[good_traces_flag]
  starttimes <- starttimes[good_traces_flag]
  endtimes <- endtimes[good_traces_flag]
  num_headers <- length(headerList)

  # Return if there is only one trace
  if (num_headers == 1) {
    up_down_times <- c(starttimes[1], endtimes[1])
    return(up_down_times)
  }

  # Create a dummy vector of POSIXct values of length 2 X num_headers.
  up_down_times <- c(starttimes, endtimes)

  # Now go through and insert appropriate times, checking for any overlaps or
  # spaces < min_gap.
  # NOTE:  up_down_times[1] and up_down_times[num_headers*2] are already in
  # place.
  for (i in seq(num_headers - 1)) {
    # get the difference in seconds between the current end time and the next
    # start time.
    delta <- difftime(starttimes[i + 1], endtimes[i], units = "secs")
    if (delta < 0 || delta < min_gap) {
      # insert NA empty space to indicate continuity
      up_down_times[(2 * i)] <- NA
      up_down_times[(2 * i) + 1] <- NA
    } else {
      # insert our current end time and write the next start time to the vector,
      # indicating discontinuity
      up_down_times[(2 * i)] <- endtimes[i]
      up_down_times[(2 * i) + 1] <- starttimes[i + 1]
    }
  }

  # remove all of the NA values that represent continuity
  # return the vector of start, end, start, end....
  return(stats::na.omit(up_down_times))
}
setMethod(
  "getUpDownTimes", signature(
    x = "Stream",
    min_signal = "numeric",
    min_gap = "numeric"
  ),
  function(x, min_signal, min_gap) {
    getUpDownTimes.Stream(x, min_signal, min_gap)
  }
)
setMethod(
  "getUpDownTimes", signature(
    x = "Stream",
    min_signal = "missing",
    min_gap = "missing"
  ),
  function(x, min_signal, min_gap) {
    getUpDownTimes.Stream(x, min_signal = 30, min_gap = 60)
  }
)


################################################################################
# Method to return a new Stream object sliced out of an existing Stream
#
# Unlike the ObsPy method, this method does no padding and does not retain empty
# Traces.  The returned Stream will always be a subset of the original Stream
#
# Slicing is rounded to the nearest second.
#
# https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.slice.html#obspy.core.stream.Stream.slice
################################################################################

if (!isGeneric("slice")) {
  setGeneric(
    "slice",
    function(x, starttime, endtime) standardGeneric("slice")
  )
}
slice.Stream <- function(x, starttime, endtime) {
  num_traces <- length(x@traces)
  stream_start <- x@traces[[1]]@stats@starttime
  stream_end <- x@traces[[num_traces]]@stats@endtime

  # Sanity check
  if (starttime >= endtime) {
    stop(paste(
      "slice.Stream: requested starttime \"",
      starttime,
      "\" >= requested endtime \"",
      endtime,
      "\""
    ))
  }
  if (starttime >= stream_end) {
    stop(paste(
      "slice.Stream: requested starttime \"",
      starttime,
      "\" >= Stream endtime \"",
      stream_end,
      "\""
    ))
  }
  if (endtime <= stream_start) {
    stop(paste(
      "slice.Stream: requested endtime \"",
      endtime,
      "\" <= Stream starttime \"",
      stream_start,
      "\""
    ))
  }

  # Set up a new list for traces we are keeping
  traces <- list()

  # TODO:  Get rid of append() in slice.Stream().

  # Check each trace to see if it should be ignored, included entirely or
  # included after slicing
  for (i in seq(num_traces)) {
    tr <- x@traces[[i]]
    if (starttime >= tr@stats@endtime || endtime <= tr@stats@starttime) {
      # skip this trace
    } else if (starttime <= tr@stats@starttime && endtime >= tr@stats@endtime) {
      # this traces does not need to be sliced
      traces <- append(traces, tr)
    } else {
      sliced_trace <- slice(tr, starttime, endtime)
      traces <- append(traces, sliced_trace)
    }
  }

  return(new("Stream",
    url = x@url,
    requestedStarttime = starttime,
    requestedEndtime = endtime,
    act_flags = x@act_flags,
    io_flags = x@io_flags,
    dq_flags = x@dq_flags,
    timing_qual = x@timing_qual,
    traces = traces
  ))
}
setMethod(
  "slice",
  signature(x = "Stream", starttime = "POSIXct", endtime = "POSIXct"),
  function(x, starttime, endtime) {
    slice.Stream(x, starttime = starttime, endtime = endtime)
  }
)

################################################################################
# Method to merge Traces in a Stream into a single Trace, replacing gaps with
# values determined by the fillMethod.
################################################################################

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

mergeTraces.Stream <- function(x, fillMethod) {
  num_traces <- length(x@traces)

  # Return immediately if is only one trace with no initial or final gap
  if (sum(getGaps(x)$nsamples) == 0 && num_traces == 1) {
    return(x)
  }

  gapInfo <- getGaps(x)
  num_gaps <- length(gapInfo$nsamples)

  # Sanity checks
  if (num_gaps != num_traces + 1) {
    stop(paste("mergeTraces.Stream: num_gaps (",
      num_gaps,
      ") should be one more than num_traces (",
      num_traces,
      ")",
      sep = ""
    ))
  }

  headers <- lapply(x@traces, slot, "stats")
  sampling_rates <- sapply(headers, slot, "sampling_rate")
  num_rates <- length(unique(round(sampling_rates, digits = 4)))
  # if (num_rates > 1) {
  if (!all(stats::dist(unique(sampling_rates)) < 0.0002)) {
    # slightly more forgiving criteria for acceptable sample rate jitter than
    # round(sampling_rates, digits = 4)
    stop(paste(
      "mergeTraces.Stream:",
      num_rates,
      "unique sampling rates encountered in Stream."
    ))
  }

  # some fdsnws/dataselect implementations cut on record instead of sample
  # like IRIS. Trace starttime can be < requestedStarttime and
  # trace endtime can be > requestedEndtime
  if (gapInfo$nsamples[1] == 0) {
    totalStart <- x@traces[[1]]@stats@starttime
  } else {
    totalStart <- x@requestedStarttime
  }

  if (gapInfo$nsamples[length(gapInfo$nsamples)] == 0) {
    totalEnd <- x@traces[[length(x@traces)]]@stats@endtime
  } else {
    totalEnd <- x@requestedEndtime
  }

  totalSecs <- as.numeric(difftime(totalEnd, totalStart, units = "secs"))
  totalPoints <- as.integer(round(totalSecs) * x@traces[[1]]@stats@sampling_rate)

  # NOTE:  Setting up a list of vectors that we will concatenate in one fell
  # NOTE:  swoop at the end. This avoids the incremental growth and copying
  # NOTE:  that is incredibly slow.
  # NOTE:  See R Inferno Chapter 2 -- Growing Objects
  # NOTE:    https://www.burns-stat.com/pages/Tutor/R_inferno.pdf

  num_vectors <- num_gaps + num_traces
  dataList <- vector("list", num_vectors)

  # Fill in the data for gaps and traces
  if (fillMethod == "fillNA") {
    for (i in seq(num_traces)) {
      dataList[[2 * i - 1]] <- rep(NA, gapInfo$nsamples[i])
      dataList[[2 * i]] <- x@traces[[i]]@data
    }
    dataList[[num_vectors]] <- rep(NA, gapInfo$nsamples[[num_gaps]])
  } else if (fillMethod == "fillZero") {
    for (i in seq(num_traces)) {
      dataList[[2 * i - 1]] <- rep(0, gapInfo$nsamples[i])
      dataList[[2 * i]] <- x@traces[[i]]@data
    }
    dataList[[num_vectors]] <- rep(0, gapInfo$nsamples[[num_gaps]])
  } else {
    stop(paste("mergeTraces.Stream: unknown fillMethod '",
      fillMethod,
      "'",
      sep = ""
    ))
  }

  # Create a single data vector
  data <- unlist(dataList)

  missing_points <- totalPoints - length(data)

  # Sanity check -- we hope to be within twice the sampling rate of a complete
  # accounting
  if (missing_points > ceiling(2 * x@traces[[1]]@stats@sampling_rate)) {
    stop(paste(
      "mergeTraces.Stream:",
      missing_points,
      "unaccounted for points after merge"
    ))
  } else if (missing_points < ceiling(-2 * x@traces[[1]]@stats@sampling_rate)) {
    stop(paste(
      "mergeTraces.Stream:",
      abs(missing_points),
      "extra points after merge"
    ))
  }

  # Guarantee that we return the correct number of points
  if (missing_points > 0) {
    if (fillMethod == "fillNA") {
      data <- c(data, rep(NA, missing_points))
    } else if (fillMethod == "fillZero") {
      data <- c(data, rep(0, missing_points))
    }
  }


  # Create a new TraceHeader
  stats <- x@traces[[1]]@stats
  stats@npts <- as.integer(totalPoints)
  if (gapInfo$nsamples[1] == 0) {
    stats@starttime <- x@traces[[1]]@stats@starttime
  } else {
    stats@starttime <- x@requestedStarttime
  }
  if (gapInfo$nsamples[length(gapInfo$nsamples)] == 0) {
    stats@endtime <- x@traces[[length(x@traces)]]@stats@endtime
  } else {
    stats@endtime <- x@requestedEndtime
  }
  stats@processing <-
    append(
      stats@processing,
      paste(num_traces,
        " traces merged into a single trace using method '",
        fillMethod,
        "'",
        sep = ""
      )
    )

  # Other Trace info
  id <- x@traces[[1]]@id
  Sensor <- x@traces[[1]]@Sensor
  InstrumentSensitivity <- x@traces[[1]]@InstrumentSensitivity
  SensitivityFrequency <- x@traces[[1]]@SensitivityFrequency
  InputUnits <- x@traces[[1]]@InputUnits

  traces <- list(new("Trace",
    id,
    stats,
    Sensor,
    InstrumentSensitivity,
    SensitivityFrequency,
    InputUnits,
    data = data[1:totalPoints]
  ))

  return(new("Stream",
    url = x@url,
    requestedStarttime = x@requestedStarttime,
    requestedEndtime = x@requestedEndtime,
    act_flags = x@act_flags,
    io_flags = x@io_flags,
    dq_flags = x@dq_flags,
    timing_qual = x@timing_qual,
    traces = traces
  ))
}

# All parameters specified
setMethod(
  "mergeTraces", signature(x = "Stream", fillMethod = "character"),
  function(x, fillMethod) mergeTraces.Stream(x, fillMethod = fillMethod)
)
# fillMethod missing
setMethod(
  "mergeTraces", signature(x = "Stream", fillMethod = "missing"),
  function(x, fillMethod) mergeTraces.Stream(x, fillMethod = "fillNA")
)


################################################################################
# Various methods and functions for upDownTimes
#
# upDownTimes is a vector of class POSIXct that contains the GMT times at which
# a channel or channels is on (up) or off (down).  The first and all odd
# elements in the list are datetimes at which data collection starts.
# Even numbered elements correspond to datetimes at which data collection stops.
################################################################################

################################################################################
# Plotting upDownTimes from Streams or pre-generated vectors of upDownTimes
################################################################################

if (!isGeneric("plotUpDownTimes")) {
  setGeneric("plotUpDownTimes", function(x, min_signal, min_gap, ...) {
    standardGeneric("plotUpDownTimes")
  })
}

# One method that gets upDOwnTimes from a Stream Object
plotUpDownTimes.Stream <- function(x, min_signal = 30, min_gap = 60, ...) {
  upDownTimes <- getUpDownTimes(x, min_signal = min_signal, min_gap = min_gap)

  # NOTE:  Extra conversion step needed to guarantee display as GMT times
  GMTTimes <- as.POSIXct(upDownTimes, "%Y-%m-%dT%H:%M:%OS", tz = "GMT")
  onOffs <- seq_along(GMTTimes) %% 2 # 1, 0, 1, 0, ...

  # Plot the transitions with appropriate axes
  # Extend the times and onOff transitions to the full time period requested
  allTimes <- c(x@requestedStarttime, GMTTimes, x@requestedEndtime)
  allOnOffs <- c(0, onOffs, 0)
  plot(allTimes, allOnOffs, type = "s", xlab = "GMT", ylab = "", yaxt = "n", ...)
  graphics::axis(2, at = c(0, 1), labels = c("Off", "On"), las = 1, tick = TRUE)

  # Cover over any transition to down at the end
  # NOTE:  Seems not to work if I use GMTTimes so revert back to upDownTimes
  graphics::abline(v = x@requestedEndtime, lwd = 3, col = "white")
  graphics::box()

  # Add the title
  id <- stringr::str_sub(
    x@traces[[1]]@id, 1,
    stringr::str_length(x@traces[[1]]@id) - 2
  )
  main <- paste("On/Off transitions for ", id)
  sensorText <- paste("(", x@traces[[1]]@Sensor, ")")
  graphics::title(main)
  graphics::mtext(sensorText, line = 0.2)
}

# Another method that accepts a pre-generated vector of plotUpDownTimes
plotUpDownTimes.POSIXct <- function(x, min_signal = 30, min_gap = 60, ...) {
  # NOTE:  min_signal and min_gap is not used in this function.
  # NOTE:  They are included for consistency.
  upDownTimes <- x

  # NOTE:  Extra step needed to guarantee display as GMT times
  GMTTimes <- as.POSIXct(upDownTimes, "%Y-%m-%dT%H:%M:%OS", tz = "GMT")
  onOff <- seq_along(GMTTimes) %% 2 # 1, 0, 1 0, ...

  # TODO:  Check for and use additional plotting parameters like "main" or "sub"
  # TODO:  Use optional "starttime", "endtime" arguments to set xlim

  # Plot the transitions with appropriate axes
  plot(GMTTimes, onOff, type = "s", xlab = "GMT", ylab = "", yaxt = "n", ...)
  graphics::axis(2, at = c(0, 1), labels = c("Off", "On"), las = 1, tick = TRUE)

  # Remove any transition to down at the end
  # NOTE:  Seems not to work if I use GMTTimes so rever back to upDownTimes
  graphics::abline(v = upDownTimes[length(upDownTimes)], lwd = 2, col = "white")
  graphics::box()
}

setMethod(
  "plotUpDownTimes",
  signature(x = "Stream", min_signal = "numeric", min_gap = "numeric"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.Stream(x, min_signal, min_gap = min_gap, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "Stream", min_signal = "numeric", min_gap = "missing"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.Stream(x, min_signal, min_gap = 60, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "Stream", min_signal = "missing", min_gap = "numeric"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.Stream(x, min_signal = 30, min_gap, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "Stream", min_signal = "missing", min_gap = "missing"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.Stream(x, min_signal = 30, min_gap = 60, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "POSIXct", min_signal = "numeric", min_gap = "numeric"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.POSIXct(x, min_signal, min_gap, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "POSIXct", min_signal = "numeric", min_gap = "missing"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.POSIXct(x, min_signal = 30, min_gap, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "POSIXct", min_signal = "missing", min_gap = "numeric"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.POSIXct(x, min_signal, min_gap = 60, ...)
  }
)
setMethod(
  "plotUpDownTimes",
  signature(x = "POSIXct", min_signal = "missing", min_gap = "missing"),
  function(x, min_signal, min_gap, ...) {
    plotUpDownTimes.POSIXct(x, min_signal = 30, min_gap = 60, ...)
  }
)

################################################################################
# Merging upDownTimes to return a new vector of datetimes that mark when
# both (either) channel is on/off
################################################################################

if (!isGeneric("mergeUpDownTimes")) {
  setGeneric(
    "mergeUpDownTimes",
    function(udt1, udt2, bothOn) standardGeneric("mergeUpDownTimes")
  )
}
mergeUpDownTimes.POSIXct <- function(udt1, udt2, bothOn) {
  # Deal with udt1 or udt2 = c()
  if (is.null(udt1)) {
    return(udt2)
  }
  if (is.null(udt2)) {
    return(udt1)
  }

  # Combine the times and sort them
  unsorted_times <- c(udt1, udt2)
  sort_indices <- order(unsorted_times)
  times <- unsorted_times[sort_indices]

  # Create a numeric vector of on/off flags with on = 1, off = -1
  onOff1 <- (seq_along(udt1) %% 2 - 0.5) * 2
  onOff2 <- (seq_along(udt2) %% 2 - 0.5) * 2
  # Combine the flags and sort them temporally
  unsorted_onOff <- c(onOff1, onOff2)
  onOff <- unsorted_onOff[sort_indices]

  # NOTE:  The times and onOffs are now sorted and we will forgout about
  # NOTE:  'sort_indices'. The indices in the next section have nothing to do
  # NOTE:  with  the original sort_indices

  # The cumulative sum now has the following meaning: 2 = both on, 0 = 1 on,
  # -1 = both off
  cumOnOff <- cumsum(onOff)

  if (bothOn) {
    # Create up_down_times associated with "both channels on/either channel off"
    both_up_indices <- which(cumOnOff == 2)
    any_down_indices <- both_up_indices + 1
    both_upDownTimes <- sort(times[c(both_up_indices, any_down_indices)])
    return(both_upDownTimes)
  } else {
    # Create up_down_times associated with "either channel on/both channels off"
    both_down_indices <- which(cumOnOff == 0)
    # up indices start with the first time and use every time just after both
    # down (except for the last one)
    either_up_indices <- c(1, both_down_indices[-length(both_down_indices)] + 1)
    either_upDownTimes <- sort(times[c(either_up_indices, both_down_indices)])
    return(either_upDownTimes)
  }
}
setMethod(
  "mergeUpDownTimes",
  signature(udt1 = "POSIXct", udt2 = "POSIXct", bothOn = "logical"),
  function(udt1, udt2, bothOn) {
    mergeUpDownTimes.POSIXct(udt1, udt2, bothOn)
  }
)
setMethod(
  "mergeUpDownTimes",
  signature(udt1 = "POSIXct", udt2 = "POSIXct", bothOn = "missing"),
  function(udt1, udt2, bothOn) {
    mergeUpDownTimes.POSIXct(udt1, udt2, bothOn = FALSE)
  }
)
setMethod(
  "mergeUpDownTimes",
  signature(udt1 = "NULL", udt2 = "POSIXct", bothOn = "logical"),
  function(udt1, udt2, bothOn) {
    mergeUpDownTimes.POSIXct(udt1, udt2, bothOn)
  }
)
setMethod(
  "mergeUpDownTimes",
  signature(udt1 = "NULL", udt2 = "POSIXct", bothOn = "missing"),
  function(udt1, udt2, bothOn) {
    mergeUpDownTimes.POSIXct(udt1, udt2, bothOn = FALSE)
  }
)
setMethod(
  "mergeUpDownTimes",
  signature(udt1 = "POSIXct", udt2 = "NULL", bothOn = "logical"),
  function(udt1, udt2, bothOn) {
    mergeUpDownTimes.POSIXct(udt1, udt2, bothOn)
  }
)
setMethod(
  "mergeUpDownTimes",
  signature(udt1 = "POSIXct", udt2 = "NULL", bothOn = "missing"),
  function(udt1, udt2, bothOn) {
    mergeUpDownTimes.POSIXct(udt1, udt2, bothOn = FALSE)
  }
)


################################################################################
# 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.Stream <- function(x, ...) {
  tr <- mergeTraces(x)@traces[[1]]

  plot(tr, ...)
  if (length(x@traces) == 1) {
    graphics::mtext(paste(length(x@traces), "trace"),
      side = 3, line = 0.2, adj = 0.95
    )
  } else {
    graphics::mtext(paste(length(x@traces), "traces"),
      side = 3, line = 0.2, adj = 0.95
    )
  }
}

# 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 = "Stream"), function(x, ...) plot.Stream(x, ...))