PySEP v0.6.0 (#143)

* added read function for asdfdatasets and allow RecSec to plot windows on top of record section (#136)

* added test to check zero division error, (#137)

added logger warning about zero amplitude scaling
feature: user can now set kwargs for obs_color and syn_color

* Improve SAC header append functions (#138)

* extending read_sem_cartesian test to check for evdp related to #132
adding warning messages to sac header append when evdp or mag not
present

* extended cartesian sac headers to match all normal sac header values except for cmpinc and cmpaz
added better warning message for missing sac header values in normal sac header reading
extended tests to cover new functionalities

* Improve RecSec preprocessing architecture (#139)

* reorganizing docstring, restructuring preprocess flag
default value for overwrite is now set True so User does not have to unset
changed default preprocess flag to True, removed 'both' option

* fixed up check function to reflect new preprocess flag

* restructuring preprocessing function to put each individual feature behind a boolean flag so that they can be turned on/off at will
also expose some key arguments as keyword arguments so the User has more control

* all preprocessing now behind tunable knobs

* moved kwargs to args of process_st
reorganized init parameter input to organize a bit better
added new parameters to main docstring

* allow preprocess=True to be more flexible and just take st and st_syn if available, removed hard check on preprocess=True requireing st_syn to be more intutive by users
bugfix move fill value to after assignment of trace in loop

* moved parameters back to internal kwargs to match recsec structure and allow passing them in through command line

* bump version 0.6.0 and bump copyright year in docs

* update changelog

* slim down changelog to include only most recent

Author: bch0w

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/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_recsec.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/tests/test_utils.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/cap_sac.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