seismic.receiver_fn package

Submodules

seismic.receiver_fn.bulk_rf_report module

Produce PDF report of network stations showing RF waveforms

seismic.receiver_fn.generate_rf module

Generate Receiver Functions (RF) from collection of 3-channel seismic traces.

seismic.receiver_fn.generate_rf.event_waveforms_to_rf(input_file, output_file, config)[source]

Main entry point for generating RFs from event traces.

Config file consists of 3 sub-dictionaries. One named “filtering” for input stream filtering settings, one named “processing” for RF processing settings, and one named “system” for options on how the system will run the job. Each of these sub-dicts is described below:

"filtering":  # Filtering settings
{
  "resample_rate": float # Resampling rate in Hz
  "taper_limit": float   # Fraction of signal to taper at end, between 0 and 0.5
  "filter_band": (float, float) # Filter pass band (Hz). Not required for freq-domain deconvolution.
  "channel_pattern": # Ordered list of preferred channels, e.g. 'HH*,BH*',
                     # where channel selection is ambiguous.
  "baz_range": (float, float) or [(float, float), ...] # Discrete ranges of source back azimuth to use (degrees).
      # Each value must be between 0 and 360. May be a pair or a list of pairs for multiple ranges.
}

"processing":  # RF processing settings
{
  "custom_preproc":
  {
    "import": 'import custom symbols',  # statement to import required symbols
    "func": 'preproc functor'  # expression to get handle to custom preprocessing functor
    "args": {}  # additional kwargs to pass to func
  }
  "trim_start_time": float # Trace trim start time in sec, relative to onset
  "trim_end_time": float # Trace trim end time in sec, relative to onset
  "rotation_type": str # Choice of ['zrt', 'lqt']. Rotational coordinate system
                       # for aligning ZNE trace components with incident wave direction
  "deconv_domain": str # Choice of ['time', 'freq', 'iter']. Whether to perform deconvolution
                       # in time or freq domain, or iterative technique
  "gauss_width": float # Gaussian freq domain filter width. Only required for freq-domain deconvolution
  "water_level": float # Water-level for freq domain spectrum. Only required for freq-domain deconvolution
  "spiking": float # Spiking factor (noise suppression), only required for time-domain deconvolution
  "normalize": bool # Whether to normalize RF amplitude
}

"system":  # job run settings
{
  "parallel": bool # Use parallel execution
  "memmap": bool # Memmap input file for improved performance in data reading thread.
                 # Useful when data input is bottleneck, if system memory permits.
  "temp_dir": str or path # Temporary directory to use for best performance
  "aggressive_dispatch": bool # Dispatch all worker jobs as aggressively as possible to minimize
                              # chance of worker being starved of work. Uses more memory.
}
Parameters
  • input_file (str or pathlib.Path) – Event waveform source file for seismograms, generated using extract_event_traces.py script

  • output_file (str or pathlib.Path) – Name of hdf5 file to produce containing RFs

  • config (dict) – Dictionary containing job configuration parameters

Returns

None

seismic.receiver_fn.generate_rf.transform_stream_to_rf(ev_id, stream3c, config_filtering, config_processing, **kwargs)[source]

Generate P-phase receiver functions for a single 3-channel stream. See documentation for function event_waveforms_to_rf for details of config dictionary contents.

Parameters
  • ev_id (int or str) – The event id

  • stream3c (rf.RFStream) – Stream with 3 components of trace data

  • config_filtering (dict) – Dictionary containing stream filtering settings

  • config_processing (dict) – Dictionary containing RF processing settings

  • kwargs (dict) – Keyword arguments that will be passed to filtering and deconvolution functions.

Returns

RFstream containing receiver function if successful, None otherwise

Return type

rf.RFStream or NoneType

seismic.receiver_fn.generate_rf.transform_stream_to_rf_queue(oqueue, ev_id, stream3c, config_filtering, config_processing, **kwargs)[source]

Process using transform_stream_to_rf and queue results for further handling.

Parameters
  • oqueue (queue or multiprocessing.Manager.Queue) – Output queue where filtered streams are queued

  • ev_id – Event ID to pass on to transform_stream_to_rf

  • stream3c – 3-channel stream to process

  • config_filtering – Filtering settings

  • config_processing – Processing settings

  • kwargs (dict) – Keyword arguments that will be passed to filtering and deconvolution functions.

Returns

True if stream containing receiver function is pushed into output queue oqueue, False otherwise

Return type

bool

seismic.receiver_fn.plot_ccp module

Generate common conversion point (CCP) plot as per C.Sippl, “Moho geometry along a north-south passive seismic transect through Central Australia”, Technophysics 676 (2016), pp.56-69, DOI https://doi.org/10.1016/j.tecto.2016.03.031

This code adapted from Christian Sippl’s original code.

Workflow:

extract_event_traces.py –> generate_rf.py –> rf_quality_filter.py –> plot_ccp.py (this script)

Example usage:

python seismic/receiver_fn/plot_ccp.py --start-latlon -19.5 133.0 --end-latlon -19.5 140.0 --width 120         --channels T --stacked-scale 0.3 --title "Network OA CCP T-stacking (profile BS24-CF24)"         /software/hiperseis/seismic/receiver_fn/DATA/OA-ZRT-cleaned.h5         /software/hiperseis/seismic/receiver_fn/DATA/OA-ZRT-T_CCP_stack_BS24-CF24_2km_spacing.png
seismic.receiver_fn.plot_ccp.add_ccp_trace(trace, inc_p, matrx, matrx_entry, vmod, depstep, lenstep, sta_offset, az)[source]

project amplitudes from all RFs onto the profile…2D rot:

seismic.receiver_fn.plot_ccp.angular_distance(p1, p2)[source]

Compute the angular distance (in degrees) between two points p1 and p2 using Haversine formula.

Math reference: https://www.movable-type.co.uk/scripts/latlong.html

Parameters
seismic.receiver_fn.plot_ccp.bearing(p1, p2)[source]

Compute bearing (forward azimuth) in degrees from p1 to p2.

Math reference: https://www.movable-type.co.uk/scripts/latlong.html

Parameters
seismic.receiver_fn.plot_ccp.bounding_box(startpoint, endpoint)[source]

Compute a bounding box from start and end points.

Parameters
  • startpoint (pair of float) – Coordinates of starting point

  • endpoint (pair of float) – Coordinates of end point

Returns

Bounding box (left, bottom, right, top)

Return type

tuple(float, float, float, float)

seismic.receiver_fn.plot_ccp.ccp_compute_station_params(rf_stream, startpoint, endpoint, width, bm=None)[source]
Determines which stations are between startpoint and endpoint great circle profile line,

and within width distance of that profile line. Generates a dictionary of distance along profile line (sta_offset) and orthogonal distance from profile line (dist). For stations outside the region of interest, dictionary has the value None.

Parameters
  • rf_stream (rf.rfstream.RFStream) – RFStreams whose stations will be processed.

  • startpoint (tuple(float, float)) – (latitude, longitude) in degrees of the profile line start point.

  • endpoint (tuple(float, float)) – (latitude, longitude) in degrees of the profile line end point.

  • width (float) – Width of ROI in km to either side of the profile line.

  • bm (mpl_toolkits.basemap.Basemap) – Optional basemap upon which to plot stations coded according to whether they are within the ROI.

Returns

Dictionary keyed by station code containing distances relative to profile line, or None if station is outside the ROI.

Return type

dict

seismic.receiver_fn.plot_ccp.ccp_generate(rf_stream, startpoint, endpoint, width, spacing, max_depth, channels=None, v_background='ak135', station_map=True)[source]
Main function for processing RF collection and plotting common conversion point (CCP) stack of receiver

functions (RFs) along a specific line between startpoint and endpoint.

Parameters
  • rf_stream (rf.rfstream.RFStream) – Sequence of receiver function traces

  • startpoint (tuple(float, float)) – Starting point for the transect line in latitude, longitude degrees

  • endpoint (tuple(float, float)) – End point for the transect line in latitude, longitude degrees

  • width (float) – Width (km) of region about the transect line from which to project RFs to the slice.

  • spacing (float) – Size (km) of each discrete sample cell in the spatial slice beneath the transect.

  • max_depth (float) – Maximum depth (km) of the slice beneath the transect

  • channels (list(str), optional) – Filter for channels, defaults to None in which case only ‘R’ channel is used. Allowed values are in [‘Z’, ‘R’, ‘T’].

  • v_background (str, optional) – Assumed background 1D velocity model, defaults to ‘ak135’

  • station_map (bool) – Whether to generate figure for station map, defaults to True

Returns

Normalized stack matrix, normalization factors, profile length, station metadata, map figure

Return type

numpy.array, numpy.array, float, dict, matplotlib.pyplot.Figure

seismic.receiver_fn.plot_ccp.cross_along_track_distance(p1, p2, p3)[source]

Compute the cross-track distance and the along-track distance of p3 in relation to great circle from p1 to p2.

Math reference: https://www.movable-type.co.uk/scripts/latlong.html

Parameters
Returns

Cross-track distance (orthogonal to arc from p1 to p2) and along-track distance (along arc from p1 to p2) of the location p3.

Return type

tuple(float, float)

seismic.receiver_fn.plot_ccp.equirectangular_projection(x0, y0, x1, y1)[source]

Perform equirectangular projection of a pair of latitude, longitude coordinates to cartesian coordinates.

This length calculation uses the forward equirectangular projection (https://en.wikipedia.org/wiki/Equirectangular_projection). (See also https://www.movable-type.co.uk/scripts/latlong.html)

Parameters
  • x0 (float) – Point 0 longitude (deg)

  • y0 (float) – Point 0 latitude (deg)

  • x1 (float) – Point 1 longitude (deg)

  • y1 (float) – Point 1 latitude (deg)

Returns

Lengths of sides of rectangle and the diagonal. The diagonal is the distance between points 0 and 1.

Return type

float, float, float

seismic.receiver_fn.plot_ccp.get_amplitude(trace, time)[source]

retrieve amplitude value

seismic.receiver_fn.plot_ccp.matrx_lookup(xsz, sta_offset, h, depstep, lenstep)[source]

return index values for amplitude contrbution in profile matrix

seismic.receiver_fn.plot_ccp.plot_ccp(matrx, length, max_depth, spacing, vlims=None, metadata=None, title=None, colormap='seismic')[source]

Plot results of CCP stacking procedure.

Parameters
  • matrx (numpy.array) – [description]

  • length (float) – Length (km) of the transect line

  • max_depth (float) – Maximum depth (km) of the slice to plot

  • spacing (float) – Grid spacing (km)

  • vlims (tuple(float, float), optional) – Min and max values for the color scale, defaults to None

  • metadata (dict, optional) – Metadata generated by ccp_compute_station_params(), defaults to None

  • title (str, optional) – Title text to add to the plot, defaults to None

  • colormap (str) – Color map to use for the CCP intensity shading.

Returns

Figure handle

Return type

matplotlib.pyplot.Figure

seismic.receiver_fn.plot_ccp.run(rf_stream, start_latlon, end_latlon, width, spacing, max_depth, channels, background_model='ak135', stacked_scale=None, title=None, colormap=None)[source]

Run CCP generation on a given dataset of RFs.

Parameters
  • rf_stream (rf.RFStream) – Set of RFs to use for CCP plot

  • start_latlon (tuple(float, float)) – Starting (latitude, longitude) coordinates of line transect

  • end_latlon (tuple(float, float)) – End (latitude, longitude) coordinates of line transect

  • width (float) – Width of transect (km)

  • spacing (float) – Discretization size (km) for RF ray sampling

  • max_depth (float) – Maximum depth of slice below the transect line (km)

  • channels (str, comma separated) – String of comma-separated component IDs to source for the RF amplitude

  • background_model (str) – 1D background model to assume

  • stacked_scale (float) – Max value to represent on color scale of CCP plot

  • title (str) – Title to place at top of CCP plot

  • colormap (str) – Color map to use for the CCP intensity shading. Suggest ‘seismic’, ‘coolwarm’ or ‘jet’.

Returns

Figure handles: Main figure, map figure, dict of station parameters

Return type

(matplotlib.pyplot.Figure, matplotlib.pyplot.Figure, dict)

seismic.receiver_fn.plot_ccp.setup_ccp_profile(length, spacing, maxdep)[source]

Construct the grid for a CCP stacking profile

Parameters
  • length (float) – Length (km) of the profile

  • spacing (float) – Grid spacing (km)

  • maxdep (float) – Maximum depth (km) of the slice to plot

Returns

Zeroed matrix and mesh coordinates

Return type

numpy.array, numpy.array, numpy.array

seismic.receiver_fn.plot_ccp_batch module

seismic.receiver_fn.plot_ccp_batch.gravity_subplot(hf, metadata, grav_map)[source]

Custom plot modifier to add a gravity subplot along the CCP transect line.

Parameters
  • hf (matplotlib.pyplot.Figure) – Figure containing main CCP plot

  • metadata (dict) – Dictionary from CCP plotting code containing metadata of the transect and station data

  • grav_map (Instance of 2D interpolation class from scipy.interpolator) – Interpolator callable that takes a iterable of 2D coordinates and interpolates a gravity value at each, based on external data source.

seismic.receiver_fn.plot_ccp_batch.moho_annotator(hf, metadata)[source]

Custom plot annotator to add markers to the main figure showing locations of Moho estimates from other sources.

Parameters
  • hf (matplotlib.pyplot.Figure) – Figure containing main CCP plot

  • metadata (dict) – Dictionary from CCP plotting code containing metadata of the transect and station data

seismic.receiver_fn.plot_ccp_batch.run_batch(transect_file, rf_waveform_file, fed_db_file, amplitude_filter=False, similarity_filter=False, stack_scale=0.4, width=30.0, spacing=2.0, max_depth=200.0, channel='R', output_folder='', colormap='seismic', annotators=None)[source]

Run CCP generation in batch mode along a series of transects.

Parameters
  • transect_file (str or Path) – File containing specification of network and station locations of ends of transects

  • rf_waveform_file (str or Path) – HDF5 file of QA’d receiver functions for the network matching the transect file

  • fed_db_file (str or Path) – Name of file with which to initialize FederatedASDFDataBase

  • amplitude_filter (bool) – Whether to use amplitude-based filtering of waveforms beform plotting.

  • similarity_filter (bool) – Whether to use RF waveform similarity filtering of waveforms beform plotting.

  • stack_scale (float) – Max value to represent on color scale of CCP plot

  • width (float) – Width of transect (km)

  • spacing (float) – Discretization size (km) for RF ray sampling

  • max_depth (float) – Maximum depth of slice below the transect line (km)

  • channel (str length 1) – Channel component ID to source for the RF amplitude

Returns

None

seismic.receiver_fn.plot_spatial_map module

Use cartopy to plot point dataset onto map.

seismic.receiver_fn.plot_spatial_map.plot_spatial_map(point_dataset, projection_code, title=None, feature_label=None)[source]

Make spatial plot of point dataset with filled contours overlaid on map.

Parameters
  • point_dataset – Name of point dataset file. Should be in format produced by script pointsets2grid.py

  • projection_code – EPSG projection code, e.g. 3577 for Australia

  • title – Title string for top of plot

  • feature_label – Label for the color bar of the plotted feature

Returns

seismic.receiver_fn.pointsets2grid module

Blend multiple datasets on different lon/lat point sets together onto a common grid.

Uses weighted Gaussian interpolation method of Kennett, but simplified to ignore the individual point weighting.

Multiple datasets with per-dataset settings are passed in using a JSON configuration file with layout illustrated by the following example:

{
    "source_files":
    {
        "file1":
        {
            "weighting": 2.5,
            "scale_length_km": 10.0,
            "scale_length_cutoff": 3.5
        },
        "file2":
        {
            "weighting": 0.5,
            "scale_length_km": 20.0
        },
        "file3":
        {
            "weighting": 1.0,
            "scale_length_km": 10.0,
            "scale_length_cutoff": 3
        }
    },
    "proj_crs_code": 3577,
    "output_spacing_km": 10.0
}

Reference: B. L. N. Kennett 2019, “Areal parameter estimates from multiple datasets”, Proc. R. Soc. A. 475:20190352, http://dx.doi.org/10.1098/rspa.2019.0352

Requires:

  • pyepsg

  • cartopy

seismic.receiver_fn.rf_3dmigrate module

Description:
Implements migration algorithm as described in Frassetto et al. (2010):

Improved imaging with phase-weighted common conversion point stacks of receiver functions (GJI)

References:

CreationDate: 3/15/18 Developer: rakib.hassan@ga.gov.au

Revision History:

LastUpdate: 3/15/18 RH LastUpdate: dd/mm/yyyy Who Optional description

class seismic.receiver_fn.rf_3dmigrate.Geometry(start_lat_lon, azimuth, lengthkm, nx, widthkm, ny, depthkm, nz, debug=False)[source]

Bases: object

generateGrids()[source]
class seismic.receiver_fn.rf_3dmigrate.Migrate(geometry, stream, debug=False, output_folder='/tmp')[source]

Bases: object

execute()[source]
seismic.receiver_fn.rf_3dmigrate.rtp2xyz(r, theta, phi)[source]

Convert spherical to cartesian coordinates

Parameters
  • r ([type]) – [description]

  • theta ([type]) – [description]

  • phi ([type]) – [description]

Returns

[description]

Return type

[type]

seismic.receiver_fn.rf_3dmigrate.xyz2rtp(x, y, z)[source]

Convert cartesian to spherical coordinates

Parameters
  • x ([type]) – [description]

  • y ([type]) – [description]

  • z ([type]) – [description]

Returns

[description]

Return type

[type]

seismic.receiver_fn.rf_deconvolution module

Custom receiver function deconvolution methods and support functions.

seismic.receiver_fn.rf_deconvolution.iter_deconv_pulsetrain(numerator, denominator, sampling_rate, time_shift, max_pulses=1000, tol=0.001, gwidth=2.5, only_positive=False, log=None)[source]

Iterative deconvolution of source and response signal to generate seismic receiver function. Adapted to Python by Andrew Medlin, Geoscience Australia (2019), from Chuck Ammon’s (Saint Louis University) iterdeconfd Fortran code, version 1.04.

Note this is not really a frequency-domain deconvolution method, since the deconvolution is based on generating a time-domain pulse train which is filtered and convolved with source signal in time domain in order to try to replicate observation. Results should be the same even if some of the spectral techniques used in functions (such as gauss_filter) were replaced by non-spectral equivalents.

Parameters
  • numerator (numpy.array(float)) – The observed response signal (e.g. R, Q or T component)

  • denominator (numpy.array(float)) – The source signal (e.g. L or Z component)

  • sampling_rate (float) – The sampling rate in Hz of the numerator and denominator signals

  • time_shift (float) – Time shift (sec) from start of input signals until expected P-wave arrival onset.

  • max_pulses (int) – Maximum number of delta function pulses to synthesize in the unfiltered RF (up to 1000)

  • tol (float) – Convergence tolerance, iteration stops if change in error falls below this value

  • gwidth (float) – Gaussian filter width in the normalized frequency domain.

  • only_positive (bool) – If true, only use positive pulses in the RF.

  • log (logging.Logger) – Log instance to log messages to.

Returns

RF trace, pulse train, expected response signal, predicted response signal, quality of fit statistic

Return type

numpy.array(float), numpy.array(float), numpy.array(float), numpy.array(float), float

seismic.receiver_fn.rf_deconvolution.rf_iter_deconv(response_data, source_data, sr, tshift, min_fit_threshold=80.0, normalize=0, **kwargs)[source]

Adapter function to rf library. To use, add arguments deconvolve=’func’, func=rf_iter_deconv to rf.RFStream.rf() function call.

Parameters
  • response_data (list of numpy.array(float)) – List of response signals for which to compute receiver function

  • source_data (numpy.array(float)) – Source signal to use for computing the receiver functions

  • sampling_rate (float) – Sampling rate of the input signals (Hz)

  • time_shift (float) – Time shift (seconds) from start of signal to onset

  • min_fit_threshold (float) – Minimum percentage of fit to include trace in results, otherwise will be returned as empty numpy array.

  • normalize (int or None) – Component of stream to use for normalization, usually component 0 (the vertical component). Set to None to disable normalization.

Returns

Receiver functions corresponding to the list of input response signals.

Return type

list of numpy.array(float)

seismic.receiver_fn.rf_h5_file_event_iterator module

Helper class to iterate over 3-channel event traces in h5 file generated by rf library, without loading all the traces into memory. This is a scalable solution for very large files.

class seismic.receiver_fn.rf_h5_file_event_iterator.IterRfH5FileEvents(h5_filename, memmap=False, channel_pattern=None)[source]

Bases: object

Helper class to iterate over events in h5 file generated by extract_event_traces.py and pass them to RF generator. This class avoids having to load the whole file up front via obspy which is slow and not scalable.

Yields (station_id, event_id, event_time, stream)

Due to the expected hierarchy structure of the input H5 file, yielded event traces are grouped by station ID.

Data yielded per event can easily be hundreds of kB in size, depending on the length of the event traces. rf library defaults to window of (-50, 150) sec about the P wave arrival time.

Initializer

Parameters
  • h5_filename (str or pathlib.Path) – Name of file containing event seismograms in HDF5 format, indexed in seismic.stream_io.EVENTIO_H5INDEX format

  • memmap (bool) – If True, memmap the open file. Can improve tractability of handling very large files.

  • channel_pattern (str) – [OPTIONAL] Ordered list of preferred channels, e.g. ‘HH*,BH*’. Channel mask to use to select channels returned.

seismic.receiver_fn.rf_h5_file_station_iterator module

Helper class to iterate over station events in h5 file generated by rf library, without loading all the traces into memory. This is a scalable solution for very large files.

class seismic.receiver_fn.rf_h5_file_station_iterator.IterRfH5StationEvents(h5_filename, memmap=False)[source]

Bases: object

Helper class to iterate over stations in h5 file generated by extract_event_traces.py and pass them to RF generator. This class avoids having to load the whole file up front via obspy which is slow and not scalable.

This class yields 3-channel traces per station per event.

Data yielded per station can easily be many MB in size.

Initializer

Parameters
  • h5_filename (str or pathlib.Path) – Name of file containing event seismograms in HDF5 format, indexed in seismic.stream_io.EVENTIO_H5INDEX format

  • memmap (bool) – If True, memmap the open file. Can improve tractability of handling very large files.

seismic.receiver_fn.rf_handpick_tool module

Simple tool for hand picking functions from RF plots. Selected plots event IDs are exported to text file for later use as RF filter in other tools and workflows.

seismic.receiver_fn.rf_handpick_tool.main()[source]

Main entry function for RF picking tool.

seismic.receiver_fn.rf_handpick_tool.on_release(event, target_axes, select_mask, background, rect_selector)[source]

Event handler when mouse button released.

Parameters
  • target_axes (matplotlib.axes.Axes) – The original axes in which the RectangleSelector began. May differ from event.inaxes.

  • event (matplotlib.backend_bases.MouseEvent) – Button up event for end of rectangle area selection

  • select_mask (numpy.array(bool)) – Boolean mask of selected state of each receiver function in the plot

  • background (Return type from fig.canvas.copy_from_bbox) – Background raster from initial render of RFs

  • rect_selector (matplotlib.widgets.RectangleSelector) – Widget for selecting rectangular region to toggle RF selection

seismic.receiver_fn.rf_handpick_tool.on_select(e_down, e_up, select_mask)[source]

Event handler for RectangleSelector

Parameters

seismic.receiver_fn.rf_inversion_export module

Export RFs to file format for external inversion code to run on.

seismic.receiver_fn.rf_inversion_export.rf_inversion_export(input_h5_file, output_folder, network_code, component='R', resample_freq=6.25, trim_window=- 5.0, 20.0, moveout=True)[source]

Export receiver function to text format for ingestion into Fortran RF inversion code.

Parameters
  • input_h5_file (str or Path) – Input hdf5 file containing receiver function data

  • output_folder (str or Path) – Folder in which to export text files, one per channel per station. Will be appended with network code.

  • network_code (str) – Network to which this RF data belongs, used to disambiguate and track folders.

  • component (str, optional) – The channel component to export, defaults to ‘R’

  • resample_freq (float, optional) – Sampling rate (Hz) of the output files, defaults to 6.25 Hz

  • trim_window (tuple, optional) – Time window to export relative to onset, defaults to (-5.0, 20.0). If data needs to be resampled, the samples are anchored to the start of this time window.

  • moveout (bool, optional) – Whether to apply moveout correction prior to exporting, defaults to True

seismic.receiver_fn.rf_network_dict module

Class encapsulating a collection of receiver functions for stations of one network.

class seismic.receiver_fn.rf_network_dict.NetworkRFDict(rf_stream)[source]

Bases: object

Collection of RFs for a given network indexed by station code, channel code.

Initialize from rf.RFStream

Parameters

rf_stream (rf.RFStream) – RFStream data

keys()[source]

Accessor for the top level keys (station codes) of the network in an iterable container.

Returns

Iterable of top level keys to the dictionary

Return type

Python iterable

seismic.receiver_fn.rf_plot_utils module

Utility plotting functions for consistent and convenient plotting of RFs.

seismic.receiver_fn.rf_plot_utils.plot_hk_stack(k_grid, h_grid, hk_stack, title=None, save_file=None, num=None, clip_negative=True)[source]

Plot H-k stack using data generated by function seismic.receiver_fn.rf_stacking.computed_weighted_stack().

Parameters
  • k_grid (Two-dimensional numpy.array) – Grid of k-values

  • h_grid (Two-dimensional numpy.array) – Grid of H-values

  • hk_stack (Two-dimensional numpy.array) – Grid of stacked RF sample values produced by function seismic.receiver_fn.rf_stacking.computed_weighted_stack()

  • title (str, optional) – Title to add to the plot, defaults to None

  • save_file (str or Path, optional) – File name in which to save the plot, defaults to None

  • num (int, optional) – Number of RFs used to produce the stack, defaults to None

  • clip_negative (bool, optional) – Clip negative stack regions to zero, defaults to True

Returns

Handle to the figure created for the plot

Return type

matplotlib.figure.Figure

seismic.receiver_fn.rf_plot_utils.plot_iir_filter_response(filter_band_hz, sampling_rate_hz, corners)[source]

Plot one-way bandpass filter response in the frequency domain. If filter is used as zero-phase, the attenuation will be twice what is computed here.

Parameters
  • filter_band_hz (tuple(float) of length 2 (i.e. pair)) – Pair of frequencies corresponding to low cutoff and high cutoff freqs

  • sampling_rate_hz (float) – The sampling rate in Hz

  • corners (int) – The order of the filter

Returns

Figure object

Return type

matplotlib.figure.Figure

seismic.receiver_fn.rf_plot_utils.plot_iir_impulse_response(filter_band_hz, sampling_rate_hz, corners, zero_phase=False, N=1000, blip_period=1.0)[source]

Plot bandpass filter response to standard waveforms in the time domain - impulse (delta function, step function, square wave pulse. By default filter is applied one-way. Set zero_phase=True to plot two-way filter response.

Parameters
  • filter_band_hz (tuple(float) of length 2 (i.e. pair)) – Pair of frequencies corresponding to low cutoff and high cutoff freqs

  • sampling_rate_hz (float) – The sampling rate in Hz

  • zero_phase (bool) – If True, plot two-way signal response (zero phase), otherwise plot one-way signal response.

  • N (int) – Number of samples in the input test signals.

  • blip_period (float) – Period of the ‘blip’ test signal (square wave pulse).

Returns

Figure object

Return type

matplotlib.figure.Figure

seismic.receiver_fn.rf_plot_utils.plot_rf_stack(rf_stream, time_window=- 10.0, 25.0, trace_height=0.2, stack_height=0.8, save_file=None, **kwargs)[source]

Wrapper function of rf.RFStream.plot_rf() to help do RF plotting with consistent formatting and layout.

Parameters
  • rf_stream (rf.RFStream) – RFStream to plot

  • time_window (tuple, optional) – Time window to plot, defaults to (-10.0, 25.0)

  • trace_height (float, optional) – Height of a single trace (reduce to cram RFs closer together), defaults to 0.2

  • stack_height (float, optional) – Height of mean (stacked) RF at top of plot, defaults to 0.8

  • save_file (str to valid file path, optional) – File to save resulting image into, defaults to None

Returns

Figure handle to the stack plot

Return type

matplotlib.figure.Figure

seismic.receiver_fn.rf_plot_utils.plot_rf_wheel(rf_stream, max_time=15.0, deg_per_unit_amplitude=45.0, plt_col='C0', title='', figsize=10, 10, cluster=True, cluster_col='#ff4000', layout=None, fontscaling=1.0)[source]

Plot receiver functions around a polar plot with source direction used to position radial RF plot.

Parameters
  • rf_stream (rf.RFStream or list(rf.RFStream)) – Collection of RFs to plot. If passed as a list, then each stream in the list will be plotted on separate polar axes.

  • max_time (float, optional) – maximum time relative to onset, defaults to 25.0

  • deg_per_unit_amplitude (float, optional) – Azimuthal scaling factor for RF amplitude, defaults to 20

  • plt_col (str, optional) – Plot color for line and positive signal areas, defaults to ‘C0’

  • title (str, optional) – Title for the overall plot, defaults to ‘’

  • figsize (tuple, optional) – Size of figure area, defaults to (12, 12)

  • cluster (bool) – Whether to add overlaid mean RF where there are many RFs close together.

  • cluster_col (matplotlib color specification) – Color of clustered stacked overlay plots.

  • layout (tuple(int, int)) – Arrangement of polar plots in grid. If None, then arranged in a column.

Returns

Figure object

Return type

matplotlib.figure.Figure

seismic.receiver_fn.rf_plot_utils.plot_station_rf_overlays(db_station, title=None, time_range=None)[source]

Plot translucent overlaid RF traces for all traces in each channel, and overplot the mean signal of all the traces per channel.

Parameters
  • db_station (dict({str, list(RFTrace)})) – Dictionary with list of traces per channel for a given station.

  • title (str, optional) – Plot title, defaults to None

  • time_range (Pair of float) – Min and max times for the horizontal axis

Returns

Mean trace signal per channel

Return type

list(numpy.array)

seismic.receiver_fn.rf_plot_vespagram module

Plot vespagrams from receiver function data

seismic.receiver_fn.rf_process_io module

Helper for asynchronous writing of RF data to file.

seismic.receiver_fn.rf_process_io.async_write(rfstream_queue, outfile_name, max_buffered=100, metadata='')[source]

Monitors asynchronous queue for data, removes from queue to buffer, then flushes buffer intermittently and when queue termination signal is put.

When None is received on the queue, this is taken as the signal to terminate monitoring the queue.

Parameters
  • rfstream_queue (multiprocessing.Manager.Queue) – Queue into which RFStreams are pushed for writing to file.

  • outfile_name (str or Path) – Name of file into which queued RFStream results are periodically written.

  • max_buffered (int, optional) – Maximum number of RFStreams to buffer before flushing to file, defaults to 100

  • metadata (str, optional) – Metadata string to write to root attribute

seismic.receiver_fn.rf_quality_filter module

Filter out invalid RFs, and compute a variety of quality metrics for the remaining RFs. These metrics can be combined in various ways downstream to perform different kinds of filtering. Quality metrics are stored in the stats of each trace.

seismic.receiver_fn.rf_quality_filter.compute_max_coherence(orig_stream, f1, f2)[source]

Finding coherence between two signals in frequency domain f1 and f2 - normalised min and max frequencies, f1 < f2 <= ~0.5 returns array of indexes for coherent traces with median

Suggested minimum level of coherence for good results: 0.6

seismic.receiver_fn.rf_quality_filter.compute_rf_quality_metrics(station_id, station_stream3c, similarity_eps)[source]

Top level function for adding quality metrics to trace metadata.

Parameters
  • station_id (str) – Station ID

  • station_stream3c (list(rf.RFStream) with 3 components) – 3-channel stream

  • similarity_eps (float) – Distance threshold used for DBSCAN clustering

Returns

Triplet of RF streams with Z, R or Q, and T components with populated quality metrics. Otherwise return None in case of failure.

seismic.receiver_fn.rf_quality_filter.get_rf_stream_components(stream)[source]

Identify the RF component types and return them.

Parameters

stream (rf.RFStream) – Stream containing mixed RF components.

Returns

(RF component type, primary RF component (R or Q), transverse RF component (T), source component (Z or L))

Return type

(str, rf.RFStream, rf.RFStream, rf.RFStream)

seismic.receiver_fn.rf_quality_filter.rf_group_by_similarity(swipe, similarity_eps)[source]

Cluster waveforms by similarity

Parameters
  • swipe (numpy.array) – Numpy array of RF rowwise

  • similarity_eps (float) – Tolerance on similarity between traced to be considered in the same group.

Returns

Index of the group for each trace. -1 if no group is found for a given trace.

Return type

numpy.array

seismic.receiver_fn.rf_quality_filter.rf_quality_metrics_queue(oqueue, station_id, station_stream3c, similarity_eps, drop_z=True)[source]

Produce RF quality metrics in a stream and queue the QC’d components for downstream processing.

Parameters
  • oqueue (queue or multiprocessing.Manager.Queue) – Output queue where filtered streams are queued

  • station_id (str) – Station ID

  • station_stream3c (list(rf.RFStream) with 3 components) – 3-channel stream

  • similarity_eps (float) – Distance threshold used for DBSCAN clustering

seismic.receiver_fn.rf_quality_filter.spectral_entropy(stream)[source]

Compute the spectral entropy of a trace

Parameters

trace (rf.RFTrace) – Single channel seismic trace

Returns

Spectral entropy of the trace waveform

Return type

float

seismic.receiver_fn.rf_stacking module

seismic.receiver_fn.rf_stacking.compute_hk_stack(cha_data, V_p=None, h_range=None, k_range=None, root_order=1, include_t3=True)[source]

Compute H-k stacking array on a dataset of receiver functions.

Parameters
  • cha_data (Iterable(rf.RFTrace)) – List or iterable of RF traces to use for H-k stacking.

  • V_p (float, optional) – P-wave velocity in crustal layer, defaults to None in which case it is inferred from trace metadata

  • h_range (numpy.array [1D], optional) – Range of h values (Moho depth) values to cover, defaults to np.linspace(20.0, 70.0, 251)

  • k_range (numpy.array [1D], optional) – Range of k values to cover, defaults to np.linspace(1.4, 2.0, 301)

  • root_order (int, optional) – Exponent for nth root stacking as per K.J.Muirhead (1968), defaults to 1

  • include_t3 (bool, optional) – If True, include the t3 (PpSs+PsPs) multiple in the stacking, defaults to True

Returns

k-grid values [2D], h-grid values [2D], H-k stack in series of 2D layers having one layer per multiple

Return type

numpy.array [2D], numpy.array [2D], numpy.array [3D]

seismic.receiver_fn.rf_stacking.compute_theoretical_phase_times(tr, H, k, V_p, include_t3=True)[source]

Compute arrival times of Ps, PpPs and (PpSs + PsPs) phases relative to primary P-wave arrival time for an assume Moho depth H, velocity ratio k = V<sub>p</sub>/V<sub>s</sub>, and p-wave velocity V<sub>p</sub>.

Parameters
  • tr (obspy.Trace or rf.RFTrace) – Trace for which to compute the theoretical phase arrival times

  • H (float) – Presumed Moho depth (km)

  • k (float) – Presumed velocity ratio k (dimensionless)

  • V_p (float) – Presumed p-wave velocity (km/sec)

  • include_t3 (bool) – Flag whether or not to compute t<sub>3</sub> value (i.e. (PpSs + PsPs) arrival time)

Returns

Triplet of phase arrival times (t1, t2, t3). t3 is None if include_t3 is False.

Return type

tuple(float, float, float)

seismic.receiver_fn.rf_stacking.compute_weighted_stack(hk_components, weighting=0.5, 0.5, 0.0)[source]

Given stack components from function compute_hk_stack, compute the overall weighted stack.

Parameters
  • hk_components (numpy.array) – H-k stack layers returned from compute_hk_stack

  • weighting (tuple, optional) – Weightings for (t1, t2, t3) layers respectively, defaults to (0.5, 0.5, 0.0)

Returns

Weighted stack in H-k space

Return type

numpy.array

seismic.receiver_fn.rf_stacking.find_global_hk_maximum(k_grid, h_grid, hk_weighted_stack)[source]

Given the weighted stack computed from function compute_weighted_stack and the corresponding k-grid and h-grid, find the location in H-k space of the global maximum.

Parameters
  • k_grid (Two-dimensional numpy.array) – Grid of k-values

  • h_grid (Two-dimensional numpy.array) – Grid of H-values

  • hk_weighted_stack (Two-dimensional numpy.array) – Grid of stacked RF sample values produced by function rf_stacking.computed_weighted_stack()

Returns

Location of global maximum on the H-k grid of the maximum stack value.

Return type

tuple(float, float)

seismic.receiver_fn.rf_stacking.find_local_hk_maxima(k_grid, h_grid, hk_stack, min_rel_value=0.5)[source]

Given the weighted stack computed from function compute_weighted_stack and the corresponding k-grid and h-grid, find the locations in H-k space of all local maxima above a certain threshold.

Parameters
  • k_grid (Two-dimensional numpy.array) – Grid of k-values

  • h_grid (Two-dimensional numpy.array) – Grid of H-values

  • hk_stack (Two-dimensional numpy.array) – Grid of stacked RF sample values produced by function rf_stacking.computed_weighted_stack()

  • min_rel_value (float, optional) – Minimum value required relative to the largest value in the stack, defaults to 0.5

Returns

List of tuples containing parameters of local maxima solutions, with values in the following order: (H, k, stack_value, row_index, col_index)

Return type

list(tuple(float, float, float, int, int))

seismic.receiver_fn.rf_stacking.infer_Vp_from_traces(cha_data, log=None)[source]

Infer the Vp value used in earth model for computing trace stats.

Parameters
  • cha_data (Iterable(obspy.Stream) or Iterable(rf.RFStream)) – Iterable of traces for a given event (e.g. obspy.Stream or rf.RFStream)

  • log (logging.Logger) – Logging instance to log messages

Returns

Vp value

Return type

float

seismic.receiver_fn.rf_synthetic module

Helper functions for producing synthetic pseudo-Receiver function traces

seismic.receiver_fn.rf_synthetic.convert_inclination_to_distance(inclinations, model='iasp91', nominal_source_depth_km=10.0)[source]

Helper function to convert range of inclinations to teleseismic distance in degrees.

Parameters
  • inclinations (numpy.array(float)) – Array of inclination angles in degrees

  • model (str, optional) – Name of model to use for ray tracing, defaults to “iasp91”

  • nominal_source_depth_km (float, optional) – Assumed depth of source events, defaults to 10.0

Returns

Array of teleseismic distances in degrees corresponding to input inclination angles.

Return type

numpy.array(float)

seismic.receiver_fn.rf_synthetic.generate_synth_rf(arrival_times, arrival_amplitudes, fs_hz=100.0, window_sec=- 10, 30, f_cutoff_hz=2.0)[source]

Simple generator of synthetic R component receiver function with pulses at given arrival times.

Parameters
  • arrival_times (iterable of float) – Iterable of arrival times as numerical values in seconds

  • arrival_amplitudes (iterable of float) – Iterable of arrival amplitudes

  • fs_hz (float, optional) – Sampling rate (Hz) of output signal, defaults to 100.0

  • window_sec (tuple, optional) – Time window over which to create signal (sec), defaults to (-10, 30)

  • f_cutoff_hz (float, optional) – Cutoff frequency (Hz) for low-pass filtering to generate realistic result, defaults to 2.0

Returns

Array of times and corresponding signal amplitudes

Return type

numpy.array, numpy.array

seismic.receiver_fn.rf_synthetic.synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, include_t3=False, amplitudes=None, baz=0.0)[source]

Synthesize RF R-component data set over range of inclinations and distances and get result as a rf.RFStream instance.

Parameters
  • H (float) – Moho depth (km)

  • V_p (float) – P body wave velocity in uppermost layer

  • V_s (float) – S body wave velocity in uppermost layer

  • inclinations (numpy.array(float)) – Array of inclinations for which to create RFs

  • distances (numpy.array(float)) – Array of teleseismic distances corresponding to inclinations

  • ds (float) – Final sampling rate (Hz) for the downsampled output signal

  • log (logger, optional) – Logger to send output to, defaults to None

  • include_t3 (bool, optional) – If True, include the third expected multiple PpSs+PsPs

  • amplitudes (list(float), optional) – Custom amplitudes to apply to the multiples

  • baz (float, optional) – Back azimuth for metadata

Returns

Stream containing synthetic RFs

Return type

rf.RFStream

seismic.receiver_fn.rf_util module

Utility functions to help with RF processing and analysis.

seismic.receiver_fn.rf_util.choose_rf_source_channel(rf_type, db_station)[source]

Choose source channel for RF analysis.

Parameters
  • rf_type (str) – The RF rotation type, should be either ‘ZRT’ or ‘LQT’

  • db_station (dict(str, list(rf.RFTrace))) – Dict of traces for a given station keyed by channel code.

Returns

Channel code of the primary RF source channel

Return type

str

seismic.receiver_fn.rf_util.compute_extra_rf_stats(stream)[source]

Compute extra statistics for each trace and add it to the RFTrace.stats structure.

Parameters

stream (rf.RFStream) – RFStream to augment with additional metadata statistics.

seismic.receiver_fn.rf_util.compute_rf_snr(rf_stream)[source]

Compute signal to noise (S/N) ratio of the RF itself about the onset pulse (key ‘snr’). This SNR is a ratio of RMS amplitudes. Stores results in metadata of input stream traces.

In the LQT rotation case when rotation is working ideally, the onset pulse of the rotated transverse signals should be minimal, and a large pulse at t = 0 indicates lack of effective rotation of coordinate system, so for ‘snr’ we use a long time window after onset pulse, deliberately excluding the onset pulse, to maximize contribution to the SNR from the multiples after the onset pulse.

Parameters

rf_stream (rf.RFStream) – R or Q component of Receiver Function

Returns

SNR for each trace in the input stream

Return type

numpy.array

seismic.receiver_fn.rf_util.compute_vertical_snr(src_stream)[source]

Compute the SNR of the Z component (Z before deconvolution) including the onset pulse (key ‘snr_prior’). Stores results in metadata of input stream traces. This SNR is a ratio of max envelopes.

Some authors compute this prior SNR on signal after rotation but before deconvolution, however that doesn’t make sense for LQT rotation where the optimal rotation will result in the least energy in the L component. For simplicity we compute it on Z-component only which is a reasonable estimate for teleseismic events.

Parameters

src_stream (rf.RFStream or obspy.Stream) – Seismic traces before RF deconvolution of raw stream.

seismic.receiver_fn.rf_util.filter_crosscorr_coeff(rf_stream, time_window=- 2, 25, threshold_cc=0.7, min_fraction=0.15, apply_moveout=False)[source]

For each trace in the stream, compute its correlation coefficient with the other traces. Return only traces matching cross correlation coefficient criteria based on C.Sippl (2016) [see http://dx.doi.org/10.1016/j.tecto.2016.03.031]

Parameters
  • rf_stream (rf.RFStream) – Stream of RF traces to filter, should be for a single component of a single station

  • time_window (tuple, optional) – Time window to filter by, defaults to (-2, 25)

  • threshold_cc (float, optional) – Threshold cross-correlation coefficient, defaults to 0.70. Denoted Xi in Sippl, who used value 0.80.

  • min_fraction (float, optional) – Minimum fraction of coefficients above threshold_cc, defaults to 0.15. Denoted tau in Sippl, who used value 0.15.

  • apply_moveout (bool) – Whether to apply moveout correction to Ps phase prior to computing correlation coefficients.

Returns

Filtered stream of RF traces

Return type

rf.RFStream

seismic.receiver_fn.rf_util.filter_station_streams(db_station, freq_band=None, None)[source]

Perform frequency filtering on all channels’ traces for a given station. Returns a copy of db_station with streams containing filtered results.

seismic.receiver_fn.rf_util.filter_station_to_mean_signal(db_station, min_correlation=1.0)[source]

Filter out streams which are not ‘close enough’ to the mean signal, based on simple correlation score. The term “correlation” here really just means a similarity dot product (projection of individual trace onto the mean).

seismic.receiver_fn.rf_util.find_rf_group_ids(stream)[source]

For the given stream, which is expected to have an rf_group attribute in its traces’ metadata, determine the unique set of group ids that the traces contain.

Parameters

stream (obspy.Stream) – Stream containing traces with rf_group ids associated with them.

Returns

Set of rf_group ids found in the traces

Return type

set(int)

seismic.receiver_fn.rf_util.label_rf_quality_simple_amplitude(rf_type, traces, snr_cutoff=2.0, rms_amp_cutoff=0.2, max_amp_cutoff=1.0)[source]

Add RF quality label for a collection of RFs based on simple amplitude criteria computed by quality filter script. Adds quality label in-place.

Parameters
  • rf_type (str) – The RF rotation type, should be either ‘ZRT’ or ‘LQT’

  • traces (Iterable collection of rf.RFTrace) – Iterable collection of rf.RFTrace

  • snr_cutoff (float, optional) – Minimum signal SNR, defaults to 2.0

  • rms_amp_cutoff (float, optional) – Maximum accepted RMS amplitude of signal, defaults to 0.2

  • max_amp_cutoff (float, optional) – Maximum accepted amplitude of signal, defaults to 1.0

seismic.receiver_fn.rf_util.phase_weights(stream)[source]

Phase weighting takes all the traces in a stream and computes a weighting for each sample in the stream between 0 and 1. The weighting represents how consistent is the phase angle of the signal at the same point in the time series across all streams.

If phase weighting to accentuate higher multiples than Ps, then moveout should be applied first before calling this function.

See https://doi.org/10.1111/j.1365-246X.1997.tb05664.x

Note: this function should not be applied to streams with mixed components.

Parameters

stream (Iterable container of obspy.Trace) – Stream containing one or more traces from which phase coherence weightings will be generated.

Returns

Array of normalized weighting factors with same length as traces in stream.

Return type

numpy.array

seismic.receiver_fn.rf_util.read_h5_rf(src_file, network=None, station=None, loc='', root='/waveforms')[source]

Helper function to load data from hdf5 file generated by rf library or script rf_quality_filter.py. For faster loading time, a particular network and station may be specified.

Parameters
  • src_file (str or Path) – File from which to load data

  • network (str, optional) – Specific network to load, defaults to None

  • station (str, optional) – Specific station to load, defaults to None

  • root (str, optional) – Root path in hdf5 file where to start looking for data, defaults to ‘/waveforms’

Returns

All the loaded data in a rf.RFStream container.

Return type

rf.RFStream

seismic.receiver_fn.rf_util.rf_to_dict(rf_data)[source]

Convert RF data loaded from function read_h5_rf() into a dict format for easier addressing of selected station and channel RF traces.

Parameters

rf_data (rf.RFStream) – RFStream data

Returns

Nested dicts to find traces by station then channel code, with attached metadata.

Return type

seismic.receiver_fn.rf_network_dict.NetworkRFDict

seismic.receiver_fn.rf_util.signed_nth_power(arr, order)[source]

As per DOI https://doi.org/10.1038/217533a0. Muirhead, K.J. “Eliminating False Alarms when detecting Seismic Events Automatically”

Parameters
  • arr (numpy.array) – Compute n-th power of input array, preserving sign of original data.

  • order (float or int) – Order of the power to compute

Returns

Input array raised to nth power.

Return type

numpy.array

seismic.receiver_fn.rf_util.signed_nth_root(arr, order)[source]

As per DOI https://doi.org/10.1038/217533a0. Muirhead, K.J. “Eliminating False Alarms when detecting Seismic Events Automatically”

Parameters
  • arr (numpy.array) – Compute n-th root of input array, preserving sign of original data.

  • order (float or int) – Order of the root to compute

Returns

Input array raised to 1/nth power.

Return type

numpy.array

Module contents