From 9e1204588c060f71982b4d4bece6382908b88c20 Mon Sep 17 00:00:00 2001
From: Bryant Chow <bryant.h.c@gmail.com>
Date: Fri, 16 Sep 2022 14:52:00 -0800
Subject: [PATCH] Merge Devel (#34)

* bugfix plotting synthetic record section wasnt working due to premature assertion

* bugfix dont let lack of phase arrivals from TauP break SAC header appending, skip over and warn user instead.
added recsec kwarg plotting to readme

* added debug message for phase arrival list that does not return any phase arrivals during TauP search

* feature: recsec can now read specfem3D cmtsolutions that are defined in Cartesian coordinate system, allowing for record sections of specfem3d waveforms generated in cartesian domains
---
 pysep/utils/io.py | 66 ++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 63 insertions(+), 3 deletions(-)

diff --git a/pysep/utils/io.py b/pysep/utils/io.py
index ec9819c..164c061 100644
--- a/pysep/utils/io.py
+++ b/pysep/utils/io.py
@@ -74,7 +74,7 @@ def read_synthetics(fid, cmtsolution, stations, location="", precision=4):
     delta = round(times[1] - times[0], precision)
 
     # Get metadata information from CMTSOLUTION and STATIONS files
-    event = read_events(cmtsolution)[0]
+    event = read_events(cmtsolution, format="CMTSOLUTION")[0]
     inv = read_stations(stations)
 
     origintime = event.preferred_origin().time
@@ -158,8 +158,9 @@ def read_synthetics_cartesian(fid, source, stations, location="", precision=4):
 
     # Get metadata information from CMTSOLUTION and STATIONS files
     try:
-        event = read_events(source)[0]
-    except TypeError:
+        event = read_specfem3d_cmtsolution_cartesian(source)
+    # Specfem2D and 3D source/cmtsolution files have different formats
+    except ValueError:
         event = read_specfem2d_source(source)
 
     # Generate a dictionary object to store station information
@@ -224,6 +225,65 @@ def read_synthetics_cartesian(fid, source, stations, location="", precision=4):
     return st
 
 
+def read_specfem3d_cmtsolution_cartesian(path_to_cmtsolution):
+    """
+    Create a barebones ObsPy Event object from a SPECFEM3D CMTSOLUTION file with
+    coordinates defined in cartesian. Required because ObsPy read_events() will
+    throw a ValueError when CMTSOLUTIONS do not have geographically defined
+    coordinates.
+    """
+    def _get_resource_id(name, res_type, tag=None):
+        """
+        Helper function to create consistent resource ids. From ObsPy
+        """
+        res_id = f"smi:local/source/{name:s}/{res_type:s}"
+        if tag is not None:
+            res_id += "#" + tag
+        return res_id
+
+    with open(path_to_cmtsolution, "r") as f:
+        lines = f.readlines()
+
+    # First line contains meta informatino about event
+    _, year, month, day, hour, minute, sec, lat, lon, depth, mb, ms, _ = \
+        lines[0].strip().split()
+
+    origin_time = UTCDateTime(f"{year}-{month}-{day}T{hour}:{minute}:{sec}")
+
+    # Remaining lines contain information on event, some useful some not
+    source_dict = {}
+    for line in lines[1:]:
+        # Skip comments and newlines
+        if line.startswith("#") or line == "\n":
+            continue
+        key, val = line.strip().split(":")
+        # Strip trailing comments from values
+        source_dict[key.strip()] = val.strip()
+
+    origin = Origin(
+        resource_id=_get_resource_id(source_dict["event name"],
+                                     "origin", tag="source"),
+        time=origin_time, longitude=source_dict["longorUTM"],
+        latitude=source_dict["latorUTM"],
+        depth=abs(float(source_dict["depth"]) * 1E3)  # units: m
+    )
+
+    magnitude = Magnitude(
+        resource_id=_get_resource_id(source_dict["event name"], "magnitude"),
+        mag=ms, magnitude_type="Ms", origin_id=origin.resource_id.id
+    )
+
+    event = Event(resource_id=_get_resource_id(name=source_dict["event name"],
+                                               res_type="event"))
+    event.origins.append(origin)
+    event.magnitudes.append(magnitude)
+
+    event.preferred_origin_id = origin.resource_id.id
+    event.preferred_magnitude_id = magnitude.resource_id.id
+
+    return event
+
+
 def read_specfem2d_source(path_to_source, origin_time=None):
     """
     Create a barebones ObsPy Event object from a SPECFEM2D Source file, which