seisobs-webapp/scripts/mkGeolocationTable.py

#! /usr/bin/env python3
# -*- coding: utf-8 -*-

'''
    Script to lookup city names of events with Nominatim service

    The input should be an valid quakeML file passed to stdin.
    The output will will be a javascript structure to be included in the
    SeisObs map service.

    The script should be updated regularly keep the total number of all
    AJAX calls to the Nominatim service small, e. g. :
    curl -s "https://fdsnws.geophysik.ruhr-uni-bochum.de/fdsnws/event/1/query?minlat=50&maxlat=54&minlon=3&maxlon=10&minmag=1" | mkGeolocationTable.py > geolocationTable.js

    Version: v1.3 (2023-04-17)

    License
    Copyright 2020-2021 Kasper D. Fischer <kasper.fischer@rub.de>

    This program is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the Free
    Software Foundation, either version 3 of the License, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License along
    with this program. If not, see http://www.gnu.org/licenses/.
'''

def mkGeolocationTable(file=''):
    ## imports

    # XML ElementTree
    try:
        import xml.etree.cElementTree as ET
    except ImportError:
        import xml.etree.ElementTree as ET

    # json
    import json as JSON

    # sys stdin
    from sys import stdin

    # geopy
    from geopy.geocoders import Photon
    from geopy.extra.rate_limiter import RateLimiter
    from geopy.exc import GeocoderServiceError

    ## constants
    URL = 'https://photon.komoot.io/reverse?lon={lng:.3f}&lat={lat:.3f}&limit=10'
    NAMESPACES = {'sc3': 'http://geofon.gfz-potsdam.de/ns/seiscomp3-schema/0.7',
        'qml': 'http://quakeml.org/xmlns/bed/1.2'}

    # try loading the file
    geolocationTable = {}
    if file :
        try:
            jsonfile = open(file)
            jsonfileContent = jsonfile.read().split('=')[1].replace(';', '')
            geolocationTable = JSON.loads(jsonfileContent)
        except:
            geolocationTable = {}
            logging.warning('Could not parse file %s' %file)

    # parse event.xml
    DOM = ET.parse(stdin).getroot()
    geolocator = Photon()
    reverse_geolocate = RateLimiter(geolocator.reverse, min_delay_seconds=2)

    # iterate over all events
    count = 0
    for event in DOM.iterfind('qml:eventParameters/qml:event', NAMESPACES):
        count += 1
        publicID = event.attrib['publicID'].split('/')[2]
        lat = float(event.find('./qml:origin/qml:latitude/qml:value', NAMESPACES).text)
        lng = float(event.find('./qml:origin/qml:longitude/qml:value', NAMESPACES).text)
        evaluationMode = event.find('./qml:origin/qml:evaluationMode', NAMESPACES).text

        if publicID in geolocationTable:
            logging.warning('Skipping cached event {id}'.format(id=publicID))
        elif evaluationMode == 'automatic':
            logging.warning('Skipping automatic event {id}'.format(id=publicID))
        else:
            logging.info('Processing event {id}'.format(id=publicID))
            try:
                locations = reverse_geolocate("{lat:.3f}, {lng:.3f}".format(lat=lat, lng=lng),exactly_one=False,limit=10)
            except GeocoderServiceError:
                logging.warning('Reverse Geolocation failed. Skipping event.')
                continue
            place = []
            for location in locations:
                try:
                    place = location.raw['properties']['city']
                except KeyError:
                    try:
                        place = location.raw['properties']['town']
                    except KeyError:
                        try:
                            place = location.raw['properties']['village']
                        except KeyError:
                            try:
                                place = location.raw['properties']['county']
                            except KeyError:
                                logging.debug('Could not extract city for event {id} at {lat:.3f} N / {lng:.3f} E (Service: {url}), trying next returned result'
                                    .format(id=publicID, lat=lat, lng=lng, url=URL.format(lat=lat,lng=lng)))
                                logging.debug(location.raw)
            if not place:
                logging.critical('Could not extract city for event {id} at {lat:.3f} N / {lng:.3f} E (Service: {url})'
                    .format(id=publicID, lat=lat, lng=lng, url=URL.format(lat=lat,lng=lng)))
            geolocationTable[publicID] = place

    # dump json
    print('var geolocationTable = {0};'.format(JSON.dumps(geolocationTable, sort_keys=True)))
    logging.info("processed %d events", count)

# __main__
if __name__ == "__main__":
    # use module logging
    import logging

    # parse arguments
    import argparse
    versionText = 'v1.3 (2023-04-17)'
    parser = argparse.ArgumentParser(
        description='Reverse geocoding lookup of events in xml format (stdin).',
        epilog=versionText)
    parser.add_argument('-V', '--version', action='version', version=versionText,
        help="show version")
    parser.add_argument('-f', '--file', action='store', dest='file',
        help='read in JSON file containing old output.')
    group = parser.add_mutually_exclusive_group()
    group.add_argument("-v", "--verbose", action="store_true",
        help="increase output verbosity")
    group.add_argument("-q", "--quiet", action="store_true",
        help="disable output")
    group.add_argument("-d", "--debug", action="store_true",
        help="show debugging output",
                    )
    cla = parser.parse_args()

    # set logging level
    if cla.quiet:
       logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.ERROR)
    elif cla.verbose:
       logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
    elif cla.debug:
       logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.DEBUG)
    else:
       logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.WARN)
    
    # call mkGeolocationTable(...)
    mkGeolocationTable(file=cla.file)