diff --git a/setup.py b/setup.py index b5bf3cbc36d0a8805e22925f47a4ba2f1e9d9e7f..f6b48ffe7140f4055e0d20fd6e30ca1cf2b7ef6e 100644 --- a/setup.py +++ b/setup.py @@ -13,6 +13,7 @@ setup( install_requires=[ "openquake.engine", "geopandas", + "pyyaml", "Rtree", "rasterio", ], diff --git a/shakyground2/import_fdsnws_eq.py b/shakyground2/import_fdsnws_eq.py new file mode 100644 index 0000000000000000000000000000000000000000..e7e9b8d576ef2544442ce70827f4edad51170012 --- /dev/null +++ b/shakyground2/import_fdsnws_eq.py @@ -0,0 +1,280 @@ +# +# Grab an event's XML file from GEOFON and convert it +# to what we need. +# +# Author: Peter L. Evans +# +# Copyright (C) 2021 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 os +import sys +from typing import List, Union +import urllib.request as ul +from xml.etree import ElementTree as ET +import yaml + + +FDSNWS_ENDPOINT: str = "http://geofon.gfz-potsdam.de/fdsnws/event/1/" +QUAKEML_NS: dict = {"ns": "http://quakeml.org/xmlns/bed/1.2"} + + +def _perror(msg: str) -> None: + print(msg, file=sys.stderr) + + +def fetch_fm(root: ET.Element, ns: dict, preferredfmID: str) -> list: + """ + Extract interesting properties from the focal mechanism with + the given publicID. + + root: ElementTree root element. + ns (dict): its namespace. + preferredfmID (string): the focal mechanism to hunt for. + + Returns a list of one dictionary filled from each of the + and components of the + element. + + """ + 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 ! + if not np: + _perror("Oops: no object seen") + break + + d: List[dict] = [dict(), dict()] + for k in range(2): + plane = np.find("ns:nodalPlane" + str(k + 1), ns) + if plane is None: + continue + for child in plane: + found = child.find("ns:value", ns) + if found is None: + continue + v = found.text + tag = child.tag + tag_idx = tag.index("}") + 1 + tag = tag[tag_idx:] + try: + d[k][tag] = float(str(v)) + except ValueError: + pass + return d + + +def fetch_magnitude( + root: ET.Element, ns: dict, preferredmagnitudeID: str +) -> Union[float, None]: + for m in root[0][0].findall("ns:magnitude", ns): + if m.attrib["publicID"] != preferredmagnitudeID: + continue + + child = m.find("ns:mag", ns) + if child is None: + return None + value = child.find("ns:value", ns) + if value is None: + return None + v = value.text + try: + mag = float(str(v)) + except ValueError: + return None + return mag + + +def fetch_origin(root, ns: dict, preferredoriginID: str) -> dict: + """ + Extract interesting properties from the origin with given publicID. + + root: ElementTree root element. + ns (dict): the XML object's 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_idx = tag.index("}") + 1 + tag = tag[tag_idx:] + + 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 + + # QuakeML depths are in metres. + d["depth"] = d["depth"] / 1000.0 + return d + + +def fetch_quakeml_ws(evid: str) -> str: + """ + Query fdsnws-event web service, and return string + to fetch_quakeml() for parsing. + + evid (string): the event ID to query for. + + """ + url = ( + FDSNWS_ENDPOINT + + "query?" + + "&".join( + ( + "includeallmagnitudes=true", + "includefocalmechanism=true", + "eventid={}".format(evid), + ) + ) + ) + + req = ul.Request(url) + u = ul.urlopen(req) + buf = u.read().decode("utf8") + + print("Got", len(buf), "char(s).") + return buf + + +def fetch_quakeml(path: str) -> Union[dict, None]: + """ + Prepare a dictionary holding 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. + + path (string): Either a local file name, or a GEOFON event ID + such as 'gfz2044sqpr'. + + """ + if os.path.exists(path): + with open(path) as fid: + buf = fid.read() + elif path.startswith("gfz") and len(path) == len("gfz2044sqpr"): + buf = fetch_quakeml_ws(path) + else: + raise IOError("Not a local path or a GEOFON event ID") + + root = ET.fromstring(buf) + ns = QUAKEML_NS + + event = root[0][0] + if not event: + print("Couldn't get an event!") + return None + if "{" + ns["ns"] + "}event" != event.tag: + print("Got a", event.tag, "but expected an event") + return None + try: + # e.g. "smi:org.gfz-potsdam.de/geofon/gfz2021ekhv" + evid = event.attrib["publicID"].split("/")[-1] + except AttributeError: + _perror("Oops, couldn't get event id from " + event.attrib["publicID"]) + + preferredoriginIDElem = root[0][0].find("ns:preferredOriginID", ns) + if preferredoriginIDElem is None: + _perror("Oops, couldn't find the preferred origin ID") + return None + preferredoriginID = preferredoriginIDElem.text + + preferredmagID = None + preferredMagnitudeIDElem = root[0][0].find("ns:preferredMagnitudeID", ns) + if preferredMagnitudeIDElem is not None: + try: + preferredmagID = preferredMagnitudeIDElem.text + except AttributeError: + pass + preferredfmID = None + try: + preferredfmIDElem = root[0][0].find("ns:preferredFocalMechanismID", ns) + if preferredfmIDElem is not None: + preferredfmID = preferredfmIDElem.text + except AttributeError: + pass + + if not preferredoriginID: + print("Oops, no preferredOriginID was found") + return 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) + + focalmech = None + if preferredfmID: + focalmech = fetch_fm(root, ns, preferredfmID) + + # Should this be event's preferred magnitude, + # or the preferred focal mechanism's ? + # Probably they are the same thing... + if not preferredmagID: + print("Oops, no preferredMagnitudeID was found") + return None + 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) + if ev_info is None: + print("Got no dictionary from QuakeML") + outfile = evid + ".yaml" + with open(outfile, "w") as fid: + fid.write(yaml.safe_dump(ev_info, default_flow_style=False)) + print("Wrote to", outfile) diff --git a/tests/data/gfz2021ekhv.xml b/tests/data/gfz2021ekhv.xml new file mode 100644 index 0000000000000000000000000000000000000000..989e2172d72156175076dcbdafdb65a234443174 --- /dev/null +++ b/tests/data/gfz2021ekhv.xml @@ -0,0 +1 @@ +Off E. Coast of N. Island, N.Z.region nameGFZ2021-03-04T13:33:25.96342Z7.705753735e+191.744746227e+19-7.84343467e+196.098688444e+19-1.125507602e+192.367758042e+19-1.275457308e+190.93600977450.81729414130.1827058587GFZ2021-03-04T15:29:24.368901Zsmi:org.gfz-potsdam.de/geofon/Origin/20210304152924.369602.1998993smi:org.gfz-potsdam.de/geofon/Magnitude/20210304152924.369708.1998994smi:org.gfz-potsdam.de/geofon/saul_/gemini-premsmi:org.gfz-potsdam.de/geofon/BP_200s-600s307.177228877.29754332158.586775242.1057512269.1356842213.6104104263.297946124.036887887.311853502e+19355.77250325.529375306-8.046967561e+1997.8975685765.255882317.351140591e+1826.472750830.0639902255GFZ2021-03-04T15:29:24.368901Zsmi:org.gfz-potsdam.de/geofon/Origin/20210304151039.157235.246693310626.47275083GFZ2021-03-04T15:29:24.368901Z7.191210084Mwsmi:org.gfz-potsdam.de/geofon/MT179.36047361.763208609-37.376113891.97801631556418656118601.59272408823.19470215114.48722081.9876948677.20787048manualGFZ2021-03-04T15:56:27.282369Z34801.784525773.602023743.9510824381.956121.41514397horizontal uncertaintysmi:org.gfz-potsdam.de/geofon/LOCSATsmi:org.gfz-potsdam.de/geofon/iasp91confirmedsmi:org.gfz-potsdam.de/geofon/Origin/20210304155627.269946.2475276smi:org.gfz-potsdam.de/geofon/Magnitude/20210304152924.369708.1998994smi:org.gfz-potsdam.de/geofon/FocalMechanism/20210304152924.36933.1998992earthquake diff --git a/tests/test_quakexml.py b/tests/test_quakexml.py new file mode 100644 index 0000000000000000000000000000000000000000..fa12436d03f812711114f14e709cee2bed8699b5 --- /dev/null +++ b/tests/test_quakexml.py @@ -0,0 +1,84 @@ +# +# Author: Peter L. Evans +# +# Copyright (C) 2021 Helmholtz-Zentrum Potsdam - Deutsches GeoForschungsZentrum GFZ +# +# This software is free software and comes with ABSOLUTELY NO WARRANTY. +# +# ---------------------------------------------------------------------- +import os +import unittest +from shakyground2.import_fdsnws_eq import fetch_quakeml + + +# Path for storing the cache file that will be generated +DATA_PATH = os.path.join(os.path.dirname(__file__), "data") + + +class QuakeMLReadTestCase(unittest.TestCase): + """Test on digesting QuakeML from GEOFON. + + Beyond the basic test (read in QuakeML, export YAML), this would + also test the tricker input files: no focal mech available, focal + mech missing, multiple FMs, bad or incomplete XML, non-existent or + unreadable input files, and so on. + + """ + + outfile = "data/testoutput.yaml" + + def test_read_good_file(self): + """ + This GEOFON event has a focal mechanism, and a depth of 34 km. + Check that a value in km, not metres, is produced. + """ + infile = os.path.join(DATA_PATH, "gfz2021ekhv.xml") + result = fetch_quakeml(infile) + self.assertTrue(isinstance(result, dict)) + self.assertEqual(result["id"], "gfz2021ekhv") + self.assertLess(result["origin"]["depth"], 100.0) + + def test_read_good_eventid(self): + """ + Might fail if the GEOFON fdsnws-event webservice is unavailable. + """ + evid = "gfz2021ekhv" + result = fetch_quakeml(evid) + self.assertTrue(isinstance(result, dict)) + self.assertGreater(len(result), 3) + self.assertEqual(result["id"], evid) + + def test_read_no_eventid(self): + """ + The function is given a string which isn't a file and isn't an + event id other. So raise an error. + """ + infile = "/no/such/file" + with self.assertRaises(IOError): + fetch_quakeml(infile) + + def test_read_good_eventid2(self): + """ + Might fail if the GEOFON fdsnws-event webservice is unavailable. + Does the dictionary have the values expected for this event? + """ + evid = "gfz2021htpl" + result = fetch_quakeml(evid) + self.assertTrue(isinstance(result, dict)) + self.assertGreater(len(result), 3) + self.assertEqual(result["date"], "2021-04-21") + self.assertEqual(len(result["focalmechanism"]), 2) + self.assertEqual(result["id"], "gfz2021htpl") + self.assertLess(abs(float(result["magnitude"]) - 4.465), 0.01) + self.assertEqual(len(result["origin"]), 3) + self.assertTrue(result["time"].startswith("08:08:23")) + + def tearDown(self): + try: + os.unlink(self.outfile) + except FileNotFoundError: + pass + + +if __name__ == "__main__": + unittest.main()