ATaCR

Introduction

atacr is a module for the correction of vertical component data from OBS stations from tilt and compliance noise. This module is a translation of the Matlab code ATaCR and the acronym stands for Automatic Tilt and Compliance Removal. For more details on the theory and methodology, we refer the interested reader to the following papers:

  • Bell, S. W., D. W. Forsyth, and Y. Ruan (2014), Removing noise from the vertical component records of ocean-bottom seismometers: Results from year one of the Cascadia Initiative, Bull. Seismol. Soc. Am., 105, 300-313, https://doi.org/10.1785/0120140054

  • Crawford, W.C., Webb, S.C., (2000). Identifying and removing tilt noise from low-frequency (0.1 Hz) seafloor vertical seismic data, Bull. seism. Soc. Am., 90, 952-963, https://doi.org/10.1785/0119990121

  • Janiszewski, H A, J B Gaherty, G A Abers, H Gao, Z C Eilon, Amphibious surface-wave phase-velocity measurements of the Cascadia subduction zone, Geophysical Journal International, Volume 217, Issue 3, June 2019, Pages 1929-1948, https://doi.org/10.1093/gji/ggz051

The analysis can be carried out for either one (or both) compliance or tilt corrections. In all cases the analysis requires at least vertical component data. Additional data required depend on the type of analysis. The software will automatically calculate all possible corrections depending on the available channels.

Noise Corrections

Compliance

Compliance is defined as the spectral ratio between pressure and vertical displacement data. Compliance noise arises from seafloor deformation due to seafloor and water wave effects (including infragravity waves). This is likely the main source of noise in vertical component OBS data. This analysis therefore requires both vertical (?HZ) and pressure (?XH) data.

Tilt

Tilt noise arises from OBS stations that are not perfectly leveled, and therefore the horizontal seafloor deformation leaks onto the vertical component. This effect can be removed by calculating the spectral ratio between horizontal and vertical displacement data. In most cases, however, the tilt direction (measured on a compass - as opposed to tilt angle, measured from the vertical axis) is unknown and must be determined from the coherence between rotated horizontal components and the vertical component. This analysis therefore requires vertical (?HZ) and the two horizontal (?H1,2) component data.

Compliance + Tilt

It is of course possible to combine both corrections and apply them sequentially. In this case the tilt noise is removed first, followed by compliance. This analysis requires all four components: three-component seismic (?HZ,1,2) and pressure (?DH) data.

API documentation

Base Classes

atacr defines the following base classes:

The class DayNoise contains attributes and methods for the analysis of two- to four-component day-long time-series (3-component seismograms and pressure data). Objects created with this class are required in any subsequent analysis. The available methods calculate the power-spectral density (PSD) functions of sub-windows (default is 2-hour windows) and identifies windows with anomalous PSD properties. These windows are flagged and are excluded from the final averages of all possible PSD and cross-spectral density functions between all available components.

The class StaNoise contains attributes and methods for the aggregation of averaged daily spectra into a station average. An object created with this class requires that at least two DayNoise objects are available in memory. Methods available for this class are similar to those defined in the DayNoise class, but are applied to daily spectral averages, as opposed to sub-daily averages. The result is a spectral average that represents all available data for the specific station.

The class TFNoise contains attributes and methods for the calculation of transfer functions from noise traces used to correct the vertical component. A TFNoise object works with either one of DayNoise and StaNoise objects to calculate all possible transfer functions across all available components. These transfer functions are saved as attributes of the object in a Dictionary.

The class EventStream contains attributes and methods for the application of the transfer functions to the event traces for the correction (cleaning) of vertical component seismograms. An EventStream object is initialized with raw (or pre-processed) seismic and/or pressure data and needs to be processed using the same (sub) window properties as the DayNoise objects. This ensures that the component corrections are safely applied to produce corrected (cleaned) vertical components.

atacr further defines the following container classes:

These classes are used as containers for individual traces/objects that are used as attributes of the base classes.

Note

In the examples below, the SAC data were obtained and pre-processed using the accompanying scripts atacr_download_data and atacr_download_event. See the script and tutorial for details.

DayNoise

class obstools.atacr.classes.DayNoise(tr1=None, tr2=None, trZ=None, trP=None, window=7200.0, overlap=0.3, key='')

A DayNoise object contains attributes that associate three-component raw (or deconvolved) traces, metadata information and window parameters. The available methods carry out the quality control steps and the average daily spectra for windows flagged as “good”.

Note

The object is initialized with Trace objects for H1, H2, HZ and P components. Traces can be empty if data are not available. Upon saving, those traces are discarded to save disk space.

tr1, tr2, trZ, trP

Corresponding trace objects for components H1, H2, HZ and HP. Traces can be empty (i.e., Trace()) for missing components.

Type:

Trace object

window

Length of time window in seconds

Type:

float

overlap

Fraction of overlap between adjacent windows

Type:

float

key

Station key for current object

Type:

str

dt

Sampling distance in seconds. Obtained from trZ object

Type:

float

npts

Number of points in time series. Obtained from trZ object

Type:

int

fs

Sampling frequency (in Hz). Obtained from trZ object

Type:

float

year

Year for current object (obtained from UTCDateTime). Obtained from trZ object

Type:

str

julday

Julian day for current object (obtained from UTCDateTime). Obtained from trZ object

Type:

str

ncomp

Number of available components (either 2, 3 or 4). Obtained from non-empty Trace objects

Type:

int

tf_list

Dictionary of possible transfer functions given the available components.

Type:

Dict

Examples

Get demo noise data as DayNoise object

>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A

Now check its main attributes

>>> print(*[daynoise.tr1, daynoise.tr2, daynoise.trZ, daynoise.trP], sep="\n")
7D.M08A..BH1 | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
7D.M08A..BH2 | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
7D.M08A..BHZ | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
7D.M08A..BDH | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
>>> daynoise.window
7200.0
>>> daynoise.overlap
0.3
>>> daynoise.key
'7D.M08A'
>>> daynoise.ncomp
4
>>> daynoise.tf_list
{'ZP': True, 'Z1': True, 'Z2-1': True, 'ZP-21': True, 'ZH': True, 'ZP-H': True}
QC_daily_spectra(pd=[0.004, 0.2], tol=1.5, alpha=0.05, smooth=True, fig_QC=False, debug=False, save=None, form='png')

Method to determine daily time windows for which the spectra are anomalous and should be discarded in the calculation of the transfer functions.

Parameters:

  • pd (list) – Frequency corners of passband for calculating the spectra

  • tol (float) – Tolerance threshold. If spectrum > std*tol, window is flagged as bad

  • alpha (float) – Confidence interval for f-test

  • smooth (boolean) – Determines if the smoothed (True) or raw (False) spectra are used

  • fig_QC (boolean) – Whether or not to produce a figure showing the results of the quality control

  • debug (boolean) – Whether or not to plot intermediate steps in the QC procedure for debugging

  • save (Path object) – Relative path to figures folder

  • form (str) – File format (e.g., ‘png’, ‘jpg’, ‘eps’)

ftX

Windowed Fourier transform for the X component (can be either 1, 2, Z or P)

Type:

ndarray

f

Full frequency axis (Hz)

Type:

ndarray

goodwins

List of booleans representing whether a window is good (True) or not (False)

Type:

list

Examples

Perform QC on DayNoise object using default values and plot final figure

>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra(fig_QC=True)
_images/Figure_3a.png

Print out new attribute of DayNoise object

>>> daynoise.goodwins
array([False,  True,  True,  True,  True,  True,  True,  True, False,
   False,  True,  True,  True,  True,  True,  True], dtype=bool)
average_daily_spectra(calc_rotation=True, fig_average=False, fig_coh_ph=False, save=None, form='png')

Method to average the daily spectra for good windows. By default, the method will attempt to calculate the azimuth of maximum coherence between horizontal components and the vertical component (for maximum tilt direction), and use the rotated horizontals in the transfer function calculations.

Parameters:

  • calc_rotation (boolean) – Whether or not to calculate the tilt direction

  • fig_average (boolean) – Whether or not to produce a figure showing the average daily spectra

  • fig_coh_ph (boolean) – Whether or not to produce a figure showing the maximum coherence between H and Z

  • save (Path object) – Relative path to figures folder

  • form (str) – File format (e.g., ‘png’, ‘jpg’, ‘eps’)

power

Container for the Power spectra

Type:

Power

cross

Container for the Cross power spectra

Type:

Cross

rotation

Container for the Rotated power and cross spectra

Type:

Rotation, optional

Examples

Average spectra for good windows using default values and plot final figure

>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra(fig_average=True)
_images/Figure_3b.png

Print out available attributes of DayNoise object

>>> daynoise.__dict__.keys()
dict_keys(['tr1', 'tr2', 'trZ', 'trP', 'window', 'overlap', 'key',
'dt', 'npts', 'fs', 'year', 'julday', 'tkey', 'ncomp', 'tf_list',
'QC', 'av', 'f', 'goodwins', 'power', 'cross', 'rotation'])
save(filename)

Method to save the object to file using ~Pickle.

Parameters:

filename (str) – File name

Examples

Run demo through all methods

>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()

Save object

>>> daynoise.save('daynoise_demo.pkl')

Check that it has been saved

>>> import glob
>>> glob.glob("./daynoise_demo.pkl")
['./daynoise_demo.pkl']

StaNoise

class obstools.atacr.classes.StaNoise(daylist=None)

A StaNoise object contains attributes that associate three-component raw (or deconvolved) traces, metadata information and window parameters.

Note

The object is initially a container for DayNoise objects. Once the StaNoise object is initialized (using the method init() or by calling the QC_sta_spectra() method), each individual spectral quantity is unpacked as an object attribute and the original DayNoise objects are removed from memory. DayNoise objects cannot be added or appended to the StaNoise object once this is done. In addition, all spectral quantities associated with the original DayNoise objects (now stored as attributes) are discarded as the object is saved to disk and new container objects are defined and saved.

daylist

A list of DayNoise objects to process and produce a station average

Type:

list

initialized

Whether or not the object has been initialized - False unless one of the methods have been called. When True, the daylist attribute is deleted from memory

Type:

bool

Examples

Initialize empty object

>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise()

Initialize with DayNoise object

>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> stanoise = StaNoise(daylist=[daynoise])

Add or append DayNoise object to StaNoise

>>> stanoise = StaNoise()
>>> stanoise += daynoise

>>> stanoise = StaNoise()
>>> stanoise.append(daynoise)

Import demo noise data with 4 DayNoise objects

>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.daylist
[<obstools.atacr.classes.DayNoise at 0x11e3ce8d0>,
 <obstools.atacr.classes.DayNoise at 0x121c7ae10>,
 <obstools.atacr.classes.DayNoise at 0x121ca5940>,
 <obstools.atacr.classes.DayNoise at 0x121e7dd30>]
 >>> stanoise.initialized
 False
init()

Method to initialize the StaNoise object. This method is used to unpack the spectral quantities from the original DayNoise objects and allow the methods to proceed. The original DayNoise objects are deleted from memory during this process.

Note

If the original DayNoise objects have not been processed using their QC and averaging methods, these will be called first before unpacking into the object attributes.

f

Frequency axis for corresponding time sampling parameters

Type:

ndarray

nwins

Number of good windows from the DayNoise object

Type:

int

key

Station key for current object

Type:

str

ncomp

Number of available components (either 2, 3 or 4)

Type:

int

tf_list

Dictionary of possible transfer functions given the available components.

Type:

Dict

c11

Power spectra for component H1. Other identical attributes are available for the power, cross and rotated spectra: [11, 12, 1Z, 1P, 22, 2Z, 2P, ZZ, ZP, PP, HH, HZ, HP]

Type:

numpy.ndarray

direc

Array of azimuths used in determining the tilt direction

Type:

numpy.ndarray

tilt

Tilt direction from maximum coherence between rotated H1 and HZ components

Type:

float

QC

Whether or not the method QC_sta_spectra() has been called.

Type:

bool

av

Whether or not the method average_sta_spectra() has been called.

Type:

bool

Examples

Initialize demo data

>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.init()

Check that daylist attribute has been deleted

>>> stanoise.daylist
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'StaNoise' object has no attribute 'daylist'
>>> stanoise.__dict__.keys()
dict_keys(['initialized', 'c11', 'c22', 'cZZ', 'cPP', 'c12', 'c1Z',
'c1P', 'c2Z', 'c2P', 'cZP', 'cHH', 'cHZ', 'cHP', 'direc', 'tilt', 'f',
'nwins', 'ncomp', 'key', 'tf_list', 'QC', 'av'])
QC_sta_spectra(pd=[0.004, 0.2], tol=2.0, alpha=0.05, fig_QC=False, debug=False, save=None, form='png')

Method to determine the days (for given time window) for which the spectra are anomalous and should be discarded in the calculation of the long-term transfer functions.

Parameters:

  • pd (list) – Frequency corners of passband for calculating the spectra

  • tol (float) – Tolerance threshold. If spectrum > std*tol, window is flagged as bad

  • alpha (float) – Confidence interval for f-test

  • fig_QC (boolean) – Whether or not to produce a figure showing the results of the quality control

  • debug (boolean) – Whether or not to plot intermediate steps in the QC procedure for debugging

  • save (Path object) – Relative path to figures folder

  • form (str) – File format (e.g., ‘png’, ‘jpg’, ‘eps’)

gooddays

List of booleans representing whether a day is good (True) or not (False)

Type:

list

Examples

Import demo data, call method and generate final figure

>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra(fig_QC=True)
>>> stanoise.QC
True
average_sta_spectra(fig_average=False, save=None, form='png')

Method to average the daily station spectra for good windows.

Parameters:

  • fig_average (boolean) – Whether or not to produce a figure showing the average daily spectra

  • save (Path object) – Relative path to figures folder

  • form (str) – File format (e.g., ‘png’, ‘jpg’, ‘eps’)

power

Container for the Power spectra

Type:

Power

cross

Container for the Cross power spectra

Type:

Cross

rotation

Container for the Rotated power and cross spectra

Type:

Cross, optional

Examples

Average daily spectra for good days using default values and produce final figure

>>> obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
save(filename)

Method to save the object to file using ~Pickle.

Parameters:

filename (str) – File name

Examples

Run demo through all methods

>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()

Save object

>>> stanoise.save('stanoise_demo.pkl')

Check that it has been saved

>>> import glob
>>> glob.glob("./stanoise_demo.pkl")
['./stanoise_demo.pkl']

TFNoise

class obstools.atacr.classes.TFNoise(objnoise=None)

A TFNoise object contains attributes that store the transfer function information from multiple components (and component combinations).

Note

The object is initialized with either a processed DayNoise or StaNoise object. Each individual spectral quantity is unpacked as an object attribute, but all of them are discarded as the object is saved to disk and new container objects are defined and saved.

f

Frequency axis for corresponding time sampling parameters

Type:

ndarray

c11

Power spectra for component H1. Other identical attributes are available for the power, cross and rotated spectra: [11, 12, 1Z, 1P, 22, 2Z, 2P, ZZ, ZP, PP, HH, HZ, HP]

Type:

numpy.ndarray

tilt

Tilt direction from maximum coherence between rotated H1 and HZ components

Type:

float

tf_list

Dictionary of possible transfer functions given the available components.

Type:

Dict

Examples

Initialize a TFNoise object with a DayNoise object. The DayNoise object must be processed for QC and averaging, otherwise the TFNoise object will not initialize.

>>> from obstools.atacr import DayNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> tfnoise = TFNoise(daynoise)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/pascalaudet/Softwares/Python/Projects/dev/OBStools/obstools/atacr/classes.py", line 1215, in __init__
Exception: Error: Noise object has not been processed (QC and averaging) - aborting

Now re-initialized with a processed DayNoise object

>>> from obstools.atacr import DayNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise = TFNoise(daynoise)

Initialize a TFNoise object with a processed StaNoise object

>>> from obstools.atacr import StaNoise, TFNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise = TFNoise(stanoise)
class TfDict
transfer_func()

Method to calculate transfer functions between multiple components (and component combinations) from the averaged (daily or station-averaged) noise spectra.

transfunc

Container Dictionary for all possible transfer functions

Type:

TfDict

Examples

Calculate transfer functions for a DayNoise object

>>> from obstools.atacr import DayNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise = TFNoise(daynoise)
>>> tfnoise.transfer_func()
>>> tfnoise.transfunc.keys()
dict_keys(['ZP', 'Z1', 'Z2-1', 'ZP-21', 'ZH', 'ZP-H'])

Calculate transfer functions for a StaNoise object

>>> from obstools.atacr import StaNoise, TFNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise = TFNoise(stanoise)
>>> tfnoise.transfer_func()
>>> tfnoise.transfunc.keys()
dict_keys(['ZP', 'Z1', 'Z2-1', 'ZP-21'])
save(filename)

Method to save the object to file using ~Pickle.

Parameters:

filename (str) – File name

Examples

Run demo through all methods

>>> from obstools.atacr import DayNoise, StaNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise_day = TFNoise(daynoise)
>>> tfnoise_day.transfer_func()
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise_sta = TFNoise(stanoise)
>>> tfnoise_sta.transfer_func()

Save object

>>> tfnoise_day.save('tf_daynoise_demo.pkl')
>>> tfnoise_sta.save('tf_stanoise_demo.pkl')

Check that everything has been saved

>>> import glob
>>> glob.glob("./tf_daynoise_demo.pkl")
['./tf_daynoise_demo.pkl']
>>> glob.glob("./tf_stanoise_demo.pkl")
['./tf_stanoise_demo.pkl']

EventStream

class obstools.atacr.classes.EventStream(tr1=<obspy.core.trace.Trace object>, tr2=<obspy.core.trace.Trace object>, trZ=<obspy.core.trace.Trace object>, trP=<obspy.core.trace.Trace object>, correct=False)

An EventStream object contains attributes that store station-event metadata and methods for applying the transfer functions to the various components and produce corrected/cleaned vertical components.

Note

The object is initialized with Trace objects for H1, H2, HZ and P components. Traces can be empty if data are not available. Based on the available components, a list of possible corrections is determined automatically.

tr1, tr2, trZ, trP

Corresponding trace objects for components H1, H2, HZ and HP. Traces can be empty (i.e., Trace()) for missing components.

Type:

Trace object

key

Station key for current object

Type:

str

evtime

Origin time of seismic event.

Type:

UTCDateTime

tstamp

Time stamp for event

Type:

str

prefix

Associated prefix of SAC files

Type:

str

npts

Number of points in time series.

Type:

int

fs

Sampling frequency (in Hz).

Type:

float

dt

Sampling distance in seconds.

Type:

float

ncomp

Number of available components (either 2, 3 or 4). Obtained from non-empty Trace objects

Type:

int

ev_list

Dictionary of possible transfer functions given the available components. This is determined during initialization.

Type:

Dict

Examples

Get demo earthquake data as EventStream object

>>> from obstools.atacr import EventStream
>>> evstream = EventStream('demo')
Uploading demo earthquake data - March 09, 2012, station 7D.M08A
>>> evstream.__dict__.keys()
dict_keys(['tr1', 'tr2', 'trZ', 'trP, 'key', 'evtime',
'tstamp', 'prefix', 'npts', fs', 'dt', 'ncomp', 'ev_list'])

Plot the raw traces

>>> import obstools.atacr.plotting as atplot
>>> figure = atplot.fig_event_raw(evstream, fmin=1./150., fmax=2.)
>>> figure.show()
_images/Figure_11.png
class CorrectDict
correct_data(tfnoise)

Method to apply transfer functions between multiple components (and component combinations) to produce corrected/cleaned vertical components.

Parameters:

tfnoise (TFNoise) – Object that contains the noise transfer functions used in the correction

correct

Container Dictionary for all possible corrections from the transfer functions

Type:

CorrectDict

Examples

Let’s carry through the correction of the vertical component for a single day of noise, say corresponding to the noise recorded on March 04, 2012. In practice, the DayNoise object should correspond to the same day at that of the recorded earthquake to avoid bias in the correction.

>>> from obstools.atacr import DayNoise, TFNoise, EventStream
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise_day = TFNoise(daynoise)
>>> tfnoise_day.transfer_func()
>>> evstream = EventStream('demo')
Uploading demo earthquake data - March 09, 2012, station 7D.M08A
>>> evstream.correct_data(tfnoise_day)

Plot the corrected traces

>>> import obstools.atacr.plotting as atplot
>>> figure = atplot.fig_event_corrected(evstream, tfnoise_day.tf_list)
>>> figure.show()
_images/Figure_corrected_march04.png

Carry out the same exercise but this time using a StaNoise object

>>> from obstools.atacr import StaNoise, TFNoise, EventStream
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise_sta = TFNoise(stanoise)
>>> tfnoise_sta.transfer_func()
>>> evstream = EventStream('demo')
Uploading demo earthquake data - March 09, 2012, station 7D.M08A
>>> evstream.correct_data(tfnoise_sta)

Plot the corrected traces

>>> import obstools.atacr.plot as plot
>>> plot.fig_event_corrected(evstream, tfnoise_sta.tf_list)
_images/Figure_corrected_sta.png

Warning

If the noise window and event window are not identical, they cannot be compared on the same frequency axis and the code will exit. Make sure you are using identical time samples in both the noise and event windows.

save(filename)

Method to save the object to file using ~Pickle.

Parameters:

filename (str) – File name

Examples

Following from the example outlined in method correct_data(), we simply save the final object

>>> evstream.save('evstream_demo.pkl')

Check that object has been saved

>>> import glob
>>> glob.glob("./evstream_demo.pkl")
['./evstream_demo.pkl']

Container Classes

Power

class obstools.atacr.classes.Power(c11=None, c22=None, cZZ=None, cPP=None)

Container for power spectra for each component, with any shape

c11

Power spectral density for component 1 (any shape)

Type:

ndarray

c22

Power spectral density for component 2 (any shape)

Type:

ndarray

cZZ

Power spectral density for component Z (any shape)

Type:

ndarray

cPP

Power spectral density for component P (any shape)

Type:

ndarray

Cross

class obstools.atacr.classes.Cross(c12=None, c1Z=None, c1P=None, c2Z=None, c2P=None, cZP=None)

Container for cross-power spectra for each component pairs, with any shape

c12

Cross-power spectral density for components 1 and 2 (any shape)

Type:

ndarray

c1Z

Cross-power spectral density for components 1 and Z (any shape)

Type:

ndarray

c1P

Cross-power spectral density for components 1 and P (any shape)

Type:

ndarray

c2Z

Cross-power spectral density for components 2 and Z (any shape)

Type:

ndarray

c2P

Cross-power spectral density for components 2 and P (any shape)

Type:

ndarray

cZP

Cross-power spectral density for components Z and P (any shape)

Type:

ndarray

Rotation

class obstools.atacr.classes.Rotation(cHH=None, cHZ=None, cHP=None, coh=None, ph=None, tilt=None, coh_value=None, phase_value=None, direc=None)

Container for rotated spectra, with any shape

cHH

Power spectral density for rotated horizontal component H (any shape)

Type:

ndarray

cHZ

Cross-power spectral density for components H and Z (any shape)

Type:

ndarray

cHP

Cross-power spectral density for components H and P (any shape)

Type:

ndarray

coh

Coherence between horizontal components

Type:

ndarray

ph

Phase of cross-power spectrum between horizontal components

Type:

ndarray

tilt

Angle (azimuth) of tilt axis

Type:

float

coh_value

Maximum coherence

Type:

float

phase_value

Phase at maximum coherence

Type:

float

direc

Directions for which the coherence is calculated

Type:

ndarray

Utility functions

utils contains several functions that are used in the class methods of ~obstools.atacr.classes.

obstools.atacr.utils.traceshift(trace, tt)

Function to shift traces in time given travel time

Parameters:

  • trace (Trace object) – Trace object to update

  • tt (float) – Time shift in seconds

Returns:

rtrace – Updated trace object

Return type:

Trace object

obstools.atacr.utils.QC_streams(start, end, st)

Function for quality control of traces, which compares the start and end times that were requested, as well as the total n length of the traces.

Parameters:

  • start (UTCDateTime object) – Start time of requested stream

  • end (UTCDateTime object) – End time of requested stream

  • st (Stream object) – Stream object with all trace data

Returns:

  • (pass) (bool) – Whether the QC test has passed

  • st (Stream object) – Updated stream object

obstools.atacr.utils.update_stats(tr, stla, stlo, stel, cha, evla=None, evlo=None)

Function to include SAC metadata to Trace objects

Parameters:

  • tr (Trace object) – Trace object to update

  • stla (float) – Latitude of station

  • stlo (float) – Longitude of station

  • stel (float) – Station elevation (m)

  • cha (str) – Channel for component

  • evla (float, optional) – Latitude of event

  • evlo (float, optional) – Longitute of event

Returns:

tr – Updated trace object

Return type:

Trace object

obstools.atacr.utils.get_data(datapath, tstart, tend)

Function to grab all available noise data given a path and data time range

Parameters:

  • datapath (str) – Path to noise data folder

  • tstart (UTCDateTime) – Start time for query

  • tend (UTCDateTime) – End time for query

Returns:

tr1, tr2, trZ, trP – Corresponding trace objects for components H1, H2, HZ and HP. Returns empty traces for missing components.

Return type:

Trace object

obstools.atacr.utils.get_event(eventpath, tstart, tend)

Function to grab all available earthquake data given a path and data time range

Parameters:

  • eventpath (str) – Path to earthquake data folder

  • tstart (UTCDateTime) – Start time for query

  • tend (UTCDateTime) – End time for query

Returns:

tr1, tr2, trZ, trP – Corresponding trace objects for components H1, H2, HZ and HP. Returns empty traces for missing components.

Return type:

Trace object

obstools.atacr.utils.calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreq=[0.005, 0.035])

Determines tilt direction from maximum coherence between rotated H1 and Z.

Parameters:

  • ft1 (ndarray) – Fourier transform of corresponding H1, H2, HZ and HP components

  • ft2 (ndarray) – Fourier transform of corresponding H1, H2, HZ and HP components

  • ftZ (ndarray) – Fourier transform of corresponding H1, H2, HZ and HP components

  • ftP (ndarray) – Fourier transform of corresponding H1, H2, HZ and HP components

  • f (ndarray) – Frequency axis in Hz

  • goodwins (list) – List of booleans representing whether a window is good (True) or not (False). This attribute is returned from the method QC_daily_spectra()

  • tiltfreq (list, optional) – Two floats representing the frequency band at which the tilt is calculated

Returns:

  • cHH, cHZ, cHP (ndarray) – Arrays of power and cross-spectral density functions of components HH (rotated H1 in direction of maximum tilt), HZ, and HP

  • coh (ndarray) – Coherence value between rotated H and Z components, as a function of directions (azimuths)

  • ph (ndarray) – Phase value between rotated H and Z components, as a function of directions (azimuths)

  • direc (ndarray) – Array of directions (azimuths) considered

  • tilt (float) – Direction (azimuth) of maximum coherence between rotated H1 and Z

  • coh_value (float) – Coherence value at tilt direction

  • phase_value (float) – Phase value at tilt direction

obstools.atacr.utils.smooth(data, nd, axis=0)

Function to smooth power spectral density functions from the convolution of a boxcar function with the PSD

Parameters:

  • data (ndarray) – Real-valued array to smooth (PSD)

  • nd (int) – Number of samples over which to smooth

  • axis (int, optional) – axis over which to perform the smoothing

Returns:

filt – Filtered data

Return type:

ndarray, optional

obstools.atacr.utils.admittance(Gxy, Gxx)

Calculates admittance between two components

Parameters:

  • Gxy (ndarray) – Cross spectral density function of x and y

  • Gxx (ndarray) – Power spectral density function of x

Returns:

Admittance between x and y

Return type:

ndarray, optional

obstools.atacr.utils.coherence(Gxy, Gxx, Gyy)

Calculates coherence between two components

Parameters:

  • Gxy (ndarray) – Cross spectral density function of x and y

  • Gxx (ndarray) – Power spectral density function of x

  • Gyy (ndarray) – Power spectral density function of y

Returns:

Coherence between x and y

Return type:

ndarray, optional

obstools.atacr.utils.phase(Gxy)

Calculates phase angle between two components

Parameters:

Gxy (ndarray) – Cross spectral density function of x and y

Returns:

Phase angle between x and y

Return type:

ndarray, optional

Plotting functions

plot contains several functions for plotting the results of the analysis at various final and intermediate steps.

obstools.atacr.plotting.fig_QC(f, power, gooddays, ncomp, key='')

Function to plot the Quality-Control step of the analysis. This function is used in both the atacr_daily_spectra or atacr_clean_spectra scripts.

Parameters:

  • f (ndarray) – Frequency axis (in Hz)

  • power (Power) – Container for the Power spectra

  • gooddays (List) – List of booleans representing whether a window is good (True) or not (False)

  • ncomp (int) – Number of components used in analysis (can be 2, 3 or 4)

  • key (str) – String corresponding to the station key under analysis

obstools.atacr.plotting.fig_average(f, power, bad, gooddays, ncomp, key='')

Function to plot the averaged spectra (those qualified as ‘good’ in the QC step). This function is used in both the atacr_daily_spectra or atacr_clean_spectra scripts.

Parameters:

  • f (ndarray) – Frequency axis (in Hz)

  • power (Power) – Container for the Power spectra

  • bad (Power) – Container for the bad Power spectra

  • gooddays (List) – List of booleans representing whether a window is good (True) or not (False)

  • ncomp (int) – Number of components used in analysis (can be 2, 3 or 4)

  • key (str) – String corresponding to the station key under analysis

obstools.atacr.plotting.fig_av_cross(f, field, gooddays, ftype, ncomp, key='', save=False, fname='', form='png', **kwargs)

Function to plot the averaged cross-spectra (those qualified as ‘good’ in the QC step). This function is used in the atacr_daily_spectra script.

Parameters:

  • f (ndarray) – Frequency axis (in Hz)

  • field (Rotation) – Container for the Power spectra

  • gooddays (List) – List of booleans representing whether a window is good (True) or not (False)

  • ftype (str) – Type of plot to be displayed. If ftype is Admittance, plot is loglog. Otherwise semilogx

  • key (str) – String corresponding to the station key under analysis

  • **kwargs (None) – Keyword arguments to modify plot

obstools.atacr.plotting.fig_coh_ph(coh, ph, direc)

Function to plot the coherence and phase between the rotated H and Z components, used to characterize the tilt direction.

Parameters:

  • coh (ndarray) – Coherence between rotated H and Z components

  • ph (ndarray) – Phase between rotated H and Z components

  • direc (ndarray) – Directions considered in maximizing coherence between H and Z

obstools.atacr.plotting.fig_TF(f, day_trfs, day_list, sta_trfs, sta_list, skey='')

Function to plot the transfer functions available.

Parameters:

  • f (ndarray) – Frequency axis (in Hz)

  • day_trfs (Dict) – Dictionary containing the transfer functions for the daily averages

  • day_list (Dict) – Dictionary containing the list of daily transfer functions

  • sta_trfs (Dict) – Dictionary containing the transfer functions for the station averages

  • sta_list (Dict) – Dictionary containing the list of average transfer functions

  • skey (str) – String corresponding to the station key under analysis

obstools.atacr.plotting.fig_comply(f, day_comps, day_list, sta_comps, sta_list, skey=None, elev=-1000.0, f_0=None)

Function to plot the transfer functions available.

Parameters:

  • f (ndarray) – Frequency axis (in Hz)

  • day_comps (Dict) – Dictionary containing the compliance functions for the daily averages

  • day_list (Dict) – Dictionary containing the list of daily transfer functions

  • sta_comps (Dict) – Dictionary containing the compliance functions for the station averages

  • sta_list (Dict) – Dictionary containing the list of average transfer functions

  • skey (str) – String corresponding to the station key under analysis

  • elev (float) – Station elevation in meters (OBS stations have negative elevations)

  • f_0 (float) – Lowest frequency to consider in plot (Hz)

obstools.atacr.plotting.fig_event_raw(evstream, fmin=0.006666666666666667, fmax=2.0)

Function to plot the raw (although bandpassed) seismograms.

Parameters:

  • evstream (EventStream) – Container for the event stream data

  • fmin (float) – Low frequency corner (in Hz)

  • fmax (float) – High frequency corner (in Hz)

obstools.atacr.plotting.fig_event_corrected(evstream, TF_list, fmin=0.006666666666666667, fmax=2.0)

Function to plot the corrected vertical component seismograms.

Parameters:

  • evstream (EventStream) – Container for the event stream data

  • Tf_list (list) – List of Dictionary elements of transfer functions used for plotting the corrected vertical component.

Scripts

There are several Python scripts that accompany ~obstools.atacr. These can be used in bash scripts to automate data processing. These include scripts to download noise and event data, and perform tilt and compliance noise removal using either the default program values or by refining parameters. All of them use a station database provided as a StDb dictionary. These scripts are:

  • atacr_download_data

  • atacr_download_event

  • atacr_daily_spectra

  • atacr_clean_spectra

  • atacr_transfer_functions

  • atacr_correct_event

atacr_download_data

Description

Downloads up to four-component (H1, H2, Z and P), day-long seismograms to use in noise corrections of vertical component data. Station selection is specified by a network and station code. The database is provided as a StDb dictionary.

Usage

$ atacr_download_data -h
usage: atacr_download_data [options] <indb>

Script used to download and pre-process up to four-component (H1, H2, Z and
P), day-long seismograms to use in noise corrections of vertical component of
OBS data. Data are requested from the internet using the client services
framework for a given date range. The stations are processed one by one and
the data are stored to disk.

positional arguments:
  indb                  Station Database to process from.

optional arguments:
  -h, --help            show this help message and exit
  --keys STKEYS         Specify a comma-separated list of station keys for
                        which to perform the analysis. These must be contained
                        within the station database. Partial keys will be used
                        to match against those in the dictionary. For
                        instance, providing IU will match with all stations in
                        the IU network. [Default processes all stations in the
                        database]
  -C CHANNELS, --channels CHANNELS
                        Specify a comma-separated list of channels for which
                        to perform the transfer function analysis. Possible
                        options are '12' (for horizontal channels) or 'P' (for
                        pressure channel). Specifying '12' allows for tilt
                        correction. Specifying 'P' allows for compliance
                        correction. [Default looks for both horizontal and
                        pressure and allows for both tilt AND compliance
                        corrections]
  --zcomp ZCOMP         Specify the Vertical Component Channel Identifier.
                        [Default Z].
  -O, --overwrite       Force the overwriting of pre-existing data. [Default
                        False]

Server Settings:
  Settings associated with which datacenter to log into.

  --server SERVER       Base URL of FDSN web service compatible server (e.g.
                        “http://service.iris.edu”) or key string for
                        recognized server (one of 'AUSPASS', 'BGR',
                        'EARTHSCOPE', 'EIDA', 'EMSC', 'ETH', 'GEOFON',
                        'GEONET', 'GFZ', 'ICGC', 'IESDMC', 'INGV', 'IPGP',
                        'IRIS', 'IRISPH5', 'ISC', 'KNMI', 'KOERI', 'LMU',
                        'NCEDC', 'NIEP', 'NOA', 'NRCAN', 'ODC', 'ORFEUS',
                        'RASPISHAKE', 'RESIF', 'RESIFPH5', 'SCEDC', 'TEXNET',
                        'UIB-NORSAR', 'USGS', 'USP'). [Default 'IRIS']
  --user-auth USERAUTH  Authentification Username and Password for the
                        waveform server (--user-auth='username:authpassword')
                        to access and download restricted data. [Default no
                        user and password]
  --eida-token TOKENFILE
                        Token for EIDA authentication mechanism, see
                        http://geofon.gfz-
                        potsdam.de/waveform/archive/auth/index.php. If a token
                        is provided, argument --user-auth will be ignored.
                        This mechanism is only available on select EIDA nodes.
                        The token can be provided in form of the PGP message
                        as a string, or the filename of a local file with the
                        PGP message in it. [Default None]

Frequency Settings:
  Miscellaneous frequency settings

  --sampling-rate NEW_SAMPLING_RATE
                        Specify new sampling rate (float, in Hz). [Default 5.]
  --units UNITS         Choose the output seismogram units. Options are:
                        'DISP', 'VEL', 'ACC'. [Default 'DISP']
  --pre-filt PRE_FILT   Specify four comma-separated corner frequencies
                        (float, in Hz) for deconvolution pre-filter. [Default
                        0.001,0.005,45.,50.]

Time Search Settings:
  Time settings associated with searching for day-long seismograms

  --start STARTT        Specify a UTCDateTime compatible string representing
                        the start day for the data search. This will override
                        any station start times. [Default start date for each
                        station in database]
  --end ENDT            Specify a UTCDateTime compatible string representing
                        the start time for the event search. This will
                        override any station end times [Default end date for
                        each station in database]

atacr_daily_spectra

Description

Extracts two-hour-long windows from the day-long data, calculates power-spectral densities and flags windows for outlier from the PSD properties. Station selection is specified by a network and station code. The database is provided as a StDb dictionary.

Usage

$ atacr_daily_spectra -h
usage: atacr_daily_spectra [options] <indb>

Script used to extract shorter windows from the day-long seismograms,
calculate the power-spectral properties, flag windows for outlier PSDs and
calculate daily averages of the corresponding Fourier transforms. The stations
are processed one by one and the data are stored to disk. The program will
look for data saved in the previous steps and use all available components.

positional arguments:
  indb               Station Database to process from

optional arguments:
  -h, --help         show this help message and exit
  --keys STKEYS      Specify a comma separated list of station keys for which
                     to perform the analysis. These must be contained within
                     the station database. Partial keys will be used to match
                     against those in the dictionary. For instance, providing
                     IU will match with all stations in the IU network.
                     [Default processes all stations in the database]
  -O, --overwrite    Force the overwriting of pre-existing data. [Default
                     False]

Time Search Settings:
  Time settings associated with searching for day-long seismograms

  --start STARTT     Specify a UTCDateTime compatible string representing the
                     start day for the data search. This will override any
                     station start times. [Default start date of each station
                     in database]
  --end ENDT         Specify a UTCDateTime compatible string representing the
                     start time for the data search. This will override any
                     station end times. [Default end date of each station n
                     database]

Parameter Settings:
  Miscellaneous default values and settings

  --window WINDOW    Specify window length in seconds. Default value is highly
                     recommended. Program may not be stable for large
                     deviations from default value. [Default 7200. (or 2
                     hours)]
  --overlap OVERLAP  Specify fraction of overlap between windows. [Default 0.3
                     (or 30 percent)]
  --minwin MINWIN    Specify minimum number of 'good' windows in any given day
                     to continue with analysis. [Default 10]
  --freq-band PD     Specify comma-separated frequency limits (float, in Hz)
                     over which to calculate spectral features used in
                     flagging the bad windows. [Default 0.004,2.0]
  --tolerance TOL    Specify parameter for tolerance threshold. If spectrum >
                     std*tol, window is flagged as bad. [Default 2.0]
  --alpha ALPHA      Specify confidence level for f-test, for iterative
                     flagging of windows. [Default 0.05, or 95 percent
                     confidence]
  --raw              Raw spectra will be used in calculating spectral features
                     for flagging. [Default uses smoothed spectra]
  --no-rotation      Do not rotate horizontal components to tilt direction.
                     [Default calculates rotation]

Figure Settings:
  Flags for plotting figures

  --figQC            Plot Quality-Control figure. [Default does not plot
                     figure]
  --debug            Plot intermediate steps for debugging. [Default does not
                     plot figure]
  --figAverage       Plot daily average figure. [Default does not plot figure]
  --figCoh           Plot Coherence and Phase figure. [Default does not plot
                     figure]
  --save-fig         Set this option if you wish to save the figure(s).
                     [Default does not save figure]
  --format FORM      Specify format of figure. Can be any one of the
                     validmatplotlib formats: 'png', 'jpg', 'eps', 'pdf'.
                     [Default 'png']

atacr_clean_spectra

Description

Extracts daily spectra calculated from atacr_daily_spectra and flags days for which the daily averages are outliers from the PSD properties. Further averages the spectra over the whole period specified by --start and --end. Station selection is specified by a network and station code. The database is provided as a StDb dictionary.

Usage

$ atacr_clean_spectra -h
usage: atacr_clean_spectra [options] <indb>

Script used to extract daily spectra calculated from `obs_daily_spectra.py`
and flag days for outlier PSDs and calculate spectral averages of the
corresponding Fourier transforms over the entire time period specified. The
stations are processed one by one and the data are stored to disk.

positional arguments:
  indb             Station Database to process from.

optional arguments:
  -h, --help       show this help message and exit
  --keys STKEYS    Specify a comma separated list of station keys for which to
                   perform the analysis. These must be contained within the
                   station database. Partial keys will be used to match
                   against those in the dictionary. For instance, providing IU
                   will match with all stations in the IU network. [Default
                   processes all stations in the database]
  -O, --overwrite  Force the overwriting of pre-existing data. [Default False]

Time Search Settings:
  Time settings associated with searching for day-long seismograms

  --start STARTT   Specify a UTCDateTime compatible string representing the
                   start day for the data search. This will override any
                   station start times. [Default start date of each station in
                   database]
  --end ENDT       Specify a UTCDateTime compatible string representing the
                   start time for the data search. This will override any
                   station end times. [Default end date of each station in
                   database]

Parameter Settings:
  Miscellaneous default values and settings

  --freq-band PD   Specify comma-separated frequency limits (float, in Hz)
                   over which to calculate spectral features used in flagging
                   the days/windows. [Default 0.004,2.0]
  --tolerance TOL  Specify parameter for tolerance threshold. If spectrum >
                   std*tol, window is flagged as bad. [Default 1.5]
  --alpha ALPHA    Confidence level for f-test, for iterative flagging of
                   windows. [Default 0.05, or 95 percent confidence]

Figure Settings:
  Flags for plotting figures

  --figQC          Plot Quality-Control figure. [Default does not plot figure]
  --debug          Plot intermediate steps for debugging. [Default does not
                   plot figure]
  --figAverage     Plot daily average figure. [Default does not plot figure]
  --figCoh         Plot Coherence and Phase figure. [Default does not plot
                   figure]
  --figCross       Plot cross-spectra figure. [Default does not plot figure]
  --save-fig       Set this option if you wish to save the figure(s). [Default
                   does not save figure]
  --format FORM    Specify format of figure. Can be any one of the
                   validmatplotlib formats: 'png', 'jpg', 'eps', 'pdf'.
                   [Default 'png']

atacr_transfer functions

Description

Calculates transfer functions using the noise windows flagged as good, for either a single day (from atacr_daily_spectra) or for those averaged over several days (from atacr_clean_spectra), if available. The transfer functions are stored to disk. Station selection is specified by a network and station code. The database is provided as a StDb dictionary.

Usage

$ atacr_transfer_functions -h
usage: atacr_transfer_functions [options] <indb>

Script used to calculate transfer functions between various components, to be
used in cleaning vertical component of OBS data. The noise data can be those
obtained from the daily spectra (i.e., from `obs_daily_spectra.py`) or those
obtained from the averaged noise spectra (i.e., from `obs_clean_spectra.py`).
Flags are available to specify the source of data to use as well as the time
range over which to calculate the transfer functions. The stations are
processed one by one and the data are stored to disk.

positional arguments:
  indb             Station Database to process from.

optional arguments:
  -h, --help       show this help message and exit
  --keys STKEYS    Specify a comma separated list of station keys for which to
                   perform the analysis. These must be contained within the
                   station database. Partial keys will be used to match
                   against those in the dictionary. For instance, providing IU
                   will match with all stations in the IU network. [Default
                   processes all stations in the database]
  -O, --overwrite  Force the overwriting of pre-existing data. [Default False]

Time Search Settings:
  Time settings associated with searching for day-long seismograms

  --start STARTT   Specify a UTCDateTime compatible string representing the
                   start day for the data search. This will override any
                   station start times. [Default start date of each station in
                   database]
  --end ENDT       Specify a UTCDateTime compatible string representing the
                   start time for the data search. This will override any
                   station end times. [Default end date of each station in
                   database]

Parameter Settings:
  Miscellaneous default values and settings

  --skip-daily     Skip daily spectral averages in construction of transfer
                   functions. [Default False]
  --skip-clean     Skip cleaned spectral averages in construction of transfer
                   functions. Defaults to True if data cannot be found in
                   default directory. [Default False]

Figure Settings:
  Flags for plotting figures

  --figTF          Plot transfer function figure. [Default does not plot
                   figure]
  --save-fig       Set this option if you wish to save the figure(s). [Default
                   does not save figure]
  --format FORM    Specify format of figure. Can be any one of the
                   validmatplotlib formats: 'png', 'jpg', 'eps', 'pdf'.
                   [Default 'png']

atacr_download_event

Description

Downloads up to four-component (H1, H2, Z and P), two-hour-long seismograms for individual seismic events to use in noise corrections of vertical component data. Station selection is specified by a network and station code. The database is provided as a StDb dictionary.

Usage

$ atacr_download_event -h
usage: atacr_download_event [options] <indb>

Script used to download and pre-process four-component (H1, H2, Z and P), two-
hour-long seismograms for individual events on which to apply the de-noising
algorithms. Data are requested from the internet using the client services
framework for a given date range. The stations are processed one by one and
the data are stored to disk.

positional arguments:
  indb                  Station Database to process from.

optional arguments:
  -h, --help            show this help message and exit
  --keys STKEYS         Specify a comma separated list of station keys for
                        which to perform the analysis. These must be contained
                        within the station database. Partial keys will be used
                        to match against those in the dictionary. For
                        instance, providing IU will match with all stations in
                        the IU network [Default processes all stations in the
                        database]
  -C CHANNELS, --channels CHANNELS
                        Specify a comma-separated list of channels for which
                        to perform the transfer function analysis. Possible
                        options are '12' (for horizontal channels) or 'P' (for
                        pressure channel). Specifying '12' allows for tilt
                        correction. Specifying 'P' allows for compliance
                        correction. [Default looks for both horizontal and
                        pressure and allows for both tilt AND compliance
                        corrections]
  --zcomp ZCOMP         Specify the Vertical Component Channel Identifier.
                        [Default Z].
  -O, --overwrite       Force the overwriting of pre-existing data. [Default
                        False]

Server Settings:
  Settings associated with which datacenter to log into.

  --server SERVER       Base URL of FDSN web service compatible server (e.g.
                        “http://service.iris.edu”) or key string for
                        recognized server (one of 'AUSPASS', 'BGR',
                        'EARTHSCOPE', 'EIDA', 'EMSC', 'ETH', 'GEOFON',
                        'GEONET', 'GFZ', 'ICGC', 'IESDMC', 'INGV', 'IPGP',
                        'IRIS', 'IRISPH5', 'ISC', 'KNMI', 'KOERI', 'LMU',
                        'NCEDC', 'NIEP', 'NOA', 'NRCAN', 'ODC', 'ORFEUS',
                        'RASPISHAKE', 'RESIF', 'RESIFPH5', 'SCEDC', 'TEXNET',
                        'UIB-NORSAR', 'USGS', 'USP'). [Default 'IRIS']
  --user-auth USERAUTH  Authentification Username and Password for the
                        waveform server (--user-auth='username:authpassword')
                        to access and download restricted data. [Default no
                        user and password]
  --eida-token TOKENFILE
                        Token for EIDA authentication mechanism, see
                        http://geofon.gfz-
                        potsdam.de/waveform/archive/auth/index.php. If a token
                        is provided, argument --user-auth will be ignored.
                        This mechanism is only available on select EIDA nodes.
                        The token can be provided in form of the PGP message
                        as a string, or the filename of a local file with the
                        PGP message in it. [Default None]

Frequency Settings:
  Miscellaneous frequency settings

  --sampling-rate NEW_SAMPLING_RATE
                        Specify new sampling rate (float, in Hz). [Default 5.]
  --units UNITS         Choose the output seismogram units. Options are:
                        'DISP', 'VEL', 'ACC'. [Default 'DISP']
  --pre-filt PRE_FILT   Specify four comma-separated corner frequencies
                        (float, in Hz) for deconvolution pre-filter. [Default
                        0.001,0.005,45.,50.]
  --window WINDOW       Specify window length in seconds. Default value is
                        highly recommended. Program may not be stable for
                        large deviations from default value. [Default 7200.
                        (or 2 hours)]

Event Settings:
  Settings associated with refining the events to include in matching
  station pairs

  --start STARTT        Specify a UTCDateTime compatible string representing
                        the start time for the event search. This will
                        override any station start times. [Default start date
                        of each station in database]
  --end ENDT            Specify a UTCDateTime compatible string representing
                        the start time for the event search. This will
                        override any station end times [Default end date of
                        each station in database]
  --reverse-order, -R   Reverse order of events. Default behaviour starts at
                        oldest event and works towards most recent. Specify
                        reverse order and instead the program will start with
                        the most recent events and work towards older
  --min-mag MINMAG      Specify the minimum magnitude of event for which to
                        search. [Default 5.5]
  --max-mag MAXMAG      Specify the maximum magnitude of event for which to
                        search. [Default None, i.e. no limit]

Geometry Settings:
  Settings associatd with the event-station geometries

  --min-dist MINDIST    Specify the minimum great circle distance (degrees)
                        between the station and event. [Default 30]
  --max-dist MAXDIST    Specify the maximum great circle distance (degrees)
                        between the station and event. [Default 120]

atacr_correct_event

Description

Loads the transfer functions previously calculated and performs the noise corrections, for either a single day (from atacr_daily_spectra) or for those averaged over several days (from atacr_clean_spectra), if available. The database is provided as a StDb dictionary.

Usage

$ atacr_correct_event -h
usage: atacr_correct_event [options] <indb>

Script used to extract transfer functions between various components, and use
them to clean vertical component of OBS data for selected events. The noise
data can be those obtained from the daily spectra (i.e., from
`obs_daily_spectra.py`) or those obtained from the averaged noise spectra
(i.e., from `obs_clean_spectra.py`). Flags are available to specify the source
of data to use as well as the time range for given events. The stations are
processed one by one and the data are stored to disk in a new 'CORRECTED'
folder.

positional arguments:
  indb             Station Database to process from.

optional arguments:
  -h, --help       show this help message and exit
  --keys STKEYS    Specify a comma separated list of station keys for which to
                   perform the analysis. These must be contained within the
                   station database. Partial keys will be used to match
                   against those in the dictionary. For instance, providing IU
                   will match with all stations in the IU network. [Default
                   processes all stations in the database]
  -O, --overwrite  Force the overwriting of pre-existing data. [Default False]

Time Search Settings:
  Time settings associated with searching for specific event-related
  seismograms

  --start STARTT   Specify a UTCDateTime compatible string representing the
                   start day for the event search. This will override any
                   station start times. [Default start date of each station in
                   database]
  --end ENDT       Specify a UTCDateTime compatible string representing the
                   start time for the event search. This will override any
                   station end times. [Default end date of each station in
                   database]

Parameter Settings:
  Miscellaneous default values and settings

  --skip-daily     Skip daily spectral averages in application of transfer
                   functions. [Default False]
  --skip-clean     Skip cleaned spectral averages in application of transfer
                   functions. [Default False]
  --fmin FMIN      Low frequency corner (in Hz) for plotting the raw (un-
                   corrected) seismograms. Filter is a 2nd order, zero phase
                   butterworth filter. [Default 1./150.]
  --fmax FMAX      High frequency corner (in Hz) for plotting the raw (un-
                   corrected) seismograms. Filter is a 2nd order, zero phase
                   butterworth filter. [Default 1./10.]

Figure Settings:
  Flags for plotting figures

  --figRaw         Plot raw seismogram figure. [Default does not plot figure]
  --figClean       Plot cleaned vertical seismogram figure. [Default does not
                   plot figure]
  --save-fig       Set this option if you wish to save the figure(s). [Default
                   does not save figure]
  --format FORM    Specify format of figure. Can be any one of the
                   validmatplotlib formats: 'png', 'jpg', 'eps', 'pdf'.
                   [Default 'png']

Tutorial

Note

Here we roughly follow the steps highlighted in the Matlab tutorial for this code and reproduce the various figures. The examples provided below are for one month of data (March 2012) recorded at station M08A of the Cascadia Initiative Experiment. Corrections are applied to a magnitude 6.6 earthquake that occurred near Vanuatu on March 9, 2012.

0. Creating the StDb Database

All the scripts provided require a StDb database containing station information and metadata. Let’s first create this database for station M08A and send the prompt to a logfile

$ query_fdsn_stdb -N 7D -C ?H? -S M08A M08A > logfile

Note

If you are using a Z shell (as opposed to a Bourne Shell or other types of shells) on the terminal, this command above will fail due to the presence of the question marks for pattern matching. To find out, just type ps -p $$, which will return zsh under the CMD field if you are using a Z shell. In this case, just enclose the ?H? in single or double quotes (e.g., query_fdsn_stdb.py -N 7D -C “?H?” -S M08A M08A > logfile)

To check the station info for M08A, use the program ls_stdb:

$ ls_stdb M08A.pkl
 Listing Station Pickle: M08A.pkl
 7D.M08A
 --------------------------------------------------------------------------
 1) 7D.M08A
      Station: 7D M08A
       Alternate Networks: None
       Channel: BH ;  Location: --
       Lon, Lat, Elev:  44.11870, -124.89530,  -0.126
       StartTime: 2011-10-20 00:00:00
       EndTime:   2012-07-18 23:59:59
       Status:    open
       Polarity: 1
       Azimuth Correction: 0.000000

1. Download noise data

We wish to download one month of data for the station M08A for March 2012. The various options above allow us to select the additional channels to use (e.g., -C 12,P for both horizontal and pressure data - which is the default setting). Default frequency settings for data pre-processing match those of the Matlab ATaCR software and can therefore be ignored when calling the program. Since the file M08A.pkl contains only one station, it is not necessary to specify a key. This option would be useful if the database contained several stations and we were only interested in downloading data for M08A. In this case, we would specify --keys=M08A or --keys=7D.M08A. The only required options at this point are the --start and --end options to specify the dates for which data will be downloaded.

If you change your mind about the pre-processing options, you can always re-run the following line with the option -O, which will overwrite the data saved to disk.

To download all broadband seismic and pressure data, simply type in a terminal:

$ atacr_download_data --start=2012-03-01 --end=2012-04-01 M08A.pkl

An example log printed on the terminal will look like:

Path to DATA/7D.M08A/ doesn`t exist - creating it

|===============================================|
|===============================================|
|                       M08A                    |
|===============================================|
|===============================================|
|  Station: 7D.M08A                             |
|      Channel: BH; Locations: --               |
|      Lon: -124.90; Lat:  44.12                |
|      Start time: 2011-10-20                   |
|      End time:   2012-07-18                   |
|-----------------------------------------------|
| Searching day-long files:                     |
|   Start: 2012-03-01                           |
|   End:   2012-04-01                           |

***********************************************************
* Downloading day-long data for key 7D.M08A and day 2012.61
*
* Channels selected: ['12', 'P'] and vertical
*   2012.061.*SAC
*   -> Downloading Seismic data...
*      ...done
*   -> Downloading Pressure data...
*      ...done
*   -> Removing responses - Seismic data
 WARNING: FIR normalized: sum[coef]=9.767192E-01;
 WARNING: FIR normalized: sum[coef]=9.767192E-01;
 WARNING: FIR normalized: sum[coef]=9.767192E-01;
*   -> Removing responses - Pressure data
 WARNING: FIR normalized: sum[coef]=9.767192E-01;

***********************************************************
* Downloading day-long data for key 7D.M08A and day 2012.62
*
* Channels selected: ['12', 'P'] and vertical
*   2012.062.*SAC
*   -> Downloading Seismic data...

...

And so on until all day-long files have been downloaded. You will notice that a folder called DATA/7D.M08A/ has been created. This is where all day-long files will be stored on disk.

2. QC for daily spectral averages

For this step, there are several Parameter Settings that can be tuned. Once again, the default values are the ones used to reproduce the results of the Matlab ATaCR software and can be left un-changed. The Time Search Settings can be used to look at a subset of the available day-long data files. Here these options can be ignored since we wish to look at all the availble data that we just downloaded. We can therefore type in a terminal:

$ atacr_daily_spectra M08A.pkl

Path to SPECTRA/7D.M08A/ doesn`t exist - creating it

|===============================================|
|===============================================|
|                       M08A                    |
|===============================================|
|===============================================|
|  Station: 7D.M08A                             |
|      Channel: BH; Locations: --               |
|      Lon: -124.90; Lat:  44.12                |
|      Start time: 2011-10-20 00:00:00          |
|      End time:   2012-07-18 23:59:59          |
|-----------------------------------------------|

**********************************************************************
* Calculating noise spectra for key 7D.M08A and day 2012.061
*   12 good windows. Proceeding...

**********************************************************************
* Calculating noise spectra for key 7D.M08A and day 2012.062
*   14 good windows. Proceeding...

**********************************************************************
* Calculating noise spectra for key 7D.M08A and day 2012.063
*   16 good windows. Proceeding...

...

And so on until all available data have been processed. The software stores the obstools.atacr.classes.DayNoise objects to a newly created folder called SPECTRA/7D.M08A/. To produce figures for visualization, we can re-run the above script but now use the plotting options to look at one day of the month (March 04, 2012). In this case we need to overwrite the previous results (option -O) and specify the date of interest:

$ atacr_daily_spectra -O --figQC --figAverage --start=2012-03-04 --end=2012-03-05 M08A.pkl > logfile

The script will produce several figures, including Figures 2 and 3 (separated into 3a and 3b below). Several intermediate steps are also produces, which show all the raw data and the window classification into good and bad windows for subsequent analysis.

_images/Figure_2.png

Figure 2: Daily spectrogram for the vertical (Z), horizontals (H1, H2), and pressure (P) components.

_images/Figure_3a.png

Figure 3a: Power spectral density (PSD) functions for the Z, H1, H2, and P compo- nents from a single day of data (M08A, March 4, 2012, same as in Figure 2). The left column shows PSDs for each individual window; PSDs from windows that did not pass the quality control are colored red.

_images/Figure_3b.png

Figure 3b: Daily average PSD of bad (red) and good (black) windows.

3. QC for clean station averages

Now that we have processed daily spectra for all available components, it is possible to further average the spectra over multiple days to produce a cleaned station average. It is still possible to specify a date range over which to average the spectra, thus giving flexibility in the production of the station averages. Parameter settings are similar to those used in atacr_daily_spectra but further include the option of plotting the averaged cross-spectral properties. To calculate a single station average for the entire month of March 2012 (and therefore using all available data) and plot the results, we can type in a terminal:

$ atacr_clean_spectra --figQC --figAverage --figCoh --figCross M08A.pkl

Path to AVG_STA/7D.M08A/ doesn`t exist - creating it

|===============================================|
|===============================================|
|                       M08A                    |
|===============================================|
|===============================================|
|  Station: 7D.M08A                             |
|      Channel: BH; Locations: --               |
|      Lon: -124.90; Lat:  44.12                |
|      Start time: 2011-10-20 00:00:00          |
|      End time:   2012-07-18 23:59:59          |
|-----------------------------------------------|

**********************************************************************
* Calculating noise spectra for key 7D.M08A and day 2012.061
*   -> file SPECTRA/7D.M08A/2012.061.spectra.pkl found - loading

**********************************************************************
* Calculating noise spectra for key 7D.M08A and day 2012.062
*   -> file SPECTRA/7D.M08A/2012.062.spectra.pkl found - loading

...

And so on until all DayNoise objects are averaged into a StaNoise object, which is saved to a newly created folder called AVG_STA/7D.M08A/.

Note

If you don’t specify the options --start and --end, the object will be saved with a filename that corresponds to the entire deployment time of the station, but in fact the object contains the average spectra of all daily spectra available on disk, and not necessarily the average over the entire deployment time. We recommend using the --start and --end options if you want to produce time-limited spectral averages (e.g., an average per week or per month, etc.). For example:

$ atacr_clean_spectra --start=2012-03-01 --end=2012-03-08 -O M08A.pkl

Several figures are also produced, including Figures 4, 6-9.

_images/Figure_4.png

Figure 4: The orientation of maximum coherence between the vertical and the two horizontal components for M08A during March 2012. (Left) Coherence as a function of angle from the H1 component. (Right) Phase as a function of the angle. In this example, the coherence is low indicating the absence of dominant, uni-directional tilt noise.

_images/Figure_6.png

Figure 6: The daily PSDs plotted for the vertical (Z), horizontal (H1, H2), and pressure (P) components for March 2012 at station M08A. Each line is a daily PSD. Gray colors indicate days that were accepted by the second quality control step, while the red colors indicate days that were discarded.

_images/Figure_7.png

Figure 7: The daily coherences between pairs of components as indicated above each subplot for March 2012 at station M08A (e.g. 1Z - coherence between Horizontal 1 and Vertical). Each line represents average coherence for a single day in the data set. Gray colors indicate days that were accepted by the second quality control step, while the red colors indicate days that were discarded. Note that some definitions are different (e.g., ZP vs PZ) than those used in the Matlab Tutorial.

_images/Figure_8.png

Figure 8: The daily admittances between pairs of components as indicated above (see Figure 7) each subplot for March 2012 at station M08A. Each line is a daily admittance. Gray colors indicate days that were accepted by the second quality control step, while the red colors indicate days that were discarded. Note that some definitions are different (e.g., ZP vs PZ) than those used in the Matlab Tutorial.

_images/Figure_9.png

Figure 9: The daily phases between pairs of components as indicated above (see Figure 7) each subplot for March 2012 at station M08A. Each line is a daily phase. Gray colors indicate days that were accepted by the second quality control step, while the red colors indicate days that were discarded. Note that some definitions are different (e.g., ZP vs PZ) than those used in the Matlab Tutorial.

4. Transfer function calculation

Once the StaNoise objects have been produced and saved to disk, the transfer functions across all available components can be calculated. By default the software will calculate the ones for which the spectral averages are available.

For compliance only (i.e., only ?HZ and ?DH? components are available), the only transfer function possible is:

  • ZP

For tilt only (i.e., all of ?HZ,1,2 components are available, but not ?DH), the transfer functions are:

  • Z1

  • Z2-1

For both tilt and compliance (i.e., all four components are available), the following transfer functions are possible:

  • Z1

  • Z2-1

  • ZP

  • ZP-21

If you are using a DayNoise object to calculate the transfer functions, the following may also be possible (if all components are available):

  • ZH

  • ZP-H

In this example we calculate all available transfer functions for all available data. In this case we do not need to specify any option and type in a terminal:

$ atacr_transfer_functions M08A.pkl

Path to TF_STA/7D.M08A/ doesn't exist - creating it

|===============================================|
|===============================================|
|                       M08A                    |
|===============================================|
|===============================================|
|  Station: 7D.M08A                             |
|      Channel: BH; Locations: --               |
|      Lon: -124.90; Lat:  44.12                |
|      Start time: 2011-10-20 00:00:00          |
|      End time:   2012-07-18 23:59:59          |
|-----------------------------------------------|

**********************************************************************
* Calculating transfer functions for key 7D.M08A and day 2012.088

**********************************************************************
* Calculating transfer functions for key 7D.M08A and day 2012.075

...

**********************************************************************
* Calculating transfer functions for key 7D.M08A and range 2011.293-2012.200.

Note how the DayNoise objects are read randomly from disk, followed by the StaNoise object. The result is a TFNoise object that is saved to a newly created folder called TF_STA/7D.M08A/.

We can produce Figure 10 by re-running the previous command with the options -O --figTF.

_images/Figure_10.png

Figure 10: Transfer function amplitudes for the component combinations of interest, as indicated in the title of each subplot. This example is for the month of March 2012 for station M08A. The daily transfer functions are shown in grey and the average calculated for the whole month is shown in black.

5. Download earthquake data

Now we need to download the earthquake data, for which we wish to clean the vertical component using the transfer functions just calculated. This script atacr_download_event is very similar to atacr_download_data, with the addition of the Event and Geometry Settings.

Warning

Be careful with the Frequency Settings, as these values need to be exactly the same as those used in atacr_download_data, but won’t be checked against.

To download the seismograms that recorded the March 9, 2012, magnitude 6.6 Vanuatu earthquake (be conservative with the options), type in a terminal:

$ atacr_download_event --min-mag=6.3 --max-mag=6.7 --start=2012-03-08 --end=2012-03-10 M08A.pkl

Path to EVENTS/7D.M08A/ doesn`t exist - creating it

|===============================================|
|===============================================|
|                       M08A                    |
|===============================================|
|===============================================|
|  Station: 7D.M08A                             |
|      Channel: BH; Locations: --               |
|      Lon: -124.90; Lat:  44.12                |
|      Start time: 2011-10-20 00:00:00          |
|      End time:   2012-07-18 23:59:59          |
|-----------------------------------------------|
| Searching Possible events:                    |
|   Start: 2012-03-08 00:00:00                  |
|   End:   2012-03-10 00:00:00                  |
|   Mag:   6.3 - 6.7                            |
| ...                                           |
|  Found     1 possible events                  |

****************************************************
* (1/1):  20120309_070953
*   Origin Time: 2012-03-09 07:09:53
*   Lat: -19.22; Lon:  169.75
*   Dep:  33.70; Mag: 6.6
*     M08A  -> Ev: 9651.91 km;   86.80 deg; 239.43;  40.95

* Channels selected: ['12', 'P'] and vertical
*   2012.069.07.09
*   -> Downloading Seismic data...
*      ...done
*   -> Downloading Pressure data...
     ...done
*   -> Removing responses - Seismic data
 WARNING: FIR normalized: sum[coef]=9.767192E-01;
 WARNING: FIR normalized: sum[coef]=9.767192E-01;
 WARNING: FIR normalized: sum[coef]=9.767192E-01;
*   -> Removing responses - Pressure data
 WARNING: FIR normalized: sum[coef]=9.767192E-01;

The data are stored as an EventStream object, saved to disk in the newly created folder EVENTS/7D.M08A/.

6. Correct/clean earthquake data

The final step in the analysis is the application of the transfer functions to the raw earthquake seismograms to clean up the vertical component. Once again, the default settings can be used. The corrected seismograms will be saved to disk in a new folder called EVENTS/7D.M08A/CORRECTED. To make the final Figures 11 and 12, specify the --figRaw and --figClean options:

$ atacr_correct_event --figRaw --figClean M08A.pkl


|===============================================|
|===============================================|
|                       M08A                    |
|===============================================|
|===============================================|
|  Station: 7D.M08A                             |
|      Channel: BH; Locations: --               |
|      Lon: -124.90; Lat:  44.12                |
|      Start time: 2011-10-20 00:00:00          |
|      End time:   2012-07-18 23:59:59          |
|-----------------------------------------------|
TF_STA/7D.M08A/2011.293-2012.200.transfunc.pkl file found - applying transfer functions
TF_STA/7D.M08A/2012.069.transfunc.pkl file found - applying transfer functions

Results are saved both as updated EventStream objects and as .SAC files that now contain the corrected vertical components.

Note

The EventStream object contains the corrected traces in the attribute correct, which is a dictionary. If you wish to extract the corrected vertical seismogram obtained from the ZP transfer function, with a corresponding key 'ZP', you can do this by loading the .pkl file and typing:

>>> # First load the EventStream object
>>> import pickle as pkl
>>> evstream = pkl.load(open('2012.069.transfunc.pkl', 'rb'))
>>> # Extract the 'ZP' corrected seismogram
>>> corrected = evstream.correct['ZP']

You can check which transfer functions should be available based on the number of components in the EventStream object from the attribute:

>>> print(evstream.ev_list)

Note that the corrected seismogram will only be available if the corresponding transfer function has been calculated previously and applied during the correction.

Figures 11 and 12 should look like:

_images/Figure_11.png

Figure 11: Event time series for the vertical (Z), horizontal 1 (H1), horizontal 2 (H2), and pressure (P) components. No corrections have been applied. The data is for station M08A for the Mw 6.6 earthquake that occurred near Vanuatu on March 9, 2012 and has been bandpass filtered from 10 - 150 s.

_images/Figure_12a.png

Figure 12a: Event time series for the vertical (Z) components after each of the transfer functions of interest have been applied. The corrections are specified in the titles of each subplot. The data is for station M08A for the Mw 6.6 earthquake that occurred near Vanuatu on March 9, 2012 (same as Figure 11) and has been bandpass filtered from 10 - 150 s. Traces in grey show the original (raw, un-corrected) vertical component.

_images/Figure_12b.png

Figure 12b: Same as Figure 12a but using the station averaged transfer functions. In this case the ZH and ZP-H transfer functions are not available.