Merge 896ea84 into 601750b

seismic/receiver_fn/notebooks/hk_stacking.ipynb

@@ -326,7 +326,7 @@
     "\n",
     "    # Raise the final sum over phases to power >1 to increase contrast\n",
     "    hk_stack_sum = rf_util.signed_nth_power(hk_stack_sum, 2)\n",
-    "    hk_stack_sum = hk_stack_sum/np.max(hk_stack_sum[:])\n",
+    "    hk_stack_sum = hk_stack_sum/np.nanmax(hk_stack_sum[:])\n",
     "    \n",
     "    sta = station_data[channel][0].stats.station\n",
     "    num = len(station_data[channel])\n",

seismic/receiver_fn/rf_plot_utils.py

@@ -129,10 +129,11 @@ def plot_hk_stack(k_grid, h_grid, hk_stack, title=None, save_file=None, show=Tru
     colmap = 'plasma'
     plt.figure(figsize=(16, 12))
     if clip_negative:
+        hk_stack = hk_stack.copy()
         hk_stack[hk_stack < 0] = 0
     plt.contourf(k_grid, h_grid, hk_stack, levels=50, cmap=colmap)
     cb = plt.colorbar()
-    cb.mappable.set_clim(0, np.max(hk_stack[:]))
+    cb.mappable.set_clim(0, np.nanmax(hk_stack[:]))
     cb.ax.set_ylabel('Stack sum')
     plt.contour(k_grid, h_grid, hk_stack, levels=10, colors='k', linewidths=1)
     plt.xlabel(r'$\kappa = \frac{V_p}{V_s}$ (ratio)', fontsize=14)

seismic/receiver_fn/rf_synthetic.py

@@ -0,0 +1,109 @@
+#!/usr/bin/env python
+"""Helper functions for producing synthetic pseudo-Receiver function traces
+"""
+
+import numpy as np
+from scipy import signal
+
+import obspy
+import rf
+
+import seismic.receiver_fn.rf_util as rf_util
+
+# pylint: disable=invalid-name
+
+def generate_synth_rf(arrival_times, arrival_amplitudes, fs_hz=100.0, window_sec=(-10, 30), f_cutoff_hz=2.0):
+    """Simple generator of synthetic R component receiver function with pulses at given arrival times.
+
+    :param arrival_times: Iterable of arrival times as numerical values in seconds
+    :type arrival_times: iterable of float
+    :param arrival_amplitudes: Iterable of arrival amplitudes
+    :type arrival_amplitudes: iterable of float
+    :param fs_hz: Sampling rate (Hz) of output signal, defaults to 100.0
+    :type fs_hz: float, optional
+    :param window_sec: Time window over which to create signal (sec), defaults to (-10, 30)
+    :type window_sec: tuple, optional
+    :param f_cutoff_hz: Cutoff frequency (Hz) for low-pass filtering to generate realistic result, defaults to 2.0
+    :type f_cutoff_hz: float, optional
+    :return: Array of times and corresponding signal amplitudes
+    :rtype: numpy.array, numpy.array
+    """
+    # Compute array of time values and indexes of arrivals
+    duration = window_sec[1] - window_sec[0]
+    N = int(fs_hz*duration)
+    times = np.linspace(window_sec[0], window_sec[1], N)
+    arrivals_index = np.round((np.array(arrival_times) - times[0])*fs_hz).astype(int)
+
+    # Generate kernel of delta functions at specified arrival times
+    kern = np.zeros_like(times)
+    kern[arrivals_index] = np.array(arrival_amplitudes)  # pylint: disable=unsupported-assignment-operation
+
+    # Filter to pass low frequencies
+    waveform = signal.butter(4, f_cutoff_hz/fs_hz)
+    signal_filt = signal.filtfilt(waveform[0], waveform[1], kern)
+
+    # Normalize signal so max positive amplitude is 1.
+    signal_filt = signal_filt/np.max(signal_filt)
+
+    return times, signal_filt
+# end func
+
+
+def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None):
+    """Synthesize RF R-component data set over range of inclinations and distances
+    and get result as a rf.RFStream instance.
+
+    :param H: Moho depth (km)
+    :type H: float
+    :param V_p: P body wave velocity in uppermost layer
+    :type V_p: float
+    :param V_s: S body wave velocity in uppermost layer
+    :type V_s: float
+    :param inclinations: Array of inclinations for which to create RFs
+    :type inclinations: numpy.array(float)
+    :param distances: Array of teleseismic distances corresponding to inclinations
+    :type distances: numpy.array(float)
+    :param ds: Final sampling rate (Hz) for the downsampled output signal
+    :type ds: float
+    :param log: Logger to send output to, defaults to None
+    :type log: logger, optional
+    :return: Stream containing synthetic RFs
+    :rtype: rf.RFStream
+    """
+    assert len(inclinations) == len(distances), "Must provide 1:1 inclination and distance pairs"
+
+    k = V_p/V_s
+    traces = []
+    for i, inc_deg in enumerate(inclinations):
+        theta_p = np.deg2rad(inc_deg)
+        p = np.sin(theta_p)/V_p
+
+        t1 = H*(np.sqrt((k*k/V_p/V_p) - p*p) - np.sqrt(1.0/V_p/V_p - p*p))
+        t2 = H*(np.sqrt((k*k/V_p/V_p) - p*p) + np.sqrt(1.0/V_p/V_p - p*p))
+        if log is not None:
+            log.info("Inclination {:3g} arrival times: {}".format(inc_deg, [t1, t2]))
+
+        arrivals = [0, t1, t2]
+        amplitudes = [1, 0.5, 0.4]
+        window = (-10.0, 50.0)  # sec
+        fs = 100.0  # Hz
+        _, synth_signal = generate_synth_rf(arrivals, amplitudes, fs_hz=fs, window_sec=window)
+
+        now = obspy.UTCDateTime.now()
+        dt = float(window[1] - window[0])
+        end = now + dt
+        onset = now - window[0]
+        header = {'network': 'SY', 'station': 'TST', 'location': 'GA', 'channel': 'HHR', 'sampling_rate': fs,
+                  'starttime': now, 'endtime': end, 'onset': onset,
+                  'station_latitude': -19.0, 'station_longitude': 137.0,
+                  'slowness': p*rf_util.KM_PER_DEG, 'inclination': inc_deg,
+                  'back_azimuth': 0, 'distance': float(distances[i])}
+        tr = rf.rfstream.RFTrace(data=synth_signal, header=header)
+        tr = tr.decimate(int(np.round(fs/ds)), no_filter=True)
+        traces.append(tr)
+    # end for
+
+    stream = rf.RFStream(traces)
+
+    return stream
+# end func

seismic/receiver_fn/validate_hk_stacking.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python
+"""Use a basic planar, 2-layer model of only the crust and the Moho to generate
+synthetic arrival traces for known model characteristics. Intended to be used
+for model validation.
+"""
+
+import logging
+
+import numpy as np
+from scipy.interpolate import interp1d
+import matplotlib.pyplot as plt
+
+import obspy
+
+import seismic.receiver_fn.rf_plot_utils as rf_plot_utils
+import seismic.receiver_fn.rf_stacking as rf_stacking
+from seismic.receiver_fn.rf_synthetic import synthesize_rf_dataset
+
+# pylint: disable=invalid-name
+
+logging.basicConfig()
+
+
+def main():
+    H = 50  # km
+
+    V_p = 6.4  # km/s, top (crust) layer
+    V_s = 4.0  # km/s, top (crust) layer
+
+    k = V_p/V_s
+    print("H = {:3g}".format(H))
+    print("k = {:3g}".format(k))
+
+    # Generate function mapping ray parameter to teleseismic distance
+    ts_distance = np.linspace(25, 95, 71)
+    inc = np.zeros_like(ts_distance)
+    model = obspy.taup.TauPyModel(model="iasp91")
+    source_depth_km = 10.0
+    for i, d in enumerate(ts_distance):
+        ray = model.get_ray_paths(source_depth_km, d, phase_list=['P'])
+        inc[i] = ray[0].incident_angle  # pylint: disable=unsupported-assignment-operation
+    # end for
+    plt.plot(inc, ts_distance, '-+')
+    plt.xlabel('P-arrival inclination (deg)')
+    plt.ylabel('Teleseismic distance (deg)')
+    plt.title("Relationship between incident angle and teleseismic distance\n(source depth = {:2g} km)".format(
+        source_depth_km))
+    plt.grid("#80808080", linestyle=":")
+    plt.savefig("c:\\temp\\inclination_dist.png", dpi=300)
+
+    # Create interpolator to convert inclination to teleseismic distance
+    interp_dist = interp1d(inc, ts_distance, bounds_error=False, fill_value=(np.max(ts_distance), np.min(ts_distance)))
+
+    # Loop over valid range of inclinations and generate synthetic RFs
+    inclinations = np.linspace(np.min(inc), np.max(inc), 10)
+    distances = interp_dist(inclinations)
+    final_sample_rate_hz = 10.0
+    stream = synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, final_sample_rate_hz, log)
+    stream.write("c:\\temp\\synth_rf_data.h5", format='h5')
+    _ = rf_plot_utils.plot_rf_stack(stream, time_window=(-10, 30), save_file="c:\\temp\\synth_rf.png", dpi=300)
+
+    # Run H-k stacking on synthetic data
+    station_db = {'HHR': stream}
+
+    k, h, hk = rf_stacking.compute_hk_stack(station_db, 'HHR',
+                                            h_range=np.linspace(40.0, 60.0, 101),
+                                            k_range=np.linspace(1.5, 1.8, 151),
+                                            include_t3=False)
+    rf_plot_utils.plot_hk_stack(k, h, hk[0], title="Synthetic Ps component",
+                                save_file="c:\\temp\\Ps.png", show=False)
+    rf_plot_utils.plot_hk_stack(k, h, hk[1], title="Synthetic Ppps component",
+                                save_file="c:\\temp\\Ppps.png", show=False)
+    w = (0.5, 0.5)
+    w_str = "w={:2g},{:2g}".format(*w)
+    stack = rf_stacking.compute_weighted_stack(hk, weighting=w)
+    rf_plot_utils.plot_hk_stack(k, h, stack, num=len(stream), title="Synthetic H-k stack ({})".format(w_str),
+                                save_file="c:\\temp\\synth_stack_{}.png".format(w_str), show=False)
+    # stack = rf_stacking.compute_weighted_stack(hk, weighting=(0.8, 0.2))
+    # rf_plot_utils.plot_hk_stack(k, h, stack, save_file="c:\\temp\\synth_stack_0.8_0.2.png", show=False)
+
+    max_loc = np.unravel_index(np.argmax(stack), stack.shape)
+    h_max = h[max_loc[0], 0]
+    k_max = k[0, max_loc[1]]
+    print("Numerical solution (H, k) = ({:3g}, {:3g})".format(h_max, k_max))
+
+
+if __name__ == "__main__":
+    log = logging.getLogger(__name__)
+    log.setLevel(logging.INFO)
+    main()