RF Power Spectral Density (#242)

* Adding a new script for producing PSD plots.

* Adding PSD plots to bulk_rf_report.py

seismic/receiver_fn/bulk_rf_report.py

@@ -376,6 +376,16 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
                 t_stream = rf.RFStream(
                     [tr for tr in transverse_data if tr.stats.event_id in events]).sort(['back_azimuth'])
 
+                # Plot psd
+                fig, ax = plt.subplots()
+                fig.set_size_inches(paper_size_A4[1], paper_size_A4[0])
+                fig.suptitle('.'.join([nsl, channel]))
+                ax.set_rasterized(True)
+
+                rf_plot_utils.plot_rf_psd(rf_stream, ax, min_slope_ratio=min_slope_ratio)
+                pdf.savefig(dpi=300, orientation='landscape')
+                plt.close()
+
                 # Plot pinwheel of primary and transverse components
                 fig = rf_plot_utils.plot_rf_wheel([rf_stream, t_stream], fontscaling=0.8)
                 fig.set_size_inches(*paper_size_A4)

seismic/receiver_fn/rf_plot_utils.py

@@ -13,6 +13,7 @@
 from scipy import stats
 import matplotlib.pyplot as plt
 import matplotlib.patches as patches
+from scipy import signal
 import rf
 
 # pylint: disable=invalid-name, logging-format-interpolation
@@ -40,6 +41,37 @@ def apply(st):
     return result
 # end func
 
+def plot_rf_psd(rf_stream, ax, time_window=(-10.0, 25.0), min_slope_ratio=-1):
+    if(time_window): rf_stream = rf_stream.copy().slice2(*time_window, reftime='onset')
+
+    all_trace_lens = np.array([len(tr) for tr in rf_stream])
+    most_common_len, _ = stats.mode(all_trace_lens, axis=None)
+    psd_stream = rf.RFStream([tr for tr in rf_stream if len(tr) == most_common_len])
+
+    # plot psd
+    psd_array = []
+    fbins = None
+    for trace in psd_stream:
+        fbins, psd = signal.welch(trace.data, fs=trace.stats.sampling_rate,
+                                  detrend='linear')
+
+        psd_array.append(psd)
+        ax.loglog(fbins, psd, alpha=0.05, c='m')
+    # end for
+
+    if (len(psd_array)):
+        psd_array = np.array(psd_array)
+        psd_mean = np.nanmean(psd_array, axis=0)
+
+        ax.loglog(fbins, psd_mean, alpha=1, c='m', lw=2, label='Mean PSD')
+        ax.set_xlabel('Freq. [Hz]')
+        ax.set_ylabel('Power Spectral Density [arb. units]')
+        ax.text(x=0.9, y=0.85, s='{} Traces'.format(len(psd_array)), transform=ax.transAxes)
+        ax.legend()
+        ax.grid()
+    # end if
+# end func
+
 def plot_rf_stack(rf_stream, time_window=(-10.0, 25.0), trace_height=0.2, stack_height=0.8, save_file=None, **kwargs):
     """Wrapper function of rf.RFStream.plot_rf() to help do RF plotting with consistent formatting and layout.
 

seismic/receiver_fn/rf_stacking.py

@@ -182,7 +182,7 @@ def obj_func(x0, amps):
             # end for
         # end for
 
-        tio = interp1d(trc.times() - lead_time, trc.data)
+        tio = interp1d(trc.times() - lead_time, trc.data, bounds_error=False, fill_value=0)
 
         a, b, c = tio(t4), tio(t2), -tio(t3)
         tphase_amps.append([np.sign(a) * np.power(np.fabs(a), 1. / root_order),