Inherit pyatoa gathering and I/O routines (#58)

* moved all pyatoa moment tensor related functionality into a new utility 'mt'

* changed name of util function 'read_synthetics' to 'read_sem' to match Pyatoa original naming (and keep things short for quicker imports)

* added in or fixed up I/O routines from Pyatoa which have been fully migrated here

* added tests and other routines from pyatoa that need to be refactored and organized

* Update README.md

* Update README.md

* fixed failing tests and added a few tests for new mt capabilities

* small formatting changes

* fixed typo in util tests, all tests passing

pysep/pysep.py

@@ -30,7 +30,7 @@
                                  quality_check_waveforms_before_processing,
                                  quality_check_waveforms_after_processing)
 from pysep.utils.fmt import format_event_tag, format_event_tag_legacy, get_codes
-from pysep.utils.io import read_yaml, read_event_file, write_stations_file
+from pysep.utils.io import read_yaml, read_event_file, write_pysep_stations_file
 from pysep.utils.llnl import scale_llnl_waveform_amplitudes
 from pysep.utils.process import (merge_and_trim_start_end_times, resample_data,
                                  format_streams_for_rotation, rotate_to_uvw,
@@ -906,7 +906,8 @@ def rotate_streams(self):
             for sta in stations:
                 _st = st_rtz.select(station=sta)
                 _st.rotate(method="NE->RT")  # in place rot.
-                logger.debug(f"{sta}: BAz={_st[0].stats.back_azimuth}")
+                if hasattr(_st[0].stats, "back_azimuth"):
+                    logger.debug(f"{sta}: BAz={_st[0].stats.back_azimuth}")
             st_out += st_rtz
 
         st_out = format_sac_headers_post_rotation(st_out)
@@ -994,7 +995,7 @@ def write(self, write_files=None, _return_filenames=False,
                                stations_fid or "stations_list.txt")
             logger.info("writing stations file")
             logger.debug(fid)
-            write_stations_file(self.inv, self.event, fid)
+            write_pysep_stations_file(self.inv, self.event, fid)
 
         if "inv" in write_files or "all" in write_files:
             fid = os.path.join(self.output_dir, inv_fid or f"inv.xml")

pysep/recsec.py

@@ -97,7 +97,7 @@
 
 from pysep import logger
 from pysep.utils.cap_sac import origin_time_from_sac_header
-from pysep.utils.io import read_synthetics, read_synthetics_cartesian
+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
 
@@ -491,15 +491,13 @@ def _generate_synthetic_stream(self, syn_path, source, stations,
             for fid in fids:
                 logger.debug(fid)
                 if not cartesian:
-                    st_syn += read_synthetics(fid=fid,
-                                              cmtsolution=source,
-                                              stations=stations)
+                    st_syn += read_sem(fid=fid, source=source,
+                                       stations=stations)
                 else:
                     # If we are using SPECFEM2D synthetics, trying to read
                     # the SOURCE file will
-                    st_syn += read_synthetics_cartesian(fid=fid,
-                                                        source=source,
-                                                        stations=stations)
+                    st_syn += read_sem_cartesian(fid=fid, source=source,
+                                                 stations=stations)
         return st_syn
 
     def check_parameters(self):
@@ -1481,7 +1479,7 @@ def _plot_trace(self, idx, y_index, choice="st"):
             x -= self.zero_pad_s[0]  # index 0 is start, index 1 is end
 
         # Synthetics will already have a time offset stat from the 
-        # 'read_synthetics' function which grabs T0 value from ASCII
+        # 'read_sem' function which grabs T0 value from ASCII
         # Data will have a time offset relative to event origin time 
         if hasattr(tr.stats, "time_offset"):
             x += tr.stats.time_offset  

pysep/tests/test_process.py

@@ -143,6 +143,7 @@ def test_rotate_streams_fail(test_st, test_inv, test_event):
     sep.event = test_event
     sep.st = sep.preprocess()  # make sure that streams are formatted correctly
 
+    # No back aizmuth attribute found
     with pytest.raises(TypeError):
         sep.rotate_streams()
 

pysep/tests/test_utils.py

@@ -18,7 +18,9 @@
 from pysep.utils.fmt import format_event_tag, format_event_tag_legacy
 from pysep.utils.process import estimate_prefilter_corners
 from pysep.utils.plot import plot_source_receiver_map
-from pysep.utils.io import read_synthetics
+from pysep.utils.io import read_sem
+from pysep.utils.mt import get_gcmt_moment_tensor, get_usgs_moment_tensor
+
 
 
 @pytest.fixture
@@ -137,7 +139,7 @@ def test_plot_map(test_event, test_inv):
                              show=False)
 
 
-def test_read_synthetics():
+def test_read_sem():
     """
     Test reading SPECFEM-generated synthetics in as SAC formatted files
     """
@@ -147,8 +149,8 @@ def test_read_synthetics():
 
     st = Stream()
     for test_synthetic in test_synthetics:
-        st += read_synthetics(fid=test_synthetic, cmtsolution=test_cmtsolution,
-                              stations=test_stations)
+        st += read_sem(fid=test_synthetic, source=test_cmtsolution,
+                       stations=test_stations)
 
     assert(st[0].stats.sac.evla == -40.5405)
 
@@ -195,3 +197,39 @@ def test_subset_streams(test_st):
     # all station ids should be the same
     for tr_out, tr_smaller_out in zip(test_st_out, test_st_smaller_out):
         assert(tr_out.get_id() == tr_smaller_out.get_id())
+
+
+def test_get_usgs_moment_tensor():
+    """
+    Just ensure that getting via USGS works as intended using an example event
+    """
+    event = test_get_gcmt_moment_tensor()
+    del event.focal_mechanisms
+
+    cat = get_usgs_moment_tensor(event=event)
+    assert(len(cat) == 1)
+    event = cat[0]
+    assert hasattr(event, "focal_mechanisms")
+
+    m_rr = event.preferred_focal_mechanism().moment_tensor.tensor.m_rr
+    assert(pytest.approx(m_rr, 1E-20) == 4.81E20)
+
+
+@pytest.mark.skip(reason="test already called by 'test_get_usgs_moment_tensor'")
+def test_get_gcmt_moment_tensor():
+    """
+    Just ensure that getting via GCMT works as intended using an example event
+    """
+    # Kaikoura Earthquake
+    origintime = "2016-11-13T11:02:00"
+    magnitude = 7.8
+
+    cat = get_gcmt_moment_tensor(event=None, origintime=origintime,
+                                 magnitude=magnitude)
+    assert(len(cat) == 1)
+    event = cat[0]
+    assert hasattr(event, "focal_mechanisms")
+    m_rr = event.preferred_focal_mechanism().moment_tensor.tensor.m_rr
+    assert(pytest.approx(m_rr, 1E-20) == 3.56E20)
+
+    return event

pysep/utils/fmt.py

@@ -4,6 +4,35 @@
 from obspy.clients.iris import Client
 
 
+def channel_code(dt):
+    """
+    Specfem outputs seismograms with channel band codes set by IRIS. Instrument
+    codes are always X for synthetics, but band code will vary with the sampling
+    rate of the data, return the correct code given a sampling rate.
+    Taken from Appenix B of the Specfem3D cartesian manual (June 15, 2018)
+
+    :type dt: float
+    :param dt: sampling rate of the data in seconds
+    :rtype: str
+    :return: band code as specified by Iris
+    :raises KeyError: when dt is specified incorrectly
+    """
+    if dt >= 1:
+        return "L"  # long period
+    elif 0.1 < dt < 1:
+        return "M"  # mid period
+    elif 0.0125 < dt <= 0.1:
+        return "B"  # broad band
+    elif 0.001 <= dt <= 0.0125:
+        return "H"  # high broad band
+    elif 0.004 <= dt < 0.001:
+        return "C"
+    elif 0.001 <= dt < 0.0002:
+        return "F"
+    else:
+        raise KeyError("Channel code does not exist for this value of 'dt'")
+
+
 def get_codes(st=None, choice=None, suffix=None, up_to=True):
     """
     Get station codes from the internal stream attribute, where station