1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import time
8import calendar
9import re
10import logging
12from pyrocko import model, util
13from pyrocko.moment_tensor import MomentTensor
14from pyrocko.util import Request, urlopen
15from .base_catalog import EarthquakeCatalog
18import numpy as num
21logger = logging.getLogger('pyrocko.client.globalcmt')
23km = 1000.
26class Anon(object):
27 pass
30class GlobalCMT(EarthquakeCatalog):
32 def __init__(self):
33 self.events = {}
35 def flush(self):
36 self.events = {}
38 def iter_event_names(
39 self,
40 time_range=None,
41 magmin=0.,
42 magmax=10.,
43 latmin=-90.,
44 latmax=90.,
45 lonmin=-180.,
46 lonmax=180.,
47 depthmin=0.,
48 depthmax=1000*km):
50 yearbeg, monbeg, daybeg = time.gmtime(time_range[0])[:3]
51 yearend, monend, dayend = time.gmtime(time_range[1])[:3]
53 url = 'http://www.globalcmt.org/cgi-bin/globalcmt-cgi-bin/CMT5/form?' \
54 + '&'.join([
55 'itype=ymd',
56 'yr=%i' % yearbeg, 'mo=%i' % monbeg, 'day=%i' % daybeg,
57 'otype=ymd',
58 'oyr=%i' % yearend, 'omo=%i' % monend, 'oday=%i' % dayend,
59 'jyr=1976', 'jday=1', 'ojyr=1976', 'ojday=1', 'nday=1',
60 'lmw=%g' % magmin, 'umw=%g' % magmax,
61 'lms=0', 'ums=10',
62 'lmb=0', 'umb=10',
63 'llat=%g' % latmin, 'ulat=%g' % latmax,
64 'llon=%g' % lonmin, 'ulon=%g' % lonmax,
65 'lhd=%g' % (depthmin/km), 'uhd=%g' % (depthmax/km),
66 'lts=-9999', 'uts=9999',
67 'lpe1=0', 'upe1=90',
68 'lpe2=0', 'upe2=90',
69 'list=5'])
71 while True:
72 logger.debug('Opening URL: %s' % url)
73 req = Request(url)
74 page = urlopen(req).read()
75 logger.debug('Received page (%i bytes)' % len(page))
77 events, more = self._parse_events_page(page)
79 for ev in events:
80 self.events[ev.name] = ev
82 for ev in events:
83 if time_range[0] <= ev.time and ev.time <= time_range[1]:
84 yield ev.name
86 if more:
87 url = more.decode('ascii')
88 else:
89 break
91 def get_event(self, name):
92 if name not in self.events:
93 t = self._name_to_date(name)
94 for name2 in self.iter_event_names(
95 time_range=(t-24*60*60, t+2*24*60*60)):
97 if name2 == name:
98 break
100 return self.events[name]
102 def _parse_events_page(self, page):
104 lines = page.splitlines()
105 state = 0
107 events = []
109 def complete(data):
110 try:
111 t = util.to_time_float(
112 calendar.timegm((
113 data.year, data.month, data.day,
114 data.hour, data.minute, data.seconds)))
116 m = num.array(
117 [data.mrr, data.mrt, data.mrp,
118 data.mrt, data.mtt, data.mtp,
119 data.mrp, data.mtp, data.mpp],
120 dtype=float).reshape(3, 3)
122 m *= 10.0**(data.exponent-7)
123 mt = MomentTensor(m_up_south_east=m)
124 ev = model.Event(
125 lat=data.lat,
126 lon=data.lon,
127 time=t,
128 name=data.eventname,
129 depth=data.depth_km*1000.,
130 magnitude=float(mt.moment_magnitude()),
131 duration=data.half_duration * 2.,
132 region=data.region.rstrip(),
133 catalog=data.catalog)
135 ev.moment_tensor = mt
136 events.append(ev)
138 except AttributeError:
139 pass
141 catalog = 'gCMT'
143 data = None
144 more = None
145 for line in lines:
146 if state == 0:
148 m = re.search(br'<a href="([^"]+)">More solutions', line)
149 if m:
150 more = m.group(1)
152 m = re.search(br'From Quick CMT catalog', line)
153 if m:
154 catalog = 'gCMT-Q'
156 m = re.search(br'Event name:\s+(\S+)', line)
157 if m:
158 if data:
159 complete(data)
161 data = Anon()
162 data.eventname = str(m.group(1).decode('ascii'))
163 data.catalog = catalog
165 if data:
166 m = re.search(br'Region name:\s+([^<]+)', line)
167 if m:
168 data.region = str(m.group(1).decode('ascii'))
170 m = re.search(
171 br'Date \(y/m/d\): (\d\d\d\d)/(\d+)/(\d+)', line)
173 if m:
174 data.year, data.month, data.day = (
175 int(m.group(1)), int(m.group(2)), int(m.group(3)))
177 m = re.search(br'Timing and location information', line)
178 if m:
179 state = 1
181 if state == 1:
182 toks = line.split()
183 if toks and toks[0] == b'CMT':
184 data.hour, data.minute = [int(x) for x in toks[1:3]]
185 data.seconds, data.lat, data.lon, data.depth_km = [
186 float(x) for x in toks[3:]]
188 m = re.search(br'Assumed half duration:\s+(\S+)', line)
189 if m:
190 data.half_duration = float(m.group(1))
192 m = re.search(br'Mechanism information', line)
193 if m:
194 state = 2
196 if state == 2:
197 m = re.search(br'Exponent for moment tensor:\s+(\d+)', line)
198 if m:
199 data.exponent = int(m.group(1))
201 toks = line.split()
202 if toks and toks[0] == b'CMT':
203 data.mrr, data.mtt, data.mpp, \
204 data.mrt, data.mrp, data.mtp = [
205 float(x) for x in toks[1:]]
207 m = re.search(br'^Eigenvector:', line)
208 if m:
209 state = 0
211 if data is not None:
212 complete(data)
214 return events, more
216 def _name_to_date(self, name):
218 if len(name) == 7:
219 y, m, d = time.strptime(name[:6], '%m%d%y')[:3]
220 if y > 2005:
221 y -= 100
223 else:
224 y, m, d = time.strptime(name[:8], '%Y%m%d')[:3]
226 t = util.to_time_float(
227 calendar.timegm((y, m, d, 0, 0, 0)))
228 return t