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
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
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
-
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
-
seismic.receiver_fn.plot_ccp.bounding_box(startpoint, endpoint)[source]¶ Compute a bounding box from start and end points.
-
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
-
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
-
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)
-
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_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
-
class
seismic.receiver_fn.rf_3dmigrate.Migrate(geometry, stream, debug=False, output_folder='/tmp')[source]¶ Bases:
object
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, **kwargs)[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:
objectHelper 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.
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.
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:
objectHelper 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.
Data yielded per station can easily be many MB in size.
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.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
e_down (matplotlib.backend_bases.MouseEvent) – Button down event for start of rectangle
e_up (matplotlib.backend_bases.MouseEvent) – Button up event for end of rectangle
select_mask (numpy.array(bool)) – Boolean mask of selected state of each receiver function in the plot
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.
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 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 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.pyplot.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
- Returns
Figure object
- Return type
-
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).
-
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
-
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
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
- 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_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]¶
-
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 (np.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
-
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
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
- 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.
-
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.core.trace.Trace) – Stream containing traces with rf_group ids associated with them.
- Returns
Set of rf_group ids found in the traces
- Return type
-
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.Traces) – 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
- 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_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”
-
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”