trace

This module provides basic signal processing for seismic traces.

class Trace(network='', station='STA', location='', channel='', tmin=0.0, tmax=None, deltat=1.0, ydata=None, mtime=None, meta=None)

Create new trace object.

A Trace object represents a single continuous strip of evenly sampled time series data. It is built from a 1D NumPy array containing the data samples and some attributes describing its beginning and ending time, its sampling rate and four string identifiers (its network, station, location and channel code).

Parameters:

• network – network code
• station – station code
• location – location code
• channel – channel code
• tmin – system time of first sample in [s]
• tmax – system time of last sample in [s] (if set to None it is computed from length of ydata)
• deltat – sampling interval in [s]
• ydata – 1D numpy array with data samples (can be None when tmax is not None)
• mtime – optional modification time
• meta – additional meta information (not used, but maintained by the library)

The length of the network, station, location and channel codes is not resricted by this software, but data formats like SAC, Mini-SEED or GSE have different limits on the lengths of these codes. The codes set here are silently truncated when the trace is stored

♦ network

str, default: ''

♦ station

str, default: 'STA'

♦ location

str, default: ''

♦ channel

str, default: ''

♦ tmin

builtins.float (pyrocko.guts.Timestamp), default: 0.0

♦ tmax

builtins.float (pyrocko.guts.Timestamp)

♦ deltat

float, default: 1.0

♦ ydata

numpy.ndarray (pyrocko.guts_array.Array), optional

♦ mtime

builtins.float (pyrocko.guts.Timestamp), optional

name()

Get a short string description.

interpolate(t, clip=False)

Value of trace between supporting points through linear interpolation.

Parameters:

• t – time instant
• clip – whether to clip indices to trace ends

index_clip(i)

Clip index to valid range.

add(other, interpolate=True, left=0.0, right=0.0)

Add values of other trace (self += other).

Add values of other trace to the values of self, where it intersects with other. This method does not change the extent of self. If interpolate is True (the default), the values of other to be added are interpolated at sampling instants of self. Linear interpolation is performed. In this case the sampling rate of other must be equal to or lower than that of self. If interpolate is False, the sampling rates of the two traces must match.

mult(other, interpolate=True)

Muliply with values of other trace (self *= other).

Multiply values of other trace to the values of self, where it intersects with other. This method does not change the extent of self. If interpolate is True (the default), the values of other to be multiplied are interpolated at sampling instants of self. Linear interpolation is performed. In this case the sampling rate of other must be equal to or lower than that of self. If interpolate is False, the sampling rates of the two traces must match.

max()

Get time and value of data maximum.

min()

Get time and value of data minimum.

absmax()

Get time and value of maximum of the absolute of data.

set_codes(network=None, station=None, location=None, channel=None)

Set network, station, location, and channel codes.

overlaps(tmin, tmax)

Check if trace has overlap with a given time span.

is_relevant(tmin, tmax, selector=None)

Check if trace has overlap with a given time span and matches a condition callback. (internal use)

set_mtime(mtime)

Set modification time of the trace.

get_xdata()

Create array for time axis.

get_ydata()

Get data array.

set_ydata(new_ydata)

Replace data array.

drop_data()

Forget data, make dataless trace.

drop_growbuffer()

Detach the traces grow buffer.

copy(data=True)

Make a deep copy of the trace.

crop_zeros()

Remove any zeros at beginning and end.

append(data)

Append data to the end of the trace.

To make this method efficient when successively very few or even single samples are appended, a larger grow buffer is allocated upon first invocation. The traces data is then changed to be a view into the currently filled portion of the grow buffer array.

chop(tmin, tmax, inplace=True, include_last=False, snap=(<built-in function round>, <built-in function round>), want_incomplete=True)

Cut the trace to given time span.

If the inplace argument is True (the default) the trace is cut in place, otherwise a new trace with the cut part is returned. By default, the indices where to start and end the trace data array are determined by rounding of tmin and tmax to sampling instances using Python’s round() function. This behaviour can be changed with the snap argument, which takes a tuple of two functions (one for the lower and one for the upper end) to be used instead of round(). The last sample is by default not included unless include_last is set to True. If the given time span exceeds the available time span of the trace, the available part is returned, unless want_incomplete is set to False - in that case, a NoData exception is raised. This exception is always raised, when the requested time span does dot overlap with the trace’s time span.

downsample(ndecimate, snap=False, initials=None, demean=False)

Downsample trace by a given integer factor.

Parameters:

• ndecimate – decimation factor, avoid values larger than 8
• snap – whether to put the new sampling instances closest to multiples of the sampling rate.
• initials – None, True, or initial conditions for the anti-aliasing filter, obtained from a previous run. In the latter two cases the final state of the filter is returned instead of None.
• demean – whether to demean the signal before filtering.

downsample_to(deltat, snap=False, allow_upsample_max=1, initials=None, demean=False)

Downsample to given sampling rate.

Tries to downsample the trace to a target sampling interval of deltat. This runs the Trace.downsample() one or several times. If allow_upsample_max is set to a value larger than 1, intermediate upsampling steps are allowed, in order to increase the number of possible downsampling ratios.

If the requested ratio is not supported, an exception of type pyrocko.util.UnavailableDecimation is raised.

resample(deltat)

Resample to given sampling rate deltat.

Resampling is performed in the frequency domain.

stretch(tmin_new, tmax_new)

Stretch signal while preserving sample rate using sinc interpolation.

Parameters:

• tmin_new – new time of first sample
• tmax_new – new time of last sample

This method can be used to correct for a small linear time drift or to introduce sub-sample time shifts. The amount of stretching is limited to 10% by the implementation and is expected to be much smaller than that by the approximations used.

nyquist_check(frequency, intro='Corner frequency', warn=True, raise_exception=False)

Check if a given frequency is above the Nyquist frequency of the trace.

Parameters:

• intro – string used to introduce the warning/error message
• warn – whether to emit a warning
• raise_exception – whether to raise an AboveNyquist exception.

lowpass(order, corner, nyquist_warn=True, nyquist_exception=False, demean=True)

Apply Butterworth lowpass to the trace.

Parameters:

• order – order of the filter
• corner – corner frequency of the filter

Mean is removed before filtering.

highpass(order, corner, nyquist_warn=True, nyquist_exception=False, demean=True)

Apply butterworth highpass to the trace.

Parameters:

• order – order of the filter
• corner – corner frequency of the filter

Mean is removed before filtering.

bandpass(order, corner_hp, corner_lp, demean=True)

Apply butterworth bandpass to the trace.

Parameters:

• order – order of the filter
• corner_hp – lower corner frequency of the filter
• corner_lp – upper corner frequency of the filter

Mean is removed before filtering.

envelope(inplace=True)

Calculate the envelope of the trace.

Parameters:inplace – calculate envelope in place

The calculation follows:

where H is the Hilbert-Transform of the signal Y.

taper(taperer, inplace=True, chop=False)

Apply a Taper to the trace.

Parameters:

• taperer – instance of Taper subclass
• inplace – apply taper inplace
• chop – if True: exclude tapered parts from the resulting trace

whiten(order=6)

Whiten signal in time domain using autoregression and recursive filter.

Parameters:order – order of the autoregression process

ampspec_whiten(width, td_taper='auto', fd_taper='auto', pad_to_pow2=True, demean=True)

Whiten signal via frequency domain using moving average on amplitude spectra.

Parameters:

• width – width of smoothing kernel [Hz]
• td_taper – time domain taper, object of type Taper or None or 'auto'.
• fd_taper – frequency domain taper, object of type Taper or None or 'auto'.
• pad_to_pow2 – whether to pad the signal with zeros up to a length of 2^n
• demean – whether to demean the signal before tapering

The signal is first demeaned and then tapered using td_taper. Then, the spectrum is calculated and inversely weighted with a smoothed version of its amplitude spectrum. A moving average is used for the smoothing. The smoothed spectrum is then tapered using fd_taper. Finally, the smoothed and tapered spectrum is back-transformed into the time domain.

If td_taper is set to 'auto', CosFader(1.0/width) is used. If fd_taper is set to 'auto', CosFader(width) is used.

bandpass_fft(corner_hp, corner_lp)

Apply boxcar bandbpass to trace (in spectral domain).

shift(tshift)

Time shift the trace.

snap(inplace=True, interpolate=False)

Shift trace samples to nearest even multiples of the sampling rate.

Parameters:inplace – (boolean) snap traces inplace

If inplace is False and the difference of tmin and tmax of both, the snapped and the original trace is smaller than 0.01 x deltat, snap() returns the unsnapped instance of the original trace.

fix_deltat_rounding_errors()

Try to undo sampling rate rounding errors.

See fix_deltat_rounding_errors().

sta_lta_centered(tshort, tlong, quad=True, scalingmethod=1)

Run special STA/LTA filter where the short time window is centered on the long time window.

Parameters:

• tshort – length of short time window in [s]
• tlong – length of long time window in [s]
• quad – whether to square the data prior to applying the STA/LTA filter
• scalingmethod – integer key to select how output values are scaled / normalized (1, 2, or 3)

Scalingmethod	Implementation	Range
1	As/Al* Tl/Ts	[0,1]
2	(As/Al - 1) / (Tl/Ts - 1)	[-Ts/Tl,1]
3	Like 2 but clipping range at zero	[0,1]

sta_lta_right(tshort, tlong, quad=True, scalingmethod=1)

Run special STA/LTA filter where the short time window is overlapping with the last part of the long time window.

Parameters:

• tshort – length of short time window in [s]
• tlong – length of long time window in [s]
• quad – whether to square the data prior to applying the STA/LTA filter
• scalingmethod – integer key to select how output values are scaled / normalized (1, 2, or 3)

Scalingmethod	Implementation	Range
1	As/Al* Tl/Ts	[0,1]
2	(As/Al - 1) / (Tl/Ts - 1)	[-Ts/Tl,1]
3	Like 2 but clipping range at zero	[0,1]

With scalingmethod=1, the values produced by this variant of the STA/LTA are equivalent to

where are the input samples, are the number of samples in the short time window and are the number of samples in the long time window.

peaks(threshold, tsearch, deadtime=False, nblock_duration_detection=100)

Detect peaks above a given threshold (method 1).

From every instant, where the signal rises above threshold, a time length of tsearch seconds is searched for a maximum. A list with tuples (time, value) for each detected peak is returned. The deadtime argument turns on a special deadtime duration detection algorithm useful in combination with recursive STA/LTA filters.

peaks2(threshold, tsearch)

Detect peaks above a given threshold (method 2).

This variant of peak detection is a bit more robust (and slower) than the one implemented in Trace.peaks(). First all samples with a[i-1] < a[i] > a[i+1] are masked as potential peaks. From these, iteratively the one with the maximum amplitude a[j] and time t[j] is choosen and potential peaks within t[j] - tsearch, t[j] + tsearch are discarded. The algorithm stops, when a[j] < threshold or when no more potential peaks are left.

extend(tmin=None, tmax=None, fillmethod='zeros')

Extend trace to given span.

Parameters:

• tmin – begin time of new span
• tmax – end time of new span
• fillmethod – 'zeros', 'repeat', 'mean', or 'median'

transfer(tfade=0.0, freqlimits=None, transfer_function=None, cut_off_fading=True, demean=False, invert=False)

Return new trace with transfer function applied (convolution).

Parameters:

• tfade – rise/fall time in seconds of taper applied in timedomain at both ends of trace.
• freqlimits – 4-tuple with corner frequencies in Hz.
• transfer_function – FrequencyResponse object; must provide a method ‘evaluate(freqs)’, which returns the transfer function coefficients at the frequencies ‘freqs’.
• cut_off_fading – whether to cut off rise/fall interval in output trace.
• demean – remove mean before applying transfer function
• invert – set to True to do a deconvolution

differentiate(n=1, order=4, inplace=True)

Approximate first or second derivative of the trace.

Parameters:

• n – 1 for first derivative, 2 for second
• order – order of the approximation 2 and 4 are supported
• inplace – if True the trace is differentiated in place, otherwise a new trace object with the derivative is returned.

Raises ValueError for unsupported n or order.

See diff_fd() for implementation details.

misfit(candidate, setup, nocache=False, debug=False)

Calculate misfit and normalization factor against candidate trace.

Parameters:

• candidate – Trace object
• setup – MisfitSetup object

Returns:

tuple (m, n), where m is the misfit value and n is the normalization divisor

If the sampling rates of self and candidate differ, the trace with the higher sampling rate will be downsampled.

spectrum(pad_to_pow2=False, tfade=None)

Get FFT spectrum of trace.

Parameters:

• pad_to_pow2 – whether to zero-pad the data to next larger power-of-two length
• tfade – None or a time length in seconds, to apply cosine shaped tapers to both

Returns:

a tuple with (frequencies, values)

fill_template(template, **additional)

Fill string template with trace metadata.

Uses normal python ‘%(placeholder)s’ string templates. The following placeholders are considered: network, station, location, channel, tmin (time of first sample), tmax (time of last sample), tmin_ms, tmax_ms, tmin_us, tmax_us, tmin_year, tmax_year, julianday. The variants ending with '_ms' include milliseconds, those with '_us' include microseconds, those with '_year' contain only the year.

plot()

Show trace with matplotlib.

See also: Trace.snuffle().

snuffle(**kwargs)

Show trace in a snuffler window.

Parameters:

• stations – list of pyrocko.model.Station objects or None
• events – list of pyrocko.model.Event objects or None
• markers – list of pyrocko.gui.util.Marker objects or None
• ntracks – float, number of tracks to be shown initially (default: 12)
• follow – time interval (in seconds) for real time follow mode or None
• controls – bool, whether to show the main controls (default: True)
• opengl – bool, whether to use opengl (default: False)

snuffle(traces, **kwargs)

Show traces in a snuffler window.

Parameters:

• stations – list of pyrocko.model.Station objects or None
• events – list of pyrocko.model.Event objects or None
• markers – list of pyrocko.gui.util.Marker objects or None
• ntracks – float, number of tracks to be shown initially (default: 12)
• follow – time interval (in seconds) for real time follow mode or None
• controls – bool, whether to show the main controls (default: True)
• opengl – bool, whether to use opengl (default: False)

exception InfiniteResponse

This exception is raised by Trace operations when deconvolution of a frequency response (instrument response transfer function) would result in a division by zero.

exception MisalignedTraces

This exception is raised by some Trace operations when tmin, tmax or number of samples do not match.

exception NoData

This exception is raised by some Trace operations when no or not enough data is available.

exception AboveNyquist

This exception is raised by some Trace operations when given frequencies are above the Nyquist frequency.

exception TraceTooShort

This exception is raised by some Trace operations when the trace is too short.

exception ResamplingFailed

minmax(traces, key=None, mode='minmax')

Get data range given traces grouped by selected pattern.

Parameters:

• key – a callable which takes as single argument a trace and returns a key for the grouping of the results. If this is None, the default, lambda tr: (tr.network, tr.station, tr.location, tr.channel) is used.
• mode – ‘minmax’ or floating point number. If this is ‘minmax’, minimum and maximum of the traces are used, if it is a number, mean +- standard deviation times mode is used.

Returns:

a dict with the combined data ranges.

Examples:

ranges = minmax(traces, lambda tr: tr.channel)
print ranges['N']   # print min & max of all traces with channel == 'N'
print ranges['E']   # print min & max of all traces with channel == 'E'

ranges = minmax(traces, lambda tr: (tr.network, tr.station))
print ranges['GR', 'HAM3']  # print min & max of all traces with
                            # network == 'GR' and station == 'HAM3'

ranges = minmax(traces, lambda tr: None)
print ranges[None]  # prints min & max of all traces

minmaxtime(traces, key=None)

Get time range given traces grouped by selected pattern.

Parameters:key – a callable which takes as single argument a trace and returns a key for the grouping of the results. If this is None, the default, lambda tr: (tr.network, tr.station, tr.location, tr.channel) is used. Returns:a dict with the combined data ranges.

degapper(traces, maxgap=5, fillmethod='interpolate', deoverlap='use_second', maxlap=None)

Try to connect traces and remove gaps.

This method will combine adjacent traces, which match in their network, station, location and channel attributes. Overlapping parts are handled according to the deoverlap argument.

Parameters:

• traces – input traces, must be sorted by their full_id attribute.
• maxgap – maximum number of samples to interpolate.
• fillmethod – what to put into the gaps: ‘interpolate’ or ‘zeros’.
• deoverlap – how to handle overlaps: ‘use_second’ to use data from second trace (default), ‘use_first’ to use data from first trace, ‘crossfade_cos’ to crossfade with cosine taper, ‘add’ to add amplitude values.
• maxlap – maximum number of samples of overlap which are removed

Returns:

list of traces

rotate(traces, azimuth, in_channels, out_channels)

2D rotation of traces.

Parameters:

• traces – list of input traces
• azimuth – difference of the azimuths of the component directions (azimuth of out_channels[0]) - (azimuth of in_channels[0])
• in_channels – names of the input channels (e.g. ‘N’, ‘E’)
• out_channels – names of the output channels (e.g. ‘R’, ‘T’)

Returns:

list of rotated traces

rotate_to_lqt(traces, backazimuth, incidence, in_channels, out_channels=('L', 'Q', 'T'))

Rotate traces from ZNE to LQT system.

Parameters:

• traces – list of traces in arbitrary order
• backazimuth – backazimuth in degrees clockwise from north
• incidence – incidence angle in degrees from vertical
• in_channels – input channel names
• out_channels – output channel names (default: (‘L’, ‘Q’, ‘T’))

Returns:

list of transformed traces

project(traces, matrix, in_channels, out_channels)

Affine transform of three-component traces.

Compute matrix-vector product of three-component traces, to e.g. rotate traces into a different basis. The traces are distinguished and ordered by their channel attribute. The tranform is applied to overlapping parts of any appropriate combinations of the input traces. This should allow this function to be robust with data gaps. It also tries to apply the tranformation to subsets of the channels, if this is possible, so that, if for example a vertical compontent is missing, horizontal components can still be rotated.

Parameters:

• traces – list of traces in arbitrary order
• matrix – tranformation matrix
• in_channels – input channel names
• out_channels – output channel names

Returns:

list of transformed traces

project_dependencies(matrix, in_channels, out_channels)

Figure out what dependencies project() would produce.

correlate(a, b, mode='valid', normalization=None, use_fft=False)

Cross correlation of two traces.

Parameters:

• a,b – input traces
• mode – 'valid', 'full', or 'same'
• normalization – 'normal', 'gliding', or None
• use_fft – bool, whether to do cross correlation in spectral domain

Returns:

trace containing cross correlation coefficients

This function computes the cross correlation between two traces. It evaluates the discrete equivalent of

where the star denotes complex conjugate. Note, that the arguments here are swapped when compared with the numpy.correlate() function, which is internally called. This function should be safe even with older versions of NumPy, where the correlate function has some problems.

A trace containing the cross correlation coefficients is returned. The time information of the output trace is set so that the returned cross correlation can be viewed directly as a function of time lag.

Example:

# align two traces a and b containing a time shifted similar signal:
c = pyrocko.trace.correlate(a,b)
t, coef = c.max()  # get time and value of maximum
b.shift(-t)        # align b with a

same_sampling_rate(a, b, eps=1e-06)

Check if two traces have the same sampling rate.

Parameters:

• a,b – input traces
• eps – relative tolerance

fix_deltat_rounding_errors(deltat)

Try to undo sampling rate rounding errors.

Fix rounding errors of sampling intervals when these are read from single precision floating point values.

Assumes that the true sampling rate or sampling interval was an integer value. No correction will be applied if this would change the sampling rate by more than 0.001%.

merge_codes(a, b, sep='-')

Merge network-station-location-channel codes of a pair of traces.

class Taper(**kwargs)

Base class for tapers.

Does nothing by default.

class CosTaper(a, b, c, d)

Cosine Taper.

Parameters:

• a – start of fading in
• b – end of fading in
• c – start of fading out
• d – end of fading out

♦ a

float

♦ b

float

♦ c

float

♦ d

float

class CosFader(xfade=None, xfrac=None)

Cosine Fader.

Parameters:

• xfade – fade in and fade out time in seconds (optional)
• xfrac – fade in and fade out as fraction between 0. and 1. (optional)

Only one argument can be set. The other should to be None.

♦ xfade

float, optional

♦ xfrac

float, optional

class MultiplyTaper(tapers=None)

Multiplication of several tapers.

♦ tapers

list of Taper objects, default: []

class GaussTaper(alpha)

Frequency domain Gaussian filter.

♦ alpha

float

class FrequencyResponse(**kwargs)

Evaluates frequency response at given frequencies.

is_scalar()

Check if this is a flat response.

class Evalresp(respfile, trace=None, target='dis', nslc_id=None, time=None)

Calls evalresp and generates values of the instrument response transfer function.

Parameters:

• respfile – response file in evalresp format
• trace – trace for which the response is to be extracted from the file
• target – 'dis' for displacement or 'vel' for velocity

♦ respfile

str

♦ nslc_id

tuple of 4 str objects, default: (None, None, None, None)

♦ target

str, default: 'dis'

♦ instant

float

class InverseEvalresp(respfile, trace, target='dis')

Calls evalresp and generates values of the inverse instrument response for deconvolution of instrument response.

Parameters:

• respfile – response file in evalresp format
• trace – trace for which the response is to be extracted from the file
• target – 'dis' for displacement or 'vel' for velocity

♦ respfile

str

♦ nslc_id

tuple of 4 str objects, default: (None, None, None, None)

♦ target

str, default: 'dis'

♦ instant

float

class PoleZeroResponse(zeros=None, poles=None, constant=(1+0j))

Evaluates frequency response from pole-zero representation.

Parameters:

• zeros – numpy.array containing complex positions of zeros
• poles – numpy.array containing complex positions of poles
• constant – gain as floating point number

                  (j*2*pi*f - zeros[0]) * (j*2*pi*f - zeros[1]) * ...
T(f) = constant * ----------------------------------------------------
                  (j*2*pi*f - poles[0]) * (j*2*pi*f - poles[1]) * ...

The poles and zeros should be given as angular frequencies, not in Hz.

♦ zeros

list of complex objects, default: []

♦ poles

list of complex objects, default: []

♦ constant

complex, default: (1+0j)

is_scalar()

Check if this is a flat response.

class ButterworthResponse(**kwargs)

Butterworth frequency response.

Parameters:

• corner – corner frequency of the response
• order – order of the response
• type – either high or low

♦ corner

float, default: 1.0

♦ order

int, default: 4

♦ type

builtins.str (pyrocko.guts.StringChoice), default: 'low'

class SampledResponse(frequencies, values, left=None, right=None)

Interpolates frequency response given at a set of sampled frequencies.

Parameters:

• frequencies,values – frequencies and values of the sampled response function.
• left,right – values to return when input is out of range. If set to None (the default) the endpoints are returned.

♦ frequencies

numpy.ndarray (pyrocko.guts_array.Array)

♦ values

numpy.ndarray (pyrocko.guts_array.Array)

♦ left

complex, optional

♦ right

complex, optional

inverse()

Get inverse as a new SampledResponse object.

class IntegrationResponse(n=1, gain=1.0)

The integration response, optionally multiplied by a constant gain.

Parameters:

• n – exponent (integer)
• gain – gain factor (float)

            gain
T(f) = --------------
       (j*2*pi * f)^n

♦ n

int, optional, default: 1

♦ gain

float, optional, default: 1.0

class DifferentiationResponse(n=1, gain=1.0)

The differentiation response, optionally multiplied by a constant gain.

Parameters:

• n – exponent (integer)
• gain – gain factor (float)

T(f) = gain * (j*2*pi * f)^n

♦ n

int, optional, default: 1

♦ gain

float, optional, default: 1.0

class AnalogFilterResponse(b, a)

Frequency response of an analog filter.

(see scipy.signal.freqs()).

♦ b

list of float objects, default: []

♦ a

list of float objects, default: []

class MultiplyResponse(responses=None)

Multiplication of several FrequencyResponse objects.

♦ responses

list of FrequencyResponse objects, default: []

is_scalar()

Check if this is a flat response.

numpy_has_correlate_flip_bug()

Check if NumPy’s correlate function reveals old behaviour

numpy_correlate_fixed(a, b, mode='valid', use_fft=False)

Call numpy.correlate() with fixes.

c[k] = sum_i a[i+k] * conj(b[i])

Note that the result produced by newer numpy.correlate is always flipped with respect to the formula given in its documentation (if ascending k assumed for the output).

numpy_correlate_emulate(a, b, mode='valid')

Slow version of numpy.correlate() for comparison.

numpy_correlate_lag_range(a, b, mode='valid', use_fft=False)

Get range of lags for which numpy.correlate() produces values.

autocorr(x, nshifts)

Compute biased estimate of the first autocorrelation coefficients.

Parameters:

• x – input array
• nshifts – number of coefficients to calculate

yulewalker(x, order)

Compute autoregression coefficients using Yule-Walker method.

Parameters:

• x – input array
• order – number of coefficients to produce

A biased estimate of the autocorrelation is used. The Yule-Walker equations are solved by numpy.linalg.inv() instead of Levinson-Durbin recursion which is normally used.

hilbert(x, N=None)

Return the hilbert transform of x of length N.

(from scipy.signal, but changed to use fft and ifft from numpy.fft)

class States

Utility to store channel-specific state in coroutines.

exception ScipyBug

co_lfilter(*args, **kwargs)

Successively filter broken continuous trace data (coroutine).

Create coroutine which takes Trace objects, filters their data through scipy.signal.lfilter() and sends new Trace objects containing the filtered data to target. This is useful, if one wants to filter a long continuous time series, which is split into many successive traces without producing filter artifacts at trace boundaries.

Filter states are kept per channel, specifically, for each (network, station, location, channel) combination occuring in the input traces, a separate state is created and maintained. This makes it possible to filter multichannel or multistation data with only one co_lfilter() instance.

Filter state is reset, when gaps occur.

Use it like this:

from pyrocko.trace import co_lfilter, co_list_append

filtered_traces = []
pipe = co_lfilter(co_list_append(filtered_traces), a, b)
for trace in traces:
     pipe.send(trace)

pipe.close()

co_downsample(target, q, n=None, ftype='fir')

Successively downsample broken continuous trace data (coroutine).

Create coroutine which takes Trace objects, downsamples their data and sends new Trace objects containing the downsampled data to target. This is useful, if one wants to downsample a long continuous time series, which is split into many successive traces without producing filter artifacts and gaps at trace boundaries.

Filter states are kept per channel, specifically, for each (network, station, location, channel) combination occuring in the input traces, a separate state is created and maintained. This makes it possible to filter multichannel or multistation data with only one co_lfilter() instance.

Filter state is reset, when gaps occur. The sampling instances are choosen so that they occur at (or as close as possible) to even multiples of the sampling interval of the downsampled trace (based on system time).

class DomainChoice(dummy) → str

Any str out of ['time_domain', 'frequency_domain', 'envelope', 'absolute', 'cc_max_norm'].

class MisfitSetup(**kwargs)

Contains misfit setup to be used in trace.misfit()

Parameters:

• description – Description of the setup
• norm – L-norm classifier
• taper – Object of Taper
• filter – Object of FrequencyResponse
• domain – [‘time_domain’, ‘frequency_domain’, ‘envelope’, ‘absolute’, ‘cc_max_norm’]

Can be dumped to a yaml file.

♦ description

str, optional

♦ norm

int

♦ taper

Taper

♦ filter

FrequencyResponse, optional

♦ domain

builtins.str (DomainChoice), default: 'time_domain'

equalize_sampling_rates(trace_1, trace_2)

Equalize sampling rates of two traces (reduce higher sampling rate to lower).

Parameters:

• trace_1 – Trace object
• trace_2 – Trace object

Returns a copy of the resampled trace if resampling is needed.

Lx_norm(u, v, norm=2)

Calculate the misfit denominator m and the normalization devisor n according to norm.

The normalization divisor n is calculated from v.

Parameters:

• u – numpy.array
• v – numpy.array
• norm – (default = 2)

u and v must be of same size.