feat: create catalog.to_quakeml function including necessary changes …

…in other direction on quakeml parser
swiss-seismological-service · Sep 22, 2023 · f4dac0f · f4dac0f
1 parent 4e8f228
commit f4dac0f
Show file tree

Hide file tree

Showing 6 changed files with 156 additions and 13 deletions.
diff --git a/catalog_tools/io/parser.py b/catalog_tools/io/parser.py
@@ -18,25 +18,27 @@ def _get_realvalue(key: str, value: str) -> dict:
     **_get_realvalue('origintime', 'time'),
     **_get_realvalue('originlatitude', 'latitude'),
     **_get_realvalue('originlongitude', 'longitude'),
-    **_get_realvalue('origindepth', 'depth')
+    **_get_realvalue('origindepth', 'depth'),
+    'originevaluationMode': 'evaluationmode',
+    'originpublicID': 'originid',
 }
 
 MAGNITUDE_MAPPINGS = {
     **_get_realvalue('magnitudemag', 'magnitude'),
     'magnitudetype': 'magnitude_type',
-    'magnitudeevaluationMode': 'evaluationMode',
+    'magnitudepublicID': 'magnitudeid',
 }
 
 DUMMY_MAGNITUDE = {
     'magnitudemagvalue': None,
-    'magnitudetype': None,
-    'magnitudeevaluationMode': None}
+    'magnitudetype': None}
 
 DUMMY_ORIGIN = {
     'origintimevalue': None,
     'originlatitudevalue': None,
     'originlongitudevalue': None,
-    'origindepthvalue': None
+    'origindepthvalue': None,
+    'originevaluationMode': None
 }
 
 
@@ -114,8 +116,10 @@ def _extract_magnitude(magnitude: dict) -> dict:
 def _extract_secondary_magnitudes(magnitudes: list) -> dict:
     magnitude_dict = {}
     for magnitude in magnitudes:
-        mappings = _get_realvalue(
-            'magnitudemag', f'magnitude_{magnitude["magnitudetype"]}')
+        mappings = {**_get_realvalue(
+            'magnitudemag', f'magnitude_{magnitude["magnitudetype"]}'),
+            'magnitudepublicID':
+            f'magnitude_{magnitude["magnitudetype"]}_magnitudeid'}
         for key, value in mappings.items():
             if key in magnitude:
                 magnitude_dict[value] = magnitude[key]

diff --git a/catalog_tools/seismicity/catalog.py b/catalog_tools/seismicity/catalog.py
@@ -1,13 +1,21 @@
 from __future__ import annotations
 
+import os
+import uuid
+from collections import defaultdict
+
 import pandas as pd
 
-from catalog_tools.utils import _check_required_cols, require_cols
+from catalog_tools.utils import (_check_required_cols, _render_template,
+                                 require_cols)
 from catalog_tools.utils.binning import bin_to_precision
 
 REQUIRED_COLS_CATALOG = ['longitude', 'latitude', 'depth',
                          'time', 'magnitude']
 
+QML_TEMPLATE = os.path.join(os.path.dirname(os.path.abspath(__file__)),
+                            'catalog_templates', 'quakeml.j2')
+
 
 def _catalog_constructor_with_fallback(*args, **kwargs):
     df = Catalog(*args, **kwargs)
@@ -95,8 +103,43 @@ def bin_magnitudes(
             return df
 
     @require_cols(require=_required_cols)
-    def to_quakeml(self) -> str:
-        raise NotImplementedError
+    def to_quakeml(self, agencyID=' ', author=' ') -> str:
+        df = self.copy()
+        if 'eventid' not in df.columns:
+            df['eventid'] = uuid.uuid4()
+        if 'originid' not in df.columns:
+            df['originid'] = uuid.uuid4()
+        if 'magnitudeid' not in df.columns:
+            df['magnitudeid'] = uuid.uuid4()
+
+        vals = ['_uncertainty',
+                '_lowerUncertainty',
+                '_upperUncertainty',
+                '_confidenceLevel']
+
+        secondary_mags = [mag for mag in df.columns if
+                          'magnitude_' in mag
+                          and not mag == 'magnitude_type'
+                          and not any(['magnitude' + val
+                                       in mag for val in vals])]
+
+        data = dict(events=df.to_dict(orient='records'),
+                    agencyID=agencyID, author=author)
+
+        for event in data['events']:
+            event['sec_mags'] = defaultdict(dict)
+            for mag in secondary_mags:
+                if pd.notna(event[mag]) \
+                        and event['magnitude_type'] not in mag:
+
+                    mag_type = mag.split('_')[1]
+                    mag_key = mag.replace('_' + mag_type, '')
+
+                    event['sec_mags'][mag_type][mag_key] = \
+                        event[mag]
+                del event[mag]
+        print(data)
+        return _render_template(data, QML_TEMPLATE)
 
 
 class ForecastCatalog(Catalog):

diff --git a/catalog_tools/seismicity/catalog_templates/quakeml.j2 b/catalog_tools/seismicity/catalog_templates/quakeml.j2
@@ -0,0 +1,69 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<q:quakeml xmlns="http://quakeml.org/xmlns/bed/1.2" xmlns:q="http://quakeml.org/xmlns/quakeml/1.2">
+    <eventParameters publicID="smi:org.gfz-potsdam.de/geofon/EventParameters"> 
+        {% for event in events %} 
+        <event publicID="{{ event.eventid }}">
+            <creationInfo>
+                <agencyID>{{ agencyID }}</agencyID>
+                <author>{{ author }}</author>
+            </creationInfo>
+
+            <magnitude publicID="{{ event.magnitudeid }}">
+                <creationInfo>
+                    <agencyID>{{ agencyID }}</agencyID>
+                    <author>{{ author }}</author>
+                </creationInfo>
+                <mag>
+                    <value>{{ event.magnitude }}</value>
+                    <uncertainty>{{ event.magnitude_uncertainty }}</uncertainty>
+                </mag>
+                <type>{{ event.magnitude_type }}</type>
+                <originID>{{ event.originid }}</originID>
+            </magnitude>
+
+            {% for type, magnitude in event.sec_mags.items() %}
+            <magnitude publicID="{{ magnitude.magnitude_magnitudeid }}">
+                <creationInfo>
+                    <agencyID>{{ agencyID }}</agencyID>
+                    <author>{{ author }}</author>
+                </creationInfo>
+                <mag>
+                    <value>{{ magnitude.magnitude }}</value>
+                    <uncertainty>{{ magnitude.magnitude_uncertainty }}</uncertainty>
+                </mag>
+                <type>{{ type }}</type>
+                <originID>{{ event.originid }}</originID>
+            </magnitude>
+            {% endfor %}
+
+            <origin publicID="{{ event.originid }}">
+                <time>
+                    <value>{{ event.time }}</value>
+                </time>
+                <longitude>
+                    <value>{{ event.longitude }}</value>
+                    <uncertainty> {{ event.longitude_uncertainty }}</uncertainty>
+                </longitude>
+                <latitude>
+                    <value>{{ event.latitude }}</value>
+                    <uncertainty>{{ event.latitude_uncertainty }}</uncertainty>
+                </latitude>
+                <evaluationMode>{{ event.evaluationmode }}</evaluationMode>
+                <creationInfo>
+                    <agencyID>{{ agencyID }}</agencyID>
+                    <author>{{ author }}</author>
+                </creationInfo>
+                <depth>
+                    <value>{{ event.depth }}</value>
+                    <uncertainty>{{ event.depth_uncertainty }}</uncertainty>
+                </depth>
+            </origin>
+
+            <preferredOriginID>{{ event.originid }}</preferredOriginID>
+            <preferredMagnitudeID>{{ event.magnitudeid }}</preferredMagnitudeID>
+            <type>{{ event.event_type }}</type>
+
+        </event> 
+        {% endfor %} 
+        </eventParameters>
+</q:quakeml>
diff --git a/catalog_tools/seismicity/tests/test_catalog.py b/catalog_tools/seismicity/tests/test_catalog.py
@@ -23,6 +23,7 @@
 CATALOG_TEST_DATA = [
     {'depth': '1181.640625',
      'depth_uncertainty': '274.9552879',
+     'evaluationmode': 'manual',
      'event_type': 'earthquake',
      'eventid': 'smi:ch.ethz.sed/sc20a/Event/2021zqxyri',
      'latitude': '46.05144527',
@@ -31,14 +32,23 @@
      'longitude_uncertainty': '0.1007121534',
      'magnitude': '2.510115344',
      'magnitude_MLhc': '2.510115344',
+     'magnitude_MLhc_magnitudeid':
+     'smi:ch.ethz.sed/sc20ag/Magnitude/20220103070310.700951.80206',
      'magnitude_MLhc_uncertainty': '0.23854491',
      'magnitude_MLv': '2.301758471',
+     'magnitude_MLv_magnitudeid':
+     'smi:ch.ethz.sed/sc20ag/Magnitude/20220103070310.752473.80241',
      'magnitude_MLv_uncertainty': '0.2729312832',
      'magnitude_type': 'MLhc',
      'magnitude_uncertainty': '0.23854491',
+     'magnitudeid':
+     'smi:ch.ethz.sed/sc20ag/Magnitude/20220103070310.700951.80206',
+     'originid':
+     'smi:ch.ethz.sed/sc20ag/Origin/NLL.20220103070248.816904.80080',
      'time': '2021-12-30T07:43:14.681975Z'},
     {'depth': '3364.257812',
      'depth_uncertainty': '1036.395075',
+     'evaluationmode': 'manual',
      'event_type': 'earthquake',
      'eventid': 'smi:ch.ethz.sed/sc20a/Event/2021zihlix',
      'latitude': '47.37175484',
@@ -47,9 +57,15 @@
      'longitude_uncertainty': '0.1277685645',
      'magnitude': '3.539687307',
      'magnitude_MLhc': '3.539687307',
+     'magnitude_MLhc_magnitudeid':
+     'smi:ch.ethz.sed/sc20ag/Magnitude/20211228194308.87278.210164',
      'magnitude_MLhc_uncertainty': '0.272435385',
      'magnitude_type': 'MLhc',
      'magnitude_uncertainty': '0.272435385',
+     'magnitudeid':
+     'smi:ch.ethz.sed/sc20ag/Magnitude/20211228194308.87278.210164',
+     'originid':
+     'smi:ch.ethz.sed/sc20ag/Origin/NLL.20211228194249.917108.210045',
      'time': '2021-12-25T14:49:40.125942Z'}]
 
 PATH_RESOURCES = os.path.join(os.path.dirname(os.path.abspath(__file__)),
@@ -127,11 +143,12 @@ def test_to_quakeml():
     xml_file = os.path.join(PATH_RESOURCES, 'quakeml_data.xml')
 
     catalog = Catalog(CATALOG_TEST_DATA)
-    catalog_xml = catalog.to_quakeml()
-    catalog_xml = re.sub(r"[\n\t\s]*", "", catalog_xml)
+    catalog_xml = catalog.to_quakeml(agencyID='SED', author='catalog-tools')
 
+    catalog_xml = re.sub(r"[\n\t\s]*", "", catalog_xml)
+    print('\n', catalog_xml)
     with open(xml_file, 'r') as file:
         xml = file.read()
     xml = re.sub(r"[\n\t\s]*", "", xml)
-
+    print('\n', xml)
     assert catalog_xml == xml
diff --git a/catalog_tools/utils/__init__.py b/catalog_tools/utils/__init__.py
@@ -1,6 +1,7 @@
 import functools
 
 import pandas as pd
+from jinja2 import Template, select_autoescape
 
 
 def _check_required_cols(df: pd.DataFrame,
@@ -57,3 +58,11 @@ def wrapper_require(self, *args, **kwargs):
         return decorator_require
     else:
         return decorator_require(_func)
+
+
+def _render_template(data: dict, template_path: str) -> str:
+    with open(template_path) as t:
+        template = Template(t.read(), autoescape=select_autoescape())
+
+    qml = template.render(**data)
+    return qml
diff --git a/setup.cfg b/setup.cfg
@@ -10,6 +10,7 @@ packages = find:
 install_requires =
     cartopy
     geopandas
+    jinja2
     matplotlib
     numpy
     pandas