To start, we load example data into an RFStream instance by calling read_rf()
. Calling this function without other arguments loads the sample dataset provided with the rf package. The stream consists of 3-component data from 3 earthquakes, 9 traces in total.
import matplotlib.pyplot as plt
from rf import read_rf, rfstats
stream = read_rf()
# plot the first 3 traces/data from the first event
stream[:3].plot()
<Figure size 800x750 with 3 Axes>
<Figure size 800x750 with 3 Axes>
Next, the rfstats function helps us fill in header info for the receiver function calculation, including onset time.
rfstats(stream)
9 Trace(s) in Stream: CX.PB01..BHE | -20.0s - 100.0s onset:2011-02-25T13:15:38.137852Z | 5.0 Hz, 601 samples | mag:6.0 dist:46.1 baz:325.0 slow:7.83 CX.PB01..BHN | -20.0s - 100.0s onset:2011-02-25T13:15:38.137852Z | 5.0 Hz, 601 samples | mag:6.0 dist:46.1 baz:325.0 slow:7.83 CX.PB01..BHZ | -20.0s - 100.0s onset:2011-02-25T13:15:38.137852Z | 5.0 Hz, 601 samples | mag:6.0 dist:46.1 baz:325.0 slow:7.83 CX.PB01..BHE | -20.1s - 99.9s onset:2011-03-06T14:40:59.799529Z | 5.0 Hz, 601 samples | mag:6.5 dist:47.1 baz:149.2 slow:7.77 CX.PB01..BHN | -20.1s - 99.9s onset:2011-03-06T14:40:59.799529Z | 5.0 Hz, 601 samples | mag:6.5 dist:47.1 baz:149.2 slow:7.77 CX.PB01..BHZ | -20.1s - 99.9s onset:2011-03-06T14:40:59.799529Z | 5.0 Hz, 601 samples | mag:6.5 dist:47.1 baz:149.2 slow:7.77 CX.PB01..BHE | -20.0s - 100.0s onset:2011-05-13T22:54:33.294357Z | 5.0 Hz, 601 samples | mag:6.0 dist:34.2 baz:333.6 slow:8.63 CX.PB01..BHN | -20.0s - 100.0s onset:2011-05-13T22:54:33.294357Z | 5.0 Hz, 601 samples | mag:6.0 dist:34.2 baz:333.6 slow:8.63 CX.PB01..BHZ | -20.0s - 100.0s onset:2011-05-13T22:54:33.294357Z | 5.0 Hz, 601 samples | mag:6.0 dist:34.2 baz:333.6 slow:8.63
The data need a little preprocessing:
stream.filter('bandpass', freqmin=0.33, freqmax=1)
stream.trim2(-25,75,'onset')
stream[:3].plot(type='relative')
Next, we make copies of the stream and do the receiver function calculation (rotation, deconvolution, moveout correction) using different deconvolution methods: time- and frequency-domain water level, and time-domain iterative.
rf_time = stream.copy().rf(deconvolve='time',rotate='NE->RT',solve_toeplitz='scipy').moveout()
rf_freq = stream.copy().rf(deconvolve='waterlevel',rotate='NE->RT').moveout()
rf_iter = stream.copy().rf(deconvolve='iterative',rotate='NE->RT').moveout()
rf_mult = stream.copy().rf(deconvolve='multitaper',rotate='NE->RT').moveout()
rf_time.trim2(-5,20,'onset')
rf_freq.trim2(-5,20,'onset')
rf_iter.trim2(-5,20,'onset')
rf_mult.trim2(-5,20,'onset')
9 Trace(s) in Stream: Prf CX.PB01..BHT | -5.0s - 20.0s onset:2011-02-25T13:15:38.169539Z | 5.0 Hz, 126 samples | mag:6.0 dist:46.1 baz:325.0 slow:6.40 (Ps moveout) Prf CX.PB01..BHR | -5.0s - 20.0s onset:2011-02-25T13:15:38.169539Z | 5.0 Hz, 126 samples | mag:6.0 dist:46.1 baz:325.0 slow:6.40 (Ps moveout) Prf CX.PB01..BHZ | -5.0s - 20.0s onset:2011-02-25T13:15:38.169539Z | 5.0 Hz, 126 samples | mag:6.0 dist:46.1 baz:325.0 slow:6.40 (Ps moveout) Prf CX.PB01..BHT | -5.0s - 20.0s onset:2011-03-06T14:40:59.719539Z | 5.0 Hz, 126 samples | mag:6.5 dist:47.1 baz:149.2 slow:6.40 (Ps moveout) Prf CX.PB01..BHR | -5.0s - 20.0s onset:2011-03-06T14:40:59.719539Z | 5.0 Hz, 126 samples | mag:6.5 dist:47.1 baz:149.2 slow:6.40 (Ps moveout) Prf CX.PB01..BHZ | -5.0s - 20.0s onset:2011-03-06T14:40:59.719539Z | 5.0 Hz, 126 samples | mag:6.5 dist:47.1 baz:149.2 slow:6.40 (Ps moveout) Prf CX.PB01..BHT | -5.0s - 20.0s onset:2011-05-13T22:54:33.319538Z | 5.0 Hz, 126 samples | mag:6.0 dist:34.2 baz:333.6 slow:6.40 (Ps moveout) Prf CX.PB01..BHR | -5.0s - 20.0s onset:2011-05-13T22:54:33.319538Z | 5.0 Hz, 126 samples | mag:6.0 dist:34.2 baz:333.6 slow:6.40 (Ps moveout) Prf CX.PB01..BHZ | -5.0s - 20.0s onset:2011-05-13T22:54:33.319538Z | 5.0 Hz, 126 samples | mag:6.0 dist:34.2 baz:333.6 slow:6.40 (Ps moveout)
We can plot and examine the radial components (since these are PRFs) for each method
rf_time.select(component='R').plot_rf()
rf_freq.select(component='R').plot_rf()
rf_iter.select(component='R').plot_rf()
rf_mult.select(component='R').plot_rf()
We can also make a comparison plot overlaying the RFs for one event from each deconvolution method
evt = 0 # choose 0, 1, or 2
rflist = [rf_time,rf_freq,rf_iter,rf_mult]; labels = ['time','waterlevel','iterative','multitaper']
plt.figure()
for i,st in enumerate(rflist):
sel = st[int(3*evt):int(3*(evt+1))]
sel = sel.select(component='R')
plt.plot(sel[0].times(),sel[0].data,label=labels[i])
plt.legend(fontsize=9)
<matplotlib.legend.Legend at 0x7fccaa567590>
A similar plot can be made with the stacked RFs for each deconvolution method
rflist = [rf_time,rf_freq,rf_iter,rf_mult]; labels = ['time','waterlevel','iterative','multitaper']
plt.figure()
for i,st in enumerate(rflist):
sel = st.copy().select(component='R').stack()
plt.plot(sel[0].times(),sel[0].data,label=labels[i])
plt.legend(fontsize=9)
<matplotlib.legend.Legend at 0x7fcc00060910>