|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
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
>>> from obspy import read >>> from rf import RFStream >>> stream = RFStream(stream=read('infile.SAC'))
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.
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 is going to provide a function rf_batch() which will run all the necessary steps.
|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 ObsPy stream use
>>> rfstream = RFStream(stream=obspy_stream)
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.
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.
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
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
Functions for massive receiver function calculation. Untestet, in development.
TODO: doc rf_batch
TODO: doc rf_client
TODO: doc rf_dmt
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’|