add source utility

Author: wjlei1990

pytomo3d/adjoint/__init__.py

@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+:copyright:
+    Wenjie Lei (lei@princeton.edu), 2016
+:license:
+    GNU General Public License, Version 3
+    (http://www.gnu.org/copyleft/gpl.html)
+"""
+
+from __future__ import (absolute_import, division, print_function)
+
+from .adjsrc import calculate_adjsrc_on_stream  # NOQA

pytomo3d/adjoint/examples/cc_traveltime.adjoint.config.yaml

@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+:copyright:
+    Wenjie Lei (lei@princeton.edu), 2016
+:license:
+    GNU General Public License, Version 3
+    (http://www.gnu.org/copyleft/gpl.html)
+"""
+
+from __future__ import (absolute_import, division, print_function)
+
+from .process import process_stream  # NOQA

pytomo3d/adjoint/examples/multitaper.adjoint.config.yaml

@@ -157,14 +157,14 @@ def interpolate_stream(stream, sampling_rate, starttime=None, npts=None):
     return st_new
 
 
-def process(st, remove_response_flag=False, inventory=None,
-            filter_flag=False, pre_filt=None,
-            starttime=None, endtime=None,
-            resample_flag=False, sampling_rate=1.0,
-            taper_type="hann", taper_percentage=0.05,
-            rotate_flag=False, event_latitude=None,
-            event_longitude=None, sanity_check=False,
-            **kwargs):
+def process_stream(st, remove_response_flag=False, inventory=None,
+                   filter_flag=False, pre_filt=None,
+                   starttime=None, endtime=None,
+                   resample_flag=False, sampling_rate=1.0,
+                   taper_type="hann", taper_percentage=0.05,
+                   rotate_flag=False, event_latitude=None,
+                   event_longitude=None, sanity_check=False,
+                   **kwargs):
     """
     Stream processing function defined for general purpose of tomography.
     The advantage of using Stream, rather than than Trace, is that rotation

pytomo3d/adjoint/examples/waveform.adjoint.config.yaml

@@ -137,13 +137,14 @@ def test_process_obsd():
     pre_filt = [1/90., 1/60., 1/27.0, 1/22.5]
     t1 = event_time
     t2 = event_time + 6000.0
-    st_new = proc.process(st, remove_response_flag=True, inventory=inv,
-                          filter_flag=True, pre_filt=pre_filt,
-                          starttime=t1, endtime=t2, resample_flag=True,
-                          sampling_rate=2.0, taper_type="hann",
-                          taper_percentage=0.05, rotate_flag=True,
-                          event_latitude=event_lat,
-                          event_longitude=event_lon)
+    st_new = proc.process_stream(
+        st, remove_response_flag=True, inventory=inv,
+        filter_flag=True, pre_filt=pre_filt,
+        starttime=t1, endtime=t2, resample_flag=True,
+        sampling_rate=2.0, taper_type="hann",
+        taper_percentage=0.05, rotate_flag=True,
+        event_latitude=event_lat,
+        event_longitude=event_lon)
     bmfile = os.path.join(DATA_DIR, "proc", "IU.KBL.obs.proc.mseed")
     st_compare = obspy.read(bmfile)
     assert compare_stream_kernel(st_new, st_compare)
@@ -161,14 +162,15 @@ def test_process_obsd_2():
     pre_filt = [1/90., 1/60., 1/27.0, 1/22.5]
     t1 = event_time
     t2 = event_time + 6000.0
-    st_new = proc.process(st, remove_response_flag=True, inventory=inv,
-                          filter_flag=True, pre_filt=pre_filt,
-                          starttime=t1, endtime=t2, resample_flag=True,
-                          sampling_rate=2.0, taper_type="hann",
-                          taper_percentage=0.05, rotate_flag=True,
-                          event_latitude=event_lat,
-                          event_longitude=event_lon,
-                          sanity_check=True)
+    st_new = proc.process_stream(
+        st, remove_response_flag=True, inventory=inv,
+        filter_flag=True, pre_filt=pre_filt,
+        starttime=t1, endtime=t2, resample_flag=True,
+        sampling_rate=2.0, taper_type="hann",
+        taper_percentage=0.05, rotate_flag=True,
+        event_latitude=event_lat,
+        event_longitude=event_lon,
+        sanity_check=True)
     bmfile = os.path.join(DATA_DIR, "proc", "IU.KBL.obs.proc.mseed")
     st_compare = obspy.read(bmfile)
     assert len(st_new) == 1
@@ -190,13 +192,14 @@ def test_process_synt():
     pre_filt = [1/90., 1/60., 1/27.0, 1/22.5]
     t1 = event_time
     t2 = event_time + 6000.0
-    st_new = proc.process(st, remove_response_flag=False, inventory=inv,
-                          filter_flag=True, pre_filt=pre_filt,
-                          starttime=t1, endtime=t2, resample_flag=True,
-                          sampling_rate=2.0, taper_type="hann",
-                          taper_percentage=0.05, rotate_flag=True,
-                          event_latitude=event_lat,
-                          event_longitude=event_lon)
+    st_new = proc.process_stream(
+        st, remove_response_flag=False, inventory=inv,
+        filter_flag=True, pre_filt=pre_filt,
+        starttime=t1, endtime=t2, resample_flag=True,
+        sampling_rate=2.0, taper_type="hann",
+        taper_percentage=0.05, rotate_flag=True,
+        event_latitude=event_lat,
+        event_longitude=event_lon)
     bmfile = os.path.join(DATA_DIR, "proc", "IU.KBL.syn.proc.mseed")
     st_compare = obspy.read(bmfile)
     assert compare_stream_kernel(st_new, st_compare)

pytomo3d/signal/__init__.py

@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+:copyright:
+    Wenjie Lei (lei@princeton.edu), 2016
+:license:
+    GNU General Public License, Version 3
+    (http://www.gnu.org/copyleft/gpl.html)
+"""
+
+from __future__ import (absolute_import, division, print_function)
+
+from .append_cmtsolution import append_cmt_to_catalog  # NOQA

pytomo3d/signal/process.py

@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+Functions that append CMTSOLUTION information into catalog
+:copyright:
+    Wenjie Lei (lei@princeton.edu), 2016
+:license:
+    GNU Lesser General Public License, version 3 (LGPLv3)
+    (http://www.gnu.org/licenses/lgpl-3.0.en.html)
+"""
+from __future__ import print_function, division
+import obspy
+from obspy.core.event.source import ResourceIdentifier
+from obspy.core.event import Catalog
+from obspy.core.event import CreationInfo
+
+
+def _validator(event, cmt_origin, cmt_mag, cmt_focal):
+    if event.preferred_origin() != cmt_origin:
+        raise ValueError("preferred_origin_id wrong, not the same as the "
+                         "new added cmt")
+    if event.preferred_magnitude() != cmt_mag:
+        raise ValueError("preferred_magnitude_id wrong, not the same as "
+                         "the new added cmt")
+    if event.preferred_focal_mechanism() != cmt_focal:
+        raise ValueError("preferred_focal_mechanism_id wrong, not the same "
+                         "as the new added cmt")
+
+
+def prepare_cmt_origin(cmt, tag, creation_info):
+    cmt_origin = None
+    for _origin in cmt.origins:
+        if str(_origin.resource_id).endswith("origin#cmt"):
+            cmt_origin = _origin
+            break
+
+    if cmt_origin is None:
+        raise ValueError("No cmt origin found")
+
+    new_id = str(cmt_origin.resource_id).strip() + "#%s" % tag
+    # new_id = str(cmt_origin.resource_id).replace("origin#cmt",
+    #                                             "origin#%s" % tag)
+    cmt_origin.resource_id = ResourceIdentifier(new_id)
+    cmt_origin.creation_info = creation_info
+    return cmt_origin
+
+
+def prepare_cmt_mag(cmt, tag, origin_id, creation_info):
+    cmt_mag = None
+    for _mag in cmt.magnitudes:
+        if _mag.magnitude_type == "mw":
+            cmt_mag = _mag
+
+    if cmt_mag is None:
+        raise ValueError("No cmt Mw mag found")
+
+    new_id = str(cmt_mag.resource_id).strip() + "#%s" % tag
+    cmt_mag.resource_id = ResourceIdentifier(new_id)
+    cmt_mag.origin_id = origin_id
+    cmt_mag.creation_info = creation_info
+    return cmt_mag
+
+
+def prepare_cmt_focal(cmt, tag, origin_id, mag_id, creation_info):
+
+    cmt_focal = None
+    for _focal in cmt.focal_mechanisms:
+        if "cmtsolution" in str(_focal.resource_id):
+            cmt_focal = _focal
+            break
+
+    if cmt_focal is None:
+        raise ValueError("no cmt focal found")
+
+    focal_id = str(cmt_focal.resource_id).strip() + "#%s" % tag
+    cmt_focal.resource_id = ResourceIdentifier(focal_id)
+    cmt_focal.creation_info = creation_info
+    tensor = cmt_focal.moment_tensor
+    tensor_id = str(tensor.resource_id).strip() + "#%s" % tag
+    tensor.resource_id = ResourceIdentifier(tensor_id)
+    tensor.derived_origin_id = origin_id
+    tensor.moment_magnitude_id = mag_id
+    tensor.creation_info = creation_info
+    return cmt_focal
+
+
+def append_cmt_to_catalog(event, cmt, tag, change_preferred_id=True):
+    """
+    Add cmt to event. The cmt.resource_id will be appened tag to avoid
+    tag duplication problem in event.
+    :param event: the event that you want to add cmt in.
+    :type event: str, obspy.core.event.Event or obspy.core.event.Catalog
+    :param cmt: the cmt that you want to add to event.
+    :type event: str, obspy.core.event.Event or obspy.core.event.Catalog
+    :param change_preferred_id: change all preferred_id to the new added cmt
+    :type change_preferred_id: bool
+    """
+    if isinstance(event, str):
+        event = obspy.read_events(event)[0]
+    elif isinstance(event, Catalog):
+        event = event[0]
+
+    if isinstance(cmt, str):
+        cmt_event = obspy.read_events(cmt)[0]
+    elif isinstance(cmt, Catalog):
+        cmt_event = cmt_event[0]
+
+    if not isinstance(tag, str):
+        raise TypeError("tag(%s) should be type of str" % type(tag))
+
+    creation_info = CreationInfo(author="GATG", version=tag)
+
+    cmt_origin = prepare_cmt_origin(cmt_event, tag, creation_info)
+    cmt_mag = prepare_cmt_mag(cmt_event, tag, cmt_origin.resource_id,
+                              creation_info)
+    cmt_focal = prepare_cmt_focal(cmt_event, tag, cmt_origin.resource_id,
+                                  cmt_mag.resource_id, creation_info)
+    event.origins.append(cmt_origin)
+    event.magnitudes.append(cmt_mag)
+    event.focal_mechanisms.append(cmt_focal)
+
+    if change_preferred_id:
+        event.preferred_origin_id = str(cmt_origin.resource_id)
+        event.preferred_magnitude_id = str(cmt_mag.resource_id)
+        event.preferred_focal_mechanism_id = str(cmt_focal.resource_id)
+        _validator(event, cmt_origin, cmt_mag, cmt_focal)
+
+    new_cat = Catalog()
+    new_cat.append(event)
+
+    return new_cat

pytomo3d/signal/tests/test_signal.py

@@ -0,0 +1,30 @@
+import os
+import inspect
+import pytomo3d.source.append_cmtsolution as app_cmt
+
+
+def _upper_level(path, nlevel=4):
+    """
+    Go the nlevel dir up
+    """
+    for i in range(nlevel):
+        path = os.path.dirname(path)
+    return path
+
+# Most generic way to get the data folder path.
+TESTBASE_DIR = _upper_level(os.path.abspath(
+    inspect.getfile(inspect.currentframe())), 4)
+DATA_DIR = os.path.join(TESTBASE_DIR, "tests", "data")
+
+testquakeml = os.path.join(DATA_DIR, "quakeml", "C201009031635A.xml")
+testcmt = os.path.join(DATA_DIR, "quakeml", "C201009031635A.inv")
+
+
+def test_append_cmt_to_catalog():
+    tag = "GATG_M15"
+    new_cat = app_cmt.append_cmt_to_catalog(
+        testquakeml, testcmt, tag, change_preferred_id=True)
+
+    assert new_cat[0].preferred_origin()
+    assert new_cat[0].preferred_magnitude()
+    assert new_cat[0].preferred_focal_mechanism()