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

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Access to some IRIS (non-FDSN) web services. 

8''' 

9 

10import logging 

11import re 

12from xml.parsers.expat import ParserCreate 

13 

14from pyrocko import util, pz, model 

15from pyrocko.util import Request, HTTPError, urlopen, urlencode 

16 

17 

18logger = logging.getLogger('pyrocko.client.iris') 

19 

20base_url = 'http://service.iris.edu/irisws' 

21 

22 

23def tdate(s): 

24 return util.str_to_time(s, '%Y-%m-%d') 

25 

26 

27def sdate(t): 

28 return util.time_to_str(t, '%Y-%m-%d') 

29 

30 

31def tdatetime(s): 

32 return util.str_to_time(s, format='%Y-%m-%dT%H:%M:%S') 

33 

34 

35def sdatetime(t): 

36 return util.time_to_str(t, format='%Y-%m-%dT%H:%M:%S') 

37 

38 

39class Element(object): 

40 

41 def __init__(self, name, depth, attrs): 

42 self.name = name 

43 self.depth = depth 

44 self.attrs = attrs 

45 

46 def __getattr__(self, k): 

47 return self.attrs[k] 

48 

49 def __str__(self): 

50 return '%s:\n ' % self.name + '\n '.join( 

51 '%s : %s' % (k, v) for (k, v) in self.attrs.items()) 

52 

53 

54class XMLZipHandler(object): 

55 

56 def __init__(self, watch): 

57 self.watch = watch 

58 self.stack = [] 

59 self.wstack = [] 

60 self.outstack = [] 

61 

62 def startElement(self, name, attrs): 

63 self.stack.append((name, [])) 

64 if len(self.wstack) < len(self.watch) and \ 

65 name == self.watch[len(self.wstack)]: 

66 

67 el = Element(name, len(self.stack), dict(attrs.items())) 

68 self.wstack.append(el) 

69 

70 def characters(self, content): 

71 if self.wstack and len(self.stack) == self.wstack[-1].depth + 1: 

72 if content.strip(): 

73 self.stack[-1][1].append(content) 

74 

75 def endElement(self, name): 

76 if self.wstack: 

77 if len(self.stack) == self.wstack[-1].depth + 1 and \ 

78 self.stack[-1][1]: 

79 

80 self.wstack[-1].attrs[name] = ''.join(self.stack[-1][1]) 

81 

82 if name == self.watch[len(self.wstack)-1]: 

83 if len(self.wstack) == len(self.watch): 

84 self.outstack.append(list(self.wstack)) 

85 

86 self.wstack.pop() 

87 

88 self.stack.pop() 

89 

90 def getQueuedElements(self): 

91 outstack = self.outstack 

92 self.outstack = [] 

93 return outstack 

94 

95 

96def xmlzip(source, watch, bufsize=10000): 

97 parser = ParserCreate() 

98 handler = XMLZipHandler(watch) 

99 

100 parser.StartElementHandler = handler.startElement 

101 parser.EndElementHandler = handler.endElement 

102 parser.CharacterDataHandler = handler.characters 

103 

104 while True: 

105 data = source.read(bufsize) 

106 if not data: 

107 break 

108 

109 parser.Parse(data, False) 

110 for elements in handler.getQueuedElements(): 

111 yield elements 

112 

113 parser.Parse('', True) 

114 for elements in handler.getQueuedElements(): 

115 yield elements 

116 

117 

118class NotFound(Exception): 

119 ''' 

120 Raised when the server sends an 404 response. 

121 ''' 

122 def __init__(self, url): 

123 Exception.__init__(self) 

124 self._url = url 

125 

126 def __str__(self): 

127 return 'No results for request %s' % self._url 

128 

129 

130def ws_request(url, post=False, **kwargs): 

131 url_values = urlencode(kwargs) 

132 url = url + '?' + url_values 

133 logger.debug('Accessing URL %s' % url) 

134 

135 req = Request(url) 

136 if post: 

137 req.add_data(post) 

138 

139 req.add_header('Accept', '*/*') 

140 

141 try: 

142 return urlopen(req) 

143 

144 except HTTPError as e: 

145 if e.code == 404: 

146 raise NotFound(url) 

147 else: 

148 raise e 

149 

150 

151def ws_virtualnetwork(**kwargs): 

152 ''' 

153 Query IRIS virtualnetwork web service. 

154 ''' 

155 

156 for k in 'starttime', 'endtime': 

157 if k in kwargs: 

158 kwargs[k] = sdate(kwargs[k]) 

159 

160 if 'timewindow' in kwargs: 

161 tmin, tmax = kwargs.pop('timewindow') 

162 kwargs['starttime'] = sdate(tmin) 

163 kwargs['endtime'] = sdate(tmax) 

164 

165 return ws_request(base_url + '/virtualnetwork/1/query', **kwargs) 

166 

167 

168def ws_sacpz( 

169 network=None, 

170 station=None, 

171 location=None, 

172 channel=None, 

173 time=None, 

174 tmin=None, 

175 tmax=None): 

176 ''' 

177 Query IRIS sacpz web service. 

178 ''' 

179 

180 d = {} 

181 if network: 

182 d['network'] = network 

183 if station: 

184 d['station'] = station 

185 if location: 

186 d['location'] = location 

187 else: 

188 d['location'] = '--' 

189 if channel: 

190 d['channel'] = channel 

191 

192 if tmin is not None and tmax is not None: 

193 d['starttime'] = sdatetime(tmin) 

194 d['endtime'] = sdatetime(tmax) 

195 elif time is not None: 

196 d['time'] = sdatetime(time) 

197 

198 return ws_request(base_url + '/sacpz/1/query', **d) 

199 

200 

201def ws_resp( 

202 network=None, 

203 station=None, 

204 location=None, 

205 channel=None, 

206 time=None, 

207 tmin=None, 

208 tmax=None): 

209 ''' 

210 Query IRIS resp web service. 

211 ''' 

212 

213 d = {} 

214 if network: 

215 d['network'] = network 

216 if station: 

217 d['station'] = station 

218 if location: 

219 d['location'] = location 

220 else: 

221 d['location'] = '--' 

222 if channel: 

223 d['channel'] = channel 

224 

225 if tmin is not None and tmax is not None: 

226 d['starttime'] = sdatetime(tmin) 

227 d['endtime'] = sdatetime(tmax) 

228 elif time is not None: 

229 d['time'] = sdatetime(time) 

230 

231 return ws_request(base_url + '/resp/1/query', **d) 

232 

233 

234class ChannelInfo(object): 

235 def __init__( 

236 self, network, station, location, channel, start, end, azimuth, 

237 dip, elevation, depth, latitude, longitude, sample, input, output, 

238 zpk): 

239 

240 self.network = network 

241 self.station = station 

242 self.location = location 

243 self.channel = channel 

244 self.start = start 

245 self.end = end 

246 self.azimuth = azimuth 

247 self.dip = dip 

248 self.elevation = elevation 

249 self.depth = depth 

250 self.latitude = latitude 

251 self.longitude = longitude 

252 self.sample = sample 

253 self.input = input 

254 self.output = output 

255 self.zpk = zpk 

256 

257 def __str__(self): 

258 return '%s.%s.%s.%s' % ( 

259 self.network, self.station, self.location, self.channel) 

260 

261 

262def nslc(x): 

263 return x.network, x.station, x.location, x.channel 

264 

265 

266def grok_sacpz(data): 

267 pzlines = [] 

268 d = {} 

269 responses = [] 

270 float_keys = ('latitude', 'longitude', 'elevation', 'depth', 'dip', 

271 'azimuth', 'sample') 

272 string_keys = ('input', 'output', 'network', 'station', 'location', 

273 'channel') 

274 time_keys = ('start', 'end') 

275 for line in data.splitlines(): 

276 line = line.strip() 

277 if line.startswith('*'): 

278 if pzlines: 

279 if any(pzlines): 

280 d['zpk'] = pz.read_sac_zpk(string='\n'.join(pzlines)) 

281 responses.append(d) 

282 d = {} 

283 pzlines = [] 

284 

285 m = re.match(r'^\* ([A-Z]+)[^:]*:(.*)$', line) 

286 if m: 

287 k, v = m.group(1).lower(), m.group(2).strip() 

288 if k in d: 

289 assert False, 'duplicate entry? %s' % k 

290 

291 if k in float_keys: 

292 d[k] = float(v) 

293 elif k in string_keys: 

294 d[k] = v 

295 elif k in time_keys: 

296 d[k] = tdatetime(v) 

297 

298 else: 

299 pzlines.append(line) 

300 

301 if pzlines and any(pzlines): 

302 d['zpk'] = pz.read_sac_zpk(string='\n'.join(pzlines)) 

303 responses.append(d) 

304 

305 cis = {} 

306 for kwargs in responses: 

307 try: 

308 for k in float_keys + string_keys + time_keys: 

309 if k not in kwargs: 

310 logger.error('Missing entry: %s' % k) 

311 raise Exception() 

312 

313 ci = ChannelInfo(**kwargs) 

314 

315 cis[nslc(ci)] = ci 

316 

317 except Exception: 

318 logger.error('Error while parsing SACPZ data') 

319 

320 return cis 

321 

322 

323def grok_station_xml(data, tmin, tmax): 

324 

325 stations = {} 

326 

327 for (sta, sta_epo, cha, cha_epo) in xmlzip(data, ( 

328 'Station', 'StationEpoch', 'Channel', 'Epoch')): 

329 

330 sta_beg, sta_end, cha_beg, cha_end = [tdatetime(x) for x in ( 

331 sta_epo.StartDate, sta_epo.EndDate, cha_epo.StartDate, 

332 cha_epo.EndDate)] 

333 

334 if not (sta_beg <= tmin and tmax <= sta_end and 

335 cha_beg <= tmin and tmax <= cha_end): 

336 

337 continue 

338 

339 nslc = tuple([str(x.strip()) for x in ( 

340 sta.net_code, sta.sta_code, cha.loc_code, cha.chan_code)]) 

341 

342 lat, lon, ele, dep, azi, dip = [ 

343 float(cha_epo.attrs[x]) 

344 for x in 'Lat Lon Elevation Depth Azimuth Dip'.split()] 

345 

346 nsl = nslc[:3] 

347 if nsl not in stations: 

348 stations[nsl] = model.Station( 

349 nsl[0], nsl[1], nsl[2], lat, lon, ele, dep) 

350 

351 stations[nsl].add_channel(model.station.Channel(nslc[-1], azi, dip)) 

352 

353 return list(stations.values()) 

354 

355 

356def grok_virtualnet_xml(data): 

357 net_sta = set() 

358 

359 for network, station in xmlzip(data, ('network', 'station')): 

360 net_sta.add((network.code, station.code)) 

361 

362 return net_sta 

363 

364 

365def data_selection(stations, tmin, tmax, channel_prio=[ 

366 ['BHZ', 'HHZ'], 

367 ['BH1', 'BHN', 'HH1', 'HHN'], 

368 ['BH2', 'BHE', 'HH2', 'HHE']]): 

369 

370 selection = [] 

371 for station in stations: 

372 wanted = [] 

373 for group in channel_prio: 

374 gchannels = [] 

375 for channel in station.get_channels(): 

376 if channel.name in group: 

377 gchannels.append(channel) 

378 if gchannels: 

379 gchannels.sort(key=lambda a: group.index(a.name)) 

380 wanted.append(gchannels[0]) 

381 

382 if wanted: 

383 for channel in wanted: 

384 selection.append(( 

385 station.network, station.station, station.location, 

386 channel.name, tmin, tmax)) 

387 

388 return selection