Add read support for HypoDD pha files

CHANGELOG.txt

@@ -92,6 +92,7 @@ master:
  - obspy.io
     * added read support for receiver gather format v. 1.6 (see #2070)
     * added read support for FOCMEC 'out' and 'lst' files (see #2156)
+    * added read support for HypoDD 'pha' files (see #2378)
  - obspy.signal.trigger:
     * fix a bug in AR picker (see #2157)
     * option to return Baer-Kradolfer characteristic function from pk_mbaer

misc/docs/source/packages/index.rst

@@ -95,6 +95,7 @@ categories.*
     obspy.io.css
     obspy.io.gcf
     obspy.io.gse2
+    obspy.io.hypodd
     obspy.io.kinemetrics
     obspy.io.mseed
     obspy.io.nied.knet

misc/docs/source/packages/obspy.io.hypodd.rst

@@ -0,0 +1,14 @@
+.. currentmodule:: obspy.io.hypodd
+.. automodule:: obspy.io.hypodd
+
+    .. comment to end block
+
+    Modules
+    -------
+    .. autosummary::
+       :toctree: autogen
+       :nosignatures:
+
+       pha
+
+    .. comment to end block

obspy/io/hypodd/__init__.py

@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+obspy.io.hypodd - HypoDD read support for ObsPy
+===============================================
+
+This module provides read support for the HypoDD PHA phase format.
+
+:copyright:
+    The ObsPy Development Team (devs@obspy.org)
+:license:
+    GNU Lesser General Public License, Version 3
+    (https://www.gnu.org/copyleft/lesser.html)
+
+
+Usage Example
+-------------
+
+The PHA reader hooks into the standard ObsPy event handling
+mechanisms including format autodetection.
+
+>>> from obspy import read_events
+>>> cat = read_events('/path/to/example.pha')
+>>> print(cat)
+2 Event(s) in Catalog:
+2025-05-14T14:35:35.510000Z | +40.225,  +10.450 | 3.5 None
+2025-05-14T15:43:05.280000Z | +40.223,  +10.450 | 1.8 None
+
+"""
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+from future.builtins import *  # NOQA
+
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod(exclude_empty=True)

obspy/io/hypodd/pha.py

@@ -0,0 +1,128 @@
+# -*- coding: utf-8 -*-
+"""
+HypoDD PHA read support.
+
+:copyright:
+    The ObsPy Development Team (devs@obspy.org)
+:license:
+    GNU Lesser General Public License, Version 3
+    (https://www.gnu.org/copyleft/lesser.html)
+"""
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+from future.builtins import *  # NOQA
+import io
+
+from obspy import UTCDateTime
+from obspy.core.event import (
+    Catalog, Event, Origin, Magnitude, Pick, WaveformStreamID, Arrival,
+    OriginQuality)
+
+
+def _seed_id_map(inventory=None, id_map=None, id_default='.{}..{}'):
+    if id_map is None:
+        id_map = {}
+    ret = id_map.copy()
+    if inventory is not None:
+        for net in inventory:
+            for sta in net:
+                if len(sta) == 0:
+                    temp = id_map.get(sta.code, id_default)
+                    temp = temp.split('.', 2)[-1]
+                else:
+                    cha = sta[0]
+                    temp = cha.location_code + '.' + cha.code[:-1] + '{}'
+                ret[sta.code] = net.code + '.{}.' + temp
+    return ret
+
+
+def _block2event(block, seed_map, id_default, ph2comp):
+    """
+    Read HypoDD event block
+    """
+    lines = block.strip().splitlines()
+    yr, mo, dy, hr, mn, sc, la, lo, dp, mg, eh, ez, rms, id_ = lines[0].split()
+    time = UTCDateTime(int(yr), int(mo), int(dy), int(hr), int(mn), float(sc))
+    picks = []
+    arrivals = []
+    for line in lines[1:]:
+        sta, reltime, weight, phase = line.split()
+        comp = ph2comp.get(phase, '')
+        wid = seed_map.get(sta, id_default)
+        _waveform_id = WaveformStreamID(seed_string=wid.format(sta, comp))
+        pick = Pick(waveform_id=_waveform_id, phase_hint=phase,
+                    time=time + float(reltime))
+        arrival = Arrival(phase=phase, pick_id=pick.resource_id,
+                          time_weight=float(weight))
+        picks.append(pick)
+        arrivals.append(arrival)
+    qu = None if rms == '0.0' else OriginQuality(standard_error=float(rms))
+    origin = Origin(arrivals=arrivals, resource_id="smi:local/origin/" + id_
+                    quality=qu,
+                    latitude=float(la),
+                    longitude=float(lo),
+                    depth=1000 * float(dp),
+                    time=time)
+    magnitude = Magnitude(mag=mg, resource_id="smi:local/magnitude/" + id_)
+    event = Event(resource_id="smi:local/event/" + id_,
+                  picks=picks,
+                  origins=[origin],
+                  magnitudes=[magnitude],
+                  preferred_origin_id=origin.resource_id,
+                  preferred_magnitude_id=magnitude.resource_id)
+    return event
+
+
+def _is_pha(filename):
+    try:
+        with open(filename, 'rb') as f:
+            temp = f.readline()
+    except Exception:
+        return False
+    return temp.startswith(b'#') and len(temp.split()) == 15
+
+
+def _read_pha(filename, inventory=None, id_map=None, id_default='.{}..{}',
+              ph2comp={'P': 'Z', 'S': 'N'}, encoding='utf-8'):
+    """
+    Read a HypoDD PHA file and returns an ObsPy Catalog object.
+
+    .. warning::
+        This function should NOT be called directly, it registers via the
+        ObsPy :func:`~obspy.core.event.read_events` function, call this
+        instead.
+
+    The optional parameters all deal with the problem, that the PHA format
+    only stores station names for the picks, but the Pick object expects
+    a SEED id.
+
+    :param str filename: File or file-like object in text mode.
+    :type inventory: :class:`~obspy.core.inventory.inventory.Inventory`
+    :param inventory: Inventory used to retrieve network code, location code
+        and channel code of stations (SEED id).
+    :param dict id_map: If channel information was not found in inventory,
+        it will be looked up in this dictionary
+        (example: `id_map={'MOX': 'GR.{}..HH{}'`).
+        The values must contain three dots and two `{}` which are
+        substituted by station code and component.
+    :param str id_default: Default SEED id expression.
+        The value must contain three dots and two `{}` which are
+        substituted by station code and component.
+    :param dict ph2comp: mapping of phases to components
+        (default: {'P': 'Z', 'S': 'N'})
+    :param str encoding: encoding used (default: utf-8)
+
+    :rtype: :class:`~obspy.core.event.Catalog`
+    :return: An ObsPy Catalog object.
+    """
+    seed_map = _seed_id_map(inventory, id_map, id_default)
+    with io.open(filename, 'r', encoding=encoding) as f:
+        text = f.read()
+    events = [_block2event(block, seed_map, id_default, ph2comp)
+              for block in text.split('#')[1:]]
+    return Catalog(events)
+
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod(exclude_empty=True)

obspy/io/hypodd/tests/__init__.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+from future.builtins import *  # NOQA
+
+import unittest
+
+from obspy.core.util import add_doctests, add_unittests
+
+
+MODULE_NAME = "obspy.io.hypodd"
+
+
+def suite():
+    suite = unittest.TestSuite()
+    add_doctests(suite, MODULE_NAME)
+    add_unittests(suite, MODULE_NAME)
+    return suite
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest='suite')

obspy/io/hypodd/tests/data/example.pha

@@ -0,0 +1,6 @@
+#  2025    5   14   14   35   35.510000   40.2254   10.4496   9.4080   3.50  0.0   0.0   0.0       20202331
+FUR   3.52199909  1.0000  P
+WET   5.86199909  1.0000  S
+#  2025    5   14   15   43    5.280000   40.2227   10.4497   9.7740   1.80  0.0   0.0   0.0       20202429
+UBR   2.46399997  1.0000  P
+WERD  4.16799997  1.0000  X

obspy/io/hypodd/tests/test_pha.py

@@ -0,0 +1,51 @@
+# -*- coding: utf-8 -*-
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+from future.builtins import *  # NOQA
+
+import os
+import unittest
+
+from obspy import read_events, read_inventory
+from obspy.io.hypodd import pha
+
+
+class PHATestCase(unittest.TestCase):
+    """
+    Test suite for obspy.io.hypodd.pha
+    """
+
+    def setUp(self):
+        path = os.path.dirname(__file__)
+        self.fname = os.path.join(path, 'data', 'example.pha')
+
+    def test_is_pha(self):
+        self.assertEqual(pha._is_pha(self.fname), True)
+
+    def test_read_pha(self):
+        cat = read_events(self.fname)
+        self.assertEqual(len(cat), 2)
+
+    def test_populate_waveform_id(self):
+        # read invenroty - FUR with channels, WET without channels
+        inv = read_inventory().select(channel='HH?')
+        inv[0][1].channels = []
+        cat = read_events(self.fname, inventory=inv,
+                          id_default='BLA.{}.11.DH{}',
+                          id_map={'UBR': 'BLB.{}.00.BH{}'})
+        self.assertEqual(len(cat), 2)
+        picks = cat[0].picks + cat[1].picks
+        self.assertEqual(len(picks), 4)
+        waveform_ids = [p.waveform_id.get_seed_string() for p in picks]
+        self.assertIn('GR.FUR..HHZ', waveform_ids)
+        self.assertIn('GR.WET.11.DHN', waveform_ids)
+        self.assertIn('BLB.UBR.00.BHZ', waveform_ids)
+        self.assertIn('BLA.WERD.11.DH', waveform_ids)
+
+
+def suite():
+    return unittest.makeSuite(PHATestCase, 'test')
+
+
+if __name__ == '__main__':
+    unittest.main(defaultTest='suite')

setup.py

@@ -305,6 +305,7 @@
         'IMS10BULLETIN = obspy.io.iaspei.core',
         'EVT = obspy.io.sh.evt',
         'FOCMEC = obspy.io.focmec.core',
+        'PHA = obspy.io.hypodd.pha'
         ],
     'obspy.plugin.event.QUAKEML': [
         'isFormat = obspy.io.quakeml.core:_is_quakeml',
@@ -383,6 +384,10 @@
         'isFormat = obspy.io.focmec.core:_is_focmec',
         'readFormat = obspy.io.focmec.core:_read_focmec',
         ],
+    'obspy.plugin.event.PHA': [
+        'isFormat = obspy.io.hypodd.pha:_is_pha',
+        'readFormat = obspy.io.hypodd.pha:_read_pha',
+        ],
     'obspy.plugin.inventory': [
         'STATIONXML = obspy.io.stationxml.core',
         'INVENTORYXML = obspy.io.arclink.inventory',