Merge c9545cb into 358ad6c

seismic/receiver_fn/moho_discontinuity_hk_study.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+"""Study behavior of H-k stacking in presence of Moho discontinuity. Ticket PST-433.
+"""
+
+import os
+import logging
+
+import numpy as np
+# import obspy
+import matplotlib.pyplot as plt
+
+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, convert_inclination_to_distance
+
+# pylint: disable=invalid-name,logging-format-interpolation
+
+logging.basicConfig()
+
+
+def main():
+    output_folder = 'moho_discontinuity'
+    if not os.path.exists(output_folder):
+        os.mkdir(output_folder)
+    # end if
+
+    # Generate data files corresponding to two different Moho depths
+    V_s = 3.8  # km/s
+    V_p = 6.4  # km/s
+    H_0 = 40  # km
+    k = V_p/V_s
+    log.info("H_0 = {:.3f}".format(H_0))
+    log.info("k = {:.3f}".format(k))
+
+    # Loop over valid range of inclinations and generate synthetic RFs
+    num_traces = 15
+    inclinations = np.linspace(14.0, 27.0, num_traces)
+    distances = convert_inclination_to_distance(inclinations)
+    final_sample_rate_hz = 10.0
+    # stream_0 = synthesize_rf_dataset(H_0, V_p, V_s, inclinations, distances, final_sample_rate_hz, log, baz=0)
+    stream_0 = synthesize_rf_dataset(H_0, V_p, V_s, inclinations, distances, final_sample_rate_hz, log,
+                                    #  amplitudes=[1, 0.5, 0.2], baz=0)
+                                     amplitudes=[1, 0.4, 0.3], baz=0)
+    stream_0.write(os.path.join(output_folder, "synth_rf_data_H={}.h5".format(H_0)), format='h5')
+
+    H_1 = 50
+    # stream_1 = synthesize_rf_dataset(H_1, V_p, V_s, inclinations, distances, final_sample_rate_hz, log, baz=90)
+    stream_1 = synthesize_rf_dataset(H_1, V_p, V_s, inclinations, distances, final_sample_rate_hz, log,
+                                    #  amplitudes=[1, 0.4, 0.3], baz=90)
+                                     amplitudes=[1, 0.5, 0.2], baz=90)
+    stream_1.write(os.path.join(output_folder, "synth_rf_data_H={}.h5".format(H_1)), format='h5')
+
+    # Plot each dataset separately, then aggregated together
+    rf_fig_0 = rf_plot_utils.plot_rf_stack(stream_0, time_window=(-5, 30))
+    plt.title("Exact H = {:.3f}, k = {:.3f}".format(H_0, k))
+    rf_fig_0.savefig(os.path.join(output_folder, "synth_rf_H={}.png".format(H_0)), dpi=300)
+    rf_fig_1 = rf_plot_utils.plot_rf_stack(stream_1, time_window=(-5, 30))
+    plt.title("Exact H = {:.3f}, k = {:.3f}".format(H_1, k))
+    rf_fig_1.savefig(os.path.join(output_folder, "synth_rf_H={}.png".format(H_1)), dpi=300)
+    stream_all = stream_0 + stream_1
+    rf_fig_all = rf_plot_utils.plot_rf_stack(stream_all, time_window=(-5, 30))
+    plt.title("Exact H = {:.3f}+{:.3f}, k = {:.3f}".format(H_0, H_1, k))
+    rf_fig_all.savefig(os.path.join(output_folder, "synth_rf_H={}+{}.png".format(H_0, H_1)), dpi=300)
+
+    # Run H-k stacking on synthetic data
+    station_db = {'HHR': stream_all}
+    k_grid, h_grid, hk = rf_stacking.compute_hk_stack(station_db, 'HHR', include_t3=False, root_order=2)
+    w = (0.5, 0.5)
+    w_str = "w={:1.2f},{:1.2f}".format(*w)
+    stack = rf_stacking.compute_weighted_stack(hk, weighting=w)
+    hk_fig = rf_plot_utils.plot_hk_stack(k_grid, h_grid, stack, num=len(stream_all),
+                                         title="Synthetic H-k stack ({})".format(w_str))
+    maxima = rf_stacking.find_local_hk_maxima(k_grid, h_grid, stack, min_rel_value=0.9)
+    log.info("Numerical solutions:{}".format(maxima))
+    plt.text(k, H_0 + 1.5, "Solution $H_0$ = {:.3f}, k = {:.3f}".format(H_0, k),
+             color="#ffffff", fontsize=16, horizontalalignment='left')
+    plt.text(k, H_1 + 1.5, "Solution $H_1$ = {:.3f}, k = {:.3f}".format(H_1, k),
+             color="#ffffff", fontsize=16, horizontalalignment='left')
+    for h_max, k_max, _, _, _ in maxima:
+        plt.scatter(k_max, h_max, marker='+', c="#000000", s=20)
+    # end for
+    # Save result
+    hk_fig.savefig(os.path.join(output_folder, "synth_stack_{}.png".format(w_str)), dpi=300)
+
+# end main
+
+
+if __name__ == "__main__":
+    log = logging.getLogger(__name__)
+    log.setLevel(logging.INFO)
+    main()

seismic/receiver_fn/rf_synthetic.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 from scipy import signal
+from scipy.interpolate import interp1d
 
 import obspy
 import rf
@@ -49,7 +50,8 @@ def generate_synth_rf(arrival_times, arrival_amplitudes, fs_hz=100.0, window_sec
 # end func
 
 
-def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, include_t3=False):
+def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, include_t3=False,
+                          amplitudes=None, baz=0.0):
     """Synthesize RF R-component data set over range of inclinations and distances
     and get result as a rf.RFStream instance.
 
@@ -69,6 +71,10 @@ def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, in
     :type log: logger, optional
     :param include_t3: If True, include the third expected multiple PpSs+PsPs
     :type include_t3: bool, optional
+    :param amplitudes: Custom amplitudes to apply to the multiples
+    :type amplitudes: list(float), optional
+    :param baz: Back azimuth for metadata
+    :type baz: float, optional
     :return: Stream containing synthetic RFs
     :rtype: rf.RFStream
     """
@@ -90,9 +96,16 @@ def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, in
             log.info("Inclination {:3g} arrival times: {}".format(inc_deg, arrivals))
 
         arrivals = [0] + arrivals
-        amplitudes = [1, 0.5, 0.4]
-        if include_t3:
-            amplitudes.append(-0.3)
+        if amplitudes is None:
+            amplitudes = [1, 0.5, 0.4]
+            if include_t3:
+                amplitudes.append(-0.3)
+            # end if
+        else:
+            assert len(amplitudes) == 3 + int(include_t3)
+            # t3 amplitude should be negative
+            assert (not include_t3) or (amplitudes[3] <= 0)
+        # end if
         window = (-5.0, 50.0)  # sec
         fs = 100.0  # Hz
         _, synth_signal = generate_synth_rf(arrivals, amplitudes, fs_hz=fs, window_sec=window)
@@ -105,7 +118,7 @@ def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, in
                   'starttime': now, 'endtime': end, 'onset': onset,
                   'station_latitude': -19.0, 'station_longitude': 137.0,  # arbitrary (approx location of OA deployment)
                   'slowness': p*rf_util.KM_PER_DEG, 'inclination': inc_deg,
-                  'back_azimuth': 0, 'distance': float(distances[i])}
+                  'back_azimuth': baz, '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)
@@ -115,3 +128,37 @@ def synthesize_rf_dataset(H, V_p, V_s, inclinations, distances, ds, log=None, in
 
     return stream
 # end func
+
+
+def convert_inclination_to_distance(inclinations, model="iasp91", nominal_source_depth_km=10.0):
+    """Helper function to convert range of inclinations to teleseismic distance in degrees.
+
+    :param inclinations: Array of inclination angles in degrees
+    :type inclinations: numpy.array(float)
+    :param model: Name of model to use for ray tracing, defaults to "iasp91"
+    :type model: str, optional
+    :param nominal_source_depth_km: Assumed depth of source events, defaults to 10.0
+    :type nominal_source_depth_km: float, optional
+    :return: Array of teleseismic distances in degrees corresponding to input inclination angles.
+    :rtype: numpy.array(float)
+    """
+    # Generate function mapping ray parameter to teleseismic distance.
+    # The distances are not strictly required for H-k stacking, but rf behaves better when they are there.
+    ts_distance = np.linspace(25, 95, 71)
+    inc = np.zeros_like(ts_distance)
+    model = obspy.taup.TauPyModel(model=model)
+    source_depth_km = nominal_source_depth_km
+    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
+
+    # 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
+    distances = interp_dist(inclinations)
+
+    return distances
+# end func

seismic/receiver_fn/validate_hk_stacking.py

@@ -1,7 +1,7 @@
 #!/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.
+synthetic arrival traces for known model characteristics, and run the RFs through
+H-k stacking routine the recover the original H-k values.
 """
 
 import os