Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/client/globalcmt.py: 95%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Client to get earthquake catalog information from 

8`GlobalCMT <http://www.globalcmt.org/>`_. 

9''' 

10 

11import time 

12import calendar 

13import re 

14import logging 

15 

16from pyrocko import model, util 

17from pyrocko.moment_tensor import MomentTensor 

18from pyrocko.util import Request, urlopen 

19from .base_catalog import EarthquakeCatalog 

20 

21 

22import numpy as num 

23 

24 

25logger = logging.getLogger('pyrocko.client.globalcmt') 

26 

27km = 1000. 

28 

29 

30class Anon(object): 

31 pass 

32 

33 

34class GlobalCMT(EarthquakeCatalog): 

35 ''' 

36 Client to get earthquake catalog information from 

37 `GlobalCMT <http://www.globalcmt.org/>`_. 

38 ''' 

39 

40 def __init__(self): 

41 self.events = {} 

42 

43 def flush(self): 

44 self.events = {} 

45 

46 def iter_event_names( 

47 self, 

48 time_range=None, 

49 magmin=0., 

50 magmax=10., 

51 latmin=-90., 

52 latmax=90., 

53 lonmin=-180., 

54 lonmax=180., 

55 depthmin=0., 

56 depthmax=1000*km): 

57 

58 yearbeg, monbeg, daybeg = time.gmtime(time_range[0])[:3] 

59 yearend, monend, dayend = time.gmtime(time_range[1])[:3] 

60 

61 url = 'http://www.globalcmt.org/cgi-bin/globalcmt-cgi-bin/CMT5/form?' \ 

62 + '&'.join([ 

63 'itype=ymd', 

64 'yr=%i' % yearbeg, 'mo=%i' % monbeg, 'day=%i' % daybeg, 

65 'otype=ymd', 

66 'oyr=%i' % yearend, 'omo=%i' % monend, 'oday=%i' % dayend, 

67 'jyr=1976', 'jday=1', 'ojyr=1976', 'ojday=1', 'nday=1', 

68 'lmw=%g' % magmin, 'umw=%g' % magmax, 

69 'lms=0', 'ums=10', 

70 'lmb=0', 'umb=10', 

71 'llat=%g' % latmin, 'ulat=%g' % latmax, 

72 'llon=%g' % lonmin, 'ulon=%g' % lonmax, 

73 'lhd=%g' % (depthmin/km), 'uhd=%g' % (depthmax/km), 

74 'lts=-9999', 'uts=9999', 

75 'lpe1=0', 'upe1=90', 

76 'lpe2=0', 'upe2=90', 

77 'list=5']) 

78 

79 while True: 

80 logger.debug('Opening URL: %s' % url) 

81 req = Request(url) 

82 page = urlopen(req).read() 

83 logger.debug('Received page (%i bytes)' % len(page)) 

84 

85 events, more = self._parse_events_page(page) 

86 

87 for ev in events: 

88 self.events[ev.name] = ev 

89 

90 for ev in events: 

91 if time_range[0] <= ev.time and ev.time <= time_range[1]: 

92 yield ev.name 

93 

94 if more: 

95 url = more.decode('ascii') 

96 else: 

97 break 

98 

99 def get_event(self, name): 

100 if name not in self.events: 

101 t = self._name_to_date(name) 

102 for name2 in self.iter_event_names( 

103 time_range=(t-24*60*60, t+2*24*60*60)): 

104 

105 if name2 == name: 

106 break 

107 

108 return self.events[name] 

109 

110 def _parse_events_page(self, page): 

111 

112 lines = page.splitlines() 

113 state = 0 

114 

115 events = [] 

116 

117 def complete(data): 

118 try: 

119 t = util.to_time_float( 

120 calendar.timegm(( 

121 data.year, data.month, data.day, 

122 data.hour, data.minute, data.seconds))) 

123 

124 m = num.array( 

125 [data.mrr, data.mrt, data.mrp, 

126 data.mrt, data.mtt, data.mtp, 

127 data.mrp, data.mtp, data.mpp], 

128 dtype=float).reshape(3, 3) 

129 

130 m *= 10.0**(data.exponent-7) 

131 mt = MomentTensor(m_up_south_east=m) 

132 ev = model.Event( 

133 lat=data.lat, 

134 lon=data.lon, 

135 time=t, 

136 name=data.eventname, 

137 depth=data.depth_km*1000., 

138 magnitude=float(mt.moment_magnitude()), 

139 duration=data.half_duration * 2., 

140 region=data.region.rstrip(), 

141 catalog=data.catalog) 

142 

143 ev.moment_tensor = mt 

144 events.append(ev) 

145 

146 except AttributeError: 

147 pass 

148 

149 catalog = 'gCMT' 

150 

151 data = None 

152 more = None 

153 for line in lines: 

154 if state == 0: 

155 

156 m = re.search(br'<a href="([^"]+)">More solutions', line) 

157 if m: 

158 more = m.group(1) 

159 

160 m = re.search(br'From Quick CMT catalog', line) 

161 if m: 

162 catalog = 'gCMT-Q' 

163 

164 m = re.search(br'Event name:\s+(\S+)', line) 

165 if m: 

166 if data: 

167 complete(data) 

168 

169 data = Anon() 

170 data.eventname = str(m.group(1).decode('ascii')) 

171 data.catalog = catalog 

172 

173 if data: 

174 m = re.search(br'Region name:\s+([^<]+)', line) 

175 if m: 

176 data.region = str(m.group(1).decode('ascii')) 

177 

178 m = re.search( 

179 br'Date \(y/m/d\): (\d\d\d\d)/(\d+)/(\d+)', line) 

180 

181 if m: 

182 data.year, data.month, data.day = ( 

183 int(m.group(1)), int(m.group(2)), int(m.group(3))) 

184 

185 m = re.search(br'Timing and location information', line) 

186 if m: 

187 state = 1 

188 

189 if state == 1: 

190 toks = line.split() 

191 if toks and toks[0] == b'CMT': 

192 data.hour, data.minute = [int(x) for x in toks[1:3]] 

193 data.seconds, data.lat, data.lon, data.depth_km = [ 

194 float(x) for x in toks[3:]] 

195 

196 m = re.search(br'Assumed half duration:\s+(\S+)', line) 

197 if m: 

198 data.half_duration = float(m.group(1)) 

199 

200 m = re.search(br'Mechanism information', line) 

201 if m: 

202 state = 2 

203 

204 if state == 2: 

205 m = re.search(br'Exponent for moment tensor:\s+(\d+)', line) 

206 if m: 

207 data.exponent = int(m.group(1)) 

208 

209 toks = line.split() 

210 if toks and toks[0] == b'CMT': 

211 data.mrr, data.mtt, data.mpp, \ 

212 data.mrt, data.mrp, data.mtp = [ 

213 float(x) for x in toks[1:]] 

214 

215 m = re.search(br'^Eigenvector:', line) 

216 if m: 

217 state = 0 

218 

219 if data is not None: 

220 complete(data) 

221 

222 return events, more 

223 

224 def _name_to_date(self, name): 

225 

226 if len(name) == 7: 

227 y, m, d = time.strptime(name[:6], '%m%d%y')[:3] 

228 if y > 2005: 

229 y -= 100 

230 

231 else: 

232 y, m, d = time.strptime(name[:8], '%Y%m%d')[:3] 

233 

234 t = util.to_time_float( 

235 calendar.timegm((y, m, d, 0, 0, 0))) 

236 return t