1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import time
7import calendar
8import re
9import logging
11from pyrocko import model, util
12from pyrocko.moment_tensor import MomentTensor
13from pyrocko.util import Request, urlopen
14from .base_catalog import EarthquakeCatalog
17import numpy as num
20logger = logging.getLogger('pyrocko.client.globalcmt')
22km = 1000.
25class Anon(object):
26 pass
29class GlobalCMT(EarthquakeCatalog):
31 def __init__(self):
32 self.events = {}
34 def flush(self):
35 self.events = {}
37 def iter_event_names(
38 self,
39 time_range=None,
40 magmin=0.,
41 magmax=10.,
42 latmin=-90.,
43 latmax=90.,
44 lonmin=-180.,
45 lonmax=180.,
46 depthmin=0.,
47 depthmax=1000*km):
49 yearbeg, monbeg, daybeg = time.gmtime(time_range[0])[:3]
50 yearend, monend, dayend = time.gmtime(time_range[1])[:3]
52 url = 'http://www.globalcmt.org/cgi-bin/globalcmt-cgi-bin/CMT5/form?' \
53 + '&'.join([
54 'itype=ymd',
55 'yr=%i' % yearbeg, 'mo=%i' % monbeg, 'day=%i' % daybeg,
56 'otype=ymd',
57 'oyr=%i' % yearend, 'omo=%i' % monend, 'oday=%i' % dayend,
58 'jyr=1976', 'jday=1', 'ojyr=1976', 'ojday=1', 'nday=1',
59 'lmw=%g' % magmin, 'umw=%g' % magmax,
60 'lms=0', 'ums=10',
61 'lmb=0', 'umb=10',
62 'llat=%g' % latmin, 'ulat=%g' % latmax,
63 'llon=%g' % lonmin, 'ulon=%g' % lonmax,
64 'lhd=%g' % (depthmin/km), 'uhd=%g' % (depthmax/km),
65 'lts=-9999', 'uts=9999',
66 'lpe1=0', 'upe1=90',
67 'lpe2=0', 'upe2=90',
68 'list=5'])
70 while True:
71 logger.debug('Opening URL: %s' % url)
72 req = Request(url)
73 page = urlopen(req).read()
74 logger.debug('Received page (%i bytes)' % len(page))
76 events, more = self._parse_events_page(page)
78 for ev in events:
79 self.events[ev.name] = ev
81 for ev in events:
82 if time_range[0] <= ev.time and ev.time <= time_range[1]:
83 yield ev.name
85 if more:
86 url = more.decode('ascii')
87 else:
88 break
90 def get_event(self, name):
91 if name not in self.events:
92 t = self._name_to_date(name)
93 for name2 in self.iter_event_names(
94 time_range=(t-24*60*60, t+2*24*60*60)):
96 if name2 == name:
97 break
99 return self.events[name]
101 def _parse_events_page(self, page):
103 lines = page.splitlines()
104 state = 0
106 events = []
108 def complete(data):
109 try:
110 t = util.to_time_float(
111 calendar.timegm((
112 data.year, data.month, data.day,
113 data.hour, data.minute, data.seconds)))
115 m = num.array(
116 [data.mrr, data.mrt, data.mrp,
117 data.mrt, data.mtt, data.mtp,
118 data.mrp, data.mtp, data.mpp],
119 dtype=float).reshape(3, 3)
121 m *= 10.0**(data.exponent-7)
122 mt = MomentTensor(m_up_south_east=m)
123 ev = model.Event(
124 lat=data.lat,
125 lon=data.lon,
126 time=t,
127 name=data.eventname,
128 depth=data.depth_km*1000.,
129 magnitude=float(mt.moment_magnitude()),
130 duration=data.half_duration * 2.,
131 region=data.region.rstrip(),
132 catalog=data.catalog)
134 ev.moment_tensor = mt
135 events.append(ev)
137 except AttributeError:
138 pass
140 catalog = 'gCMT'
142 data = None
143 more = None
144 for line in lines:
145 if state == 0:
147 m = re.search(br'<a href="([^"]+)">More solutions', line)
148 if m:
149 more = m.group(1)
151 m = re.search(br'From Quick CMT catalog', line)
152 if m:
153 catalog = 'gCMT-Q'
155 m = re.search(br'Event name:\s+(\S+)', line)
156 if m:
157 if data:
158 complete(data)
160 data = Anon()
161 data.eventname = str(m.group(1).decode('ascii'))
162 data.catalog = catalog
164 if data:
165 m = re.search(br'Region name:\s+([^<]+)', line)
166 if m:
167 data.region = str(m.group(1).decode('ascii'))
169 m = re.search(
170 br'Date \(y/m/d\): (\d\d\d\d)/(\d+)/(\d+)', line)
172 if m:
173 data.year, data.month, data.day = (
174 int(m.group(1)), int(m.group(2)), int(m.group(3)))
176 m = re.search(br'Timing and location information', line)
177 if m:
178 state = 1
180 if state == 1:
181 toks = line.split()
182 if toks and toks[0] == b'CMT':
183 data.hour, data.minute = [int(x) for x in toks[1:3]]
184 data.seconds, data.lat, data.lon, data.depth_km = [
185 float(x) for x in toks[3:]]
187 m = re.search(br'Assumed half duration:\s+(\S+)', line)
188 if m:
189 data.half_duration = float(m.group(1))
191 m = re.search(br'Mechanism information', line)
192 if m:
193 state = 2
195 if state == 2:
196 m = re.search(br'Exponent for moment tensor:\s+(\d+)', line)
197 if m:
198 data.exponent = int(m.group(1))
200 toks = line.split()
201 if toks and toks[0] == b'CMT':
202 data.mrr, data.mtt, data.mpp, \
203 data.mrt, data.mrp, data.mtp = [
204 float(x) for x in toks[1:]]
206 m = re.search(br'^Eigenvector:', line)
207 if m:
208 state = 0
210 if data is not None:
211 complete(data)
213 return events, more
215 def _name_to_date(self, name):
217 if len(name) == 7:
218 y, m, d = time.strptime(name[:6], '%m%d%y')[:3]
219 if y > 2005:
220 y -= 100
222 else:
223 y, m, d = time.strptime(name[:8], '%Y%m%d')[:3]
225 t = util.to_time_float(
226 calendar.timegm((y, m, d, 0, 0, 0)))
227 return t