|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)|
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
Hdf5 file support can be installed with the obspyh5 package.
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
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
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:
Q-file headers DCVREG DCVINCI PWDW are used for the station information, because the Q format has a shortage of predefined headers.
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):
When calling stream.rf the following operations are performed depending on the given kwargs:
- trimming data to window relative to onset
Please see RFStream.rf() for a more detailed description. RFStream provides the possibility to perform moveout correction and piercing point calculation.
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.
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
|RFStream||Class providing the necessary functions for receiver function calculation.|
|rfstats||Calculate ray specific values like slowness for given event and station.|
Classes and functions for receiver function calculation.
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 source component of stream from other components.
All args and kwargs are passed to the function deconv().
In-place moveout correction to a reference slowness.
Needs stats attributes slowness and onset.
Create receiver function plot.
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.
NumPy array with coordinates of piercing points
phase=’S’ is usually wanted for P receiver functions and ‘P’ for S receiver functions.
Calculate receiver functions in-place.
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 traces with the same id.
Traces with same id need to have the same number of datapoints.
Class providing the Trace object for receiver function calculation.
Map event and station object to stats with attributes.
stats object with station and event attributes
Read waveform files into RFStream object.
See read() in ObsPy.
Calculate ray specific values like slowness for given event and station.
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 batch command line utility
Frequency and time domain deconvolution.
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.
Frequency-domain deconvolution using waterlevel method.
Deconvolve src from arrays in rsp_list.
(list of) array(s) with deconvolution(s)
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
(list of) array(s) with deconvolution(s)
Simple move out and piercing point calculation.
This module will probably be replaced by something more sophisticated in future.
Simple 1D velocity model.
Calculate travel time delays for slowness and reference slowness.
Travel time delays are calculated between the direct wave and the converted phase or multiples.
original times, times stretched to reference slowness
In-place moveout correction.
Piercing point calculation.
Piercing point coordinates are saved in the plat and plon attributes of the stats objects.
piercing point latitude and longitude
Calculate horizontal distance of piercing point to station.
horizontal distance in km
Load model from file.
|Parameters:||fname – path to model file or ‘iasp91’|