Welcome to rf’s documentation!

rf: Receiver function calculation in seismology

This module heavily depends on ObsPy. The main functionality is provided by the class RFStream which is derived from ObsPy’s Stream class.

Author:Tom Richter
License:MIT
Documentation:http://rf.readthedocs.org/ (installation instructions, usage, examples and module reference)
Project page:https://github.com/trichter/rf (bug reports and feature requests via issue tracker)
Pypi page:https://pypi.python.org/pypi/rf
Test status:buildstatus

Installation

Install ObsPy, its dependencies and pip, e.g. by (Ubuntu)

sudo apt-get install python-obspy python-pip

rf can then be installed by

pip install rf

The tests can be run with the script

rf-runtests

Hdf5 file support can be installed with the obspyh5 package.

Manual installation of dev

Install

Then download the source code from GitHub eg. by

mkdir rf
git clone https://github.com/trichter/rf.git rf

Now install rf with

python setup.py install

Note

At the moment rf-dev is only working with the development version of obspy. Install obspy-dev by e.g. pip install –allow-external obspy –allow-insecure obspy obspy==dev

Basic Usage

The canonical way to load a waveform file into a RFStream is to use the read_rf() function.

>>> from rf import read_rf
>>> stream = read_rf('myfile.SAC')

If you already have an ObsPy Stream and you want to turn it into a RFStream use the generator of RFStream:

>>> from rf import RFStream
>>> stream = RFStream(obspy_stream)

The stream is again written to disc as usual by its write method:

>>> stream.write('outfile', 'SAC')

The RFStream object inherits a lot of useful methods from its ObsPy ancestor (e.g. filter, taper, simulate, ...).

The module automatically maps important (for rf calculation) header information from the stats object attached to every trace to the format specific headers. At the moment only SAC and SH/Q headers are supported. When initializing an RFStream the header information in the format specific headers are written to the stats object and before writing the information stored in the stats object is written back to the format specific headers. In this way the important header information is guaranteed to be saved in the waveform files. The following table reflects the mapping:

stats SH/Q SAC
station_latitude DCVREG stla
station_longitude DCVINCI stlo
station_elevation PWDW stel
event_latitude LAT evla
event_longitude LON evlo
event_depth DEPTH evdp
event_magnitude MAGNITUDE mag
event_time ORIGIN o
onset P-ONSET a
distance DISTANCE gcarc
back_azimuth AZIMUTH baz
inclination INCI user0
slowness SLOWNESS user1

Note

Q-file headers DCVREG DCVINCI PWDW are used for the station information, because the Q format has a shortage of predefined headers.

Note

Alternatively the hdf5 file format can be used. It is supported via the obspyh5 package. In this case all supported stats entries are automatically attached to the stored data.

The first task when calculating receiver functions is calculating some ray specific values like azimuth and epicentral distance. An appropriate stats dictionary can be calculated with rfstats():

>>> from rf import rfstats
>>> stats = rfstats(station=station, event=event, phase='P',dist_range=(30,90))
>>> for tr in stream:
>>>     tr.stats.update(stats)

or if the station and event information is already stored in the stats object:

>>> for tr in stream:
>>>     rfstats(stats=tr.stats)

Now P receiver functions can be calculated by

>>> stream.filter('bandpass', freqmin=0.05, freqmax=1.)
>>> stream.rf()
>>> stream.write('rf', 'Q')

rf can also calculate S receiver functions (not much tested):

>>> stream.rf(method='S')

When calling stream.rf the following operations are performed depending on the given kwargs:

  • filtering
  • trimming data to window relative to onset
  • downsampling
  • rotation
  • deconvolution

Please see RFStream.rf() for a more detailed description. RFStream provides the possibility to perform moveout correction and piercing point calculation.

Batch Usage

rf provides a command line utility ‘rf’ which runs all the necessary steps to perform receiver function calculation. All you need is a StationXML file with your stations as an inventory and a QuakeML file with the events you want to analyze.

The command

rf init test

initializes a new rf project in a new folder with an example configuration file conf.py. You will have to edit this file so that the program finds your data (events, stations and waveforms). After that you can use the various subcommands of rf to perform different tasks (e.g. receiver function calculation, plotting).

To init a project with a small included dataset and working configuration you can use

rf init -t test2

Now switch to the newly created directy and start using rf...

cd test2
rf calc P
rf moveout Prf Ps
rf plot Prf_Ps

Miscellaneous

Please feel free to request features, report bugs or contribute code on GitHub. The code is continiously tested by travis-ci.

Classes & Functions

RFStream Class providing the necessary functions for receiver function calculation.
rfstats Calculate ray specific values like slowness for given event and station.

Modules

rfstream Classes and functions for receiver function calculation.
batch rf batch command line utility
deconvolve Frequency and time domain deconvolution.
simple_model Simple move out and piercing point calculation.

Module Reference

rfstream Module

Classes and functions for receiver function calculation.

class rf.rfstream.RFStream(traces=None, warn=True)[source]

Bases: Stream

Class providing the necessary functions for receiver function calculation.

To initialize a RFStream from a Stream object use

>>> rfstream = RFStream(stream)

To initialize a RFStream from a file use

>>> rfstream = read_rf('test.SAC')

Format specific headers are loaded into the stats object of all traces.

deconvolve(method='P', deconvolve_method='time', **kwargs)[source]

Deconvolve source component of stream from other components.

All args and kwargs are passed to the function deconv().

moveout(phase='Ps', ref=6.4, model='iasp91')[source]

In-place moveout correction to a reference slowness.

Needs stats attributes slowness and onset.

Parameters:
  • phase – ‘Ps’, ‘Sp’, ‘Ppss’ or other multiples
  • ref – reference ray parameter in s/deg
  • model – path to model file or ‘iasp91’
plot_rf(fname=None, norm=1.0, fig_width=7.0, trace_height=0.5, stack_height=0.5, fill=False, window=None, downsample=None, title=True, info=[('back_azimuth', u'baz (xb0)', 'b'), ('distance', u'dist (xb0)', 'r')])[source]

Create receiver function plot.

Parameters:
  • fname – Filename to save plot to. Can be None. In this case the figure is left open.
  • fig_width – Width of figure in inches.
  • trace_height – Height of one trace in inches.
  • fill – Waether to fill receiver functions or not.
  • downsample – Downsample to frequency (in Hz) with Stream.decimate. Filtering is not performed. When saving in a vector format the plot size can be reduced in this way.
  • title – Print seed id as a title.
  • info – Plot one additional axes showing maximal two entries of the stats object. Each entry in this list is a list consisting of three entries: key, label and color. info can be None. In this case no additional axes is plotted.
ppoint(depth, phase='S', model='iasp91')[source]

Calculate coordinates of piercing point by 1D ray tracing.

The iasp91 model is used. Piercing point coordinates are stored in the stats attributes plat and plon. Needs stats attributes station_latitude, station_longitude, slowness and back_azimuth.

Parameters:
  • depth – depth of interface in km
  • phase – ‘P’ for piercing points of P wave, ‘S’ for piercing points of S wave. Multiples are possible, too.
  • model – path to model file or ‘iasp91’
Returns:

NumPy array with coordinates of piercing points

Note

phase=’S’ is usually wanted for P receiver functions and ‘P’ for S receiver functions.

rf(method='P', filter=None, window=None, downsample=None, rotate='ZNE->LQT', source_component='L', deconvolve='time', **kwargs)[source]

Calculate receiver functions in-place.

Parameters:
  • method – ‘P’ for P receiver functions, ‘S’ for S receiver functions
  • filter (dictionary) – filter stream with its filter() method and given kwargs
  • window (tuple of length 2) – trim stream relative to P- or S-onset with trim() (seconds)
  • downsample (float) – downsample stream with its decimate() method
  • rotate – ‘ZNE->LQT’ or ‘NE->RT’, rotate stream with its rotate() method with the angles given by the back_azimuth and inclination attributes of the traces stats objects. You can set these to your needs or let them be computed by rfstats().
  • source_component – character of the source component (i.e. ‘L’ or ‘Z’ depending on rotaion)
  • deconvolve – ‘time’ or ‘freq’ for time or frequency domain deconvolution by the streams deconvolve() method. See deconv(), deconvt() and deconvf() for further documentation.
  • **kwargs – all other kwargs not mentioned here are passed to deconvolve

After performing the deconvolution the Q/R and T components are multiplied by -1 to get a positive phase for a Moho-like positive velocity contrast. Furthermore for method=’S’ all components are mirrored at t=0 for a better comparison with P receiver functions. See source code of this function for the default deconvolution windows.

stack()[source]

Stack traces with the same id.

Traces with same id need to have the same number of datapoints.

write(filename, format, **kwargs)[source]

Save stream to file including format specific headers.

See Stream.write() in ObsPy.

class rf.rfstream.RFTrace(data=array([], dtype=float64), header={}, trace=None, warn=True)[source]

Bases: Trace

Class providing the Trace object for receiver function calculation.

write(filename, format, **kwargs)[source]

Save current trace into a file including format specific headers.

See Trace.write() in ObsPy.

rf.rfstream.obj2stats(event=None, station=None)[source]

Map event and station object to stats with attributes.

Parameters:
  • event – ObsPy Event object
  • station – station object with attributes latitude, longitude and elevation
Returns:

stats object with station and event attributes

rf.rfstream.read_rf(fname=None, format_=None, **kwargs)[source]

Read waveform files into RFStream object.

See read() in ObsPy.

rf.rfstream.rfstats(stats=None, event=None, station=None, stream=None, phase='P', dist_range=None)[source]

Calculate ray specific values like slowness for given event and station.

Parameters:
  • stats – stats object with event and/or station attributes. Can be None if both event and station are given.
  • event – ObsPy Event object
  • station – station object with attributes latitude, longitude and elevation
  • stream – If a stream is given, stats has to be None. In this case rfstats will be called for every stats object in the stream.
  • phase – string with phase to look for in result of getTravelTimes(). Usually this will be ‘P’ or ‘S’ for P and S receiver functions, respectively.
  • dist_range (tuple of length 2) –

    if epicentral of event is not in this intervall, None is returned by this function,

    if phase == ‘P’ defaults to (30, 90),

    if phase == ‘S’ defaults to (50, 85)

Returns:

stats object with event and station attributes, distance, back_azimuth, inclination, onset and slowness or None if epicentral distance is not in the given intervall

rf.rfstream.set_index(index='rf')[source]

batch Module

rf batch command line utility

rf.batch.batch_convert(events, inventory, path, root, newformat, format, **kwargs)[source]
rf.batch.batch_moveout(phase, events, inventory, path, root, format, kwargs_moveout, **kwargs)[source]
rf.batch.batch_plot(events, inventory, path, root, format, kwargs_plot, **kwargs)[source]
rf.batch.batch_rf(init, events, inventory, path, format='H5', request_window=None, phase='P', dist_range=None, **rf_kwargs)[source]
rf.batch.batch_stack(inventory, path, root, format, **kwargs)[source]
rf.batch.main(args=None)[source]

deconvolve Module

Frequency and time domain deconvolution.

rf.deconvolve.deconv(stream, source_component, method='time', **kwargs)[source]

Deconvolve one component of a stream from all components.

The deconvolutions are written to the data arrays of the stream. To keep the original data use the copy method of Stream. The stats dictionaries of the traces inside stream must have an ‘onset’ entry with a UTCDateTime object. This will be used for determining the data windows.

Parameters:
  • stream – Stream object including response and source
  • source_component – Name of component using for source function
  • method

    ‘time’ -> use time domain deconvolution in deconvt(),

    ‘freq’ -> use freqeuency domain deconvolution in deconvf()

  • winsrc (tuple (start, end, taper)) –

    data window for source function, in seconds relative to onset,

    default:

    (-10, 30, 5) for method=’time’,

    (-20, 80, 5) for method=’freq’

  • winrsp (tuple (start, end)) –

    data window for response functions,

    just for method=’time’, default: (-20, 80)

  • winrf (tuple (start, end)) –

    data window for results/deconvolution/receiver functions,

    just for method=’time’, default: (-20, 80)

Other optional parameters are passed to the underlying deconvolution functions deconvt() and deconvf() .

rf.deconvolve.deconvf(rsp_list, src, sampling_rate, water=0.05, gauss=2.0, tshift=10.0, pad=0, length=None, normalize=True, normalize_to_src=False, return_dict=False)[source]

Frequency-domain deconvolution using waterlevel method.

Deconvolve src from arrays in rsp_list.

Parameters:
  • rsp_list – either a list of arrays containing the response functions or a single array
  • src – array of source function
  • sampling_rate – sampling rate of the data
  • water – waterlevel to stabilize the deconvolution
  • gauss – Gauss parameter of Low-pass filter
  • tshift – delay time 0s will be at time tshift afterwards
  • pad – multiply number of samples used for fft by 2**pad
  • length – number of data points in results, optional
  • normalize – if results are normalized
  • normalize_to_src

    True -> normalized so that the maximum of a deconvolution of the source with itself is 1

    False -> normalized so that the maximum of the deconvolution of the first response array in rsp_list is 1

  • return_dict – return additionally a lot of different parameters in a dict for debugging purposes
Returns:

(list of) array(s) with deconvolution(s)

rf.deconvolve.deconvt(rsp_list, src, shift, spiking=1.0, length=None, normalize=True)[source]

Time domain deconvolution.

Deconvolve src from arrays in rsp_list. Calculate Toeplitz auto-correlation matrix of source, invert it, add noise and multiply it with cross-correlation vector of response and source.

In one formula:

RF = (STS + spiking*I)^-1 * STR

N... length
    ( S0   S-1  S-2 ... S-N+1 )
    ( S1   S0   S-1 ... S-N+2 )
S = ( S2   ...                )
    ( ...                     )
    ( SN-1 ...          S0    )
R = (R0 R1 ... RN-1)^T
RF = (RF0 RF1 ... RFN-1)^T
S... source matrix (shape N*N)
R... response vector (length N)
RF... receiver function (deconvolution) vector (length N)
STS = S^T*S = symetric Toeplitz autocorrelation matrix
STR = S^T*R = cross-correlation vector
I... Identity
Parameters:
  • rsp_list – either a list of arrays containing the response functions or a single array
  • src – array of source function
  • shift

    shift the source by that amount of samples to the left side to get onset in RF at the desired time (negative -> shift source to the right side)

    shift = (middle of rsp window - middle of src window) + (0 - middle rf window)

  • spiking – random noise added to autocorrelation (eg. 1.0, 0.1)
  • length – number of data points in results
  • normalize – normalize all results so that the maximum of the first result is 1
Returns:

(list of) array(s) with deconvolution(s)

simple_model Module

Simple move out and piercing point calculation.

This module will probably be replaced by something more sophisticated in future.

class rf.simple_model.SimpleModel(z, vp, vs, n=None)[source]

Bases: object

Simple 1D velocity model.

calculate_phase_delays(slowness, phase='Ps', ref=6.4)[source]

Calculate travel time delays for slowness and reference slowness.

Travel time delays are calculated between the direct wave and the converted phase or multiples.

Parameters:
  • slowness – ray parameter in s/deg
  • phase – ‘Ps’, ‘Sp’ or multiples
  • ref – reference ray parameter in s/deg
Returns:

original times, times stretched to reference slowness

moveout(stream, phase='Ps', ref=6.4)[source]

In-place moveout correction.

Parameters:
  • stream – stream with stats attributes onset and slowness.
  • phase – ‘Ps’, ‘Sp’, ‘Ppss’ or other multiples
  • ref – reference ray parameter in s/deg
ppoint(stats, depth, phase='S')[source]

Piercing point calculation.

Piercing point coordinates are saved in the plat and plon attributes of the stats objects.

Parameters:
  • stats – Stats object with attributes slowness, back_azimuth and station coordinates.
  • depth – depth of interface in km
  • phase – ‘P’ for piercing points of P wave, ‘S’ for piercing points of S wave. Multiples are possible, too.
Returns:

piercing point latitude and longitude

ppoint_distance(depth, slowness, phase='S')[source]

Calculate horizontal distance of piercing point to station.

Parameters:
  • depth – depth of interface in km
  • slowness – ray parameter in s/deg
  • phase – ‘P’ or ‘S’ for P wave or S wave. Multiples possible.
Returns:

horizontal distance in km

rf.simple_model.load_model(fname='iasp91')[source]

Load model from file.

Parameters:fname – path to model file or ‘iasp91’
Returns:SimpleModel instance