Package 'eseis'

Description

Environmental seismoloy studies the seismic signals, emitted by Earth surface processes. This package eseis provides all relevant functions to read/write seismic data files, prepare, analyse and visualise seismic data, and generate reports of the processing history.

Author(s)

Maintainer: Michael Dietze [email protected] [translator]

Authors:

• Clement Hibert [contributor]

Other contributors:

• Christoph Burow [contributor]

• Sophie Lagarde [contributor, translator]

Description

The function checks seismic files organised by aux_organisecubefiles or aux_organisecentaurfiles for completeness. The tests include agreement of file name and seismic file meta data.

Usage

aux_checkfiles(
  dir,
  station,
  component = "BHZ",
  method = "thorough",
  period = "weekly",
  format = "sac",
  duration_set = 3600,
  plot = TRUE
)

Arguments

Value

Data frame containing check results and meta data.

Author(s)

Michael Dietze

Examples

## Not run: 

## set seismic data directory
dir_data <- paste0(system.file("extdata", package="eseis"), "/")

## check data archive for record completeness
chk <- aux_checkfiles(dir = dir_data, 
                      station = "RUEG1", 
                      component = "BHZ",
                      plot = TRUE)


## End(Not run)

Description

The function compares the sampling intervals of a list of eseis objects and identifies the highest common sampling interval (dt) as well as the aggregation factors for each eseis object needed to reach this common sampling interval.

Usage

aux_commondt(data, dt)

Arguments

Value

list object with elements dt (highest common sampling interval) and agg (aggregation factors for each of the input data sets to reach the common sampling interval)

Author(s)

Michael Dietze

Examples

## Not run: 

## TO BE WRITTEN

## End(Not run)

Description

This is a simple wrapper for the Gipptools program cubeinfo, providing a short summary of the cube file meta data, in a coherent data frame structure.

Usage

aux_cubeinfo(file, gipptools)

Arguments

Value

data frame with cube meta data

Author(s)

Michael Dietze

Examples

## Not run: 

## get cube info
x = aux_cubeinfo(file = "data/cube/example.ATB", 
                 gipptools = "/software/gipptools-2019.332/")


## End(Not run)

Description

The function converts an eseis object to an ObsPy stream object. The functionality is mainly useful when running ObsPy through R using the package 'reticulate'. Currently, only single traces (i.e., single eseis objects) can be converted. Thus, to convert multiple traces, these need to be converted individually and added to the first trace using ObsPy functionalities.

Usage

aux_eseisobspy(data)

Arguments

Value

ObsPy stream object as defined by the architecture of package 'reticulate'.

Author(s)

Michael Dietze

Examples

## Not run: 

## load ObsPy library with package 'reticulate'
## (requires ObsPy to be installed on the computer)
obspy <- reticulate::import("obspy")

## load example data set
data(rockfall)

## convert example eseis object to ObsPy stream object
x <- aux_eseisobspy(data = rockfall_eseis)

## filter data set using ObsPy
x_filter <- obspy$traces[[1]]$filter(type = "bandpass", 
                                     freqmin = 0.5, 
                                     freqmax = 1.0)
                                     
## plot filtered trace using ObsPy plotting routine
x$traces[[1]]$plot()


## End(Not run)

Description

This function is a wrapper for the library 'dataselect' from IRIS. It reads a corrupt mseed file and saves it in fixed state. Therefore, the function requires dataselect being installed (see details).

Usage

aux_fixmseed(file, input_dir, output_dir, software)

Arguments

Details

The library 'dataselect' can be downloaded at https://github.com/iris-edu/dataselect and requires compilation (see README file in dataselect directory). The function goes back to an email discussion with Gillian Sharer (IRIS team), many thanks for pointing me at this option to process corrupt mseed files.

Value

a set of mseed files written to disk.

Author(s)

Michael Dietze

Examples

## Not run: 

aux_fixmseed(file = list.files(path = "~/data/mseed", 
                               pattern = "miniseed"), 
                    input_dir = "~/data/mseed",
                    software = "~/software/dataselect-3.17")


## End(Not run)

Description

The function loads seismic data from a data directory structure (see aux_organisecubefiles()) based on the event start time, duration, component and station ID.

Usage

aux_getevent(
  start,
  duration,
  station,
  component = "BHZ",
  format,
  dir,
  simplify = TRUE,
  eseis = TRUE,
  try = FALSE,
  verbose = FALSE
)

Arguments

Details

The data to be read needs to be adequately structured. The data directory must contain SAC files organised by year (e.g.2022) then by Julian Day in full three digits (e.g. 001) and then by a dedicated SAC file name, containing the station ID, two-digit year, three-digit Julian Day, start time hour, minute and second, three channel ID and the file extension SAC. All these items need to be separated by stops (e.g. sac/2022/001/LAU01.22.001.08.00.00. BHZ.SAC). This data structure will be most conveniently created by the functions aux_organisecubefiles() or aux_organisecentaurfiles(), or by manually written R code.

The function assumes complete data sets, i.e., not a single hourly data set must be missing. The time vector is loaded only once, from the first station and its first component. Thus, it is assumed that all loaded seismic signals are of the same sampling frequency and length.

Value

A list object containing either a set of eseis objects or a data set with the time vector ($time) and a list of seismic stations ($station_ID) with their seismic signals as data frame ($signal). If simplify = TRUE (the default option) and only one seismic station is provided, the output object containseither just one eseis object or the vectors for $time and $signal.

Author(s)

Michael Dietze

Examples

## set seismic data directory
dir_data <- paste0(system.file("extdata", package="eseis"), "/")

## load the z component data from a station
data <- aux_getevent(start = as.POSIXct(x = "2017-04-09 01:20:00", 
                                        tz = "UTC"), 
                      duration = 120,
                      station = "RUEG1",
                      component = "BHZ",
                      dir = dir_data)                       
## plot signal
plot_signal(data = data)

## load data from two stations
data <- aux_getevent(start = as.POSIXct(x = "2017-04-09 01:20:00", 
                                        tz = "UTC"), 
                     duration = 120,
                     station = c("RUEG1", "RUEG2"),
                     component = "BHZ",
                     dir = dir_data)

## plot both signals
par(mfcol = c(2, 1))
lapply(X = data, FUN = plot_signal)

Description

The function accesses the specified FDSN internet data base(s) and downloads seismic data based on the network and station IDs and time constraints.

Usage

aux_getFDSNdata(
  start,
  duration,
  channel = "BHZ",
  network,
  station,
  url,
  link_only = FALSE,
  eseis = TRUE
)

Arguments

Details

A convenient way to get all the required input data is using the function aux_getFDSNstation before. It will return all the information in a structured way.

It is possible to use the function to process more than one data set. In this case, the arguments network, station and url must match pairwise. The arguments start, duration and channel will be treated as constants if not also provided as vectors.

Value

List object with imported seismic data for each provided set of input arguments.

Author(s)

Michael Dietze

See Also

aux_get_FDSNstation, read_mseed

Examples

## Not run: 

## get stations < 0.6 degrees away from Piz Chengalo collapse
x <- aux_getFDSNstation(centre = c(46.3, 9.6),
                        radius = 0.6,
                        access = TRUE)

## sort statiions by distance
x <- x[order(x$distance),]

## download available data
d <- aux_getFDSNdata(start = as.POSIXct(x = "2017-08-23 07:30:00", 
                                        tz = "UTC"),
                     duration = 180, 
                     network = x$network_ID, 
                     station = x$station_code, 
                     url = x$network_url)

## remove stations without available data
x <- x[!unlist(lapply(d, is.null)),]
d <- d[!unlist(lapply(d, is.null))]

## generate plots of the three nearest stations
par(mfcol = c(3, 1))

for(i in 1:3) {

  plot_signal(data = d[[i]],
              main = paste(x$ID[i], 
                           " | ",
                           round(x$distance[i], 2),
                           "distance (DD)"))
} 

## End(Not run)

Description

This function queries as series of data bases for seismic stations that match a set of criteria for seismic data. The criteria include signal time stamp and location, and component. The returned data can be used to download data using the function aux_FDSNdata.

Usage

aux_getFDSNstation(centre, radius, start, access, url)

Arguments

Details

The function requires a working internet connection to perform the query. It uses the following FDSN data bases by default:

• orfeus "http://www.orfeus-eu.org"

• geofon "http://geofon.gfz-potsdam.de/"

• bgr "http://eida.bgr.de"

• sss "http://eida.ethz.ch"

Other FDSN data base addresses can be provided in the same way as the addresses in the above list. They need to be provided as character vector. For a list of addresses see "http://www.fdsn.org/webservices/datacenters/" and "http://docs.obspy.org/packages/obspy.clients.fdsn.html#module-obspy.clients.fdsn".

Value

Data frame with query results. The data frame contains information for all seismic stations fulfilling the defined criteria.

Author(s)

Michael Dietze

See Also

aux_get_FDSNdata, aux_getIRISstation

Examples

## Not run: 

x <- aux_getFDSNstation(start = as.POSIXct(x = "2010-01-01 22:22:22", 
                                           tz = "UTC"), 
                        centre = c(45, 10), 
                        radius = 1)
                           
## optionally plot station locations on a map (requires RgoogleMaps)
center <- c(mean(x$station_latitude), 
            mean(x$station_longitude))

zoom <- min(RgoogleMaps::MaxZoom(range(x$station_latitude), 
                                 range(x$station_longitude)))
                                 
Map <- RgoogleMaps::GetMap(center = center,
                           zoom = zoom, 
                           maptype = "terrain")
                           
RgoogleMaps::PlotOnStaticMap(MyMap = Map, 
                             lat = x$station_latitude, 
                             lon = x$station_longitude, 
                             pch = 15, 
                             col = 4)

## End(Not run)

Description

This function reads auxiliary information stored in Omnirecs/Digos Datacube files and extracts the temperature data that is stored along with each GPS tag. Optionally, the data is interpolated to equal intervals.

Usage

aux_gettemperature(input_dir, logger_ID, interval, cpu, gipptools)

Arguments

Details

This feature is ony available for Omnirecs/Digos Datacube that were produced since 2015, i.e., whose GPS output files also record the temperature inside the logger. Generating an ACSII GPS tag file using the gipptools software requires a few minutes time per daily file.

Value

A list of data frames with time and temperature values for each cube data logger.

Author(s)

Michael Dietze

Examples

## uncomment to use
# t <- aux_gettemperature(input_dir = "input",
#                         logger_ID = c("ANN", "ABT"),
#                         interval = 15,
#                         gipptools = "~/software/gipptools-2015.225/")

Description

This function either downloads an online station XML file and (optionally) saves it and/or reads an already existing local station XML file.

Usage

aux_getxml(
  xml,
  start,
  duration,
  network,
  station,
  component,
  url,
  level = "response"
)

Arguments

Details

Currently, the function uses Obspy python code. Hence, both python and the package 'obspy' need to be installed in order to use the function.

Value

Obspy object with station meta info inventory.

Author(s)

Michael Dietze

Examples

## Not run: 

x <- aux_getxml(start = "2010-10-10", 
                duration = 60,
                network = "GE",
                station = "BRNL",
                component = "BHZ",
                url = "http://service.iris.edu")


## End(Not run)

Description

This function cuts a three component seismic data set into time windows that may or may not overlap and calculates the spectral ratio for each of these windows. It returns a matrix with the ratios for each time slice. Thus, it is a wrapper for the function signal_hvratio. For further information about the technique and function arguments see the description of signal_hvratio.

Usage

aux_hvanalysis(
  data,
  time,
  window,
  overlap = 0,
  dt,
  method = "periodogram",
  kernel,
  order = "xyz",
  plot = FALSE,
  ...
)

Arguments

Value

A matrix with the h-v-frequency ratios for each time slice.

Author(s)

Michael Dietze

Examples

## load example data set
data(earthquake)

## ATTENTION, THIS EXAMPLE DATA SET IS FAR FROM IDEAL FOR THIS PURPOSE

## detrend data
s <- signal_detrend(data = s)

## calculate the HV ratios straightforward
HV <- aux_hvanalysis(data = s,
                     dt = 1 / 200,
                     kernel = 100)

## calculate the HV ratios with plot output, userdefined palette
plot_col <- colorRampPalette(colors = c("grey", "darkblue", "blue", "orange"))
HV <- aux_hvanalysis(data = s,
                     dt = 1 / 200,
                     kernel = 100,
                     plot = TRUE,
                     col = plot_col(100))

## calculate the HV ratios with optimised method settings
HV <- aux_hvanalysis(data = s, 
                     time = t,
                     dt = 1 / 200, 
                     window = 10, 
                     overlap = 0.9, 
                     method = "autoregressive",
                     plot = TRUE,
                     col = plot_col(100),
                     xlab = "Time (UTC)",
                     ylab = "f (Hz)")
                     
## calculate and plot stack (mean and sd) of all spectral ratios
HV_mean <- apply(X = HV, MARGIN = 1, FUN = mean)
HV_sd <- apply(X = HV, MARGIN = 1, FUN = sd)
HV_f <- as.numeric(rownames(HV))

plot(x = HV_f, y = HV_mean, type = "l", ylim = c(0, 50))
lines(x = HV_f, y = HV_mean + HV_sd, col = 2)
lines(x = HV_f, y = HV_mean - HV_sd, col = 2)

Description

The function generates an empty eseis object, starting with processing step 0. The object contains no data and the history only contains the system information.

Usage

aux_initiateeseis()

Details

The S3 object class eseis contains the data vector ($signal), a meta information list ($meta) with all essential seismic meta data - such as sampling interval, station ID, component, start time of the stream or file name of the input file - a list with header data of the seismic source file ($header), and a history list ($history), which records all data manipulation steps of an (eseis) object. The element ($meta) will be used by functions of the package to look for essential information to perform data manipulations (e.g., the sampling interval). Thus, working with (eseis) objects is convenient and less prone to user related errors/bugs, given that the meta information is correct and does not change during the processing chain; package functions will update the meta information whenever necessary (e.g., signal_aggregate). The element $header will only be present if a seismic source file has been imported.

The history element is the key feature for transparent and reproducable research using this R package. An eseis object records a history of every function it has been subject to, including the time stamp, the function call, all used function arguments and their associated values, and the overall processing duration in seconds. The history is updated whenever an eseis object is manipulated with one of the functions of this package (with a few exceptions, mainly from the aux_... category).

Value

S3 list object of class eseis.

Author(s)

Michael Dietze

Examples

## initiate eseis object
aux_initiateeseis()

Description

The function converts an ObsPy stream object to an eseis object. The functionality is mainly useful when running ObsPy through R using the package 'reticulate'.

Usage

aux_obspyeseis(data, simplify = TRUE)

Arguments

Details

Initiation of the reticulate-based python toolbox support can be cumbersome. The following suggestions from Guthub (https://github.com/rstudio/reticulate/issues/578) helped in a case study:

library(reticulate) use_condaenv("r-reticulate") py_install("obspy", pip = TRUE)

Value

eseis object.

Author(s)

Michael Dietze

Examples

## Not run: 

## load ObsPy library with package 'reticulate'
## (requires ObsPy to be installed on the computer)
obspy <- reticulate::import("obspy")

## set seismic data directory
dir_data <- paste0(system.file("extdata", package="eseis"), "/")

## read miniseed file to stream object via ObsPy
x <- obspy$read(pathname_or_url = "dir_data/2017/99/RUEG1.17.99.00.00.00.BHZ.SAC")

## convert ObsPy stream object to eseis object
y <- aux_obspyeseis(data = x)

## plot eseis object
plot_signal(y)

## End(Not run)

Description

This function optionally converts mseed files to sac files and organises these in a coherent directory structure, by year, Julian day, (station, hour and channel). It depends on the cubetools or gipptools software package (see details). The function is at an experimental stage and only used for data processing at the GFZ Geomorphology section, currently.

Usage

aux_organisecentaurfiles(
  stationfile,
  input_dir,
  output_dir,
  format = "sac",
  channel_name = "bh",
  cpu,
  gipptools,
  file_key = "miniseed",
  network
)

Arguments

Details

The function assumes that the Nanometrics Centaur data logger directory contains only hourly mseed files. These hourly files are organised in a coherent directory structure which is organised by year and Julian day. In each Julian day directory the hourly files are placed and named according to the following scheme: STATIONID.YEAR.JULIANDAY.HOUR.MINUTE.SECOND.CHANNEL.
The function requires that the software cubetools (http://www.omnirecs.de/documents.html) or gipptools (http://www.gfz-potsdam.de/en/section/geophysical-deep-sounding/infrastructure/geophysical-instrument-pool-potsdam-gipp/software/gipptools/) are installed.
Specifying an input directory (input_dir) is mandatory. This input directory must only contain the subdirectories with mseed data for each Centaur logger. The subdirectory must be named after the four digit Centaur ID and contain only mseed files, regardless if further subdirectories are used (e.g., for calendar days).

In the case a six-channel Centaur is used to record signals from two sensors, in the station info file (cf. aux_stationinfofile()) the logger ID field must contain the four digit logger ID and the channel qualifiers, e.g., "AH" (first three channels) or "BH" (last three channels), separated by an underscore.

Value

A set of hourly seismic files written to disk.

Author(s)

Michael Dietze

Examples

## Not run: 

## basic example with minimum effort
aux_organisecentaurfiles(stationfile = "output/stationinfo.txt", 
                         input_dir = "input", 
                         gipptools = "software/gipptools-2015.225/")


## End(Not run)

Description

The function converts Omnirecs/Digos Datacube files to mseed or sac files and organises these in a coherent directory structure (see details) for available structures. The conversion depends on the gipptools software package (see details) provided externally.

Usage

aux_organisecubefiles(
  station,
  input,
  output,
  gipptools,
  format = "sac",
  pattern = "eseis",
  component = "BH",
  mode = "dir-wise",
  fringe = "constant",
  cpu,
  verbose = TRUE
)

Arguments

Details

The function converts seismic data from the binary cube file format to mseed (cf. read_mseed) or sac (cf. read_sac) and organises the resulting files into a consistent structure, expected by 'eseis' for convenient data handling (cf. read_data).

Currently, there are two data structure schemes supported, "eseis" and "seiscomp". In the "eseis" case, the daily cube files are cut to hourly files and organised in directories structured by four digit year and three digit Julian day numbers. In each Julian day directory, the hourly files are placed and named after the following scheme: STATION.YEAR.JULIANDAY.HOUR.MINUTE.SECOND.COMPONENT.

The "seiscomp" case will yield daily files, which are organised by four digit year, seismic network, seismic station, and seismic component, each building a separate directory. In the deepest subdirectory, files are named by: NETWORK.STATION.LOCATION.COMPONENT.TYPE.YEAR.JULIANDAY.

The component naming scheme defines the codes for the sensor's band code (first letter) and instrument code (second letter). The third letter, defining the spatial component, will be added automatically. For definitions of channel codes see https://migg-ntu.github.io/SeisTomo_Tutorials/seismology/seismic-data/seismic-time-series-data.html.

The function requires that the software gipptools (http://www.gfz-potsdam.de/en/section/geophysical-deep-sounding/infrastructure/geophysical-instrument-pool-potsdam-gipp/software/gipptools/) is installed. Note that the gipptools are provided at regular update intervals, including an up to date GPS leap second table, essential to convert recently recorded files.

The Cube files will be imported in place but a series of temporary files will be created in a temporary directory in the specified output directory. Hence, if the routine stops due to a processing issue, one needs to delete the temporary data manually. The path to the temporary directory will be provided as screen output when the argument verbose = TRUE.

The Cube files can be converted in two modes: "file-wise" and "dir-wise". In "file-wise" mode, each Cube file will be converted individually. This option has the advantage that if one file in a month-long sequence of records is corrupt, the conversion will not stop, but only discard the part from the corrupted section until the file end. The disadvantage is however, that the data before the first and after the last GPS tags will not be converted unless the option fringe = "constant" (by default this is the case) is used.

In "dir-wise" mode, the fringe sample issue reduces to the margins of the total sequence of daily files but the corrupt file issue will become a more severe danger to the success when converting a large number of files.

Specifying an input directory (input) is mandatory. That directory should only contain the directories with the cube files to process. Files downloaded from a Cube are usually contained in one or more further directories, which should be moved into a single one before running this function.

Each set of cube files from a given logger should be located in a separate directory per logger and these directories should have the same name as the logger IDs (logger_ID). An appropriate structure for files from two loggers, A1A and A1B, would be something like:

1. input

   1. A1A

      1. file1.A1A

      2. file2.A1A

   2. A1B

      1. file1.A1B

      2. file2.A1B

The component definition can follow the typical keywords and key letters defined in seismology: https://migg-ntu.github.io/SeisTomo_Tutorials/seismology/seismic-data/seismic-time-series-data.html, hence the first letter indicating the instrument's band type and the second letter indicating the instrument code or instrument type.

Value

A set of converted and organised seismic files written to disk.

Author(s)

Michael Dietze

Examples

## Not run: 

## basic example with minimum effort
aux_organisecubefiles(stationfile = data.frame(logger = c("A1A", "A1B"),
                                               station = c("ST1", "ST2")), 
                      input = "input", 
                      gipptools = "software/gipptools-2023.352")


## End(Not run)

Description

The function allows detection of discrete seismic events using the STA-LTA picker approach and including several stations of a seismic network. It will pick events at all input stations and identify events picked by more than a user defined minimum number of stations, within a given time window. It further allows rejecting events that are too short or too long, or occur to short after a previous event. The function can be used to process long data archives organised in a consistent file structure.

Usage

aux_picknetwork(
  start,
  stop,
  res,
  buffer,
  station,
  component,
  dir,
  f,
  envelope = TRUE,
  sta,
  lta,
  on,
  off,
  freeze,
  dur_min,
  dur_max,
  n_common,
  t_common,
  t_pause,
  verbose = FALSE,
  ...
)

Arguments

Details

The data will be imported time slice by time slice. Optionally, a signal deconvolution can be done. The data will be filtered to isolate the frequency range of interest, and tapered to account for edge effects. The taper size is automatically set be the two-sided buffer (see below). Optionally, a signal envelope can be calculated as part of the preprocessing work. For all successfully imported and preprocessed signal time snippets, events will be detected by the STL-LTA approach. Events from all input signals will be checked for a minium and maximum allowed event duration, set as thresholds by the user. All remaining events are checked for the number of stations of the network that did consistently detect an event. That means an event must be found in each signal within a user defined time window, a window that corresponds to the maximum travel time of a seismic signal through the network. The minimum number of stations that detect the event can be set by the user. In addition, events can be rejected if they occur during a user defined ban time after the onset of a previous event, which assumes that no subsequent events are allowed during the implementation of long duration events happening.

The time period to process is defined by start, end and the length of time snippets res (s). In addition, a left sided and right sided buffer can and should be added buffer = c(60, 60) (s), to account for artefacts during signal preprocessing and events that occur across the time snippets. The start and end time should be provided in POSIXct format. If it is provided as text string, the function will try to convert to POSIXct assuming the time zone is UTC.

There must be at least two seismic stations. The seismic components component must be a vector of the same length as station and contain the corresponding information for each station.

Deconvolution of the input signals can be done, although usually this is not essential for just picking events if the filter frequencies are in the flat response range of the sensors. If wanted, one needs to provide the additional vectors (all of the same length as station) sensor, logger and optionally gain (otherwise assumed to be 1). See signal_deconvolve() for further information on the deconvolution routine and the supported keywords for sensors and loggers.

STA-LTA picking is usually done on signal envelopes enevelope = TRUE. However, for array seismology or other specific cases it might be useful to work with the waveforms directly, instead of with their envelopes.

The STA-LTA routine requires information on the length of the STA window, the LTA window, the on-threshold to define the start of an event, the off-threshold to define the end of an event, and the option to freeze the STA-LTA ratio at the onset of and event. See the documentation of the function pick_stalta() for further information.

Events that last shorter than dur_min or longer than dur_max are rejected. This criterion is applied before the test of how many stations have commonly detected an event. Hence, the minimum and maximum duration criteria need to be fulfilled by an event for all seismic stations, or at least as many as defined as minimum number of stations n_common.

A valid event should be detected by more than one station. At least two stations (n_common = 2) should detect it to confirm it is not just a local effect like rain drop impacts or animal activity. At least three stations (n_common = 3) are needed to locate an event. Detection at more stations increases the validity of an event as proper signal. However, the likelihood of stations not operating correctly or that some stations may be too far to detect a smaller event limits the number of stations that are be set as (n_common).

An event that happens just outside a seismic network will generate seismic waves that travel through the network and will require a given time for that. The time is set by the longest inter-station distance and the apparent seismic wave velocity. For example, an event who's signal propagated with 1000 m/s through a network with a maximum aperture of 1000 m will need up to 1 s. Hence, the parameter t_common should be set to 1, in this case. This parameter is good to reject events that show a longer travel time, e.g. because they are pressure waves that travel through the atmosphere, like helicopter signals or lightning strikes.

Value

data frame, detected events (start time, duration in seconds)

Author(s)

Michael Dietze

Examples

picks <- aux_picknetwork(start = "2017-04-09 01:00:00 UTC",
                         stop = "2017-04-09 02:00:00 UTC",
                         res = 1800,
                         buffer = c(90, 90),
                         station = c("RUEG1","RUEG2"),
                         component = rep("BHZ", 2),
                         dir = paste0(path.package("eseis"), "/extdata/"),
                         f = c(0.1, 0.2),
                         envelope = TRUE,
                         sta = 0.5,
                         lta = 300,
                         on = 3,
                         off = 1,
                         freeze = TRUE,
                         dur_min = 2,
                         dur_max = 90,
                         n_common = 2,
                         t_common = 5,
                         t_pause = 30,
                         verbose = TRUE)
                         
print(picks)

Description

The function allows detection of discrete seismic events using the STA-LTA picker approach and including several stations of a seismic network. It will pick events at all input stations and identify events picked by more than a user defined minimum number of stations, within a given time window. It further allows rejecting events that are too short or too long, or occur to short after a previous event. The function can be used to process long data archives organised in a consistent file structure. This is the parallel computation version of the routine. For the single CPU version, allowing to show verbose progress information, use aux_picknetwork().

Usage

aux_picknetworkparallel(
  start,
  stop,
  res,
  buffer,
  station,
  component,
  dir,
  f,
  envelope = TRUE,
  sta,
  lta,
  on,
  off,
  freeze,
  dur_min,
  dur_max,
  n_common,
  t_common,
  t_pause,
  cpu,
  ...
)

Arguments

Details

The data will be imported time slice by time slice. Optionally, a signal deconvolution can be done. The data will be filtered to isolate the frequency range of interest, and tapered to account for edge effects. The taper size is automatically set be the two-sided buffer (see below). Optionally, a signal envelope can be calculated as part of the preprocessing work. For all successfully imported and preprocessed signal time snippets, events will be detected by the STL-LTA approach. Events from all input signals will be checked for a minium and maximum allowed event duration, set as thresholds by the user. All remaining events are checked for the number of stations of the network that did consistently detect an event. That means an event must be found in each signal within a user defined time window, a window that corresponds to the maximum travel time of a seismic signal through the network. The minimum number of stations that detect the event can be set by the user. In addition, events can be rejected if they occur during a user defined ban time after the onset of a previous event, which assumes that no subsequent events are allowed during the implementation of long duration events happening.

The time period to process is defined by start, end and the length of time snippets res (s). In addition, a left sided and right sided buffer can and should be added buffer = c(60, 60) (s), to account for artefacts during signal preprocessing and events that occur across the time snippets. The start and end time should be provided in POSIXct format. If it is provided as text string, the function will try to convert to POSIXct assuming the time zone is UTC.

There must be at least two seismic stations. The seismic components component must be a vector of the same length as station and contain the corresponding information for each station.

Deconvolution of the input signals can be done, although usually this is not essential for just picking events if the filter frequencies are in the flat response range of the sensors. If wanted, one needs to provide the additional vectors (all of the same length as station) sensor, logger and optionally gain (otherwise assumed to be 1). See signal_deconvolve() for further information on the deconvolution routine and the supported keywords for sensors and loggers.

STA-LTA picking is usually done on signal envelopes enevelope = TRUE. However, for array seismology or other specific cases it might be useful to work with the waveforms directly, instead of with their envelopes.

The STA-LTA routine requires information on the length of the STA window, the LTA window, the on-threshold to define the start of an event, the off-threshold to define the end of an event, and the option to freeze the STA-LTA ratio at the onset of and event. See the documentation of the function pick_stalta() for further information.

Events that last shorter than dur_min or longer than dur_max are rejected. This criterion is applied before the test of how many stations have commonly detected an event. Hence, the minimum and maximum duration criteria need to be fulfilled by an event for all seismic stations, or at least as many as defined as minimum number of stations n_common.

A valid event should be detected by more than one station. At least two stations (n_common = 2) should detect it to confirm it is not just a local effect like rain drop impacts or animal activity. At least three stations (n_common = 3) are needed to locate an event. Detection at more stations increases the validity of an event as proper signal. However, the likelihood of stations not operating correctly or that some stations may be too far to detect a smaller event limits the number of stations that are be set as (n_common).

An event that happens just outside a seismic network will generate seismic waves that travel through the network and will require a given time for that. The time is set by the longest inter-station distance and the apparent seismic wave velocity. For example, an event who's signal propagated with 1000 m/s through a network with a maximum aperture of 1000 m will need up to 1 s. Hence, the parameter t_common should be set to 1, in this case. This parameter is good to reject events that show a longer travel time, e.g. because they are pressure waves that travel through the atmosphere, like helicopter signals or lightning strikes.

Value

data frame, detected events (start time, duration in seconds)

Author(s)

Michael Dietze

Examples

picks <- aux_picknetwork(start = "2017-04-09 01:00:00 UTC",
                         stop = "2017-04-09 02:00:00 UTC",
                         res = 1800,
                         buffer = c(90, 90),
                         station = c("RUEG1","RUEG2"),
                         component = rep("BHZ", 2),
                         dir = paste0(path.package("eseis"), "/extdata/"),
                         f = c(0.1, 0.2),
                         envelope = TRUE,
                         sta = 0.5,
                         lta = 300,
                         on = 3,
                         off = 1,
                         freeze = TRUE,
                         dur_min = 2,
                         dur_max = 90,
                         n_common = 2,
                         t_common = 5,
                         t_pause = 30)
                         
print(picks)

Description

The function generates a long time PSD with aggregated time and frequency resolution.

Usage

aux_psdsummary(
  start,
  stop,
  ID,
  component = "BHZ",
  dir,
  window,
  sensor,
  logger,
  gain = 1,
  hours_skip,
  res = 1000,
  n = 100,
  cpu,
  verbose = FALSE
)

Arguments

Details

The function will calculate PSDs using the Welch method (see signal_spectrogram), with no overlap of the main time windows. The sub-windows will be automatically set to 10 the overlap of sub-windows to 0.5.

Value

eseis object, a spectrogram

Author(s)

Michael Dietze

Examples

## Not run: 

p <- aux_psdsummary(start = "2017-04-15 19:00:00 UTC",
                    stop = "2017-04-15 22:00:00 UTC",
                    ID = "RUEG1",
                    component = "BHE",
                    dir = "~/data/sac/",
                    sensor = "TC120s", 
                    logger = "Cube3ext",
                    window = 600,
                    res = 1000,
                    verbose = TRUE)
                    
plot(p)


## End(Not run)

Description

The function converts a seismic signal to sound and saves it as a wav file.

Usage

aux_sonifysignal(
  data,
  file,
  aggregate = 1,
  amplification = 10^6,
  speed = 1,
  dt
)

Arguments

Value

Sound file in wav format, written to disk.

Author(s)

Michael Dietze

Examples

## Not run: 

## load example data
data(rockfall)

## deconvolve and taper signal
s <- signal_deconvolve(data = rockfall_eseis)
s <- signal_taper(data = s, p = 0.05)

## sonify as is (barely sensible, due to too low frequency)
aux_sonifysignal(data = s, 
                 file = "~/Downloads/r1.wav")

## sonify at 20-fold speed
aux_sonifysignal(data = s, 
                 file = "~/Downloads/r1.wav", 
                 speed = 20)

## End(Not run)

Description

The function changes the meta data and file names of seismic data sets (organised in an eseis-compatible data structure). Data cubes in combination with splitter break-out boxes allow to record vertical component signals from three single-channel geophones, each connected by an individual cable. The logger interprets the data as E, N and Z components, regardless of the actual sensor and component connected to the respective cable. Hence, it is possible to record data from three independent (vertical) sensors with a single Data cube. However, the channels will be named according to the internal scheme (E, N, Z). To correct for that effect, the function will look up the right combination of station ID and component, and change meta data and name of each file in the input directory accordingly, before saving the updated file in the output directory.

Usage

aux_splitcubechannels(
  input,
  output,
  ID_in,
  component_in,
  ID_out,
  component_out,
  gipptools
)

Arguments

Value

A set of converted and organised seismic files written to disk.

Author(s)

Michael Dietze

Examples

## Not run: 

## the example will convert the E, N and Z components of station RUEG1
## to the station IDs RUEGA, RUEGB, RUEGC, all with BHZ as component

aux_splitcubechannels(input =  paste0(system.file("extdata", 
                                      package="eseis"), "/"), 
                      output = tempdir(), 
                      ID_in = c("RUEG1", "RUEG1", "RUEG1"), 
                      component_in = c("BHE", "BHN", "BHZ"),
                      ID_out = c("RUEGA", "RUEGB", "RUEGC"),
                      component_out = c("BHZ", "BHZ", "BHZ"))


## End(Not run)

Description

This function reads GPS tags from Omnirecs/Digos Datacube files and creates a station info file from additional input data. It depends on the cubetools or gipptools software package (see details).

Usage

aux_stationinfofile(
  name,
  input_dir,
  output_dir,
  station_ID,
  station_name,
  station_z,
  station_d,
  sensor_type,
  logger_type,
  sensor_ID,
  logger_ID,
  unit = "dd",
  n,
  quantile = 0.95,
  random = TRUE,
  cpu,
  gipptools,
  write_file = TRUE,
  write_raw = FALSE,
  write_data = FALSE
)

Arguments

Details

A station info file is an ASCII file that contains all relevant information about the individual stations of a seismic network. The variables contain a station ID (containing not more than 5 characters), station name, latitude, longitude, elevation, deployment depth, sensor type, logger type, sensor ID and logger ID.
The function requires that the software cubetools (http://www.omnirecs.de/documents.html) or gipptools (http://www.gfz-potsdam.de/en/section/geophysical-deep-sounding/infrastructure/geophysical-instrument-pool-potsdam-gipp/software/gipptools/) are installed. Note that GPS tag extraction may take several minutes per cube file. Hence, depending on the number of files and utilised CPUs the processing may take a while.
Specifying an input directory (input_dir) is mandatory. This input directory must only contain the subdirectories with the cube files to process, each set of cube files must be located in a separate subdirectory and these subdiretories must have the same name as specified by the logger IDs (logger_ID). An appropriate structure would be something like:

1. input

   1. A1A

      1. file1.A1A

      2. file2.A1A

   2. A1B

      1. file1.A1B

      2. file2.A1B

Value

A set of files written to disk and a data frame with seismic station information.

Author(s)

Michael Dietze

Examples

## Not run: 

## basic example with minimum effort
aux_stationinfofile(name = "stationinfo.txt", 
                    input_dir = "input", 
                    logger_ID = c("864", "876", "AB1"),
                    gipptools = "software/gipptools-2015.225")

## example with more adjustments
aux_stationinfofile(name = "stationinfo.txt",
                    input_dir = "input",
                    logger_ID = c("864", "876", "AB1"),
                    station_name = c("KTZ01", "KTZ02", "KTZ03"), 
                    station_z = c(30, 28, 29), 
                    station_d = rep(0.5, 3), 
                    sensor_type = rep("TC120s", 3), 
                    logger_type = rep("Cube3ext", 3), 
                    unit = "utm", 
                    n = 1, 
                    cpu = 0.9,
                    gipptools = "software/gipptools-2015.225", 
                    write_raw = TRUE, 
                    write_data = TRUE)


## End(Not run)

Description

The dataset comprises the seismic signal (all three components) of a small earthquake. The data have been recorded at 200 Hz sampling frequency with an Omnirecs Cube ext 3 data logger.

The dataset comprises the time vector associated with the data set earthquake.

Usage

s

t

Format

The format is: List of 3 $ BHE: num [1:8001] -3.95e-07 ... $ BHN: num [1:8001] -2.02e-07 ... $ BHZ: num [1:8001] -1.65e-07 ...

The format is: POSIXct[1:98400], format: "2015-04-06 13:16:54" ...

Examples

## load example data set
data(earthquake)

## plot signal vector
plot(x = t, y = s$BHZ, type = "l")


## load example data set
data(earthquake)

## show range of time vector
range(t)

Description

The fluvial model inversion (FMI) routine uses a predefined look-up table with reference spectra to identify those spectra that fit the empirical data set best, and returns the corresponding target variables and fit errors.

Usage

fmi_inversion(reference, data, n_cores = 1)

Arguments

Details

Note that the frequencies of the empiric and modelled data sets must match.

Value

List object containing the inversion results.

Author(s)

Michael Dietze

Examples

## NOTE THAT THE EXAMPLE IS OF BAD QUALITY BECAUSE ONLY 10 REFERENCE 
## PARAMETER SETS AND SPECTRA ARE CALCULATED, DUE TO COMPUATATION TIME
## CONSTRAINTS. SET THE VALUE TO 1000 FOR MORE MEANINGFUL RESULTS.

## create 100 example reference parameter sets
ref_pars <- fmi_parameters(n = 10,
                           h_w = c(0.02, 1.20),
                           q_s = c(0.001, 8.000) / 2650,
                           d_s = 0.01,
                           s_s = 1.35,
                           r_s = 2650,
                           w_w = 6,
                           a_w = 0.0075,
                           f_min = 5,
                           f_max = 80,
                           r_0 = 6,
                           f_0 = 1,
                           q_0 = 10,
                           v_0 = 350,
                           p_0 = 0.55,
                           e_0 = 0.09,
                           n_0_a = 0.6,
                           n_0_b = 0.8,
                           res = 100)

## create corresponding reference spectra
ref_spectra <- fmi_spectra(parameters = ref_pars)

## define water level and bedload flux time series
h <- c(0.01, 1.00, 0.84, 0.60, 0.43, 0.32, 0.24, 0.18, 0.14, 0.11)
q <- c(0.05, 5.00, 4.18, 3.01, 2.16, 1.58, 1.18, 0.89, 0.69, 0.54) / 2650
hq <- as.list(as.data.frame(rbind(h, q)))

## calculate synthetic spectrogram
psd <- do.call(cbind, lapply(hq, function(hq) {

  psd_turbulence <- eseis::model_turbulence(h_w = hq[1],
                                            d_s = 0.01,
                                            s_s = 1.35,
                                            r_s = 2650,
                                            w_w = 6,
                                            a_w = 0.0075,
                                            f = c(10, 70),
                                            r_0 = 5.5,
                                            f_0 = 1,
                                            q_0 = 18,
                                            v_0 = 450,
                                            p_0 = 0.34,
                                            e_0 = 0.0,
                                            n_0 = c(0.5, 0.8),
                                            res = 100, 
                                            eseis = FALSE)$power

  psd_bedload <- eseis::model_bedload(h_w = hq[1],
                                      q_s = hq[2],
                                      d_s = 0.01,
                                      s_s = 1.35,
                                      r_s = 2650,
                                      w_w = 6,
                                      a_w = 0.0075,
                                      f = c(10, 70),
                                      r_0 = 5.5,
                                      f_0 = 1,
                                      q_0 = 18,
                                      v_0 = 450,
                                      x_0 = 0.34,
                                      e_0 = 0.0,
                                      n_0 = 0.5,
                                      res = 100,
                                      eseis = FALSE)$power

  ## combine spectra
  psd_sum <- psd_turbulence + psd_bedload

  ## return output
  return(10 * log10(psd_sum))
}))

graphics::image(t(psd))

## invert empiric data set
X <- fmi_inversion(reference = ref_spectra, 
                   data = psd)

## plot model results
plot(X$parameters$q_s * 2650, 
     type = "l")
plot(X$parameters$h_w, 
     type = "l")

Description

In order to run the fluvial model inversion (FMI) routine, a set of randomised target parameter combinations needs to be created. This function does this job.

Usage

fmi_parameters(
  n,
  d_s,
  s_s,
  r_s,
  q_s,
  h_w,
  w_w,
  a_w,
  f_min,
  f_max,
  r_0,
  f_0,
  q_0,
  v_0,
  p_0,
  e_0,
  n_0_a,
  n_0_b,
  res
)

Arguments

Details

All parameters must be provided as single values, except for those parameters that shall be randomised, which must be provided as a vector of length two. This vector defines the range within which uniformly distributed random values will be generated and assigned.

Value

List object with model reference parameters.

Author(s)

Michael Dietze

Examples

## create two parameter sets where h_w (water level) and q_s (sediment
## flux) are randomly varied.

ref_pars <- fmi_parameters(n = 2,
                           h_w = c(0.02, 2.00),
                           q_s = c(0.001, 50.000) / 2650,
                           d_s = 0.01,
                           s_s = 1.35,
                           r_s = 2650,
                           w_w = 6,
                           a_w = 0.0075,
                           f_min = 5,
                           f_max = 80,
                           r_0 = 6,
                           f_0 = 1,
                           q_0 = 10,
                           v_0 = 350,
                           p_0 = 0.55,
                           e_0 = 0.09,
                           n_0_a = 0.6,
                           n_0_b = 0.8,
                           res = 100)

Description

In order to run the fluvial model inversion (FMI) routine, a look-up table with reference spectra needs to be created (based on predefined model parameters). This function does this job.

Usage

fmi_spectra(parameters, n_cores = 1)

Arguments

Value

List object containing the calculated reference spectra and the corresponding input parameters. Note that the spectra are given in dB for a seamless comparison with the empirical PSD data, while the original output of the models are in linear scale.

Author(s)

Michael Dietze

Examples

## create 2 example reference parameter sets
ref_pars <- fmi_parameters(n = 2,
                           h_w = c(0.02, 2.00),
                           q_s = c(0.001, 50.000) / 2650,
                           d_s = 0.01,
                           s_s = 1.35,
                           r_s = 2650,
                           w_w = 6,
                           a_w = 0.0075,
                           f_min = 5,
                           f_max = 80,
                           r_0 = 6,
                           f_0 = 1,
                           q_0 = 10,
                           v_0 = 350,
                           p_0 = 0.55,
                           e_0 = 0.09,
                           n_0_a = 0.6,
                           n_0_b = 0.8,
                           res = 100)

## create corresponding reference spectra
ref_spectra <- fmi_spectra(parameters = ref_pars)

Description

This function starts a browser-based graphic user interface to explore the parameter space of seismic models that predict the spectra of turbulent water flow and bedload flux. Empirical spectra can be added if they are present in the current environment as eseis objects. Note that those spectra should be made from deconvolved input data in order to have the same units as the model spectra (dB). ATTENTION, due to an unresolved issue, there must be at least one eseis object present in the working environment.

Usage

gui_models(...)

Arguments

Author(s)

Michael Dietze

See Also

runApp

Examples

## Not run: 
# Start the GUI
gui_models()

## End(Not run)

Description

The function returns the list of supported data loggers to extract signal deconvolution parameters.

Usage

list_logger()

Details

The value AD is the analogue-digital conversion factor.

Value

List object, supported loggers with their parameters.

Author(s)

Michael Dietze

Examples

## show documented loggers
list_logger()

## show names of loggers in list
names(list_logger())

Description

The function returns a data frame with all header values of a sac file. It may be used for advanced modifications by the user.

Usage

list_sacparameters()

Value

List object, parameters supported by a sac file.

Author(s)

Michael Dietze

Examples

## show sac parameters
list_sacparameters()

Description

The function returns the list of supported sensors to extract signal deconvolution parameters.

Usage

list_sensor()

Details

Poles and zeros must be given in rad/s. Characteristics of further sensors can be added manually. See examples of signal_deconvolve for further information. The value s is the generator constant (sensitivity) given in Vs/m. The value k is the normalisation factor of the sensor.

Value

List object, supported sensors with their parameters.

Author(s)

Michael Dietze

Examples

## show sensors
list_sensor()

Description

The function fits one of several models of signal amplitude attenuation and returns a set of model parameters, including the source amplitude (a_0).

Usage

model_amplitude(
  data,
  model = "SurfSpreadAtten",
  distance,
  source,
  d_map,
  coupling,
  v,
  q,
  f,
  k,
  a_0
)

Arguments

Details

Depending on the choice of the model to fit, several parameters can (or should) be provided, e.g. f,q, v, k, and most importantly, a_0. If more signals than free parameters are available, the missing parameters may be estimated during the fit, but without any checks of quality and meaningfulness. The parameter a_0 will be defined as 100 times the maximum input amplitude, by default. The parameters f will be set to 10 Hz, q to 50, v to 1000 m/s and k to 0.5.

ISSUES: account for non-fixed parameters, especially k

The following amplitude-distance models are available:

• "SurfSpreadAtten", Surface waves including geometric spreading and unelastic attenuation

• "BodySpreadAtten", Body waves including geometric spreading and unelastic attenuation

• "SurfBodySpreadAtten", Surface and body waves including geometric spreading and unelastic attenuation

• "SurfSpread", Surface waves including geometric spreading, only

• "BodySpread", Body waves including geometric spreading, only

• "SurfBodySpread", Surface and body waves including geometric spreading, only

**SurfSpreadAtten** The model is based on Eq. 17 from Burtin et al. (2016):

$a_d = a_0 / sqrt(d) * exp(-(pi * f * d) / (q * v))$

where a_0 is the source amplitude, a_d the amplitude as recorded by a sensor at distance d, f is the center frequency of the signal, q the ground quality factor and v the seismic wave velocity.

**BodySpreadAtten** The model is based on Eq. 16 from Burtin et al. (2016):

$a_d = a_0 / d * exp(-(pi * f * d) / (q * v))$

where a_0 is the source amplitude, a_d the amplitude as recorded by a sensor at distance d, f is the center frequency of the signal, q the ground quality factor and v the seismic wave velocity.

**SurfBodySpreadAtten** The model based on Eqs. 16 and 17 from Burtin et al. (2016):

$a_d = k * a_0 / sqrt(d) * exp(-(pi * f * d) / (q * v)) + (1 - k) * a_0 / d * exp(-(pi * f * d) / (q * v))$

where a_0 is the source amplitude, a_d the amplitude as recorded by a sensor at distance d, f is the center frequency of the signal, q the ground quality factor, v the seismic wave velocity, and n and m two factors determining the relative contributions of the two wave types, thus summing to 1.

**BodySpread** The model is simply accounting for geometric spreading

$a_d = a_0 / d$

where a_0 is the source amplitude, a_d the amplitude as recorded by a sensor at distance d.

**SurfSpread** The model is simply accounting for geometric spreading

$a_d = a_0 / sqrt(d)$

where a_0 is the source amplitude, a_d the amplitude as recorded by a sensor at distance d.

**SurfBodySpread** The model is simply accounting for geometric spreading

$a_d = k * (a_0 / d) + (1 - k) * a_d / sqrt(d)$

where a_0 is the source amplitude, a_d the amplitude as recorded by a sensor at distance d, and n and m two factors determining the relative contributions of the two wave types, thus summing to 1.

**References** - Burtin, A., Hovius, N., and Turowski, J. M.: Seismic monitoring of torrential and fluvial processes, Earth Surf. Dynam., 4, 285–307, https://doi.org/10.5194/esurf-4-285-2016, 2016.

Value

List with model results, including a_0 (source amplitude), residuals (model residuals), coefficients model coefficients.

Author(s)

Michael Dietze

Examples

## Not run: 

 ## create synthetic DEM
dem <- terra::rast(nrows = 20, ncols = 20, 
                   xmin = 0, xmax = 10000, 
                   ymin= 0, ymax = 10000, 
                   vals = rep(0, 400))

## define station coordinates
sta <- data.frame(x = c(1000, 9000, 5000, 9000),
                  y = c(1000, 1000, 9000, 9000),
                  ID = c("A", "B", "C", "D"))
 
## create synthetic signal (source in towards lower left corner of the DEM)
s <- rbind(dnorm(x = 1:1000, mean = 500, sd = 50) * 50,
           dnorm(x = 1:1000, mean = 500, sd = 50) * 2,
           dnorm(x = 1:1000, mean = 500, sd = 50) * 1,
           dnorm(x = 1:1000, mean = 500, sd = 50) * 0.5)

## calculate spatial distance maps and inter-station distances
D <- eseis::spatial_distance(stations = sta[,1:2],
                             dem = dem)

model_amplitude(data = s, 
                source = c(500, 600), 
                d_map = D$maps, 
                v = 500, 
                q = 50, 
                f = 10)

model_amplitude(data = s, 
                distance = c(254, 8254, 9280, 11667),
                model = "SurfBodySpreadAtten", 
                v = 500, 
                q = 50, 
                f = 10, 
                k = 0.5)


## End(Not run)

Description

The function calculates a seismic spectrum as predicted by the model of Tsai et al. (2012) for river bedload transport. The code was written to R by Sophie Lagarde and integrated to the R package 'eseis' by Michael Dietze.

Usage

model_bedload(
  gsd,
  d_s,
  s_s,
  r_s,
  q_s,
  h_w,
  w_w,
  a_w,
  f = c(1, 100),
  r_0,
  f_0,
  q_0,
  e_0,
  v_0,
  x_0,
  n_0,
  n_c,
  res = 100,
  adjust = TRUE,
  eseis = FALSE,
  ...
)

Arguments

Details

The model uses a set of predefined constants. These can also be changed by the user, using the ... argument:

• g = 9.81, gravitational acceleration (m/s^2)

• r_w = 1000, fluid specific density (kg/m^3)

• k_s = 3 * d_50, roughness length (m)

• log_lim = c(0.0001, 100), limits of grain-size distribution function template (m)

• log_length = 10000, length of grain-size distribution function template

• nu = 10^(-6), specific density of the fluid (kg/m^3)

• power_d = 3, grain-size power exponent

• gamma = 0.9, gamma parameter, after Parker (1990)

• s_c = 0.8, drag coefficient parameter

• s_p = 3.5, drag coefficient parameter

• c_1 = 2 / 3, inter-impact time scaling, after Sklar & Dietrich (2004)

When no user defined grain-size distribution function is provided,the function calculates the raised cosine distribution function as defined in Tsai et al. (2012) using the default range and resolution as specified by log_lim and log_length (see additional arguments list above). These default values are appropriate for mean sediment sizes between 0.001 and 10 m and log standard deivations between 0.05 and 1. When more extreme distributions are to be used, it is necessary to either adjust the arguments log_lim and log_length or use a user defined distribution function.

The adjustment option (implemented with package version 0.6.0) is only relevant for wide grain-size distributions, i.e., s_s > 0.2. In such cases, the unadjusted version tends to underestimate seismic power.

Value

eseis object containing the modelled spectrum.

Author(s)

Sophie Lagarde, Michael Dietze

References

Tsai, V. C., B. Minchew, M. P. Lamb, and J.-P. Ampuero (2012), A physical model for seismic noise generation from sediment transport in rivers, Geophys. Res. Lett., 39, L02404, doi:10.1029/2011GL050255.

Examples

## calculate spectrum (i.e., fig. 1b in Tsai et al., 2012)
p_bedload <- model_bedload(d_s = 0.7,
                           s_s = 0.1,
                           r_s = 2650,
                           q_s = 0.001,
                           h_w = 4,
                           w_w = 50,
                           a_w = 0.005,
                           f = c(0.1, 20),
                           r_0 = 600,
                           f_0 = 1,
                           q_0 = 20,
                           e_0 = 0,
                           v_0 = 1295,
                           x_0 = 0.374,
                           n_0 = 1,
                           res = 100,
                           eseis = TRUE)

## plot spectrum
plot_spectrum(data = p_bedload, 
              ylim = c(-170, -110))
              
## define empiric grain-size distribution
gsd_empiric <- data.frame(d = c(0.70, 0.82, 0.94, 1.06, 1.18, 1.30),
                          p = c(0.02, 0.25, 0.45, 0.23, 0.04, 0.00))
                  
## calculate spectrum
p_bedload <- model_bedload(gsd = gsd_empiric,
                           r_s = 2650,
                           q_s = 0.001,
                           h_w = 4,
                           w_w = 50,
                           a_w = 0.005,
                           f = c(0.1, 20),
                           r_0 = 600,
                           f_0 = 1,
                           q_0 = 20,
                           e_0 = 0,
                           v_0 = 1295,
                           x_0 = 0.374,
                           n_0 = 1,
                           res = 100,
                           eseis = TRUE)
                  
## plot spectrum
plot_spectrum(data = p_bedload, 
              ylim = c(-170, -110))
              
## define mean and sigma for parametric distribution function
d_50 <- 1
sigma <- 0.1

## define raised cosine distribution function following Tsai et al. (2012)
d_1 <- 10^seq(log10(d_50 - 5 * sigma), 
              log10(d_50 + 5 * sigma), 
              length.out = 20)

sigma_star <- sigma / sqrt(1 / 3 - 2 / pi^2)

p_1 <- (1 / (2 * sigma_star) * 
          (1 + cos(pi * (log(d_1) - log(d_50)) / sigma_star))) / d_1
p_1[log(d_1) - log(d_50) > sigma_star] <- 0
p_1[log(d_1) - log(d_50) < -sigma_star] <- 0
p_1 <- p_1 / sum(p_1)

gsd_raised_cos <- data.frame(d = d_1,
                             p = p_1)

Description

The function calculates the seismic spectrum as predicted by the model of Gimbert et al. (2014) for hydraulic turbulence. The code was written to R by Sophie Lagarde and integrated to the R package 'eseis' by Michael Dietze.

Usage

model_turbulence(
  d_s,
  s_s,
  r_s = 2650,
  h_w,
  w_w,
  a_w,
  f = c(1, 100),
  r_0,
  f_0,
  q_0,
  v_0,
  p_0,
  n_0,
  res = 1000,
  eseis = FALSE,
  ...
)

Arguments

Details

The model uses a set of predefined constants. These can also be changed by the user, using the ... argument:

• c = 0.5, instantaneous fluid-grain friction coefficient (dimensionless)

• g = 9.81, gravitational acceleration (m/s^2)

• k = 0.5, Kolmogrov constant (dimensionless)

• k_s = 3 * d_s, roughness length (m)

• h = k_s / 2, reference height of the measurement (m)

• e_0 = 0, exponent of Q increase with frequency (dimensionless)

• r_w = 1000, specific density of the fluid (kg/m^3)

• c_w = 0.5, instantaneous fluid-grain friction coefficient (dimensionless)

Value

eseis object containing the modelled spectrum.

Author(s)

Sophie Lagarde, Michael Dietze

Examples

## model the turbulence-related power spectrum
P <- model_turbulence(d_s = 0.03, # 3 cm mean grain-size
                      s_s = 1.35, # 1.35 log standard deviation
                      r_s = 2650, # 2.65 g/cm^3 sediment density
                      h_w = 0.8, # 80 cm water level
                      w_w = 40, # 40 m river width
                      a_w = 0.0075, # 0.0075 rad river inclination
                      f = c(1, 200), # 1-200 Hz frequency range
                      r_0 = 10, # 10 m distance to the river
                      f_0 = 1, # 1 Hz Null frequency 
                      q_0 = 10, # 10 quality factor at f = 1 Hz
                      v_0 = 2175, # 2175 m/s phase velocity
                      p_0 = 0.48, # 0.48 power law variation coefficient
                      n_0 = c(0.6, 0.8), # Greens function estimates
                      res = 1000) # 1000 values build the output resolution

## plot the power spectrum
plot_spectrum(data = P)

Description

The function creates a cross correlation time series of two input data sets. The input data is cut into overlapping snippets and, optionally, further smaller sub snippets for averaging the results per time snippet. The data can be made subject to a series of optional preprocessing steps, i.e. be aggregated, deconvolved, filtered, cut by standard deviation, and sign-cut. the cross correlation function is calculated and returned for a defined lag time window. The output of the function is supposed to be used as input for the function ncc_process().

Usage

ncc_correlate(
  start,
  stop,
  ID,
  component,
  dir,
  window,
  overlap = 0,
  window_sub,
  lag,
  dt,
  deconvolve = FALSE,
  f,
  pick = FALSE,
  whiten = FALSE,
  sd,
  sign = FALSE,
  cpu,
  buffer = 0.05,
  eseis = TRUE,
  ...
)

Arguments

Details

The sampling interval (dt must be defined). It is wise to set it to more than twice the filter's higher corner frequency (f[2]). Aggregation is recommended to improve computational efficiency, but is mandatory if data sets of different sampling intervals are to be analysed. In that case, it must be possible to aggregate the data sets to the provided aggregation sampling interval (dt), otherwise an error will arise. As an example, if the two data sets have sampling intervals of 1/200 and 1/500, the highest possible aggregated sampling interval is 1/100. See aux_commondt() for further information.

The function supports parallel processing. However, keep in mind that calculating the cross correlation functions for large data sets and large windows will draw serious amounts of memory. For example, a 24 h window of two seismic signals recorded at 200 Hz will easily use 15 GB of RAM. Combining this with parallel processing will multiply that memory size. Therefore, it is better think before going for too high ambitions, and check how the computer system statistics evolve with increasing windows and parallel operations.

Deconvolution is recommended if different station hardware and setup is used for the stations to analyse (i.e., different sensors, loggers or gain factors).

To account for biases due to brief events in the signals, the data sets can be truncated (cut) in their amplitude. This cutting can either be done based on the data sets' standard deviations (or a multiplicator of those standard deviations, see signal_cut() for further details), using the argument sd = 1, for one standard deviation. Alternatively, the data can also be cut by their sign (i.e., positive and negative values will be converted to 1 and -1, respectively).

Value

List with spectrogram matrix, time and frequency vectors.

Author(s)

Michael Dietze

Examples

## Not run: 

## calculate correlogram
cc <- ncc_correlate(start = "2017-04-09 00:30:00", 
                     stop = "2017-04-09 01:30:00", 
                     ID = c("RUEG1", "RUEG2"), 
                     dt = 1/10,
                     component = c("Z", "Z"), 
                     dir = paste0(system.file("extdata", 
                                              package = "eseis"), "/"), 
                     window = 600, 
                     overlap = 0, 
                     lag = 20, 
                     deconvolve = TRUE, 
                     sensor = "TC120s",
                     logger = "Cube3extBOB",
                     gain = 1,
                     f = c(0.05, 0.1), 
                     sd = 1)

## plot output
plot_correlogram(cc)


## End(Not run)

Description

The function estimates the relative seismic wave velocity changes over time based on matching iteratively stretched master correlations to previously calculated correlograms (cf. ncc_correlate).

Usage

ncc_stretch(
  data,
  lag,
  master = "mean",
  normalise = TRUE,
  sides = "both",
  range = 0.01,
  steps = 100,
  method = "r",
  reject = 0,
  ...
)

Arguments

Value

A data.frame, object with the time and relative wave velocity change estimate.

Author(s)

Michael Dietze

Examples

## Not run: 

cc <- ncc_preprocess(start = "2017-04-09 00:30:00", 
                     stop = "2017-04-09 01:30:00", 
                     ID = c("RUEG1", "RUEG2"), 
                     component = c("Z", "Z"), 
                     dir = paste0(system.file("extdata", 
                                              package = "eseis"), "/"), 
                     window = 600, 
                     overlap = 0, 
                     lag = 20, 
                     deconvolve = TRUE, 
                     sensor = "TC120s",
                     logger = "Cube3extBOB",
                     gain = 1,
                     f = c(0.05, 0.1), 
                     sd = 1)
   
   ## estimate dv/v
   dv <- ncc_stretch(data = cc, range = 0.05)
   
   ## plot result
   plot(dv$time, dv$dvv, type = "l")
                     

## End(Not run)

Description

The function picks (identifies) events from continuous data by comparing the data patterns against a template signal using Pearson's correlation coefficient, defining an event when that coefficient is above a threshold value.

Usage

pick_correlation(data, on, template, dur_min, time, dt)

Arguments

Value

data.frame, picked events.

Author(s)

Michael Dietze

Examples

## create synthetic event signal
p <- sin(seq(0, 10 * pi, by = 0.35)) * 0.2 * 
  (1 + sin(seq(0, pi, length.out = 90)))^5

## show event signal
plot(p, type = "l")

## create synthetic noise signal
x <- runif(n = 1000, min = -1, max = 1)
t <- seq(from = Sys.time(), length.out = length(x), by = 1/200)
ii <- floor(runif(n = 3, min = 100, max = 900))

## add events to noise
for(k in 1:length(ii)) {
  
  nn <- ii[k]:(ii[k] + 89)
  x[nn] <- x[nn] + p
}

## show resulting time series
plot(x = t, y = x, type = "l")

## pick events based on template
picks <- eseis::pick_correlation(data = x, 
                                 on = 0.8, 
                                 template = p, 
                                 time = t, 
                                 dt = 1/200)
                                 
## show result
print(picks)

Description

The function picks (identifies) events from continuous data using the kurtosis of the signal, and when it reaches beyond a defined threshold value. The end of an event is determined by the signal-to-noise ratio (SNR)

Usage

pick_kurtosis(
  data,
  on,
  off = 1,
  dur_min = 0,
  dur_max,
  window_kurt,
  window_amp,
  time,
  dt
)

Arguments

Details

Further reading:

Baillard, C., Crawford, W.C., Ballu, V., Hibert, C., Mangeney, A., 2014. An automatic kurtosis-based p- and s-phase picker designed for local seismic networks. Bull. Seismol. Soc. Am. 104 (1), 394–409.

Hibert, C., Mangeney, A., Grandjean, G., Baillard, C., Rivet, D., Shapiro, N.M., Satriano, C., Maggi, A., Boissier, P., Ferrazzini, V., Crawford, W., 2014. Automated identification, location, and volume estimation of rockfalls at Piton de la Fournaise Volcano. J. Geophys. Res. Earth Surf. 119 (5), 1082–1105. http://dx.doi.org/10.1002/2013JF002970.

Value

data.frame, picked events.

Author(s)

Michael Dietze

Examples

## load example data set
data(rockfall)

## preprocess signal (aggregate to increase speed, filter, envelope)
s <- signal_aggregate(data = rockfall_eseis, n = 4)
s <- signal_filter(data = s, f = c(5, 20), lazy = TRUE)
e <- signal_envelope(data = s)

## pick events based on signal kurtosis
p <- eseis::pick_kurtosis(data = e, 
                          window_kurt = 200, 
                          on = 15, 
                          off = 5, 
                          dur_min = 10, 
                          dur_max = 90, 
                          window_amp = 300)

p$picks

Description

The function calculates the ratio of the short-term-average and long-term-average of the input signal.

Usage

pick_stalta(data, time, dt, sta, lta, freeze = FALSE, on, off)

Arguments

Value

data frame, detected events (ID, start time, duration in seconds, STA-LTA vaue).

Author(s)

Michael Dietze

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

Description

The function visualises the time evolution of three seismic components of the same signal against each other as line graphs. There are three different visualisation types available: 2D (a panel of three 2D plots), 3D (a perspective threedimensional plot) and scene (an interactive threedimensional plot, mainly for exploratory purpose).

Usage

plot_components(data, type = "2D", order = "xyz", ...)

Arguments

Details

The plot type type = "3D" requires the package plot3D being installed. The plot type type = "scene" requires the package rgl being installed.

Value

A plot

Author(s)

Michael Dietze

Examples

## load example data set
data(earthquake)

## filter seismic signals
s <- eseis::signal_filter(data = s, 
                          dt = 1/200, 
                          f = c(0.05, 0.1))

## integrate signals to get displacement
s_d <- eseis::signal_integrate(data = s, dt = 1/200)

## plot components in 2D
plot_components(data = s_d, 
                type = "2D")

## plot components with time colour-coded
plot_components(data = s_d, 
                type = "2D",
                col = rainbow(n = length(s$BHE)))

## plot components with used defined coulour ramp
col_user <- colorRampPalette(colors = c("grey20", "darkblue", "blue", 
                                        "green", "red", "orange"))

plot_components(data = s_d, 
                type = "2D",
                col = col_user(n = length(s$BHE)))

## plot components as 3D plot, uncomment to use
#plot_components(data = s_d, 
#                 type = "3D",
#                 col = rainbow(n = length(s$BHE)))

Description

The function uses the output of ncc_correlate() to show an image plot of a noise cross correlation analysis.

Usage

plot_correlogram(data, agg = c(1, 1), legend = TRUE, keep_par = FALSE, ...)

Arguments

Value

Graphic output of a correlogram.

Author(s)

Michael Dietze

See Also

ncc_correlate

Examples

## Not run: 
  
  ## calculate correlogram
  cc <- ncc_correlate(start = "2017-04-09 00:30:00", 
                       stop = "2017-04-09 01:30:00", 
                       ID = c("RUEG1", "RUEG2"), 
                       component = c("Z", "Z"), 
                       dir = paste0(system.file("extdata", 
                                    package = "eseis"), "/"), 
                       window = 600, 
                       overlap = 0, 
                       lag = 20, 
                       f = c(0.05, 0.1), 
                       sd = 1)
                       
   ## explicit plot function call with adjusted resolution
   plot_correlogram(data = cc, agg = c(2, 5))
   
   ## define plot colour scale
   cls <- colorRampPalette(colors = c("brown", "white", "green"))
   
   ## simple function call with user-defined colour scale
   plot(cc, col = cls(100))

## End(Not run)

Description

The function creates from an input waveform a multi panel plot, including a seismogram, spectrum and spectrogram, and additional frequency statistics information.

Usage

plot_event(data, ratio = c(0.3, 0.3), ...)

Arguments

Details

Note that plot generation time can get long when other than short events are passed to the function. The axes limits can only be changed for the spectrum and spectrogram plots, ylim affects the frequency range, zlim affects the spectral power range.

The function uses the native plot function plot_signal(), plot_spectrum() and plot_spectrogram() along with signal_stats() to build a four panel plot.

Value

Graphic output of an event waveform.

Author(s)

Michael Dietze

Examples

## Not run: 

## load and deconvolve example event
data(rockfall)
rockfall_eseis <- signal_deconvolve(rockfall_eseis)

## plot event straight away
plot_event(data = rockfall_eseis)

## plot event with adjusted parameters
plot_event(data = rockfall_eseis, 
           ratio = c(0.4, 0.3),
           method = "periodogram",
           n = 100, 
           window = 6,
           overlap = 0.8, 
           window_sub = 4, 
           overlap_sub = 0.8,
           format ="%M:%S", 
           col = "jet",
           ylim = c(5, 80),
           zlim = c(-170, -100))

## End(Not run)

Description

The function uses the output of signal_spectrogram() to plot a probabilistic power spectral density estimate.

Usage

plot_ppsd(data, res = c(500, 500), n, ...)

Arguments

Value

Graphic output of a spectrogram.

Author(s)

Michael Dietze

See Also

signal_spectrogram

Examples

## load example data set
data(rockfall)

## deconvolve data set
r <- signal_deconvolve(data = rockfall_eseis)

## calculate PSD
p <- signal_spectrogram(data = r)

## plot PPSD
plot_ppsd(data = p$PSD)

## plot PPSD with lower resolution, more smoothing and other colour
ppsd_color <- colorRampPalette(c("white", "black", "red"))

plot_ppsd(data = p$PSD, 
          res = c(200, 200), 
          n = c(15, 20), 
          col = ppsd_color(200))

Description

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.

Usage

plot_signal(data, time, n = 10000, ...)

Arguments

Details

The 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 strptime for a list of available keywords), e.g., " "

Value

A line plot of a seismic wave form.

Author(s)

Michael Dietze

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)

Description

This function plots spectrograms of seismic signals. It uses the output of signal_spectrogram.

Usage

plot_spectrogram(data, legend = TRUE, keep_par = FALSE, agg = c(1, 1), ...)

Arguments

Details

As of version 0.7.2, the value range (zlim) is no longer set to the full data range but to the range between quantiles 0.01 and 0.99. For the full value range to be plotted, use zlim = range(data$PSD$S).

As of version 0.7.2, the default plot colour has changed from the "jet" colour palette to the "Inferno" palette. This due to perception issues with the "jet" palette. If one wants to decisively use the "jet" colours, this can be done by adding the keyword col = "jet". To use other colour schemes, such as sequential HCL schemes from the colorspace package, specify them as additional argument, e.g. col = colorspace::sequential_hcl(200, palette = "Plasma"), col = colorspace::sequential_hcl(200, palette = "Inferno"), col = colorspace::sequential_hcl(200, palette = "Viridis").

Value

Graphic output of a spectrogram.

Author(s)

Michael Dietze

See Also

signal_spectrogram

Examples

## load example data set
data(rockfall)

## deconvolve signal
rockfall <- signal_deconvolve(data = rockfall_eseis)

## calculate spectrogram
PSD <- signal_spectrogram(data = rockfall)

## plot spectrogram
plot_spectrogram(data = PSD)

## plot spectrogram with legend and labels in rainbow colours
plot_spectrogram(data = PSD, 
                 xlab = "Time (min)", 
                 ylab = "f (Hz)", 
                 main = "Power spectral density estimate", 
                 legend = TRUE, 
                 zlim = c(-220, -70),
                 col = rainbow(100)) 
                 
## plot spectrogram with frequencies in log scale
plot_spectrogram(data = PSD, log = "y")

## plot spectrogram with formatted time axis (minutes and seconds)
plot_spectrogram(data = PSD, format = "%M:%S")

Description

This function plots a line graph of the spectrum of a seismic signal.

Usage

plot_spectrum(data, unit = "dB", n = 10000, ...)

Arguments

Value

A line plot.

Author(s)

Michael Dietze

See Also

signal_spectrum

Examples

## load example data set
data(rockfall)

## calculate spectrum
spectrum_rockfall <- signal_spectrum(data = rockfall_eseis)

## plot data set with lower resolution
plot_spectrum(data = spectrum_rockfall)

Description

The function loads seismic data from a data directory structure (see aux_organisecubefiles) based on the event start time, duration, component and station ID. The data to be read needs to be adequately structured. The data directory must contain mseed or SAC files. These files will either be identified automatically or can be defined explicitly by the parameter format.

Usage

read_data(
  start,
  duration,
  station,
  component = "BHZ",
  format,
  dir,
  pattern = "eseis",
  simplify = TRUE,
  interpolate = FALSE,
  eseis = TRUE,
  try = TRUE,
  ...
)

Arguments

Details

Data organisation must follow a consistent scheme. The default scheme, eseis (Dietze, 2018 ESurf) requires hourly files organised in a directory for each Julian Day, and in each calendar year. The file name must be entirely composed of station ID, 2-digit year, Julian Day, hour, minute, second and channel name. Each item must be separated by a full stop, e.g. "2013/203/IGB01.13.203.16.00.00.BHZ" for a file from 2013, Julian Day 203, from station IGB01, covering one hour from "16:00:00 UTC", and containing the BHZ component. Each Julian Day directory can contain files from different components and stations. The respective pattern string to describe that file organisation is "%Y/%j/%STA.%y.%j.%H.%M.%S.%CMP". The percent sign indicates a wild card, where %Y is the 4-digit year, %j the 3-digit Julian Julian Day, %STA the station ID, %y the 2-digit year, %H the 2-digit hour, %M the 2-digit minute, %S the 2-digit second and %CMP the component ID. The files can have a further file extension which does not need to be explicitly defined in the pattern string. The slashes in the above pattern string define subdirectories.

An alternative organisation scheme is the one used by SeisComP, indicated by the keyword "seiscomp" or the pattern string "%Y/%NET/%STA/%CMP/%NET.%STA.%LOC.%CMP.%TYP.%Y.%j". The wild card "NET" means the network ID, "LOC" the location abbreviation and "TYP" the data type. The other wild cards are as defined above. Hence, the SeisComP scheme consists of directories of the calendar year, the network to which the data belongs, the station it has been recorded by, and the component it belongs to. The files in that latter directory must be daily files.

Value

A list object containing either a set of eseis objects or a data set with the time vector ($time) and a list of seismic stations ($station_ID) with their seismic signals as data frame ($signal). If simplify = TRUE (the default option) and only one seismic station is provided, the output object containseither just one eseis object or the vectors for $time and $signal.

Author(s)

Michael Dietze

Examples

## set seismic data directory
dir_data <- paste0(system.file("extdata", package="eseis"), "/")

## load the z component data from a station
data <- read_data(start = as.POSIXct(x = "2017-04-09 01:20:00", 
                                        tz = "UTC"), 
                      duration = 120,
                      station = "RUEG1",
                      component = "BHZ",
                      dir = dir_data)
## plot signal
plot_signal(data = data)

## load data from two stations
data <- read_data(start = as.POSIXct(x = "2017-04-09 01:20:00", 
                                     tz = "UTC"), 
                  duration = 120,
                  station = c("RUEG1", "RUEG2"),
                  component = "BHZ",
                  dir = dir_data)

## plot both signals
par(mfcol = c(2, 1))
lapply(X = data, FUN = plot_signal)

Description

The function implements download and import of seismic data from FDSN data providers via the fdsnws-dataselect service (see https://www.fdsn.org/webservices/). It is basically a wrapper for the query approach.

Usage

read_fdsn(
  start,
  duration,
  station,
  network,
  component = "BHZ",
  url,
  eseis = TRUE,
  ...
)

Arguments

Details

The FDSN (International Federation of Digital Seismograph Networks) provides access to a large number of seismic stations worldwide. The data are organised by network, station, component and further arguments. In order to use the eseis function read_fdsn, one must know at least the former three criteria for the data of interest. A list of networks is available here: https://www.fdsn.org/networks/. The function expects the 2-digit network code, the 3- or 4-digit station code, a single seismic component ID, and the URL to the data archive. Additional query arguments can be added (and must be added to point at a single seismic trace to download and import). A complete list of query arguments is available here: https://www.fdsn.org/webservices/fdsnws-dataselect-1.1.pdf.

For each network listed there, one can find the URL that gives access to the data (if existing) under "Data Access". Note that the function only requires the first URL part, e.g., https://geofon.gfz-potsdam.de.

Value

An eseis object or a list with the time ($time) and $signal vectors as well as meta information.

Author(s)

Michael Dietze

Examples

## Not run: 

## read and plot 10 min of data from Ecuador, specifying the component
s <- read_fdsn(start = "2020-05-16 22:42:00",
               duration = 360, 
               station = "IMBA", 
               network = "EC", 
               component = "HHZ")
plot(s)

## read and plot 10 min of data from Germany, specifying the URL
s <- read_fdsn(start = "2017-03-21 04:38:00",
               duration = 360, 
               station = "RGN", 
               network = "GE", 
               url = "http://geofon.gfz-potsdam.de")
plot(s)


## End(Not run)