WaveLoc API

This page contains the documentation on each of the modules / classes / functions that make up the waveloc package. The modules are ordered by subject type. This is documentation is automatically generated from the doc-strings in the code itself. It is meant as an aide to developpers or advanced users who need access to the intrigacies of the code. It is not a tutorial (one will be coming soon).

Waveloc Options

Waveloc is piloted entirely through a set of options / parameters that are set in a WavelocOptions object. All options will (soon!) be explained on the tutorial pages.

options

class options.WavelocOptions[source]

The WavelocOptions class contains a single attribute opdict containing all the options and parameters to control the behaviour of waveloc. It also contains many methods to verify the consistency of the options before launching the various steps of the waveloc procedure.

set_test_options()[source]

Sets up a standard version of opdict. Used notably in the example cases.

verify_SDS_processing_options()[source]

Verify presence of all options necessary for SDS_processing.

verify_base_path()[source]

Verifies that the base_path is set. If the ‘base_path’ option is not directly set in the opdict, its value is read from the environment variable $WAVELOC_PATH.

verify_cluster_options()[source]

Verify presence of all options necessary for clustering.

verify_correlation_options()[source]

Verify presence of all options necessary for correlation.

verify_doublediff_options()[source]

Verify presence of all options necessary for double difference location.

verify_kurtogram_options()[source]

Verify presence of all options necessary for kurtogram.

verify_location_options()[source]

Verify presence of all options necessary for location.

verify_magnitude_options()[source]

Verify presence of all options necessary for magnitude.

verify_migration_options()[source]

Verify presence of all options necessary for migration.

verify_plotting_options()[source]

Verify presence of all options necessary for plotting waveloc results.

verify_probloc_plotting_options()[source]

Verify presence of all options necessary for plotting waveloc using the probability density location method.

verify_synthetic_options()[source]

Verify presence of all options necessary for synthetic migration runs.

Travel-times

Waveloc needs travel-times to migrate the kurtosis waveforms. These can in principle be obtained in many different ways. For simple Earth models, a convenient way to obtain travel-time grids is to take advantage of the NonLinLoc implementation of the eikonal solver of Podvin & Lecomte (1991).

NllGridLib

The NllGridLib module contains functions to read the output of NonLinLoc (NLL) in order to use the travel-times within waveloc. Contributed by Alessia Maggi

NllGridLib.latlon2rect(proj_name, lat, lon, proj_info={})[source]

Projects latitude and longitude coordinates into an x-y coordinate system. Takes the projection from the NLL header info as read by read_hdr_file. The inverse projection is provided by rect2latlon.

Parameters:
  • proj_name (string) – ‘TRANS_GLOBAL’ | ‘TRANS_NONE’ | ‘TRANS_SIMPLE’ Only these projection types are supported.
  • lat (float) – Latitude in decimal degrees (positive = North, negative = South)
  • lon (float) – Longitude in decimal degrees (positive = East, negative = West)
  • proj_info (dictionary, optional) – NLL header info as read by read_hdr_file.
Raises UserWarning:
 

if unsupported projection type is encountered.

Return type:

tuple of floats

Returns:

x and y coordinates in km wrt the origin of the projection

NllGridLib.qd_read_hyp_file(filename)[source]

Quick and dirty routine to read a NLL hypocenter file.

Parameters:filename (string) – File to read
Return type:tuple of floats
Returns:(otime, hypo_x, sigma_x, hypo_y, sigma_y, hypo_z, sigma_z)
NllGridLib.qd_read_picks_from_hyp_file(filename)[source]

Quick and dirty routine to read P-wave picks from NLL hypocenter file.

Parameters:filename (string) – File to read
Return type:dictionary
Returns:A dictionary where the keys are the station names and the values are the P-arrival times as UTCDateTime objects.
NllGridLib.read_hdr_file(filename)[source]

Reads a NLL .hdr file.

Parameters:filename (string) – File to read
Returns:Dictionary containing the information in the file. If filename is a grid geometry header, the second line will be the projection info. If filename is a travel-time grid header, there will also be a line with the station information. The dictionary keys for the geometry are ‘nx’, ‘ny’, ‘nz’, ‘x_orig’, ‘y_orig’, ‘z_orig’; those for the projection are ‘proj_name’, ‘orig_lat’, ‘orig_lon’, ‘map_rot’; those for the station information are ‘station’, ‘sta_x’, ‘sta_y’, ‘sta_z’.
NllGridLib.read_stations_file(filename)[source]

Reads a file formatted according the the NLL station format.

Parameters:filename (string) – File to be read
Returns:Dictionary of stations. Each station is itself a dictionary containing the following information: ‘station’ = station name; ‘loc_type’ = [‘XYZ’ | ‘LATLON’]; if ‘loc_type’ = ‘XYZ’, then the coordinates of the station are given by ‘x’ and ‘y’, otherwise they are ‘lat’, ‘lon’; the depth and elevation of the station are given respectively by ‘depth’ and ‘elevation’.
Raises UserWarning:
 if unknown ‘loc_type’ encoutered.
NllGridLib.rect2latlon(proj_name, x, y, proj_info={})[source]

Projects x and y coordinates into an latitude-longitude coordinate system. Takes the projection from the NLL header info as read by read_hdr_file. The inverse projection is provided by latlon2rect.

Parameters:
  • proj_name (string) – ‘TRANS_GLOBAL’ | ‘TRANS_NONE’ | ‘TRANS_SIMPLE’ Only these projection types are supported.
  • x (float) – x-coordinate in km wrt to the origin of the projection
  • y (float) – y-coordinate in km wrt to the origin of the projection
  • proj_info (dictionary, optional) – NLL header info as read by read_hdr_file.
Raises UserWarning:
 

if unsupported projection type is encountered.

Return type:

tuple of floats

Returns:

latitude and longitude (latitude positive North, longitude positive East)

hdf5_grids

Waveloc manipulates many large arrays, some of which are written / read from disk. On some machines, these arrays are larger than the available RAM. In order to avoid the drammatic slow-downs when swapping to disc occurs in these cases, we use hdf (hierarchical data formats) that provide memory-like access to files on disc.

All large arrays that are infrequently accessed (e.g. travel-time arrays) are written and used in this manner. When there is sufficient RAM available, large arrays that are frequently accessed (search-point travel-times, 4D-migration grids) are kept in memory.

class hdf5_grids.H5NllSingleGrid(filename, nll_filename)[source]

Class derived from H5SingleGrid that upon initialization reads a NLL pair of files (.hdr and .buf) and writes an equivalent hdf5 file.

Parameters:
  • nll_name – NLL filename without extension
  • h5_name – hdf5 filename
class hdf5_grids.H5SingleGrid(filename=None, grid_data=None, grid_info=None)[source]

Class that wraps several methods for regular grids in HDF5 format.

Attributes

grid_data

An HDF5 dataset.

grid_info

Attributes of an HDF5 dataset.

Methods

interp_to_newgrid(new_filename, new_grid_info)[source]

Intrepolates to a new grid.

Parameters:
  • new_filename – Filename for new grid
  • new_grid_info – Grid info dictionary for new grid
value_at_point(x, y, z, epsilon=0.001)[source]

Performs n-linear interpolation on the regular grid.

Parameters:
  • x (float) – x-coordinate of the point of interest
  • y (float) – y-coordinate of the point of interest
  • z (float) – z-coordinate of the point of interest
  • epsilon (float, optional) – tollerance to edges of the grid
Return type:

float

Returns:

value at point

hdf5_grids.get_interpolated_time_grids(opdict)[source]

Interpolates the NLL time grids onto the search grids. Uses options contained in a WavelocOptions.opdict.

Parameters:opdict – Dictionary of options in WavelocOptions.opdict format
hdf5_grids.nll2hdf5(nll_name, h5_name)[source]

Translates NLL files to hdf5 format.

Parameters:
  • nll_name – NLL filename without extension
  • h5_name – hdf5 filename

Waveforms

The waveform manipulation routines in waveloc are heavily based on obspy. As waveloc development started before obspy was fully functionnal, some external functions are used where obspy equivalents were not available. Updating all the waveform manipulation routines to use the latest obspy features will be completed at some time.

OP_waveform

The OP_waveforms module provides the Waveform class through which most time-series manipulation is carried out, as well as several useful functions.

Note

Requires scikits.samplerate for full functionality (on linux systems install with easy_install scikits.samplerate after having installed the libsamplerate0 libraries using e.g. apt-get install libsamplerate0). scikits.samplerate is used to provide high-fidelity resampling when standard decimation is not possible (e.g. to pass from 125Hz sampling to 100Hz sampling).

class OP_waveforms.Waveform[source]

Wrapper class for obspy streams. Adds many convenience functions, and functionality specific to WaveLoc processing.

Attributes

stream

An obspy.core.Stream object. Can be accessed directly. Initialises to None. If you set this attribute directly, you should also set trace to stream.traces[0]

trace

An obspy.Trace object. Can be accessed directly. Initialises to None. If you set this directly, make sure it points to stream.traces[0]

proc

Text string indicating processing type. Can be accessed directly. Initialises to None.

Methods

bp_filter(freqmin, freqmax, zerophase=False, rmean=False, taper=False)[source]

Apply a band-pass filter to the data. If data are segmented into multiple traces, apply the same filter to all traces. Calls obspy.signal.filter() to do the filtering.

Parameters:
  • freqmin (float) – Low frequency corner
  • freqmax (float) – High frequency corner
  • zerophase (bool) – If True applies a non-causal bandpass filter. If False applies a causal bandpass filter.
  • rmean (bool) – If True remove mean before filtering.
  • taper (bool) – If True apply taper before filtering.
Raises UserWarning:
 

if no data in stream.

channel

Channel name from self.trace (read-only).

comp

Channel name from self.trace (read-only).

compute_signature()[source]

Computes a signature for a waveform, made up by the maximum of the trace and its sum.

Return type:tuple of floats
Returns:maximum, datasum
decimate(factor=1)[source]

Applies obspy decimation, after applying an anti-aliasing, non-causal filter of 0.4*(new sample rate).

Parameters:factor (integer) – Decimation factor.
delta

Time step between data points in trace (read-only).

display(title='', filename='')[source]

Makes a quick and dirty plot of the waveform. If filename is given (and is different from “”) then the plot of the waveform is written to file, otherwise it is shown on the screen.

Parameters:
  • title (string) – Title for the plot.
  • filename (string) – Filename for the plot (format defined by the file extension).
dt

Time step between data points in trace (read-only).

get_snr(o_time, left_time, right_time)[source]

Returns signal-to-noise ratio, where signal = max(abs(signal between left_time and right_time)) and noise = median(abs(signal between left_time- and otime)).

Parameters:
  • o_time (UTCDateTime) – center time for SNR calculation
  • left_time (UTCDateTime) – left time for SNR calculation
  • right_time (UTCDateTime) – right time for SNR calculation
npts

Number of points in trace (read-only).

process_envelope()[source]

Runs envelope processing on a waveform. Replaces self.trace and sets self.proc to “Envelope’.

process_gaussian(threshold, mu=0.0, sigma=0.1)[source]

Replace local maxima by a series of dirac and convolve them with a gaussian distribution. Overwrites the waveform. Sets self.proc to ‘Gaussian’

Parameters:
  • threshold (float) – value over which the convolution is applied
  • mu (float) – expected value of the gaussian distribution
  • sigma (float) – variance of the gaussian distribution
process_kurtosis(win, recursive=False, pre_rmean=False, pre_taper=False, post_taper=True)[source]

Processing waveform using kurtosis (from statlib package). Calls filters.sw_kurtosis1(), and overwrites the waveform. Sets self.prof to ‘Kurtosis’

Parameters:
  • win – length of the window (in seconds) on which to calculate the kurtosis
  • recursive (boolean, optional) – If True applies recursive kurtosis calculation
  • pre_rmean (boolean, optional) – If True removes mean of signal before processing.
  • pre_taper (boolean, optional) – If True applies taper to signal before processing.
  • post_taper (boolean, optional) – If True applies taper to signal after processing.
process_none()[source]

No processing on a waveform. Sets self.proc to “None”

process_sta_lta(stawin, ltawin)[source]

Runs classic short to long term average ratio processing on waveform. Replaces self.trace and sets self.proc to StaLta.

Parameters:
  • stawin – length of the short term average window in seconds
  • ltawin – length of the long term average window in seconds
read_from_SDS(sds_root, net_name, sta_name, comp_name, starttime=None, endtime=None, rmean=False, taper=False, pad_value=None)[source]

Read waveform data from an SDS structured archive. Simple overlaps and adjacent traces are merged if possile.

Parameters:
  • sds_root (string) – root of the SDS archive
  • net_name (string) – network name
  • sta_name (string) – station name
  • comp_name (string) – component name
  • starttime (obspy.core.utcdatetime.UTCDateTime object, optional) – Start time of data to be read.
  • endtime (obspy.core.utcdatetime.UTCDateTime object, optional) – End time of data to be read.
  • rmean (boolean, optional) – If True removes the mean from the data upon reading. If data are segmented, the mean will be removed from all segments individually.
  • taper (boolean, optional) – If True applies a cosine taper to the data upon reading. If data are segmented, tapers are applied to all segments individually.
  • pad_value (float, optional) – If this parameter is set, points between starttime and the first point in the file, and points between the last point in the file and endtime, will be set to pad_value. You may want to also use the rmean and taper parameters, depending on the nature of the data.
Raises UserWarning:
 

If there are no data between starttime and endtime

read_from_file(filename, format=None, starttime=None, endtime=None, rmean=False, taper=False, pad_value=None)[source]

Read waveform data from file. Multiple traces are merged if they overlap exactly or are adjacent.

Parameters:
  • filename – Waveform filename
  • format (string) – obspy format type (e.g. ‘SAC’, ‘mseed’...)
  • starttime (obspy.core.utcdatetime.UTCDateTime object) – Start time of data to be read.
  • endtime (obspy.core.utcdatetime.UTCDateTime object) – End time of data to be read.
  • rmean (boolean) – If True removes the mean from the data upon reading. If data are segmented, the mean will be removed from all segments individually.
  • taper (boolean) – If True applies a cosine taper to the data upon reading. If data are segmented, tapers are applied to all segments individually.
  • pad_value (float) – If this parameter is set, points between starttime and the first point in the file, and points between the last point in the file and endtime, will be set to pad_value. You may want to also use the rmean and taper parameters, depending on the nature of the data.
Raises UserWarning:
 

If there are no data between starttime and endtime

resample(new_samplerate, resample_type='sinc_best')[source]

Applies audio-quality resampling in place. Requires installation of scikits.samplerate.

Parameters:
  • new_samplerate (float) – New sample rate.
  • resample_type (string) – Can be 'sinc_best', ...
rmean()[source]

Removes the mean of the stream (iterates over all traces).

Raises UserWarning:
 if no data in stream.
starttime

Start time of self.trace (read-only).

station

Station name from self.trace (read-only).

t_array

A numpy.array containing times in seconds from self.trace.stats.starttime for data points in self.trace (read-only).

take_positive_derivative(pre_rmean=False, pre_taper=False, post_taper=True)[source]

Takes positive derivative of a stream (used to calculate the positive kurtosis gradient).

Parameters:
  • pre_mean – If True then remove the mean before taking positive gradient
  • pre_taper – If True then apply a taper before taking positive

gradient :param post_taper: If True then apply a taper after taking positive gradient.

Raises UserWarning:
 if no data in stream.
taper()[source]

Applies a cosine taper the stream (iterates over all traces).

Raises UserWarning:
 if no data in stream.
values

Values of self.trace.data (read-only).

write_to_file(filename, format=None)[source]

Write waveform to file.

Parameters:
  • filename (string) – Output filename.
  • format (string) – obspy format type (e.g. ‘SAC’, ‘mseed’...)
write_to_file_filled(filename, format=None, fill_value=0, rmean=False, taper=False)[source]
Write waveform to file, after merging and filling blanks with
fill_value.
Parameters:
  • filename (string) – Output filename.
  • formatobspy format type (e.g. ‘SAC’, ‘mseed’...)
  • fill_value (float, optional) – Value used to fill in gaps
  • rmean (bool, optional) – Remove the mean before merging.
  • taper (bool, optional) – Apply taper before merging.
OP_waveforms.compute_gauss(dt, mu, sig)[source]

Given a time-step, a mean and sigma, returns an array with a Gaussian signal from -4*sigma to +4*sigma around the mean.

Parameters:
  • dt (float) – time-step in seconds
  • mu (float) – mean of the Gaussian
  • sig (float) – half-width sigma of the Gaussian
Return type:

numpy array

Returns:

array with Gaussian signal

OP_waveforms.read_data_compatible_with_time_dict(filenames, time_dict, starttime, endtime)[source]

Reads data that have the same station names as those used as keys to a given time dictionary

Parameters:
  • filename (list of strings) – list of filenames
  • time_dict (dictionary) – dictionary indexed by station names
  • starttime (UTCDateTime) – start-time to read data
  • endtime (UTCDateTime) – end-time to read data
Raises UserWarning:
 

if sampling is different between the files

OP_waveforms.stream_positive_derivative(st)[source]

Takes first time derivative of a stream (iterates over all available traces) and keep only positive values (set negative values to zero).

Parameters:st (An obspy.stream object) – Input stream, processed in place
OP_waveforms.stream_rmean(st)[source]

Removes mean from a stream (iterates over all available traces).

Parameters:st (An obspy.stream object) – Input stream, processed in place
OP_waveforms.stream_taper(st)[source]

Applies cosine taper to a stream (iterates over all available traces).

Parameters:st (An obspy.stream object) – Input stream, processed in place

filters

Functions for timeseries filtering.

Created by Alessia Maggi and Alberto Michelini.

filters.allens_stalta(x, dt, STAwin, LTAwin)[source]

Determines the short to long time average ratio of a timeseries. using Allen (1978) caracteristic function (see report)

Parameters:
  • x – signal
  • dt – sampling interval
  • dt – sampling interval in seconds.
  • STAwin – short time average window in seconds.
  • LTAwin – long time average in seconds.
Returns:

processed signal

filters.envelope(x, N)[source]

Determines the envelope of a signal using the Hilbert transform. Uses scipy.signal.hilbert() to do the Hilbert transform.

Parameters:
  • x – signal
  • N
Returns:

envelope as numpy array

filters.filtfilt(b, a, x)[source]

TODO What does this function do? Alberto??

Parameters:
  • b
  • a
  • x
Raises ValueError:
 

filters.gradient(x, dt, LENwin)[source]

Determines the gradient of a timeseries.

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
LENwin:

time window (in secs) over which the gradient is determined

filters.kurto(x, dt, LENwin)[source]

Determines the kurtosis of a timeseries.

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
LENwin:

time window (in secs) over which the kurtosis is determined

filters.kurto_improved(x, dt, LENwin)[source]

Determines the kurtosis of a timeseries as described by kuperkoch et al. 2010: calculate the kurtosis recursively. Results are not satisfying, but one may want to improve this... it might save some time!

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
LENwin:

time window (in secs) over which the kurtosis is determined

filters.lfilter_zi(b, a)[source]

Computes the zi state from the filter parameters. See [Gust96].

Based on: [Gust96] Fredrik Gustafsson, Determining the initial states in forward-backward filtering, IEEE Transactions on Signal Processing, pp. 988–992, April 1996, Volume 44, Issue 4

Parameters:
  • b
  • a
Returns:

zi state as numpy array

filters.projection(x, y, theta)[source]

Project x and y onto a line of angle theta.

filters.rec_dx2(x, C1)[source]

Recursive dx2 (from Chassande-Mottin, 2002)

Parameters:
  • x – signal
  • C
Returns:

processed signal as numpy array

filters.rec_kurtosis(x, C1)[source]

Recursive Kurtosis calculated using Chassande-Mottin (2002). Beware, this function is unstable for strongly non-stationary signals.

Parameters:
  • x – signal
  • C
Returns:

processed signal as numpy array

filters.rec_kurtosis_old(x, C)[source]

Recursive pseudo-kurtosis calculation as used by Langet et al. (2014)

Parameters:
  • x – signal
  • C
Returns:

procesed signal as numpy array

filters.skewness(x, dt, LENwin)[source]

Determines the absolute value of the skewness of a timeseries.

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
  • LENwin – time window (in secs) over which the skewness is determined
filters.smooth(x, window_len=11, window='hanning')[source]

Smooth the data using a window with requested size.

This method is based on the convolution of a scaled window with the signal. The signal is prepared by introducing reflected copies of the signal (with the window size) in both ends so that transient parts are minimized in the begining and end part of the output signal.

Parameters:
  • x – the input signal
  • window_len – the dimension of the smoothing window; should be an odd integer
  • window – the type of window from ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’. Flat window will produce a moving average smoothing.
Returns:

the smoothed signal

Raises ValueError:
 

example:

t=linspace(-2,2,0.1) x=sin(t)+randn(len(t))*0.1 y=smooth(x)

see also: numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman,
numpy.convolve scipy.signal.lfilter
TODO: the window parameter could be the window itself if an array instead
of a string
NOTE: length(output) != length(input), to correct this: return
y[(window_len/2-1):-(window_len/2)] instead of just y.
filters.stalta(x, dt, STAwin, LTAwin)[source]

Determines the short to long time average ratio of a timeseries.

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
  • STAwin – short time average window in seconds.
  • LTAwin – long time average in seconds.
Returns:

processed signal

filters.sw_kurtosis1(x, n)[source]

Returns kurtosis calculated over data set x using sliding windows of length n points. Length of returned array is len(x)-n+1.

Parameters:
  • x – signal
  • n – length of sliding window in points
Returns:

procesed signal as numpy array

filters.sw_kurtosis2(x, n)[source]

Returns kurtosis calculated over data set x using sliding windows of length n points. Length of returned array is len(x)-n+1. This version uses filter.window() to slice the data array and is faster than sw_kurtosis1.

Parameters:
  • x – signal
  • n – length of sliding window in points
Returns:

procesed signal as numpy array

filters.variance(x, dt, LENwin)[source]

Determines the variance of a timeseries.

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
LENwin:

time window (in secs) over which the variance is determined

Returns:

processed signal as numpy array

filters.variance_under_rotation(x, y, dt, LENwin)[source]

Determines the variance under rotation of a timeseries.

Parameters:
  • x – signal
  • dt – sampling interval in seconds.
LENwin:

time window (in secs) over which the variance under rotation is determined

filters.varrot(x, y)[source]

TODO : figure out what this function actually does !!

filters.window(seq, n=2)[source]
Returns a sliding window (of width n) over data from the iterable
s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ...
Parameters:
  • seq – an iterable sequence
  • n – length of window (number of elements)

kurtogram

kurtogram.DBFB(x, h, g)[source]

Double-band filter-bank. [a,d] = DBFB(x,h,g) computes the approximation coefficients vector a and detail coefficients vector d, obtained by passing signal x though a two-band analysis filter-bank.

Parameters:
  • x (numpy array) – signal
  • h (numpy array) – The decomposition low-pass filter and
  • g (numpy array) – The decomposition high-pass filter.
Return type:

numpy array

Returns:

a, d

kurtogram.Fast_Kurtogram(x, nlevel, verbose=False, Fs=1, NFIR=16, fcut=0.4, opt1=1, opt2=1)[source]

Computes the fast kurtogram Kwav of signal x up to level ‘nlevel’ Maximum number of decomposition levels is log2(length(x)), but it is recommended to stay by a factor 1/8 below this. Also returns the vector of k-levels Level_w, the frequency vector freq_w, the complex envelope of the signal c and the extreme frequencies of the “best” bandpass f_lower and f_upper.

J. Antoni : 02/2005 Translation to Python: T. Lecocq 02/2012

Parameters:
  • x (numpy array) – signal to analyse
  • nlevel (integer) – number of decomposition levels
  • verbose – If True outputs debugging information
  • Fs (integer) – Sampling frequency of signal x
  • NFIR (integer) – Length of FIR filter
  • fcut (float) – Fraction of Nyquist for filter
  • opt1 – [1 | 2]: * opt1 = 1: classical kurtosis based on 4th order statistics * opt1 = 2: robust kurtosis based on 2nd order statistics of the envelope (if there is any difference in the kurtogram between the two measures, this is due to the presence of impulsive additive noise)
  • opt2 – [1 | 2]: * opt2=1: the kurtogram is computed via a fast decimated filterbank * opt2=2: the kurtogram is computed via the short-time Fourier transform (option 1 is faster and has more flexibility than option 2 in the design of the analysis filter: a short filter in option 1 gives virtually the same results as option 2)
Returns:

Kwav, Level_w, freq_w, c, f_lower, f_upper * Kwav: kurtogram * Level_w: vector of levels * freq_w: frequency vector * c: complex envelope of the signal filtered in the frequency band that maximizes the kurtogram * f_lower: lower frequency of the band pass * f_upper: upper frequency of the band pass

kurtogram.Find_wav_kurt(x, h, g, h1, h2, h3, nlevel, Sc, Fr, Fs=1, verbose=False)[source]

TODO flesh out this doc-string

J. Antoni : 12/2004 Translation to Python: T. Lecocq 02/2012

Parameters:
  • x (numpy array) – signal
  • h (numpy array) – lowpass filter
  • g (numpy array) – highpass filter
  • h1 (numpy array) – filter parameter returned by get_h_parameters
  • h2 (numpy array) – filter parameter returned by get_h_parameters
  • h3 (numpy array) – filter parameter returned by get_h_parameters
  • nlevel (integer) – number of decomposition levels
  • Sc – Sc = -log2(Bw)-1 with Bw the bandwidth of the filter
  • Fr (float) – in the range [0, 0.5]
  • Fs – Sampling frequency of signal x
  • verbose – If True outputs debugging information
Type:

Fs: integer

Returns:

c, s, threshold, Bw, fc

kurtogram.K_wpQ(x, h, g, h1, h2, h3, nlevel, verbose, opt, level=0)[source]

Calculates the kurtosis K (2-D matrix) of the complete quinte wavelet packet transform w of signal x, up to nlevel, using the lowpass and highpass filters h and g, respectively. The WP coefficients are sorted according to the frequency decomposition. This version handles both real and analytical filters, but does not yield WP coefficients suitable for signal synthesis.

J. Antoni : 12/2004 Translation to Python: T. Lecocq 02/2012

Parameters:
  • x (numpy array) – signal
  • h (numpy array) – lowpass filter
  • g (numpy array) – higpass filter
  • h1 (numpy array) – filter parameter returned by get_h_parameters
  • h2 (numpy array) – filter parameter returned by get_h_parameters
  • h3 (numpy array) – filter parameter returned by get_h_parameters
  • nlevel (integer) – number of decomposition levels
  • verbose – If True outputs debugging information
  • opt (string) – [‘kurt1’ | ‘kurt2’] * ‘kurt1’ = variance of the envelope magnitude * ‘kurt2’ = kurtosis of the complex envelope
  • level – decomposition level for this call
Returns:

kurtosis

kurtogram.K_wpQ_filt(x, h, g, h1, h2, h3, acoeff, bcoeff, level=0)[source]

Calculates the kurtosis K of the complete quinte wavelet packet transform w of signal x, up to nlevel, using the lowpass and highpass filters h and g, respectively. The WP coefficients are sorted according to the frequency decomposition. This version handles both real and analytical filters, but does not yield WP coefficients suitable for signal synthesis.

J. Antoni : 12/2004 Translation to Python: T. Lecocq 02/2012

Parameters:
  • x – signal
  • h – lowpass filter
  • g – highpass filter
  • h1 – filter parameter returned by get_h_parameters
  • h2 – filter parameter returned by get_h_parameters
  • h3 – filter parameter returned by get_h_parameters
  • acoeff
  • acoeff
  • level
kurtogram.K_wpQ_filt_local(x, h, g, h1, h2, h3, acoeff, bcoeff, level)[source]

Performs one analysis level into the analysis tree TODO : flesh out this doc-string

Parameters:
  • x – signal
  • h – lowpass filter
  • g – higpass filter
  • h1 – filter parameter returned by get_h_parameters
  • h2 – filter parameter returned by get_h_parameters
  • h3 – filter parameter returned by get_h_parameters
  • acoeff
  • acoeff
  • level
kurtogram.K_wpQ_local(x, h, g, h1, h2, h3, nlevel, verbose, opt, level)[source]

Is a recursive funtion. Computes and returns the 2-D vector K, which contains the kurtosis value of the signal as well as the 2 kurtosis values corresponding to the signal filtered into 2 different band-passes. Also returns and computes the 2-D vector KQ which contains the 3 kurtosis values corresponding to the signal filtered into 3 different band-passes.

Parameters:
  • x (numpy array) – signal
  • h (numpy array) – lowpass filter
  • g (numpy array) – highpass filter
  • h1 (numpy array) – filter parameter returned by get_h_parameters
  • h2 (numpy array) – filter parameter returned by get_h_parameters
  • h3 (numpy array) – filter parameter returned by get_h_parameters
  • nlevel (integer) – number of decomposition levels
  • verbose – If True outputs debugging information
  • opt (string) – [‘kurt1’ | ‘kurt2’] * ‘kurt1’ = variance of the envelope magnitude * ‘kurt2’ = kurtosis of the complex envelope
  • level (integer) – decomposition level for this call
Returns:

K, KQ

kurtogram.TBFB(x, h1, h2, h3)[source]

Triple-band filter-bank. [a1,a2,a3] = TBFB(x,h1,h2,h3)

Parameters:
  • x (numpy array) – signal
  • h1 (numpy array) – filter parameter
  • h2 (numpy array) – filter parameter
  • h3 (numpy array) – filter parameter
Return type:

numpy array

Returns:

a1, a2, a3

kurtogram.binary(i, k)[source]

Computes the coefficients of the binary expansion of i: i = a(1)*2^(k-1) + a(2)*2^(k-2) + ... + a(k)

Parameters:
  • i – integer to expand
  • k – nummber of coefficients
Returns:

coefficients a

kurtogram.do_kurtogram_setup_and_run(opdict)[source]

Run the kurtogram analysis using the parameters contained in the WavelocOptions.opdict.

Parameters:opdict – Dictionary containing the waveloc parameters and options.
kurtogram.eps = 2.2204460492503131e-16

Implements the kurtogram routine of Antoni (2005).

kurtogram.getBandwidthAndFrequency(nlevel, Fs, level_w, freq_w, level_index, freq_index)[source]

Gets bandwidth bw and frequency parameters knowing the level and the frequency indexes.

Parameters:
  • nlevel (integer) – number of decomposition levels
  • Fs (integer) – sampling frequency of the signal
  • level_w (numpy array) – vector of decomposition levels
  • freq_w (numpy array) – vector of frequencies
  • level_index (integer) – index of the level
  • freq_index (integer) – index of the frequency
Returns:

bw, fc, fi, l1 * bw: bandwidth * fc: central frequency * fi: index of the frequency sequence within the level l1 * l1: level

kurtogram.getFTSquaredEnvelope(c)[source]

Calculates the Fourier transform of the squared envelope

Parameters:c – signal
Returns S:FT of squared envelope
kurtogram.get_GridMax(grid)[source]

Gets maximum of a nD grid and its unraveled index

Parameters:grid – an nD-grid
Returns:
  • M : grid maximum
  • index : index of maximum in unraveled grid
kurtogram.get_h_parameters(NFIR, fcut)[source]

Calculates h-parameters used in Antoni (2005)

Parameters:
  • NFIR (integer) – length of FIR filter
  • fcut (float) – fraction of Nyquist for filter
Return type:

numpy array

Returns:

h-parameters: h, g, h1, h2, h3

kurtogram.kurt(x, opt)[source]

Calculates kurtosis of a signal according to the option chosen

Parameters:
  • x (numpy array) – signal
  • opt (string) – [‘kurt1’ | ‘kurt2’] * ‘kurt1’ = variance of the envelope magnitude * ‘kurt2’ = kurtosis of the complex envelope
Return type:

float

Returns:

Kurtosis

kurtogram.kurto(origin_time, info, opdict)[source]

Finds for each Waveloc event and for each station the best filtering parameters for kurtosis computation. Writes them into the dictionary info.

Parameters:
  • origin_time (utcdatetime) – origin time of the signal
  • info (dictionary) – dictionary of parameters
  • opdict (dictionary) – dictionary of the Waveloc parameters and options
Return type:

dictionary

Returns:

info

kurtogram.nextpow2(n)[source]

Calculates next power of 2.

Parameters:n (integer) – a number
Return type:integer
Returns:The next power of 2
kurtogram.plotKurtogram(Kwav, freq_w, nlevel, Level_w, Fs, fi, I)[source]

Plots the kurtogram.

Parameters:
  • Kwav (numpy array) – kurtogram
  • freq_w (numpy array) – frequency vector
  • nlevel (integer) – number of decomposition levels
  • level_w (numpy array) – vector of levels
  • Fs (integer) – sampling frequency of the signal
  • fi
  • I (integer) – level index
kurtogram.plot_envelope(x, Fs, c, fc, level, spec=False)[source]

Plots envelope (with or without its spectrum)

Parameters:
  • x (numpy array) – signal
  • Fs (float) – sampling frequency of signal
  • c (numpy array) – complex envelope of signal
  • fc (float) – central frequency of the bandpass in Hz
  • level (float) – index of the decomposition level
  • spec – If True also plots the envelope spectrum
kurtogram.plot_trace(fig, G, x, xfilt, kurtx, tr, info, f_lower, f_upper, snr, snr_ref, snr_kurt, kmax, kmax_ref, tstack)[source]

Plots both signal and kurtosis for comparison. TODO : flesh out this doc-string.

Parameters:
  • fig (matplotlib.pyplot.figure) – figure
  • G (matplotlib.gridspec) – location of the subplot in the figure
  • x (numpy array) – initial signal
  • xfilt (numpy array) – signal filtered in bandpass determined after the kurtogram analysis
  • kurtx (numpy array) – kurtosis of the initial signal
  • tr (numpy array) – kurtosis of the filtered signal
  • info (dictionary) – dictionary of parameters
  • f_lower (float) – lower frequency of the filtering bandpass
  • f_upper (flaot) – upper frequency of the filtering bandpass
  • snr (float) – signal-to-noise ratio of xfilt
  • snr_ref (float) – signal-to-noise ratio of x
  • snr_kurt (float) – signal-to-noise ratio of tr
  • kmax (float) – maximum value of tr
  • kmax_ref (float) – maximum value of kurtx
  • tstack (utcdatetime) – origin time of the signal
kurtogram.raylinv(p, b)[source]

Inverse of the Rayleigh cumulative distribution function (cdf). X = RAYLINV(P,B) returns the Rayleigh cumulative distribution function with parameter B at the probabilities in P.

Parameters:
  • p – probabilities
  • b
kurtogram.read_kurtogram_frequencies(filename)[source]

Reads the binary file with kurtogram frequencies. Plots the histograms of lower and upper frequencies for each station. Aims at determining the best filtering parameters.

Parameters:filename (string) – File to read
kurtogram.waveval(xall, tstart, tend, dt, tdeb)[source]

Returns the part of signal corresponding to a given event.

Parameters:
  • xall (numpy array) – complete signal
  • tstart (utcdatetime) – start time of the signal
  • tend (utcdatetime) – end time of the signal
  • dt (float) – time sampling in seconds
  • tdeb (utcdateime) – start time of the whole data
Return type:

numpy array

Returns:

val

kurtogram.write_file(info, tstart, tend, tr)[source]

Replaces the initial kurtosis by a new one (from the signal filtered in the preferred band pass). Writes it in the dictionary info.

Parameters:
  • info (dictionary) – dictionary of parameters
  • tstart (utcdatetime) – start-time
  • tend (utcdatetime) – end-time
  • tr (numpy array) – trace to write
Return type:

dictionary

Returns:

info

Migration

migration

migration.do_migration_loop_continuous(opdict, data, delta, start_time, grid_info, time_grids, keep_grid=False, keep_stacks=True)[source]

Do continuous migration loop. TODO : flesh out this doc-string.

migration.do_migration_setup_and_run(opdict)[source]

Do setup and launch migration.

Parameters:opdict – WavelocOptions.opdict
migration.extract_max_values(stack_grid, search_info, f_stack, use_ram=False, n_max=500000000.0)[source]

Extracts maximum stack value from a 4D migrated grid. Also extracts x_max, y_max, z_max. TODO : flesh out this doc-string.

Parameters:
  • stack_grid
  • search_info
  • f_stack
  • use_fam
  • n_max
migration.migrate_4D_stack(data, delta, time_grids, stack_grid, use_ram=False)[source]

Migrate the 4D stack. Central part of the waveloc process. TODO : flesh out this doc-string

Parameters:
  • data
  • delta
  • time_grids
  • stack_grid
  • use_ram

synth_migration

synth_migration.generateSyntheticDirac(opdict, time_grids=None)[source]

Generates a synthetic test and does the migration. All options are given in the WavelocOptions.opdict.

Parameters:
  • opdict – Waveloc options / parameters
  • time_grids

integrate4D

integrate4D.compute_expected_coordinates1D(grid, x0)[source]

Computes expected value and variance of a 1D grid along its (only) axis.

Parameters:
  • grid (numpy array) – 1D grid to integrate
  • x0 (numpy array) – coordinates of all points along the integration axis
Return type:

float

Returns:

  • exp_x0 : expected value
  • var_x0 : variance

integrate4D.compute_expected_coordinates3D(grid, x0, x1, x2, return_2Dgrids=False)[source]

Computes expected coordinates of a 3D grid

Parameters:
  • grid (numpy array) – 3D grid to integrate
  • x0 (numpy array) – coordinates of all points along axis=0
  • x1 (numpy array) – coordinates of all points along axis=1
  • x2 (numpy array) – coordinates of all points along axis=2
  • return_2Dgrids (boolean) – if True returns the 1D and 2D marginals in a dictionary.
Return type:

numpy array

Returns:

  • exp_x0 : expected value along axis=0
  • exp_x1 : expected value along axis=1
  • exp_x2 : expected value along axis=2
  • cov_matrix : covariance matrix
  • grid_dict : dictionary containing 1D and 2D marginals (optional)

integrate4D.compute_expected_coordinates4D(grid, x0, x1, x2, x3, return_2Dgrids=False)[source]

Computes expected coordinates of a 4D grid

Parameters:
  • grid (numpy array) – 4D grid to integrate
  • x0 (numpy array) – coordinates of all points along axis=0
  • x1 (numpy array) – coordinates of all points along axis=1
  • x2 (numpy array) – coordinates of all points along axis=2
  • x3 (numpy array) – coordinates of all points along axis=3
  • return_2Dgrids (boolean) – if True returns the 1D and 2D marginals in a dictionary.
Return type:

numpy array

Returns:

  • exp_x0 : expected value along axis=0
  • exp_x1 : expected value along axis=1
  • exp_x2 : expected value along axis=2
  • exp_x3 : expected value along axis=2
  • cov_matrix : covariance matrix
  • grid_dict : dictionary containing 1D and 2D marginals (optional)

integrate4D.compute_integral1D(grid, x0)[source]

Computes norm of a 1D grid

Parameters:
  • grid (numpy array) – 1D grid to integrate
  • x0 (numpy array) – coordinates of all points along the integration axis
Return type:

float

Returns:

grid norm

integrate4D.compute_integral3D(grid, x0, x1, x2)[source]

Computes norm of a 3D grid

Parameters:
  • grid (numpy array) – 3D grid to integrate
  • x0 (numpy array) – coordinates of all points along axis=0
  • x1 (numpy array) – coordinates of all points along axis=1
  • x2 (numpy array) – coordinates of all points along axis=2
Return type:

float

Returns:

grid norm

integrate4D.compute_integral4D(grid, x0, x1, x2, x3)[source]

Computes norm of a 4D grid

Parameters:
  • grid (numpy array) – 4D grid to integrate
  • x0 (numpy array) – coordinates of all points along axis=0
  • x1 (numpy array) – coordinates of all points along axis=1
  • x2 (numpy array) – coordinates of all points along axis=2
  • x3 (numpy array) – coordinates of all points along axis=3
Return type:

float

Returns:

grid norm

Location

locations_trigger

locations_trigger.do_locations_trigger_setup_and_run(opdict)[source]
Run locations trigger. Takes all options and paramters from
WavelocOptions.opdict dictionary.
Parameters:opdict – Dictionary of waveloc options / parameters.
locations_trigger.number_good_kurtosis_for_location(kurt_files, data_files, loc, time_dict, snr_limit=10.0, snr_tr_limit=10.0, sn_time=10.0)[source]

Analyses the filtered data and the kurtosis time-series to determine the number of stations whose traces have sufficiently high signal-to-noise ratio (SNR) to be useful for the location. Both time-series need to satisfy the conditions for a station to be counted as contributing to the location.

Parameters:
  • kurt_files – Filenames for kurtosis files. Depending on the filenames given in this list, the function will analyse kurtosis, kurtosis-gradient or gaussian waveforms.
  • data_files – Filenames for filtered data.
  • loc – Location dictionary for the event to be analysed
  • time_dict – Dictionary of travel-times for the location to be analysed
  • snr_limit – SNR limit for kurtosis-type data
  • snr_limit_tr – SNR limit for filtered data
  • snr_time – Length of time in seconds before the event for computation of SNR.
Return type:

integer

Returns:

Numer of stations that have contributed to the location.

locations_trigger.plot_location_triggers(trace, trig_start, trig_end, trig_95_start, trig_95_end, show=True)[source]

Plot the location (in time) of the triggers. TODO : Flesh out this doc-string

Parameters:
  • trace
  • trig_start
  • trig_end
  • trig_95_start
  • trig_95_end
  • show
locations_trigger.read_header_from_file(filename, opdict)[source]

Read header information from a location file into a WavelocOptions.opdict.

Parameters:
  • filename – File to read.
  • opdict – Waveloc options dictionary to write into.
Returns:

opdict

locations_trigger.read_locs_from_file(filename)[source]

Read locations from file.

Parameters:filename – File to read.
Returns:List of locations (each location is a dictionary)
locations_trigger.trigger_locations_inner(max_val, max_x, max_y, max_z, left_trig, right_trig, start_time, delta)[source]

Inner loop of the location process.

Parameters:
  • max_val – Time-series of stack-max.
  • max_x – Time-series of the x-positions corresponding to stack-max.
  • max_y – Time-series of the y-positions corresponding to stack-max.
  • max_z – Time-series of the z-positions corresponding to stack-max.
  • left_trig – Amplitude for trigger-on.
  • right_trig – Amplitude for trigger-off.
  • start_time – UTCDateTime of the first point in the time-series.
  • delta – Sampling interval of the time-series
Returns:

List of locations. Each location is a dictionary containing all the relevant information.

locations_trigger.write_header_options(loc_file, opdict)[source]

Write the information regarding how the locations were generated to an already opened location file. This information becomes a location “header”.

Parameters:
  • loc_file – File object to which the information will be written.
  • opdict – Waveloc options dictionary

locations_prob

locations_prob.do_locations_prob_setup_and_run(opdict)[source]

Setup and run probability-based locations on migration grids. Takes all parameters from WavelocOptions.opdict.

Parameters:opdict – Parameters and options for Waveloc.
locations_prob.read_prob_locs_from_file(filename)[source]

Read file containing probability determined locations.

Parameters:filename – File to be read.
Returns:Dictionary of locations

magnitude

magnitude.bvalue(mag, r)[source]

Compute the b-value by a simple linear fitting

magnitude.do_comp_mag(opdict)[source]

Do the magnitude computation, using options in WavelocOptions.opdict.

Parameters:opdict – Dictionary containing waveloc options and parameters.
magnitude.fill_values(vals, tdeb, data_glob, data_dir, comp)[source]

Create a dictionnary with all values for each channel of each station. TODO : Flesh out this doc-string.

Parameters:
  • vals
  • tdeb
  • data_glob
  • data_dir
  • comp
magnitude.read_paz(files)[source]

Read dataless and extract poles and zeros from dataless SEED files. Based on obspy.

Parameters:files – Files to read.

plot_locations2

plot_locations2.do_plotting_setup_and_run(opdict, plot_wfm=True, plot_grid=True)[source]

Plot the results of a wavloc run (migration and location). All options and parameters are taken from an opdict.

Parameters:
  • opdict – WavlocOptions.opdict that contains the options / parameters.
  • plot_wfm (boolean) – If True plots waveforms after location (filtered data and kurtosis).
  • plot_grid (boolean) – If ``True``plots the migration grid.
plot_locations2.do_probloc_plotting_setup_and_run(opdict)[source]

Plot the results of a wavloc run (migration and location using probability density). All options and parameters are taken from an opdict.

Parameters:opdict – WavlocOptions.opdict that contains the options / parameters.

plot_mpl

plot_mpl.plotDiracTest(test_info, fig_dir, otime_window)[source]

Creates plot for synthetic test. TODO : flesh out this doc-string

Parameters:
  • test_info
  • fig_dir
  • otime_window
plot_mpl.plotLocationGrid(loc, grid_info, fig_dir, otime_window)[source]

Plots location grid. TODO : Flesh out this doc-string.

Parameters:
  • loc
  • grid_info
  • fig_dir
  • otime_window
plot_mpl.plotLocationWaveforms(loc, start_time, dt, data_dict, grad_dict, stack_wfm, fig_dir)[source]

Creates plot for located waveforms. Assumes data and grad are ready for plotting. TODO : Flesh out this doc-string.

Parameters:
  • loc
  • start_time
  • dt
  • data_dict
  • grad_dict
  • stack_wfm
  • fig_dir
plot_mpl.plotProbLoc(marginals, prob_loc, loc, fig_dir, space_only)[source]

Creates plot for probabilistic location TODO : flesh out this doc-string

Parameters:
  • marginals
  • prob_loc
  • loc
  • fig_dir
  • space_only

Clustering

This part of the waveloc code was contributed by Nadège Langet during her PhD thesis.

correlation

Provides classes and functions for cross-correlation.

class correlation.BinaryFile(filename)[source]

Class that reads, writes and stores the binary files containing the correlation matrix.

read_binary_file()[source]

Reads the binary file from disk and returns un-pickled data :rtype: numpy array ? :returns: data read from tile

write_binary_file(input)[source]

Writes a binary file from disk

Parameters:input – input data ready to be pickled
correlation.corr_freq(f, Cxy, v)[source]

Computes the time lag entirely in the spectral domain, using the slope of the phase spectrum where the amplitude is strong.

Parameters:
  • f (numpy array) – frequency vector
  • Cxy (numpy array) – frequency domain cross-correlation vector
  • v (boolean) – If True displays Cxy
Return type:

float

Returns:

time-lag

correlation.corr_time(Cxy, Cxx, Cyy, dt, v)[source]

Returns maximum of cross-correlation and time-lag given frequency domain cross- and auto- correlation vectors.

Parameters:
  • Cxy (numpy array) – frequency domain cross-correlation
  • Cxx (numpy array) – frequency domain auto-correlation
  • Cyy (numpy array) – frequency domain auto-correlation
  • dt (float) – time step
  • v (boolean) – If True plots cross-correlation in the time-domain
Return type:

float

Returns:

value_t, tau_t, respectively the cross-correlation maximum value and time-lag

correlation.correlate(x, y, dt, v, a)[source]

Performs correlation on two signals.

Parameters:
  • x (numpy array) – waveform
  • y (numpy array) – waveform
  • dt (float) – time-step
  • v (boolean) – If True makes plots
  • a (char) – [‘t’ | ‘f’] for time-domain or frequency-domain extraction of time-lag
Return type:

float

Returns:

The returned values depend on the a input parameter:

  • a=’t’ : returns time-lag and correlation value
  • a=’f’ : returns time-lag

correlation.cum(x, v)[source]

Determines the minimum and maximum indices of the most significative part of a signal. (Here, determines the part of the amplitude spectrum which is useful)

Parameters:
  • x (numpy array) – signal
  • v (boolean) – If True then plot x and cumulative sum of x
Return type:

int

Returns:

mini, maxi

correlation.display(Cxy, f, f_min_max)[source]

Displays amplitude and phase of the frequency domain cross-correlation.

Parameters:
  • Cxy (numpy array) – frequency domain cross-correlation
  • f (numpy array) – frequency vector
  • f_min_max (numpy array) –
correlation.do_correlation_setup_and_run(opdict)[source]

Runs correlation using options contained in a WavelocOptions.opdict dictionary.

Explores and cross-correlates all possible event pairs of the Waveloc location file at all stations. Writes the correlation coefficients and time delays in 2-D numpy arrays for each station and saves the final dictionaries into 2 binary files.

The correlation first takes place in the time domain. If the correlation value is over a given threshold, the correlation is also performed in the frequency domain so that a subsample precision can be obtained on the time delay.

correlation.fourier(x, y, dt)[source]

Computes cross-correlation and auto-correlation vectors of two signals in the frequency domain.

Parameters:
  • x (numpy array) – first signal
  • y (numpy array) – second signal
  • dt (float) – time-step of both signals
Return type:

numpy arrays

Returns:

Cxy, Cxx, Cyy, f

correlation.plot_waveform(x, y, dt, tau, ev1, ev2)[source]

Plots the waveforms. On the first plot, both waveforms are superimposed. On the second plot, they are plotted separately.

Parameters:
  • x (numpy array) – waveform
  • y (numpy array) – waveform
  • dt (float) – time-step
  • tau (float) – time-lag
  • ev1 (string) – name of event shown in waveform x
  • ev2 (string) – name of event shown in waveform y
correlation.waveform(filename)[source]

Convenience function to read a seismogram and return its data, time-step, station name and start-time

Parameters:filename (string) – File to read
Returns:val, dt, name, tdeb
  • val = seismogram values
  • dt = time-step
  • name = station name
  • tdeb = start-time (UTCDateTime)

clustering

Provides classes and functions for clustering of earthquakes based on the depth-first algorithm.

clustering.DFS(GRAPH, summit_first, cluster_ind)[source]

Implementation of the depth first search (DFS) algorithm. This is a recursive algorithm. See detailed description in do_clustering function.

Parameters:
  • GRAPH (Graph object) – the algorithm updates GRAPH.flag and GRAPH.cluster_index
  • summit_first (integer) – event index by which the research of neighbours begins
  • cluster_ind (integer) – cluster index
class clustering.Graph[source]

Class of Graph objects. Contains 4 attributes:

** Attributes **

flag

Indicates if an event already belongs to a cluster (1) or not (0).

cluster_index

Indicates the cluster number of all events. 0 means the event does not belong to any cluster.

neighbours

Lists the indexes of the neighbour events for each event. GRAPH.neighbours is then a list of lists.

set_cluster_index(value)[source]

Appends value to cluster_index attribute

Parameters:value (integer) – Value to append.
set_flag(value)[source]

Appends value to flag attribute.

Parameters:value (integer) – Value to append (0=False, 1=True).
set_neighbours(value)[source]

Appends value to neighbours attribute

Parameters:value (list) – Value to append (list of integers)
clustering.compute_nbsta(event, coeff, threshold)[source]

Computes the number of stations where the correlation value is greater than the threshold for every event pair

Parameters:
  • event (integer) – total number of events located by Waveloc
  • coeff (dictionary) – cross-correlation coefficients of all possible event pairs for all stations
  • threshold (float) – correlation value threshold set to form a cluster
Return type:

numpy matrix containing integers

Returns:

matrix of number of stations

clustering.do_clustering(event, nbsta, nbmin)[source]

First step: all events are examined one by one. For each of them, the indexes of the events where there is a sufficient number of stations (i.e. the number exceeds or equals nbmin) are kept in GRAPH.neighbours. All events are flagged to 0 and belong to cluster 0 by default. The indexes of the events for which neighbours were found are written in the list summits.

From this list (summits), we extract the event which has the greatest number of neighbours. The clustering process will start with that event.

Then the DFS algorithm is applied: when an event is examined, it is flagged to 1 and its cluster number is also updated (GRAPH.cluster_index). The DFS algorithm is recursive and explores each possible path until it ends: it means that when the research starts with a given event, it will search for the the first neighbour, and then the first neighbour of this first neighbour and so on... When all the first neighbours are found, the research can concentrates on second neighbours, and so on until all the neighbours are found.

Once the DFS algorithm has found all events linked to summit_first, we write the corresponding indexes into the dictionary CLUSTER (where the keys are the cluster indexes) and look for events with neighbours which are still not flagged to 1. The process keeps going until the whole events with neighbours belong to a cluster.

The function finally returns the dictionary CLUSTER.

Parameters:
  • event – total number of events in the Waveloc location file
  • nbsta – 2-D matrix containing the number of stations where the cross-correlation value

is greater than a given threshold for all possible event pairs :param nbmin: minimum number of stations where the cross-correlation value should be greater than a given threshold to form a cluster

Return type:dictionary
Returns:indexes of events forming each cluster
clustering.do_clustering_setup_and_run(opdict)[source]

Does clustering by applying the depth first search algorithm and saves the result (= a dictionary containing the event indexes forming each cluster) in a binary file. Needs to define the correlation value threshold and the minimum number of stations where this threshold should be reached to form a cluster (should be done in the options dictionary)

Parameters:opdict – Dictionary of waveloc options
clustering.plot_graphs(locs, stations, nbsta, CLUSTER, nbmin, threshold)[source]

Displays two figures. On the first plot, all events are represented. Those belonging to a cluster are color-coded and labelled with the cluster index. On the second plot, all events are also represented, but those belonging to a cluster are colored in black. The links between events of a same cluster are also plotted and color-coded in function of the number of stations where the correlation value exceeds the threshold. Uses mlab from mayavi

Parameters:
  • locs – list of the whole Waveloc locations (each element of the list is a dictionary)
  • stations – dictionary of stations

:param nbsta:2-D matrix containing the number of stations where the cross-correlation value is greater than a given threshold for all possible event pairs :param CLUSTER: dictionary containing the indexes of the events for each cluster :param nbmin: minimum number of stations where the cross-correlation value should be greater than a given threshold to form a cluster :param threshold: correlation value threshold set to form a cluster

clustering.plot_traces(CLUSTER, delay_file, coeff, locs, datadir, data_files, threshold)[source]

Plots the waveforms of all possible event pairs within a cluster. On the same figure, displays the superimposed waveforms of the event pair for all stations. Also displays the correlation value.

Parameters:
  • CLUSTER (dictionary) – dictionary containing the event indexes belonging to each cluster
  • delay_file (string) – file name of the file containing the time delays
  • coeff (dictionary) – cross-correlation values of all possible event pairs for all stations
  • locs (list) – list of the whole Waveloc locations (each element of the list is a dictionary)
  • datadir (string) – data directory path
  • data_files (list) – list of paths of data files
  • threshold (float) – correlation coefficient threshold
clustering.waveval(stack_time, t_before, t_after, dt, tdeb)[source]

Finds the indexes corresponding to t_before and t_after

Parameters:
  • stack_time (UTCDateTime) – origin time of the event
  • t_before (float) – time taken before the origin time
  • t_after (float) – time taken after the origin time
  • dt (float) – waveform sample distance in seconds
  • tdeb (UTCDateTime) – start time of the waveform
Return type:

integer

Returns:

indexes corresponding to t_before and t_after

double_diff

Routines to do double-difference relative location of events.

double_diff.centroid_constraint(G, d, W)[source]

Applies centroid constraint: sum(delta m)=0.

Parameters:
  • G – matrix of partial derivatives
  • d – vector of double-difference times
  • W – weighting matrix
Returns:

Updated versions of G, d, W

double_diff.coord_cluster(cluster, locs)[source]

Extract the coordinates of the events of a given cluster

Parameters:
  • cluster (list) – indices of events in the cluster
  • locs (list) – list of the whole locations (each element of the list is a dictionary)
Return type:

list

Returns:

xini, yini, zini, zini_ph, to_ini

double_diff.do_double_diff(x, y, z, to, stations, coeff, delay, cluster, threshold, t_th, arr_times)[source]

Do double difference location (inner routine) and return new coordinates.

Parameters:
  • x (list) – x-coordinates of events of a given cluster
  • y (list) – y-coordinates of events of a given cluster
  • z (list) – z-coordinates of events of a given cluster
  • to (list) – origin times of events of a given cluster
  • stations (dictionary) – dictionary of stations
  • coeff (dictionary) – cross-correlation coefficients between all possible pairs of events
  • delay (dictionary) – time delays measured between all possible pairs of events
  • cluster (list) – indices of events in the cluster
  • threshold (float) – minimum value of cross-correlation coefficient used to form a cluster
  • t_th (dictionary) – theoretical traveltimes
  • arr_times (dictionary) – theoretical arrival times
Return type:

list

Returns:

x, y, z, to

double_diff.do_double_diff_setup_and_run(opdict)[source]

Do double difference (outer routine). Takes options from a WavelocOptions.opdict dictionary.

Parameters:opdict – Dictionary of parameters and options
double_diff.fill_matrix(cluster, x, y, z, t_orig, stations, t_th, t_arr, coeff, delay, threshold)[source]

Fills in matrix G (partial derivatives) data-vector d (double differences) and matrix W (weights)

Parameters:
  • cluster (list) – indices of events in the cluster
  • x (list) – x-coordinates of all events
  • y (list) – y-coordinates of all events
  • z (list) – z-coordinates of all events
  • t_orig (list) – origin times of all events
  • stations (dictionary) – dictionary of station coordinates
  • t_th (dictionary) – dictionary of theoretical traveltimes
  • t_arr (dictionary) – dictionary of theoretical arrival times
  • coeff (dictionary) – cross-correlation coefficients between all possible pairs of events
  • delay (dictionary) – time delays measured between all possible pairs of events
  • threshold (float) –
double_diff.inversion(G, d, W)[source]

Inverts G, d, W using LSQR to obtain parameter vector m

Parameters:
  • G – matrix of partial derivatives
  • d – vector of double-difference times
  • W – weighting matrix
Returns:

Parameter vector m

double_diff.partial_deriv(coord, ev, tth)[source]

Computes partial derivatives

Parameters:
  • coord (list) – list of coordinates [x,y,z] of a given station
  • ev (list) – list of coordinates [x,y,z] of the event
  • tth (float) – theoretical traveltime from the event to the station
Return type:

list

Returns:

dpx, dpy, dpz

double_diff.plot_events(cluster, locs, stations, x, y, z, i, threshold, nbmin, area, nbsta)[source]

Plot old and new locations (uses mayavi)

Parameters:
  • cluster (numpy array) – indices of events composing all clusters
  • locs (list) – list of the whole locations (each element of the list is a dictionary)
  • stations (dictionary) – dictionary of stations
  • x (list) – new x-coordinates of events
  • y (list) – new y-coordinates of events
  • z (list) – new z-coordinates of events
  • i (int) – cluster index
  • threshold (float) – minimum value of cross-correlation coefficient used to form a cluster
  • nbmin (int) – minimum number of stations required to form a cluster
  • area (list) – coordinates of the study area
  • nbsta (numpy array) – number of stations where the measured correlation coefficient was greater than the given threshold for all possible event pairs
double_diff.traveltimes(x, y, z, t_orig, stations, time_grids)[source]

Compute theoretical traveltimes and arrival times for a point to all stations using pre-calculated time-grids

Parameters:
  • x (float) – x-coordinate of point
  • y (float) – y-coordinate of point
  • z (float) – z-coordinate of point
  • t_orig (UTCDateTime) – origin time of envent (used to calculate arrival times)
  • stations (string) – list of station names
  • time_grids – dictionary of time-grids
Return type:

dictionary

Returns:

t_th, arr_times

CZ_color

Provides colors for clustering routines.

CZ_color.CZ_Clust_2_color(dt)[source]

TODO: Add doc-string.

Parameters:dt (float) –
Returns:List of 3 values between 0 and 1 for color indexing.
CZ_color.CZ_W_2_color(dt)[source]

TODO: Add doc-string.

Parameters:dt (float) –
Returns:List of 3 values between 0 and 1 for color indexing.