particleman package

Submodules

particleman.core module

Core Stockwell transform and inverse transform functions.

particleman.core.get_TF_arrays(N, Fs=0, hp=0, lp=0)[source]

Make the Stockwell time, frequency arrays for plotting.

Parameters

  • N (int) – Number of samples in the time series.

  • hp (float) – high-pass point in samples (if Fs is not specified) or in Hz (if Fs is specified)

  • lp (float) – low-pass point in samples (if Fs is not specified) or in Hz (if Fs is specified)

  • Fs (float) – sampling rate in Hz

Returns

T, F

Return type

numpy.ndarray (complex, rank 2)

particleman.core.istransform(X, Fs=0, hp=0, lp=0)[source]

Perform inverse Stockwell transform

particleman.core.stransform(x, Fs=0, hp=0, lp=0, return_time_freq=False)[source]

Perform a Stockwell transform on a time-series.

Returns the transform (S), and time (T) and frequency (F) matrices suitable for use with the contour/contourf functions.

Parameters

  • x (numpy.ndarray) – array containing time-series data

  • hp (float) – high-pass point in samples (if Fs is not specified) or in Hz (if Fs is specified)

  • lp (float) – low-pass point in samples (if Fs is not specified) or in Hz (if Fs is specified)

  • Fs (float) – sampling rate in Hz

  • return_time_freq (bool) – If True, also return the correct-sized time and frequency domain tiles.

Returns

  • S (numpy.ndarray (numpy.complex128, rank 2)) – Stockwell transform (S) matrix

  • T, F (numpy.ndarray (float64, rank 2), optional) – Time (T) and frequency (F) matrices.

Examples

Transform a 100 Hz time series

>>> S, T, F = stransform(data, Fs=100, return_time_freq=True)
>>> plt.contourf(T, F, abs(S))

References

particleman.filter module

Time-frequency filtering for seismic waves.

To estimate propagation azimuth for a prograde/retrograde Rayleigh wavetrain:

  1. Make complex component s-transforms, Sn, Se, Sv

  2. Apply phase delay to Sv for normal retrograde motion, advance for prograde motion

  3. Convert complex transforms into MxNx2 real (“vector”) arrays

  4. Calculate theta from theta_r and theta_l

To make a normalized inner product (NIP) filter, given azimuth

  1. Rotate Sn, Se through theta, into radial (Sr) and transverse (St) component transforms

  2. Calculate NIP using the vector Sr and phase-shifted vector Sv

  3. Extract surface waves using NIP >= 0.8 (smoothed) as a filter on Sr, St, Sv, and invert to time domain. Rayleigh waves will be on the time-domain radial, Love on the time-domain transverse.

References

Meza-Fajardo, K. C., Papageorgiou, A. S., and Semblat, J. F. (2015). Identification and Extraction of Surface Waves from Three Component Seismograms Based on the Normalized Inner Product. Bulletin of the Seismological Society of America.

particleman.filter.NIP(Sr, Sv, polarization=None, eps=None)[source]

Get the normalized inner product of two complex MxN stockwell transforms.

Parameters

  • Sr (numpy.ndarray (complex, rank 2)) – The radial and vertical component s-transforms. If the polarization argument is omitted, Sv is assumed to be phase-shifted according to the desired polarization.

  • Sv (numpy.ndarray (complex, rank 2)) – The radial and vertical component s-transforms. If the polarization argument is omitted, Sv is assumed to be phase-shifted according to the desired polarization.

  • polarization (str, optional) – If provided, the Sv will be phase-shifted according to this string before calculating the NIP. ‘retrograde’ will apply a pi/2 phase advance (1j * Sv) ‘prograde’ or ‘linear’ will apply a pi/2 phase delay (-1j * Sv)

  • eps (float, optional) – Tolerance for small denominator values, for numerical stability. Useful for synthetic noise-free data. Authors used 0.04.

Returns

nip – MxN array of floats between -1 and 1.

Return type

numpy.ndarray (rank 2)

References

Equation (16) and (26) from Meza-Fajardo et al. (2015)

particleman.filter.NIP_filter(n, e, v, fs, xpr, polarization, threshold=0.8, width=0.1, eps=None)[source]

Filter a 3-component seismogram based on the NIP criterion.

This is a composite convenience routine that uses sane defaults. If you want to get intermediate products, call internal routines individually.

Parameters

  • n (numpy.ndarray (rank 1)) – Equal-length data arrays for North, East, and Vertical components, respectively.

  • e (numpy.ndarray (rank 1)) – Equal-length data arrays for North, East, and Vertical components, respectively.

  • v (numpy.ndarray (rank 1)) – Equal-length data arrays for North, East, and Vertical components, respectively.

  • fs (float) – Sampling frequency [Hz]

  • xpr (int) – Sense of wave propagation. -1 for westward, 1 for eastward. Try -int(np.sign(np.sin(np.radians(baz)))), unless they’re directly N-S from each other.

  • polarization (str) – ‘retrograde’ to extract normal retrograde Rayleigh waves ‘prograde’ to extract prograde Rayleigh waves ‘linear’ to extract Love waves

  • threshold (float) – The critical value (“x_r”) and width (“Delta x”) for the NIP filter (cosine) taper.

  • width (float) – The critical value (“x_r”) and width (“Delta x”) for the NIP filter (cosine) taper.

  • eps (float) – Tolerance for small NIP denominator values, for numerical stability.

Returns

  • n, e, v (numpy.ndarray (rank 1)) – Filtered north, east, and vertical components.

  • theta (float) – Average propagation azimuth for the given polarization [degrees]. Use this angle to rotate the data to radial and transverse.

Examples

# compare filtered north, east, vertical to original >>> import obspy.signal as signal >>> nf, ef, vf, theta = NIP_filter(n, e, v, fs, xpr) >>> if theta > 180:

nip_baz = theta - 180

else:

nip_baz = theta + 180

>>> rf, tf = signal.rotate_NE_RT(nf, ef, nip_baz)
>>> r, t = signal.rotate_NE_RT(n, e, baz)
particleman.filter.get_filter(nip, polarization, threshold=None, width=0.1)[source]

Get an NIP-based filter that will pass waves of the specified type.

The filter is made from the NIP and cosine taper for the specified wave type. The nip and the polarization type must match.

Parameters

  • nip (numpy.ndarray (real, rank 2)) – The NIP array [-1.0, 1.0]

  • polarization (str) – The type of polarization that was used to calculate the provided NIP. ‘retrograde’, ‘prograde’, or ‘linear’. See “NIP” function.

  • threshold (float) – The cosine taper critical/crossover value (“x_r”) and width (“Delta x”). If not supplied, the default for retrograde polarization is 0.8, and for prograde or linear polarization is 0.2.

  • width (float) – The cosine taper critical/crossover value (“x_r”) and width (“Delta x”). If not supplied, the default for retrograde polarization is 0.8, and for prograde or linear polarization is 0.2.

Returns

The NIP-based filter array [0.0, 1.0] to multiply into the complex Stockwell arrays, before inverse transforming to the time-domain.

Return type

numpy.ndarray (real, rank 2)

References

Equation (27) and (28) from Meza-Fajardo et al. (2015)

particleman.filter.get_shift(polarization)[source]

Get the appropriate pi/2 phase advance or delay for the provided polarization.

Parameters

polarization (str, {'retrograde', 'prograde', 'linear'}) – ‘retrograde’ returns i, for a pi/2 phase advance. ‘prograde’ or ‘linear’ returns -i, for a pi/2 phase delay

Returns

Multiply this value (i or -i) into a complex vertical S-transform to shift its phase.

Return type

numpy.complex128

particleman.filter.instantaneous_azimuth(Sv, Sn, Se, polarization, xpr)[source]

Get instantaneous propagation angle [degrees], under the Rayleigh wave assumption, using the Meza-Fajardo et al. Normalized Inner Product criterion.

Parameters

  • Sv (numpy.ndarray (complex, rank 2)) – The vertical, North, and East component equal-sized complex s-transforms.

  • Sn (numpy.ndarray (complex, rank 2)) – The vertical, North, and East component equal-sized complex s-transforms.

  • Se (numpy.ndarray (complex, rank 2)) – The vertical, North, and East component equal-sized complex s-transforms.

  • polarization (str, {'retrograde', 'prograde', 'linear'}) – ‘retrograde’ will apply a pi/2 phase advance. ‘prograde’ or ‘linear’ will apply a pi/2 phase delay

  • xpr (int) – Sense of propagation. 1 for eastward, -1 for westward. Try -int(np.sign(np.sin(np.radians(baz)))), unless they’re directly N-S from each other.

Returns

az – Instantaneous Rayleigh wave propagation angle [degrees]

Return type

numpy.ndarray (real, rank 2)

References

Equations (19), (20), and (21) from Meza-Fajardo et al. (2015)

“Let us emphatically note that if the sense of propagation of the phase under investigation is not established, prograde or retrograde motion cannot be defined without ambiguity.”

particleman.filter.rotate_NE_RT(Sn, Se, az)[source]

Rotate North and East s-transforms to radial and transverse, through the propagation angle.

Parameters

  • Sn (numpy.ndarray (complex, rank 2)) – Complex, equal-sized s-transform arrays, for North and East components, respectively.

  • Se (numpy.ndarray (complex, rank 2)) – Complex, equal-sized s-transform arrays, for North and East components, respectively.

  • az (float) – Rotation angle [degrees].

Returns

Sr, St – Complex s-transform arrays for radial and transverse components, respectively.

Return type

numpy.ndarray (rank 2)

References

Equation (17) from Meta-Fajardo et al. (2015)

particleman.filter.scalar_azimuth(e, n, vhat)[source]

Time domain estimation the scalar/average azimuth, in degrees.

References

Equations (10) and (12) from Meza-Fajardo et al. (2015)

“If the extracted signal is composed of more than one dispersive wave propagating in distinct, albeit similar, directions, then, equations (10)–(12) should be applied independently to each one of them. By inspecting the Stockwell transform of the signal, the analyst can observe if there are several wavetrains.”

particleman.filter.shift_phase(Sv, polarization)[source]

Phase-shift an s-transform by the appropriate phase shift for prograde/retrograde motion.

Shift is done on a complex MxN array by multiplication with either i or -i (imaginary unit). This is mostly a reference for how to do/interpret phase shifts, as it’s such a simple thing to do outside of a function.

Parameters

  • Sv (numpy.ndarray (complex, rank 2)) –

  • polarization (str, {'retrograde', 'prograde', 'linear'}) – ‘retrograde’ will apply a pi/2 phase advance (normal Rayleigh waves) ‘prograde’ or ‘linear’ will apply a pi/2 phase delay

Returns

Return type

numpy.ndarray (real, rank 2)

References

Pages 5 and 10 from Meta-Fajardo et al. (2015)

particleman.filter.xpr(az)[source]

Get the Meza-Fajardo “xpr” sense-of-propagation of wavefield. propagation azimuth.

Parameters

az (int or float) – Propagation direction in degrees.

Returns

  • 1 for eastward propagation

  • -1 for westward

Notes

If the azimuth is 0 or 180, polarization type may be ambiguous.

particleman.plotting module

Plotting functions for Stockwell transforms and normalized inner-product (NIP) filtering.

These plotting routines, particularly tile_comparison and NIP_filter_plot, are useful to see how much surface wave energy is in a packet, or how much off-great-circle propagation there is.

These functions plot a number of configurations of “tiles.” Each tile consists of the Stockwell transform on top and an aligned time-series waveform below. The transform may have a hatching overlay that would normally correspond to a NIP filter, and the time-series axis may have a reference (gray) trace overlayed, which normally corresponds to the unfiltered trace.

## Channel nomenclature:

` [nevrt][sd][f] `

  • n : north component

  • e : east

  • v : vertical

  • s : scalar rotation (great circle)

  • d : dynamic rotation (NIP estimate)

  • f : NIP filtered

particleman.plotting.NIP_filter_plots(T, F, theta, fs, Sr, St, Sv, r, t, v, rf, tf=None, vf=None, arrivals=None, flim=None, hatch=None, hatchlim=None, fig=None)[source]
def NIP_filter_plots(T, F, theta, fs, Sr, St, Sv, rf, r, vf, v, t, tf=None,

Quad plot of NIP, and 3 tiles of Stockwell transform with NIP filter hatch and filtered+unfiltered time-series for each component.

T, F, theta: numpy.ndarray (ndim 2)

Time, frequency, instantaneous azimuth tiles.

fsfloat

Sampling rate of underlying time-series data.

Sr, St, Svnumpy.ndarray (ndim 2, complex)

Stockwell transform of the radial, transverse, and vertical component data.

r, t, v: numpy.ndarray (ndim 1)

Unfiltered radial, transverse, and vertical component time-series.

rf, tf, vfnumpy.ndarray (ndim 1)

NIP-filtered radial, transverse, and vertical component time-series.

arrivalssequence of (str, float) 2-tuples

Sequence of arrivals to plot, of the form (label, time_in_seconds)

flimtuple

Frequency limits as (fmin, fmax) 2-tuple, in Hz.

hatchnumpy.ndarray (ndim 2)

Optional tile used for cross-hatch visual mask.

hatchlimtuple

Hatch range used to display mask. 2-tuple of floats (hmin, hmax).

fig : matplotlib.Figure

tileslist

List of two-tuples of axes objects (axis_top, axis_bottom) of each tile, clockwise from top-left. Can be unpacked to get all axes as follows:

>>> tiles = NIP_filter_plots(...)
>>> (ax11, ax12), (ax21, ax22), (ax31, ax32), (ax41, ax42) = tiles
particleman.plotting.check_filters(T, F, Sv, Srs, Sts, vsf, rsf, ts, arrivals, flim, clim, dlim, xlim, hatch=None, hatchlim=None, fig=None)[source]
Parameters

  • T (numpy.ndarray (ndim 2)) – The time and freqency domain tiles/grids.

  • F (numpy.ndarray (ndim 2)) – The time and freqency domain tiles/grids.

  • Sv (numpy.ndarray (ndim 2)) – The vertical, scalar-rotated radial, and scalar-rotated transverse Stockwell transform tiles.

  • Srs (numpy.ndarray (ndim 2)) – The vertical, scalar-rotated radial, and scalar-rotated transverse Stockwell transform tiles.

  • Sts (numpy.ndarray (ndim 2)) – The vertical, scalar-rotated radial, and scalar-rotated transverse Stockwell transform tiles.

  • vsf (numpy.ndarray (ndim 1)) – The Stockwell NIP-filtered vertical and radial, and transverse time-series vectors.

  • rsf (numpy.ndarray (ndim 1)) – The Stockwell NIP-filtered vertical and radial, and transverse time-series vectors.

  • ts (numpy.ndarray (ndim 1)) – The Stockwell NIP-filtered vertical and radial, and transverse time-series vectors.

  • arrivals (sequence of (str, float) 2-tuples) – Sequence of arrivals to plot, of the form (label, time_in_seconds)

  • dlim (2-tuple of floats) – Limits on the time-series amplitudes (y axis limits).

  • hatch (numpy.ndarray (ndim 2)) – Optional tile used for hatch mask.

  • hatchlim (tuple) – Hatch range used to display mask. 2-tuple of floats (hmin, hmax).

  • fig (matplotlib.Figure) –

Returns

Return type

matplotlib.Figure

particleman.plotting.compare_waveforms(v, vsf, rs, rsf, ts, arrivals)[source]

Compare the static and dynamically filtered waveforms.

A 3-panel waveform plot of 3 components. Unfiltered waves are in gray, filtered are overplotted in black.

Parameters

  • v (numpy.ndarray (rank 1)) – Unfiltered and static-rotated filtered vertical waveform.

  • vsf (numpy.ndarray (rank 1)) – Unfiltered and static-rotated filtered vertical waveform.

  • rs (numpy.ndarray (rank 1)) – Unfiltered and static-rotated filtered radial waveform.

  • rsf (numpy.ndarray (rank 1)) – Unfiltered and static-rotated filtered radial waveform.

  • ts (numpy.ndarray (rank 1)) – Unfiltered transverse waveform.

particleman.plotting.make_tiles(fig, gs0, full=None)[source]

Give a list of (ax_top, ax_bottom) axis tuples for each SubPlotSpec in gs0.

fulllist

Integer subplotspec numbers for which the ax_top is to take up the whole tile, so no ax_bottom is to be created. Returns these tiles’ axis handles as (ax_top, None).

particleman.plotting.plot_NIP(T, F, nips, fs=1.0, flim=None, fig=None, ax=None)[source]

Plot the normalized inner product tile.

Parameters

  • T (numpy.ndarray (ndim 2)) – Time, frequency, normalized inner-product tiles.

  • F (numpy.ndarray (ndim 2)) – Time, frequency, normalized inner-product tiles.

  • nips (numpy.ndarray (ndim 2)) – Time, frequency, normalized inner-product tiles.

  • fs (float) – Sampling frequency of the underlying time-series data.

  • flim (tuple) – Frequency limits as (fmin, fmax) 2-tuple, in Hz.

Returns

Return type

matplotlib.Figure

particleman.plotting.plot_arrivals(ax, arrivals, dmin, dmax)[source]
arrivalsdict

are a dict of {name: seconds}.

dmin, dmaxfloat

y value in axis “ax” at which labels are plotted.

particleman.plotting.plot_image(T, F, C, hatch=None, hatchlim=None, flim=None, clim=None, fig=None, ax=None, cmap=None, alpha=None)[source]

Plot a fime-frequency image, optionally with a hatched mask.

T, F, Cnumpy.ndarray (ndim 2)

Time and frequency domain arrays, and data arrays.

hatchnumpy.ndarray (rank 2)

Optional array mask for hatching or for opacity.

hatchlimfloat

Optional hatch value cutoff for masking, above which masking is not done. If hatch is provided but hatchlim is not, hatch will be used as an opacity mask, and it values are assumed to be between 0 and 1.

Returns

Return type

matplotlib.collections.QuadMesh from pcolormesh

particleman.plotting.plot_instantaneous_azimuth(T, F, theta, fs=1.0, flim=None, dlim=None, clim=None, fig=None, ax=None)[source]

Plot the instantanous azimuth TF tile using imshow.

Parameters

  • T (numpy.ndarray (ndim 2)) – Time [sec], frequency [Hz], and instantaneous azimuth calculated by stockwell.filter.instantanous_azimuth

  • F (numpy.ndarray (ndim 2)) – Time [sec], frequency [Hz], and instantaneous azimuth calculated by stockwell.filter.instantanous_azimuth

  • theta (numpy.ndarray (ndim 2)) – Time [sec], frequency [Hz], and instantaneous azimuth calculated by stockwell.filter.instantanous_azimuth

  • fs (float) – Sampling rate of data used.

  • flim (tuple (min, max)) – Optional frequency [Hz], time [sec], or color min/max limits for plots.

  • dlim (tuple (min, max)) – Optional frequency [Hz], time [sec], or color min/max limits for plots.

  • clim (tuple (min, max)) – Optional frequency [Hz], time [sec], or color min/max limits for plots.

  • fig (matplotlib.Figure instance) –

  • ax (matplotlib.axis.Axis instance) –

Returns

Return type

matplotlib.Figure

particleman.plotting.plot_tile(fig, ax1, T, F, S, ax2, d1, label1, color1='k', d2=None, label2=None, arrivals=None, flim=None, clim=None, hatch=None, hatchlim=None, dlim=None, amp_fmt='%.2e', cmap=None, alpha=None)[source]

Plot time-frequency pcolormesh tiles above time-aligned aligned time-series.

Parameters

  • fig (matplotlib.Figure) –

  • ax1 (matplotlib.axis.Axis) – Axis for time-frequency pcolormesh tile and optional hatching.

  • ax2 (matplotlib.axis.Axis) – Axis for time-series plot.

  • T (numpy.ndarray (ndim 2)) – Time, frequency, S-transform tiles from stockwell.stransform

  • F (numpy.ndarray (ndim 2)) – Time, frequency, S-transform tiles from stockwell.stransform

  • S (numpy.ndarray (ndim 2)) – Time, frequency, S-transform tiles from stockwell.stransform

  • d1 (numpy.ndarray (ndim 1)) – Time-series, plotted black. Optional d2 plotted gray. These need to be registered in time to T.

  • d2 (numpy.ndarray (ndim 1)) – Time-series, plotted black. Optional d2 plotted gray. These need to be registered in time to T.

  • color1 (str) – Color of plotted d1 line.

  • label1 (str) – Time-series legend label strings.

  • label2 (str) – Time-series legend label strings.

  • arrivals (dict) – Sequence of arrivals to plot, of the form {label: time_in_seconds, …}

  • dlim (2-tuple of floats) – Limits on the time-series amplitudes (y axis limits).

  • hatch (numpy.ndarray (ndim 2)) – Optional tile used for hatch mask.

  • hatchlim (tuple) – Hatch range used to display mask. 2-tuple of floats (hmin, hmax).

  • amp_fmt (str) – Matplotlib format string used for amplitudes in colorbars and axes.

  • cmap (matplotlib.colors.Colormap) –

  • alpha (float) – Optionally, use a white transparency mask for overlying the hatch, with the given alpha value.

Returns

The ax1 image from pcolormesh.

Return type

matplotlib.collections.QuadMesh

Examples

# filtered versus unfiltered radial, and set color limits >>> plot_tile(fig, ax21, T, F, Srs, ax22, rs, ‘unfiltered’, rsf, ‘NIP filtered’,

arrivals=arrivals, flim=(0.0, fmax), clim=(0.0, 5e-5), hatch=sfilt, hatchlim=(0.0, 0.8))

# scalar versus dynamic rotated radial >>> plot_tile(fig, ax21, T, F, Srs, ax22, rs, ‘scalar’, rd, ‘dynamic’, arrivals

flim=(0.0, fmax), clim=(0.0, 5e-5), hatch=dfilt, hatchlim=(0.0, 0.8))

particleman.plotting.tile_comparison(T, F, Sv, Srs, Srd, Sts, Std, v, rs, rd, ts, td, arrivals, flim, clim, dlim, hatch=None, hatchlim=None, fig=None, xlim=None)[source]

Make a 6-panel side-by-side comparison of tiles, such as scalar versus dynamic rotations.

Sv, Srs, Srd, Sts, Stdnumpy.ndarray (ndim 2)

The vertical, radial-scalar, radial-dynamic, transverse-scalar, and transverse-dynamic stockwell transforms.

v, rs, rd, ts, tdnumpy.ndarray (ndim 1)

The corresponding vertical, radial-scalar, radial-dynamic, transverse-scalar, and transverse-dynamic time-series vectors.

arrivalssequence of (str, float) 2-tuples

Sequence of arrivals to plot, of the form (label, time_in_seconds)

flim, clim, dlim, xlimtuple

Frequency, stockwell amplitude, time-series amplitude, and time-series time limits of display. 2-tuples of (min, max) floats.

hatchnumpy.ndarray (ndim 2)

Optional tile used for hatch mask.

hatchlimtuple

Hatch range used to display mask. 2-tuple of floats (hmin, hmax).

Returns

Return type

matplotlib.Figure

particleman.st module

ctypes interface to st.c

particleman.st.ist(X, lo=None, hi=None)[source]
particleman.st.st(data, lo=None, hi=None)[source]

st(x[, lo, hi]) returns the 2d, complex Stockwell transform of the real array x. If lo and hi are specified, only those frequencies (rows) are returned; lo and hi default to 0 and n/2, resp., where n is the length of x.

Stockwell transform of the real array data. The number of time points need not be a power of two. The lo and hi arguments specify the range of frequencies to return, in Hz. If they are both zero, they default to lo = 0 and hi = len / 2. The result is returned in the complex array result, which must be preallocated, with n rows and len columns, where n is hi - lo + 1. For the default values of lo and hi, n is len / 2 + 1.

Module contents

Stockwell transform and its inverse.