This notebook demonstrates the most simple use case of the rf package.
If you want to run this notebook locally, have a look at the repository readme. Dependencies of this notebook are ObsPy and rf.
We load example data into a RFStream instance by simply calling read_rf()
. The stream consists of 9 traces with three component data from three different earthquakes. The header values of the first trace in the stream is printed and we can observe that the SAC headers in stats.sac
are correctly mapped to corresponding entries in stats
. Raw waveforms of the first event are plotted with ObsPy.
import matplotlib.pyplot as plt
from rf import read_rf, rfstats
stream = read_rf()
print(stream)
print('\nStats:\n', stream[0].stats)
stream[:3].plot()
9 Trace(s) in Stream: CX.PB01..BHE | 2011-02-25T13:15:18.169539Z - 2011-02-25T13:17:18.169539Z | 5.0 Hz, 601 samples CX.PB01..BHN | 2011-02-25T13:15:18.169539Z - 2011-02-25T13:17:18.169539Z | 5.0 Hz, 601 samples CX.PB01..BHZ | 2011-02-25T13:15:18.169539Z - 2011-02-25T13:17:18.169539Z | 5.0 Hz, 601 samples CX.PB01..BHE | 2011-03-06T14:40:39.719539Z - 2011-03-06T14:42:39.719539Z | 5.0 Hz, 601 samples CX.PB01..BHN | 2011-03-06T14:40:39.719538Z - 2011-03-06T14:42:39.719538Z | 5.0 Hz, 601 samples CX.PB01..BHZ | 2011-03-06T14:40:39.719539Z - 2011-03-06T14:42:39.719539Z | 5.0 Hz, 601 samples CX.PB01..BHE | 2011-05-13T22:54:13.319538Z - 2011-05-13T22:56:13.319538Z | 5.0 Hz, 601 samples CX.PB01..BHN | 2011-05-13T22:54:13.319536Z - 2011-05-13T22:56:13.319536Z | 5.0 Hz, 601 samples CX.PB01..BHZ | 2011-05-13T22:54:13.319538Z - 2011-05-13T22:56:13.319538Z | 5.0 Hz, 601 samples Stats: network: CX station: PB01 location: channel: BHE starttime: 2011-02-25T13:15:18.169539Z endtime: 2011-02-25T13:17:18.169539Z sampling_rate: 5.0 delta: 0.2 npts: 601 calib: 1.0 _format: SAC back_azimuth: 325.033 distance: 46.1505 event_depth: 130.6 event_latitude: 17.8214 event_longitude: -95.1708 event_magnitude: 6.0 event_time: 2011-02-25T13:07:26.980005Z sac: AttribDict({'delta': 0.2, 'depmin': -198.0, 'depmax': 1325.0, 'scale': 1.0, 'b': 0.00053899997, 'e': 120.00054, 'o': -471.189, 'stla': -21.04323, 'stlo': -69.487396, 'stel': 900.0, 'evla': 17.8214, 'evlo': -95.170799, 'evdp': 130.60001, 'mag': 6.0, 'dist': 5131.6958, 'az': 145.81221, 'baz': 325.03323, 'gcarc': 46.150452, 'depmen': 579.33942, 'nzyear': 2011, 'nzjday': 56, 'nzhour': 13, 'nzmin': 15, 'nzsec': 18, 'nzmsec': 169, 'nvhdr': 6, 'npts': 601, 'iftype': 1, 'iztype': 9, 'leven': 1, 'lpspol': 0, 'lovrok': 1, 'lcalda': 1, 'kstnm': 'PB01', 'kcmpnm': 'BHE', 'knetwk': 'CX'}) station_elevation: 900.0 station_latitude: -21.0432 station_longitude: -69.4874
Now, we fill the stats dictionary with entries which we need for receiver function calculation with the rfstats
function. The print statement automatically displays some important header values. Data is filtered, trimmed and the preprocessed data of the first event is again plotted with ObsPy.
rfstats(stream)
stream.filter('bandpass', freqmin=0.4, freqmax=1)
stream.trim2(5, 95, 'starttime')
print(stream)
stream[:3].plot(type='relative', reftime=stream[0].stats.onset)
9 Trace(s) in Stream: CX.PB01..BHE | -15.0s - 75.0s onset:2011-02-25T13:15:38.137852Z | 5.0 Hz, 451 samples | mag:6.0 dist:46.1 baz:325.0 slow:7.83 CX.PB01..BHN | -15.0s - 75.0s onset:2011-02-25T13:15:38.137852Z | 5.0 Hz, 451 samples | mag:6.0 dist:46.1 baz:325.0 slow:7.83 CX.PB01..BHZ | -15.0s - 75.0s onset:2011-02-25T13:15:38.137852Z | 5.0 Hz, 451 samples | mag:6.0 dist:46.1 baz:325.0 slow:7.83 CX.PB01..BHE | -15.1s - 74.9s onset:2011-03-06T14:40:59.799529Z | 5.0 Hz, 451 samples | mag:6.5 dist:47.1 baz:149.2 slow:7.77 CX.PB01..BHN | -15.1s - 74.9s onset:2011-03-06T14:40:59.799529Z | 5.0 Hz, 451 samples | mag:6.5 dist:47.1 baz:149.2 slow:7.77 CX.PB01..BHZ | -15.1s - 74.9s onset:2011-03-06T14:40:59.799529Z | 5.0 Hz, 451 samples | mag:6.5 dist:47.1 baz:149.2 slow:7.77 CX.PB01..BHE | -15.0s - 75.0s onset:2011-05-13T22:54:33.294357Z | 5.0 Hz, 451 samples | mag:6.0 dist:34.2 baz:333.6 slow:8.63 CX.PB01..BHN | -15.0s - 75.0s onset:2011-05-13T22:54:33.294357Z | 5.0 Hz, 451 samples | mag:6.0 dist:34.2 baz:333.6 slow:8.63 CX.PB01..BHZ | -15.0s - 75.0s onset:2011-05-13T22:54:33.294357Z | 5.0 Hz, 451 samples | mag:6.0 dist:34.2 baz:333.6 slow:8.63
Finally, we perform the receiver function calculation, i.e. rotation and deconvolution, by calling RFStream.rf()
. The L components of the receiver functions have the expected peak at 0s. The Q component stack shows a distinct phase at around 9s. That's it.
stream.rf()
stream.moveout()
stream.trim2(-5, 22, 'onset')
print(stream)
stream.select(component='L').plot_rf()
stream.select(component='Q').plot_rf()
plt.show()
9 Trace(s) in Stream: Prf CX.PB01..BHT | -5.0s - 22.0s onset:2011-02-25T13:15:38.169539Z | 5.0 Hz, 136 samples | mag:6.0 dist:46.1 baz:325.0 slow:6.40 (Ps moveout) Prf CX.PB01..BHQ | -5.0s - 22.0s onset:2011-02-25T13:15:38.169539Z | 5.0 Hz, 136 samples | mag:6.0 dist:46.1 baz:325.0 slow:6.40 (Ps moveout) Prf CX.PB01..BHL | -5.0s - 22.0s onset:2011-02-25T13:15:38.169539Z | 5.0 Hz, 136 samples | mag:6.0 dist:46.1 baz:325.0 slow:6.40 (Ps moveout) Prf CX.PB01..BHT | -5.0s - 22.0s onset:2011-03-06T14:40:59.719539Z | 5.0 Hz, 136 samples | mag:6.5 dist:47.1 baz:149.2 slow:6.40 (Ps moveout) Prf CX.PB01..BHQ | -5.0s - 22.0s onset:2011-03-06T14:40:59.719539Z | 5.0 Hz, 136 samples | mag:6.5 dist:47.1 baz:149.2 slow:6.40 (Ps moveout) Prf CX.PB01..BHL | -5.0s - 22.0s onset:2011-03-06T14:40:59.719539Z | 5.0 Hz, 136 samples | mag:6.5 dist:47.1 baz:149.2 slow:6.40 (Ps moveout) Prf CX.PB01..BHT | -5.0s - 22.0s onset:2011-05-13T22:54:33.319538Z | 5.0 Hz, 136 samples | mag:6.0 dist:34.2 baz:333.6 slow:6.40 (Ps moveout) Prf CX.PB01..BHQ | -5.0s - 22.0s onset:2011-05-13T22:54:33.319538Z | 5.0 Hz, 136 samples | mag:6.0 dist:34.2 baz:333.6 slow:6.40 (Ps moveout) Prf CX.PB01..BHL | -5.0s - 22.0s onset:2011-05-13T22:54:33.319538Z | 5.0 Hz, 136 samples | mag:6.0 dist:34.2 baz:333.6 slow:6.40 (Ps moveout)