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

Revision 318, 8.8 KB checked in by marcus, 12 years ago (diff)

r157 | walther | 2011-02-10 12:46:56 +0100 (Do, 10 Feb 2011) | 4 lines

Subversion export from BGR/SZO/SZGRF internal used utilities at revision

  1. Since this software is only used in Seismic Handler context,

future enhancements will be done here ("the software got a new home").

  • Property svn:eol-style set to native
  • Property svn:executable set to *
Line 
1#! /usr/bin/python
2# -*- coding: UTF8 -*-
3
4# requires correctly set PYTHONPATH
5try:
6    from SzoDc.quakeml.parse import QuakeML
7except ImportError:
8    from parse import QuakeML
9
10import urllib2 as urllib
11from datetime import datetime, date, timedelta
12import sys
13from StringIO import StringIO
14import decimal
15
16DATEFORMAT = "%Y-%m-%d"
17
18class Emsc(object):
19    """
20    Handles ESMC data request via seismicportal.eu
21    """
22
23    BASEURI = "http://www.seismicportal.eu/services/event/search?"
24
25    def __init__(self, **kwargs):
26        """
27        Query web service
28        """
29
30        # dict of events
31        self.events = {}
32        # list of all origins
33        self.origins = []
34
35        # build query string
36        agencies = kwargs.get("agency", None)
37        if agencies:
38            agencies = agencies.upper().split(",")
39
40            # If only one agency is given, include in query. Otherwise it will be
41            # filtered later.
42            if len(agencies) == 1:
43                kwargs["auth"] = kwargs["agency"]
44                agencies = None
45        self.agencies = agencies
46       
47        # Also supply default values here, since they might not defined yet.
48        self.out = kwargs.get("outfile", sys.stdout)
49        self.debugLevel = kwargs.get("debug", False)
50        self.nomagnitude = kwargs.get("nomagnitude", False)
51        kwargs["limit"] = kwargs.get("limit", "2000")
52
53        if self.debugLevel:
54            print "Query EMSC from %s to %s (agency: %s)" % (
55                        kwargs["dateMin"], kwargs["dateMax"], kwargs["agency"])
56
57        try:
58            self.out.closed
59        except:
60            self.out = open(self.out, "w")
61
62        # clean up kwargs
63        discard = ["agency", "outfile", "debug", "nomagnitude",]
64        for i in discard:
65            if kwargs.has_key(i):
66                del kwargs[i]
67
68        # build uri
69        uri = self.BASEURI + "&".join(map(lambda x: "=".join(x), \
70                                           zip(kwargs.keys(), kwargs.values())))
71
72        self.__load(uri)
73
74    def __load(self, uri):
75        """
76        Load data from ESMC web service.
77        """
78
79#        data = StringIO(open("example.xml").read())
80        data = StringIO(urllib.urlopen(uri).read())
81
82        # Some data handling
83        for e in QuakeML(data, self.debugLevel).events:
84            # two decimal places for locations
85            e.longitude = e.longitude.quantize(decimal.Decimal('.01'), \
86                                                    rounding=decimal.ROUND_DOWN)
87            e.latitude = e.latitude.quantize(decimal.Decimal('.01'), \
88                                                    rounding=decimal.ROUND_DOWN)
89
90            # add south/north/east/west indicator (separates by two spaces)
91            if e.longitude < 0:
92                e.longitudeEW = str(-1 * e.longitude) + "  W"
93            else:
94                e.longitudeEW = str(e.longitude) + "  E"
95
96            if e.latitude < 0:
97                e.latitudeNS = str(-1*e.latitude) + "  S"
98            else:
99                e.latitudeNS = str(e.latitude) + "  N"
100
101            e.agencyOrigin = e.agencyOrigin.split("/")[-1]
102
103            # Check for multiple agency criteria.
104            if self.agencies and e.agencyOrigin not in self.agencies:
105                continue
106
107            # some localisations do not include a magnitude
108            try:
109                e.agencyMagnitude = e.agencyMagnitude.split("/")[-1]
110            except AttributeError:
111                # skip output if option not set
112                if not self.nomagnitude:
113                    continue
114
115                e.agencyMagnitude = "NOTSET"
116                e.magnitude = Decimal("NaN")
117                e.magnitudeType = "no"
118
119            self.origins.append(e)
120            try:
121                self.events[e.eventID].append(e)
122            except KeyError:
123                self.events[e.eventID] = [e]
124
125    def Ascii(self):
126        """
127        Ascii output of data.
128        """
129        for e in self.origins:
130            print >> self.out, "  ".join(e.origin.split("T"))[:-6],
131            print >> self.out, e.latitudeNS, " ", e.longitudeEW, " ", e.depth,
132            print >> self.out, " ", e.magnitude, e.magnitudeType,
133            print >> self.out, e.mode[0].upper(), e.description, " ",
134            print >> self.out, e.agencyOrigin
135
136    def Oed(self):
137        """
138        Returns origin-epicentre-depth list. There will be only one origin per
139        event. If no trusted agencies are set, the first origin found is used,
140        otherwise the first origin of a trusted agencies is returned.
141
142        In any case manual solutions are preferred against automatic ones.
143        """
144        oed = []
145
146        for eventID in sorted(self.events.keys()):
147            # classify by mode (manual, automatic) solutions
148            evts = {}
149            for e in self.events[eventID]:
150                try:
151                    evts[e.mode].append(e)
152                except KeyError:
153                    evts[e.mode] = [e]
154
155            # look into manual locations
156            solution = evts.get("manual", None)
157
158            # only automatic solution found:
159            if not solution:
160                solution = evts["automatic"]
161
162            # reformat to meet oed format
163            try:
164                res = (
165                    datetime.strptime(solution[0].origin, \
166                                                       "%Y-%m-%dT%H:%M:%S%f%z"),
167                    float(solution[0].latitude),
168                    float(solution[0].longitude),
169                    float(solution[0].depth),
170                    float(solution[0].magnitude),
171                )
172            except ValueError:
173                # %f and %z are only supported from python 2.6+
174                res = (
175                    datetime.strptime(solution[0].origin[:-8], \
176                                                           "%Y-%m-%dT%H:%M:%S"),
177                    float(solution[0].latitude),
178                    float(solution[0].longitude),
179                    float(solution[0].depth),
180                    float(solution[0].magnitude),
181                )
182
183            oed.append(res)
184
185        return oed
186
187def main():
188    from optparse import OptionParser
189    usage = """
190    %prog [options] [dateMin] [dateMax]
191
192    If dateMin/dateMax is not given a query for the last 5 days is issued. Date
193    must be formatted as YYYY-MM-DD (e.g. 2010-06-24).
194
195    If you enter only one datum, the query will cover data from this date to
196    today. If you want to get all events from today, simply enter today's date."""
197
198    parser = OptionParser(usage=usage)
199    parser.add_option("-a", "--agency", dest="agency",
200                  help="Query only for given agency (e.g. NEIR). Separate "
201                       "multiple agencies by comma: NEIR,GSRC.")
202    parser.add_option("-o", "--outfile", dest="outfile", default=sys.stdout,
203                  help="Write result to FILE [default stdout]", metavar="FILE")
204    parser.add_option("-l", "--limit", dest="limit", default='2000',
205                  help="Maximum number of events returned [default 2.000]. "
206                       "Please note: There might be multiple origin locations "
207                       "per event.")
208    parser.add_option("-d", "--debug", action="count", dest="debug",
209                  default=False,
210                  help="Switch on debugging [default off]. If this parameter "
211                       "is used twice ('-dd' or '-d -d') the input XML data "
212                       "will be dumped to 'debug-dump.xml' in the current "
213                       "directory.")
214    parser.add_option("--include-no-magnitude", action="store_true",
215                  dest="nomagnitude", default=False,
216                  help="Also output localisations that include no magnitude "
217                       "information [default off]. If this option is set, these "
218                       "events will have 'NaN' as magnitude and 'no' as type of "
219                       "magnitude.")
220
221    (options, args) = parser.parse_args()
222    query = vars(options)
223
224    query["dateMin"] = None
225    query["dateMax"] = None
226
227    if len(args) == 0:
228        # no arguments: today and 5 days back
229        query["dateMax"] = date.today().strftime(DATEFORMAT)
230        query["dateMin"] = (date.today() - timedelta(5)).strftime(DATEFORMAT)
231    elif len(args) == 1:
232        # one argument: from date till today
233        query["dateMin"] = args[0]
234        query["dateMax"] = date.today().strftime(DATEFORMAT)
235    elif len(args) == 2:
236        # timespan given: sort it
237        if args[0] > args[1]:
238            query["dateMin"] = args[1]
239            query["dateMax"] = args[0]
240        else:
241            query["dateMin"] = args[0]
242            query["dateMax"] = args[1]
243    else:
244        parser.error("Only two arguments can be processed.")
245
246    # add time information to query the hole day
247    query["dateMin"] += "T00:00:00"
248    query["dateMax"] += "T23:59:59"
249
250    _ = Emsc(**query).Ascii()
251
252if __name__ == "__main__":
253    main()
Note: See TracBrowser for help on using the repository browser.