Commit d4652715 authored by Peter Evans's avatar Peter Evans
Browse files

import_fdsnws_eq: Get the XML, dump some YAML

Prepares a dictionary with the interesting things found
by reading the QuakeML file served for a given event.

- If no focal mechanism is present, return None in
  those components of the dictionary; YAML serialises
  these as 'null'.
- Parse FM XML object for all six nodal plane angles.
- Catch float conversion problems by omitting output of these.
- Check event id looks okay.
- If called from command line, report outfile name.
parent e61ff9d3
#
# Grab an event's XML file from GEOFON and convert it
# to what we need.
#
# Author: Peter L. Evans <pevans@gfz-potsdam.de>
#
# Copyright (C) 2018 Helmholtz-Zentrum Potsdam - Deutsches GeoForschungsZentrum GFZ
#
# This software is free software and comes with ABSOLUTELY NO WARRANTY.
#
# ----------------------------------------------------------------------
"""
Can run standalone:
>>> python import_fdsnws_eq.py gfz2021ekhv
"""
import datetime
import sys
import urllib.request as ul
from xml.etree import ElementTree as ET
import yaml
FDSNWS_ENDPOINT = "http://geofon.gfz-potsdam.de/fdsnws/event/1/"
QUAKEML_NS = {"ns": "http://quakeml.org/xmlns/bed/1.2"}
def fetch_fm(root, ns, preferredfmID):
"""
Extract interesting properties from the focal mechanism with
the given publicID.
root - ElementTree root element.
ns - string, its namespace.
preferredfmID - string, the focal mechanism to hunt for.
"""
for fm in root[0][0].findall("ns:focalMechanism", ns):
if fm.attrib["publicID"] != preferredfmID:
continue
# We've found our guy!
np = fm.find("ns:nodalPlanes", ns)
# Expect there's only one <nodalPlanes>!
d = [None, None]
for k in range(2):
plane = np.find("ns:nodalPlane" + str(k + 1), ns)
d[k] = dict()
for child in plane:
v = child.find("ns:value", ns).text
tag = child.tag
tag = tag[tag.index("}") + 1 :]
try:
d[k][tag] = float(v)
except ValueError:
pass
return d
def fetch_magnitude(root, ns, preferredmagnitudeID):
for m in root[0][0].findall("ns:magnitude", ns):
if m.attrib["publicID"] != preferredmagnitudeID:
continue
child = m.find("ns:mag", ns)
if not child:
return None
v = child.find("ns:value", ns).text
try:
mag = float(v)
except ValueError:
mag = None
return mag
def fetch_origin(root, ns, preferredoriginID):
"""
Extract interesting properties from the origin with given publicID.
root - ElementTree root element.
ns - string, its namespace.
preferredoriginID - string, the origin to hunt for.
"""
for o in root[0][0].findall("ns:origin", ns):
if o.attrib["publicID"] != preferredoriginID:
continue
d = dict()
for child in o:
tag = child.tag
tag = tag[tag.index("}") + 1 :]
if tag in ("depth", "latitude", "longitude", "time"):
v = child.find("ns:value", ns).text
if tag == "time":
d[tag] = v
else:
try:
d[tag] = float(v)
except ValueError:
pass
return d
def fetch_quakeml(evid):
"""
Query fdsnws-event web service, and prepare a dictionary with the
interesting things found by reading the QuakeML file served for a given event.
If there is no focal mechanism, return None where that would otherwise go.
evid - string, the event ID to query for.
"""
url = (
FDSNWS_ENDPOINT
+ "query?"
+ "&".join(
(
"includeallmagnitudes=true",
"includefocalmechanism=true",
"eventid={}".format(evid),
)
)
)
# print(url)
req = ul.Request(url)
u = ul.urlopen(req)
buf = u.read().decode("utf8")
print("Got", len(buf), "char(s).")
with open("tmp.xml", "w") as fid:
fid.write(str(buf))
root = ET.fromstring(buf)
ns = QUAKEML_NS
preferredoriginID = root[0][0].find("ns:preferredOriginID", ns).text
preferredmagID = root[0][0].find("ns:preferredMagnitudeID", ns).text
try:
preferredfmID = root[0][0].find("ns:preferredFocalMechanismID", ns).text
except AttributeError:
preferredfmID = None
origin = fetch_origin(root, ns, preferredoriginID)
(d, t) = origin.pop("time").split("T", 2)
# There's little point in forcing these into datetimes,
# since writing to YAML requires they be converted back
# to strings. :(
# d = datetime.date(2020, 1, 1)
# t = datetime.time(12, 0, 30, 123000)
if preferredfmID:
focalmech = fetch_fm(root, ns, preferredfmID)
else:
focalmech = None
# Should this be event's preferred magnitude,
# or the preferred focal mechanism's <momentMagnitudeID>?
# Probably they are the same thing...
mag = fetch_magnitude(root, ns, preferredmagID)
return {
"id": evid,
"date": str(d),
"time": str(t),
"origin": origin,
"magnitude": mag,
"focalmechanism": focalmech,
"preferred_origin_id": preferredoriginID,
"preferred_magnitude_id": preferredmagID,
"preferred_focalmechanism_id": preferredfmID,
}
if __name__ == "__main__":
try:
evid = sys.argv[1]
except IndexError:
print("Usage:", sys.argv[0], "EVENTID")
print(
"""
where EVENTID is a GEOFON event ID. These are listed at
http://geofon.gfz-potsdam.de/eqinfo/list.php
"""
)
sys.exit(1)
if not evid.startswith("gfz20"):
print("GEOFON event IDs are generally of the form 'gfz20xxnnnn'. Are you sure?")
sys.exit(1)
ev_info = fetch_quakeml(evid)
outfile = evid + ".yaml"
with open(outfile, "w") as fid:
fid.write(yaml.safe_dump(ev_info, default_flow_style=False))
print("Wrote to", outfile)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment