PySEP v0.6.0

CHANGELOG.md

@@ -69,6 +69,7 @@
 - Adds version release documentation
 - Slightly modifies pysep-docs conda environment to accomodate converted nbooks
 
+
 ## Version 0.5.0 
 - Improves functions 'read_forcesolution' and 'read_source', which now return
   `obspy.core.event.Event` objects, rather than the makeshift Source objects 
@@ -109,6 +110,15 @@
 - Removed unnused parameters 'legacy_naming' and 'log_level' from 
 	`Pysep.write_config`
 
+
 ## Version 0.5.1
 
 - Bugfix: RecSec subset streams, which checked that 'st' and 'st_syn' had the same stations, would not run for streams of the same length, leading to edge case where same length streams would plot out of order because they had not been sorted. removed the criteria and now subset streams runs at all times
+
+
+## Version 0.6.0
+
+- #136: New read function for ASDFDataSets for misfit window plotting
+- #137: More control over RecSec kwargs and better warning messages
+- #138: Improved SAC header creation for SPECFEM synthetics
+- #139: Improved RecSec preprocessing setup, more manual control for the User 

docs/conf.py

@@ -15,7 +15,7 @@ def setup(app):
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
 
 project = 'PySEP'
-copyright = '2023, adjTomo Dev Team'
+copyright = '2024, adjTomo Dev Team'
 author = 'adjTomo Dev Team'
 release = ''
 # Grab version number from 'pyproject.toml'

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pysep-adjtomo"
-version = "0.5.1"
+version = "0.6.0"
 description = "Python Seismogram Extraction and Processing"
 readme = "README.md"
 requires-python = ">=3.8"

pysep/tests/test_recsec.py

@@ -137,3 +137,24 @@ def test_recsec_calc_time_offset(recsec_w_synthetics):
     recsec_w_synthetics.get_parameters()
     for tr in recsec_w_synthetics.st:
         assert(tr.stats.time_offset == -100)
+
+def test_recsec_zero_amplitude(recsec):
+    """
+    waveforms that have zero amplitude and are normalized should be able 
+    to bypass normalizations which lead to weird plotting (see #131).
+
+    .. note::
+
+        This does not really test that the method is working correctly because
+        dividing a NumPy array by zero leads to NaNs in the array which just
+        won't plot. This is more of a visual test to make sure that the 
+        zero amplitude is plotting correctly, look for green lines
+    """
+    recsec.kwargs.scale_by = "normalize"
+    recsec.kwargs.obs_color = "green"
+    recsec.linewidth = 30
+    for tr in recsec.st:
+        tr.data *= 0
+    recsec.process_st()
+    recsec.get_parameters()
+    recsec.plot()

pysep/tests/test_utils.py

@@ -59,6 +59,16 @@ def test_append_sac_headers(test_st, test_inv, test_event):
     assert(st[0].stats.sac["evla"] == test_event.preferred_origin().latitude)
 
 
+def test_append_sac_headers_cartesian(test_st, test_inv, test_event):
+    """
+    Make sure we can write SAC headers correctly
+    """
+    st = append_sac_headers(st=test_st, inv=test_inv, event=test_event)
+    assert(not hasattr(test_st[0].stats, "sac"))
+    assert(hasattr(st[0].stats, "sac"))
+    assert(st[0].stats.sac["evla"] == test_event.preferred_origin().latitude)
+
+
 def test_event_tag_and_event_tag_legacy(test_event):
     """
     Check that event tagging works as expected
@@ -166,6 +176,14 @@ def test_read_sem():
         st += read_sem(fid=test_synthetic, source=test_cmtsolution,
                        stations=test_stations)
     assert(st)
+
+    expected_headers = ["iztype", "b", "e", "evla", "evlo", "stla", "stlo", 
+                        "stel", "kevnm", "nzyear", "nzjday", "nzhour", "nzmin", 
+                        "nzsec", "nzmsec", "dist", "az", "baz", "gcarc", 
+                        "lpspol", "lcalda", "evdp", "mag"]
+    for expected_header in expected_headers:
+        assert(expected_header in st[0].stats.sac)
+
     assert(st[0].stats.sac.evla == -40.5405)
 
 
@@ -182,7 +200,17 @@ def test_read_sem_cartesian():
         st += read_sem(fid=test_synthetic, source=test_cmtsolution,
                        stations=test_stations)
     assert(st)
+
+    expected_headers = ["iztype", "b", "e", "evla", "evlo", "stla", "stlo", 
+                        "kevnm", "nzyear", "nzjday", "nzhour", "nzmin", 
+                        "nzsec", "nzmsec", "dist", "az", "baz", "gcarc", 
+                        "lpspol", "lcalda", "evdp"]
+    for expected_header in expected_headers:
+        assert(expected_header in st[0].stats.sac)
+
     assert(st[0].stats.sac.stla == 67000.0)
+    assert(st[0].stats.sac.evdp == 30.)
+    assert(st[0].stats.sac.b == -10.)
 
 def test_estimate_prefilter_corners(test_st):
     """

pysep/utils/cap_sac.py

@@ -155,7 +155,7 @@ def _append_sac_headers_trace(tr, event, inv):
 
         We explicitely set 'iztype, 'b' and 'e' in the SAC header to tell ObsPy
         that the trace start is NOT the origin time. Otherwise all the relative
-        timing (e.g., picks) will be wrong.
+        timing (e.g., picks) in SAC will be wrong.
 
     :type tr: obspy.core.trace.Trace
     :param tr: Trace to append SAC header to
@@ -206,6 +206,7 @@ def _append_sac_headers_trace(tr, event, inv):
         "lpspol": 0,  # 1 if left-hand polarity (usually no in passive seis)
         "lcalda": 1,  # 1 if DIST, AZ, BAZ, GCARC to be calc'd from metadata
     }
+
     # Some Inventory objects will not go all the way to channel, only to station
     try:
         cha = sta[0]
@@ -227,6 +228,13 @@ def _append_sac_headers_trace(tr, event, inv):
     except Exception:  # NOQA
         pass
 
+    # Warn User that the following SAC headers could not be found
+    _warn_about = []
+    for key in ["cmpinc", "cmpaz", "evdp", "mag"]:
+        if key not in sac_header:
+            _warn_about.append(key)
+    logger.warning(f"no SAC header values found for: {_warn_about}")
+
     # Append SAC header and include back azimuth for rotation
     tr.stats.sac = sac_header
     tr.stats.back_azimuth = baz
@@ -297,24 +305,29 @@ def _append_sac_headers_cartesian_trace(tr, event, rcv_x, rcv_y):
     :rtype: obspy.core.trace.Trace
     :return: Trace with appended SAC header
     """
-    net_sta = f"{tr.stats.network}.{tr.stats.station}"
-
     src_y = event.preferred_origin().latitude
     src_x = event.preferred_origin().longitude
+    otime = event.preferred_origin().time
+    evdepth_km = event.preferred_origin().depth / 1E3  # units: m -> km
 
     # Calculate Cartesian distance and azimuth/backazimuth
     dist_m = np.sqrt(((rcv_x - src_x) ** 2) + ((rcv_y - src_y) ** 2))
+    dist_km = dist_m * 1E-3  # units: m -> km
+    dist_deg = kilometer2degrees(dist_km)  # spherical earth approximation
     azimuth = np.degrees(np.arctan2((rcv_x - src_x), (rcv_y - src_y))) % 360
     backazimuth = (azimuth - 180) % 360
-    otime = event.preferred_origin().time
-
-    # Barebones SAC header, we only append values required by RecSec
+
     sac_header = {
+        "iztype": 9,  # Ref time equivalence, IB (9): Begin time
+        "b": tr.stats.starttime - otime,  # begin time
+        "e": tr.stats.npts * tr.stats.delta,  # end time
         "stla": rcv_y,
         "stlo": rcv_x,
         "evla": src_y,
         "evlo": src_x,
-        "dist": dist_m * 1E-3,  # units: km
+        "evdp": evdepth_km,
+        "dist": dist_km,
+        "gcarc": dist_deg,  # degrees
         "az": azimuth,
         "baz": backazimuth,
         "kevnm": format_event_tag_legacy(event),  # only take date code
@@ -324,7 +337,9 @@ def _append_sac_headers_cartesian_trace(tr, event, rcv_x, rcv_y):
         "nzmin": otime.minute,
         "nzsec": otime.second,
         "nzmsec": otime.microsecond,
-        }
+        "lpspol": 0,  # 1 if left-hand polarity (usually no in passive seis)
+        "lcalda": 1,  # 1 if DIST, AZ, BAZ, GCARC to be calc'd from metadata
+    }
 
     tr.stats.sac = sac_header
 

pysep/utils/io.py

@@ -14,7 +14,7 @@
 
 from pysep import logger
 from pysep.utils.mt import moment_magnitude, seismic_moment
-from pysep.utils.fmt import format_event_tag_legacy, channel_code
+from pysep.utils.fmt import channel_code
 from pysep.utils.cap_sac import append_sac_headers, append_sac_headers_cartesian
 
 
@@ -498,6 +498,67 @@ def read_event_file(fid):
     return list_out
 
 
+def read_asdfdataset(path, evaluation):
+    """
+    Read an ASDFDataSet created by a SeisFlows Pyaflowa inversion run. 
+    The dataset should contain observed and synthetic waveforms, and 
+    optionally contains misfit windows. Will return Streams with SAC headers
+
+    .. note::
+
+        This function makes assumptions about the PyASDF data structure which
+        is defined by the external package Pyatoa. If Pyatoa changes that
+        structure, this function will break.
+
+    :type path: str
+    :param path: Path to the ASDF dataset to be read
+    :type evaluation: str
+    :param evaluation: evaluation to take synthetics from. These are saved
+        following a format specified by Pyatoa, but usually follow the form
+        iteration/step_count, e.g., i01s00 gives iteration 1, step count 0.
+        Take a look at the waveform tags in `ASDFDataSet.waveforms[<station>]`
+        for tags following the 'synthetic_' prefix
+    """
+    # PySEP, by default will not require PyASDF to be installed
+    try:
+        from pyasdf import ASDFDataSet  # NOQA
+    except ImportError:
+        logger.critical("pyasdf is not installed. Please install pyasdf "
+                        "to read ASDF datasets")
+        return None, None
+
+    with ASDFDataSet(path) as ds:
+        event = ds.events[0]
+        st_out = Stream()
+        st_syn_out = Stream()
+        for station in ds.waveforms.list():
+            inv = ds.waveforms[station].StationXML
+            st = ds.waveforms[station].observed
+            st_syn = ds.waveforms[station][f"synthetic_{evaluation}"]
+
+            st_out += append_sac_headers(st, event, inv)
+            st_syn_out += append_sac_headers(st_syn, event, inv)
+
+        # Read windows from the dataset
+        windows = {}
+        if hasattr(ds.auxiliary_data, "MisfitWindows"):
+            iter_ = evaluation[:3]  # 'i01s00' -> 'i01'
+            step = evaluation[3:]
+            for station in ds.auxiliary_data.MisfitWindows[iter_][step].list():
+                parameters = ds.auxiliary_data.MisfitWindows[iter_][step][
+                    station].parameters
+                trace_id = parameters["channel_id"]
+                starttime = parameters["relative_starttime"]
+                endtime = parameters["relative_endtime"]    
+                # initialize empty window array
+                if trace_id not in windows:
+                    windows[trace_id] = []
+
+                windows[trace_id].append((starttime, endtime))
+
+    return st_out, st_syn_out, windows
+
+
 def write_sem(st, unit, path="./", time_offset=0):
     """
     Write an ObsPy Stream in the two-column ASCII format that Specfem uses