improve TauP theoretical arrival usage (#94)

* Update README.md

Update readme with better description of PySEP

* remove hardcode phase list for getting TauP arrival times

* change phase_list default value NoneType -> ('ttall',) to match ObsPy default which gets all phase arrivals

* REMOVED function 'get_taup_arrival_w_sac_headers' which was merged into
'get_taup_arrivals' to avoid redundant code. 'get_taup_arrivals' now
either takes a Stream or Inventory + Event object"

* changed the keys of the phase dict outputted by TauP arrival getting to be more uniform, and better warning messages when phase arrivals can not be found for a station.
arrival time fetching now finds earliest arriving P and S wave, rather than just the P or S arrival

* removed 'model' from SAC header TauP arrival time labels because SAC cuts off the values for long phase names (ie pdiff_ak135 -> pdiff_ak), which seems more confusing than its worht

* allowed SAC header value 'a' to be earliest arriving phase regardless of P or S derived. P and S are still hardcoded later

* removed s and p phase names from incident angle and takeoff angle names in SAC header because long phase names were cutting them off (8 char max), hardcoded to string that sorta looks like the words'

* BUGFIX: SAC header incident angle was being incorrectly set to takeoff angle value, causing no takeoff angle name and incorrect incidient angle name to be attributed

* added 'earliest arrival' to SAC DICT for use in record section

* RecSec now accepts phase arrival times as an acceptable time shift value to align on specific phases

* add log message about how time shift is applied in RecSec

* flip tshift value for phase arrivals to align rather than space out

* fixed failing test for sac header

README.md

@@ -4,8 +4,10 @@ Python Seismogram Extraction and Processing
 [![PyPI Version](https://img.shields.io/pypi/v/pysep-adjtomo.svg)](https://pypi.python.org/pypi/pysep-adjtomo)
 [![SCOPED](https://img.shields.io/endpoint?url=https://runkit.io/wangyinz/scoped/branches/master/pysep)](https://github.com/SeisSCOPED/container/pkgs/container/pysep)
 
-`PySEP` uses ObsPy routines to request and standardize seismic data and metadata, outputting formatted SAC files that are ready for moment tensor and waveform inversion techniques. The package also contains `RecSec` a **rec**ord **sec**tion plotter for rapid visualization of
-waveform data.
+`PySEP` uses ObsPy routines to request and standardize seismic data and metadata, returning uniform, consistent, and minimally processed waveforms for use in moment tensor and waveform inversions. 
+
+The package also contains core classes `RecSec` a **rec**ord **sec**tion plotter for rapid visualization of
+waveform data, and `Declust`, for earthquake catalog declustering and geographical source-receiver weighting for waveform inversions.
 
 - Documentation can found on GitHub pages: https://adjtomo.github.io/pysep/
 - Found a bug, requesting a new feature, or wanting to know how to use`PySEP`? [Feel free to open a GitHub issue](https://github.com/adjtomo/pysep/issues).

pysep/pysep.py

@@ -63,7 +63,7 @@ def __init__(self, config_file=None, event_selection="default",
                  mindistance=0, maxdistance=20E3, minazimuth=0, maxazimuth=360,
                  minlatitude=None, minlongitude=None, maxlatitude=None,
                  maxlongitude=None, resample_freq=None, scale_factor=1,
-                 phase_list=None, seconds_before_event=20,
+                 phase_list=("ttall",), seconds_before_event=20,
                  seconds_after_event=20, seconds_before_ref=100,
                  seconds_after_ref=300, taup_model="ak135", output_unit="VEL",
                  user=None, password=None, client_debug=False,
@@ -328,9 +328,11 @@ def __init__(self, config_file=None, event_selection="default",
 
         :type phase_list: list of str
         :param phase_list: phase names to get ray information from TauP with.
-            Defaults to direct arrivals 'P' and 'S'. Must match Phases expected
-            by TauP (see ObsPy TauP documentation for acceptable phases). Phase
-            information will be added to SAC headers.
+            Defaults to 'ttall', which is ObsPy's default for getting all phase
+            arrivals. Must match Phases expected by TauP (see ObsPy TauP
+            documentation for acceptable phases). Earliest P and S phase
+            arrivals will be added to SAC headers, the remainder will be
+            discarded.
         :type taup_model: str
         :param taup_model: name of TauP model to use to calculate phase arrivals
             See also `phase_list` which defines phases to grab arrival data

pysep/recsec.py

@@ -99,7 +99,7 @@
 from obspy.geodetics import (kilometers2degrees, gps2dist_azimuth)
 
 from pysep import logger
-from pysep.utils.cap_sac import origin_time_from_sac_header
+from pysep.utils.cap_sac import origin_time_from_sac_header, SACDICT
 from pysep.utils.io import read_sem, read_sem_cartesian
 from pysep.utils.curtail import subset_streams
 from pysep.utils.plot import plot_geometric_spreading, set_plot_aesthetic
@@ -233,14 +233,23 @@ def __init__(self, pysep_path=None, syn_path=None, stations=None,
                 distnace, d. 'k' is the `geometric_spreading_k_val` and 'f' is
                 the `geometric_spreading_factor`.
                 'k' is calculated automatically if not given.
-        :type time_shift_s: float OR list of float
+        :type time_shift_s: float OR list of float OR str
         :param time_shift_s: apply static time shift to waveforms, two options:
 
             1. float (e.g., -10.2), will shift ALL waveforms by
                 that number (i.e., -10.2 second time shift applied)
             2. list (e.g., [5., -2., ... 11.2]), will apply individual time
                 shifts to EACH trace in the stream. The length of this list MUST
                 match the number of traces in your input stream.
+            3. str: apply time shift based on a theoretical TauP phase arrival
+                if available in the SAC header. These should have been appended
+                by PySEP during data download. If no value is available in the
+                SAC header, defaults to 0. This may have unintended consequences
+                so you should manually check that all values are okay.
+                Available options are:
+                - 'first_arrival_time': shift based on earliest phase arrival
+                - 'p_arrival_time': shift based on earliest P phase arrival
+                - 's_arrival_time': shift based on earliest S phase arrival
         :type zero_pad_s: list
         :param zero_pad_s: zero pad data in units of seconsd. applied after
             tapering and before filtering. Input as a tuple of floats,
@@ -453,7 +462,12 @@ def __init__(self, pysep_path=None, syn_path=None, stations=None,
 
         # Time shift parameters
         self.move_out = move_out
-        self.time_shift_s = time_shift_s
+        # float check incase command line input provides as a string, and also
+        # to ensure int -> float
+        try:
+            self.time_shift_s = float(time_shift_s)
+        except (TypeError, ValueError):
+            self.time_shift_s = time_shift_s
         self.zero_pad_s = zero_pad_s
 
         # Filtering parameters
@@ -610,12 +624,15 @@ def check_parameters(self):
                 err.scale_by = f"must be in {acceptable_scale_by}"
 
         if self.time_shift_s is not None:
-            acceptable_time_shift_s = [int, float, list]
+            acceptable_time_shift_s = [float, list, str]
             if type(self.time_shift_s) not in acceptable_time_shift_s:
                 err.time_shift_s = f"must be in {acceptable_time_shift_s}"
             if isinstance(self.time_shift_s, list) and \
                     len(self.time_shift_s) != len(self.st):
                 err.time_shift_s = f"must be list of length {len(self.st)}"
+            elif isinstance(self.time_shift_s, str):
+                if self.time_shift_s not in SACDICT:
+                    err.time_shift_s = f"must be {SACDICT} ending w/ '_time'"
 
         if self.zero_pad_s is not None:
             assert(isinstance(self.zero_pad_s, list)), (
@@ -956,11 +973,18 @@ def get_time_shifts(self):
         time_shift_arr = np.zeros(len(self.st))
         if self.time_shift_s is not None:
             # User inputs a static time shift
-            if isinstance(self.time_shift_s, (int, float)):
+            if isinstance(self.time_shift_s, float):
+                logger.info(f"apply constant time shift {self.time_shift_s}s")
                 time_shift_arr += self.time_shift_s
             # User input an array which should have already been checked for len
-            else:
+            elif isinstance(self.time_shift_s, list):
+                logger.info(f"apply user-defined array time shift values")
                 time_shift_arr = self.time_shift_s
+            # Allow shifting by a given phase arrival in the SAC header
+            elif isinstance(self.time_shift_s, str):
+                sac_key = SACDICT[self.time_shift_s]
+                logger.info(f"apply time shift by {self.time_shift_s}")
+                time_shift_arr = [-1 * tr.stats.sac[sac_key] for tr in self.st]
         time_shift_arr = np.array(time_shift_arr)
 
         # Further change the time shift if we have move out input
@@ -2103,8 +2127,9 @@ def parse_args():
                         nargs="?", type=float,
                         help="Controls the y-max value on the geometric "
                              "spreading plot.")
-    parser.add_argument("--time_shift_s", default=None, type=float, nargs="?",
-                        help="Set a constant time shift in unit: seconds")
+    parser.add_argument("--time_shift_s", default=None, nargs="?",
+                        help="Set a constant time shift in unit: seconds OR "
+                             "shift by a given phase arrival in SAC header")
     parser.add_argument("--move_out", default=None, type=float, nargs="?",
                         help="Set a constant velocity-based move out in units:"
                              "`distance_units`/s")

pysep/tests/test_utils.py

@@ -77,7 +77,7 @@ def test_format_sac_headers_w_taup_traveltimes(test_st, test_inv, test_event):
     """
     st = append_sac_headers(st=test_st, inv=test_inv, event=test_event)
     st = format_sac_header_w_taup_traveltimes(st=st, model="ak135")
-    assert(pytest.approx(st[0].stats.sac["user4"], .01) == 57.66)
+    assert(pytest.approx(st[0].stats.sac["a"], .01) == 224.535)
 
 
 def test_sac_header_correct_origin_time(tmpdir, test_st, test_inv, test_event):

pysep/utils/cap_sac.py

@@ -12,10 +12,11 @@
 
 from pysep import logger
 from pysep.utils.fmt import format_event_tag_legacy
-from pysep.utils.fetch import get_taup_arrivals_with_sac_headers
+from pysep.utils.fetch import get_taup_arrivals
 
 # SAC HEADER CONSTANTS DEFINING NON-INTUITIVE QUANTITIES
 SACDICT = {
+    "first_arrival_time": "a",
     "p_arrival_time": "t5",
     "p_incident_angle": "user1",
     "p_takeoff_angle": "user3",
@@ -32,7 +33,6 @@ def write_cap_weights_files(st, path_out="./", order_by="dist"):
 
     TODO re-add Ptime setting with event.picks
 
-
     Replaces `write_cap_weights`
 
     The weight file has columns corresponding to the following:
@@ -233,7 +233,8 @@ def _append_sac_headers_trace(tr, event, inv):
     return tr
 
 
-def format_sac_header_w_taup_traveltimes(st, model="ak135", phase_list=None):
+def format_sac_header_w_taup_traveltimes(st, model="ak135",
+                                         phase_list=("ttall",)):
     """
     Add TauP travel times to the SAC headers using information in the SAC header
     Also get some information from TauP regarding incident angle, takeoff angle
@@ -256,60 +257,66 @@ def format_sac_header_w_taup_traveltimes(st, model="ak135", phase_list=None):
         defaults to 'ak135'
     :type phase_list: list of str
     :param phase_list: phase names to get ray information from TauP with.
-        Defaults to direct arrivals 'P' and 'S'. Must match Phases expected
-        by TauP (see ObsPy TauP documentation for acceptable phases).
+        Defaults to 'ttall', which is ObsPy's default for getting all phase
+        arrivals. Must match Phases expected by TauP (see ObsPy TauP
+        documentation for acceptable phases).
     """
     st_out = st.copy()
+    phase_dict = get_taup_arrivals(st=st, model=model, phase_list=phase_list)
 
-    # Call TauP with a specific model to retrieve travel times etc.
-    if not phase_list:
-        phase_list = ["p", "P", "s", "S"]
-
-    phase_dict = get_taup_arrivals_with_sac_headers(st=st, model=model,
-                                                    phase_list=phase_list)
     # Arrivals may return multiple entires for each phase, pick earliest
     for tr in st_out[:]:
+        net_sta = ".".join(tr.get_id().split(".")[:2])  # NN.SSS
+
         # Missing SAC header values may cause certain or all stations to not 
         # be present in the `phase_dict`
-        if tr.get_id() not in phase_dict:
+        if net_sta not in phase_dict:
+            logger.warning(f"{tr.get_id()} not present in TauP arrivals, cant "
+                           f"append arrival time SAC headers")
             continue
-        arrivals = phase_dict[tr.get_id()]
+        arrivals = phase_dict[net_sta]
         # If the TauP arrival calculation fails, `arrivals` will be empty
         if not arrivals:
-            logger.warning(f"{tr.get_id()} has no TauP phase arrivals, cannot "
+            logger.warning(f"{tr.get_id()} has no TauP phase arrivals, cant "
                            f"append arrival time SAC headers")
             continue
-        # Find earliest arriving P-wave (P or p)
-        idx_times = [(i, a.time) for i, a in enumerate(arrivals) if
-                     a.name.upper() == "P"]
+
+        # Find earliest arriving phase (likely same as P phase but maybe not)
+        idx_times = [(i, a.time) for i, a in enumerate(arrivals)]
         idx, _ = min(idx_times, key=lambda x: x[1])  # find index of min time
-        p = arrivals[idx]  # Earliest P-wave Arrival object
 
-        tr.stats.sac["a"] = p.time  # relative time sec: float
-        tr.stats.sac["ka"] = f"{p.name}_{model}"  # name: str
+        phase = arrivals[idx]  # Earliest Arrival object
+        tr.stats.sac[SACDICT["first_arrival_time"]] = phase.time  # 'a'
+        tr.stats.sac[f"k{SACDICT['first_arrival_time']}"] = f"{phase.name}" #ka
 
+        # Find earliest arriving P phase (must start with letter 'P')
+        idx_times = [(i, a.time) for i, a in enumerate(arrivals) if
+                     a.name.upper().startswith("P")]
+        idx, _ = min(idx_times, key=lambda x: x[1])  # find index of min time
+
+        p = arrivals[idx]  # Earliest P-wave Arrival object
         tr.stats.sac[SACDICT["p_arrival_time"]] = p.time
-        tr.stats.sac[f"k{SACDICT['p_arrival_time']}"] = f"{p.name}_{model}"
+        tr.stats.sac[f"k{SACDICT['p_arrival_time']}"] = f"{p.name}"
 
-        # P-wave incident angle (ia) and takeoff angle (ta)
+        # P phase incident angle (ia) and takeoff angle (ta)
         tr.stats.sac[SACDICT["p_incident_angle"]] = p.incident_angle
-        tr.stats.sac[f"k{SACDICT['p_incident_angle']}"] = f"{p.name}_ia_{model}"
+        tr.stats.sac[f"k{SACDICT['p_incident_angle']}"] = f"p_IncAng"
         tr.stats.sac[SACDICT["p_takeoff_angle"]] = arrivals[idx].takeoff_angle
-        tr.stats.sac[f"k{SACDICT['p_incident_angle']}"] = f"{p.name}_ta_{model}"
+        tr.stats.sac[f"k{SACDICT['p_takeoff_angle']}"] = f"p_TkfAng"
 
-        # Find earliest arriving S-wave (S or s)
+        # Find earliest arriving S phase (must start with letter 'S')
         idx_times = [(i, a.time) for i, a in enumerate(arrivals) if
-                     a.name.upper() == "S"]
+                     a.name.upper().startswith("S")]
         idx, _ = min(idx_times, key=lambda x: x[1])
         s = arrivals[idx]  # Earliest S-wave Arrival object
 
         tr.stats.sac[SACDICT["s_arrival_time"]] = s.time
-        tr.stats.sac[f"k{SACDICT['s_arrival_time']}"] = f"{s.name}_{model}"
+        tr.stats.sac[f"k{SACDICT['s_arrival_time']}"] = f"{s.name}"
 
         tr.stats.sac[SACDICT["s_incident_angle"]] = s.incident_angle
-        tr.stats.sac[f"k{SACDICT['s_incident_angle']}"] = f"{s.name}_ia_{model}"
+        tr.stats.sac[f"k{SACDICT['s_incident_angle']}"] = f"s_IncAng"
         tr.stats.sac[SACDICT["s_takeoff_angle"]] = s.takeoff_angle
-        tr.stats.sac[f"k{SACDICT['s_incident_angle']}"] = f"{s.name}_ta_{model}"
+        tr.stats.sac[f"k{SACDICT['s_takeoff_angle']}"] = f"s_TkfAng"
 
     return st_out