source: SH_SHM/trunk/util/quakeml/emsc.py @ 318

#! /usr/bin/python
# -*- coding: UTF8 -*-

# requires correctly set PYTHONPATH
try:
    from SzoDc.quakeml.parse import QuakeML
except ImportError:
    from parse import QuakeML

import urllib2 as urllib
from datetime import datetime, date, timedelta
import sys
from StringIO import StringIO
import decimal

DATEFORMAT = "%Y-%m-%d"

class Emsc(object):
    """
    Handles ESMC data request via seismicportal.eu
    """

    BASEURI = "http://www.seismicportal.eu/services/event/search?"

    def __init__(self, **kwargs):
        """
        Query web service
        """

        # dict of events
        self.events = {}
        # list of all origins
        self.origins = []

        # build query string
        agencies = kwargs.get("agency", None)
        if agencies:
            agencies = agencies.upper().split(",")

            # If only one agency is given, include in query. Otherwise it will be
            # filtered later.
            if len(agencies) == 1:
                kwargs["auth"] = kwargs["agency"]
                agencies = None
        self.agencies = agencies
        
        # Also supply default values here, since they might not defined yet.
        self.out = kwargs.get("outfile", sys.stdout)
        self.debugLevel = kwargs.get("debug", False)
        self.nomagnitude = kwargs.get("nomagnitude", False)
        kwargs["limit"] = kwargs.get("limit", "2000")

        if self.debugLevel:
            print "Query EMSC from %s to %s (agency: %s)" % (
                        kwargs["dateMin"], kwargs["dateMax"], kwargs["agency"])

        try:
            self.out.closed
        except:
            self.out = open(self.out, "w")

        # clean up kwargs
        discard = ["agency", "outfile", "debug", "nomagnitude",]
        for i in discard:
            if kwargs.has_key(i):
                del kwargs[i]

        # build uri
        uri = self.BASEURI + "&".join(map(lambda x: "=".join(x), \
                                           zip(kwargs.keys(), kwargs.values())))

        self.__load(uri)

    def __load(self, uri):
        """
        Load data from ESMC web service.
        """

#        data = StringIO(open("example.xml").read())
        data = StringIO(urllib.urlopen(uri).read())

        # Some data handling
        for e in QuakeML(data, self.debugLevel).events:
            # two decimal places for locations
            e.longitude = e.longitude.quantize(decimal.Decimal('.01'), \
                                                    rounding=decimal.ROUND_DOWN)
            e.latitude = e.latitude.quantize(decimal.Decimal('.01'), \
                                                    rounding=decimal.ROUND_DOWN)

            # add south/north/east/west indicator (separates by two spaces)
            if e.longitude < 0:
                e.longitudeEW = str(-1 * e.longitude) + "  W"
            else:
                e.longitudeEW = str(e.longitude) + "  E"

            if e.latitude < 0:
                e.latitudeNS = str(-1*e.latitude) + "  S"
            else:
                e.latitudeNS = str(e.latitude) + "  N"

            e.agencyOrigin = e.agencyOrigin.split("/")[-1]

            # Check for multiple agency criteria.
            if self.agencies and e.agencyOrigin not in self.agencies:
                continue

            # some localisations do not include a magnitude
            try:
                e.agencyMagnitude = e.agencyMagnitude.split("/")[-1]
            except AttributeError:
                # skip output if option not set
                if not self.nomagnitude:
                    continue

                e.agencyMagnitude = "NOTSET"
                e.magnitude = Decimal("NaN")
                e.magnitudeType = "no"

            self.origins.append(e)
            try:
                self.events[e.eventID].append(e)
            except KeyError:
                self.events[e.eventID] = [e]

    def Ascii(self):
        """
        Ascii output of data.
        """
        for e in self.origins:
            print >> self.out, "  ".join(e.origin.split("T"))[:-6],
            print >> self.out, e.latitudeNS, " ", e.longitudeEW, " ", e.depth, 
            print >> self.out, " ", e.magnitude, e.magnitudeType,
            print >> self.out, e.mode[0].upper(), e.description, " ",
            print >> self.out, e.agencyOrigin

    def Oed(self):
        """
        Returns origin-epicentre-depth list. There will be only one origin per
        event. If no trusted agencies are set, the first origin found is used,
        otherwise the first origin of a trusted agencies is returned.

        In any case manual solutions are preferred against automatic ones.
        """
        oed = []

        for eventID in sorted(self.events.keys()):
            # classify by mode (manual, automatic) solutions
            evts = {}
            for e in self.events[eventID]:
                try:
                    evts[e.mode].append(e)
                except KeyError:
                    evts[e.mode] = [e]

            # look into manual locations
            solution = evts.get("manual", None)

            # only automatic solution found:
            if not solution:
                solution = evts["automatic"]

            # reformat to meet oed format
            try:
                res = (
                    datetime.strptime(solution[0].origin, \
                                                       "%Y-%m-%dT%H:%M:%S%f%z"),
                    float(solution[0].latitude),
                    float(solution[0].longitude),
                    float(solution[0].depth),
                    float(solution[0].magnitude),
                )
            except ValueError:
                # %f and %z are only supported from python 2.6+
                res = (
                    datetime.strptime(solution[0].origin[:-8], \
                                                           "%Y-%m-%dT%H:%M:%S"),
                    float(solution[0].latitude),
                    float(solution[0].longitude),
                    float(solution[0].depth),
                    float(solution[0].magnitude),
                )

            oed.append(res)

        return oed

def main():
    from optparse import OptionParser
    usage = """
    %prog [options] [dateMin] [dateMax]

    If dateMin/dateMax is not given a query for the last 5 days is issued. Date
    must be formatted as YYYY-MM-DD (e.g. 2010-06-24).

    If you enter only one datum, the query will cover data from this date to
    today. If you want to get all events from today, simply enter today's date."""

    parser = OptionParser(usage=usage)
    parser.add_option("-a", "--agency", dest="agency",
                  help="Query only for given agency (e.g. NEIR). Separate "
                       "multiple agencies by comma: NEIR,GSRC.")
    parser.add_option("-o", "--outfile", dest="outfile", default=sys.stdout,
                  help="Write result to FILE [default stdout]", metavar="FILE")
    parser.add_option("-l", "--limit", dest="limit", default='2000',
                  help="Maximum number of events returned [default 2.000]. "
                       "Please note: There might be multiple origin locations "
                       "per event.")
    parser.add_option("-d", "--debug", action="count", dest="debug",
                  default=False,
                  help="Switch on debugging [default off]. If this parameter "
                       "is used twice ('-dd' or '-d -d') the input XML data "
                       "will be dumped to 'debug-dump.xml' in the current "
                       "directory.")
    parser.add_option("--include-no-magnitude", action="store_true",
                  dest="nomagnitude", default=False,
                  help="Also output localisations that include no magnitude "
                       "information [default off]. If this option is set, these "
                       "events will have 'NaN' as magnitude and 'no' as type of "
                       "magnitude.")

    (options, args) = parser.parse_args()
    query = vars(options)

    query["dateMin"] = None
    query["dateMax"] = None

    if len(args) == 0:
        # no arguments: today and 5 days back
        query["dateMax"] = date.today().strftime(DATEFORMAT)
        query["dateMin"] = (date.today() - timedelta(5)).strftime(DATEFORMAT)
    elif len(args) == 1:
        # one argument: from date till today
        query["dateMin"] = args[0]
        query["dateMax"] = date.today().strftime(DATEFORMAT)
    elif len(args) == 2:
        # timespan given: sort it
        if args[0] > args[1]:
            query["dateMin"] = args[1]
            query["dateMax"] = args[0]
        else:
            query["dateMin"] = args[0]
            query["dateMax"] = args[1]
    else:
        parser.error("Only two arguments can be processed.")

    # add time information to query the hole day
    query["dateMin"] += "T00:00:00"
    query["dateMax"] += "T23:59:59"

    _ = Emsc(**query).Ascii()

if __name__ == "__main__":
    main()