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

    License
    Copyright 2020 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/.

    $Id$
'''

def mkGeolocationTable(file=''):
    # imports
    try:
        import xml.etree.cElementTree as ET
    except ImportError:
        import xml.etree.ElementTree as ET
    from sys import stdin
    import warnings
    from geopy.geocoders import Photon
    from geopy.extra.rate_limiter import RateLimiter
    from geopy.exc import GeocoderServiceError
    import json as JSON

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

    def simple_warning(message, category, filename, lineno, file=None, line=None):
        return 'Warning: %s\n' % (message)
    warnings.formatwarning = simple_warning

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

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

    # iterate over all events
    for event in DOM.iterfind('qml:eventParameters/qml:event', namespaces):
        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:
            warnings.warn('Skipping cached event %s' %(publicID))
        elif evaluationMode == 'automatic':
            warnings.warn('Skipping automatic event %s' %(publicID))
        else:
            try:
                locations = reverse_geolocate("{lat:.3f}, {lng:.3f}".format(lat=lat, lng=lng),exactly_one=False,limit=5)
            except GeocoderServiceError:
                warnings.warn('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:
                                pass
            if not place:
                warnings.warn('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)))

# __main__
if __name__ == "__main__":
    # parse arguments
    import argparse
    versionText = '$Revision$ ($Date$, $Author$)'.replace('$', '').replace(':','')
    parser = argparse.ArgumentParser(
        description='Reverse geocoding lookup of events in xml format (stdin).',
        epilog=versionText)
    parser.add_argument('-v', '-V', '--version', action='version',
        version=versionText)
    parser.add_argument('-f', '--file', action='store', dest='file',
        help='read in JSON file containing old output.')
    cla = parser.parse_args()

    # call mkGeolocationTable(...)
    mkGeolocationTable(file=cla.file)