Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/rdseed.py: 20%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Access to SEED volumes by interfacing `rdseed 

8<https://ds.iris.edu/ds/nodes/dmc/software/downloads/rdseed/>`_ (*deprecated*). 

9''' 

10 

11import os 

12import sys 

13import subprocess 

14import tempfile 

15import calendar 

16import time 

17import re 

18import logging 

19import shutil 

20 

21from pyrocko import ExternalProgramMissing 

22from pyrocko import pile, model, util, response 

23from pyrocko.io import eventdata 

24 

25pjoin = os.path.join 

26logger = logging.getLogger('pyrocko.io.rdseed') 

27 

28 

29def cmp(a, b): 

30 return (a > b) - (a < b) 

31 

32 

33def read_station_header_file(fn): 

34 

35 def m(i, *args): 

36 if len(ltoks) >= i + len(args) \ 

37 and (tuple(ltoks[i:i+len(args)]) == args): 

38 return True 

39 return False 

40 

41 def f(i): 

42 return float(toks[i]) 

43 

44 def s(i): 

45 if len(toks) > i: 

46 return toks[i] 

47 else: 

48 return '' 

49 

50 with open(fn, 'r') as fi: 

51 

52 stations = [] 

53 atsec, station, channel = None, None, None 

54 

55 for line in fi: 

56 toks = line.split() 

57 ltoks = line.lower().split() 

58 if m(2, 'station', 'header'): 

59 atsec = 'station' 

60 station = {'channels': []} 

61 stations.append(station) 

62 continue 

63 

64 if m(2, 'station') and m(5, 'channel'): 

65 atsec = 'channel' 

66 channel = {} 

67 station['channels'].append(channel) 

68 continue 

69 

70 if atsec == 'station': 

71 if m(1, 'station', 'code:'): 

72 station['station'] = s(3) 

73 

74 elif m(1, 'network', 'code:'): 

75 station['network'] = s(3) 

76 

77 elif m(1, 'name:'): 

78 station['name'] = ' '.join(toks[2:]) 

79 

80 if atsec == 'channel': 

81 if m(1, 'channel:'): 

82 channel['channel'] = s(2) 

83 

84 elif m(1, 'location:'): 

85 channel['location'] = s(2) 

86 

87 elif m(1, 'latitude:'): 

88 station['lat'] = f(2) 

89 

90 elif m(1, 'longitude:'): 

91 station['lon'] = f(2) 

92 

93 elif m(1, 'elevation:'): 

94 station['elevation'] = f(2) 

95 

96 elif m(1, 'local', 'depth:'): 

97 channel['depth'] = f(3) 

98 

99 elif m(1, 'azimuth:'): 

100 channel['azimuth'] = f(2) 

101 

102 elif m(1, 'dip:'): 

103 channel['dip'] = f(2) 

104 

105 nsl_stations = {} 

106 for station in stations: 

107 for channel in station['channels']: 

108 def cs(k, default=None): 

109 return channel.get(k, station.get(k, default)) 

110 

111 nsl = station['network'], station['station'], channel['location'] 

112 if nsl not in nsl_stations: 

113 nsl_stations[nsl] = model.Station( 

114 network=station['network'], 

115 station=station['station'], 

116 location=channel['location'], 

117 lat=cs('lat'), 

118 lon=cs('lon'), 

119 elevation=cs('elevation'), 

120 depth=cs('depth', None), 

121 name=station['name']) 

122 

123 nsl_stations[nsl].add_channel(model.Channel( 

124 channel['channel'], 

125 azimuth=channel['azimuth'], 

126 dip=channel['dip'])) 

127 

128 return list(nsl_stations.values()) 

129 

130 

131def cmp_version(a, b): 

132 ai = [int(x) for x in a.split('.')] 

133 bi = [int(x) for x in b.split('.')] 

134 return cmp(ai, bi) 

135 

136 

137def dumb_parser(data): 

138 

139 (in_ws, in_kw, in_str) = (1, 2, 3) 

140 

141 state = in_ws 

142 

143 rows = [] 

144 cols = [] 

145 accu = '' 

146 for c in data: 

147 if state == in_ws: 

148 if c == '"': 

149 new_state = in_str 

150 

151 elif c not in (' ', '\t', '\n', '\r'): 

152 new_state = in_kw 

153 

154 if state == in_kw: 

155 if c in (' ', '\t', '\n', '\r'): 

156 cols.append(accu) 

157 accu = '' 

158 if c in ('\n', '\r'): 

159 rows.append(cols) 

160 cols = [] 

161 new_state = in_ws 

162 

163 if state == in_str: 

164 if c == '"': 

165 accu += c 

166 cols.append(accu[1:-1]) 

167 accu = '' 

168 if c in ('\n', '\r'): 

169 rows.append(cols) 

170 cols = [] 

171 new_state = in_ws 

172 

173 state = new_state 

174 

175 if state in (in_kw, in_str): 

176 accu += c 

177 

178 if len(cols) != 0: 

179 rows.append(cols) 

180 

181 return rows 

182 

183 

184class Programs(object): 

185 rdseed = 'rdseed' 

186 avail = None 

187 

188 @staticmethod 

189 def check(): 

190 if Programs.avail is not None: 

191 return Programs.avail 

192 

193 else: 

194 try: 

195 rdseed_proc = subprocess.Popen( 

196 [Programs.rdseed], 

197 stdin=subprocess.PIPE, 

198 stdout=subprocess.PIPE, 

199 stderr=subprocess.PIPE) 

200 

201 (out, err) = rdseed_proc.communicate() 

202 

203 except OSError as e: 

204 if e.errno == 2: 

205 reason = "Could not find executable: '%s'." \ 

206 % Programs.rdseed 

207 else: 

208 reason = str(e) 

209 

210 logging.debug('Failed to run rdseed program. %s' % reason) 

211 Programs.avail = False 

212 return Programs.avail 

213 

214 ms = [re.search( 

215 r'Release (\d+(\.\d+(\.\d+)?)?)', s.decode()) 

216 for s in (err, out)] 

217 ms = list(filter(bool, ms)) 

218 if not ms: 

219 logger.error('Cannot determine rdseed version number.') 

220 else: 

221 version = ms[0].group(1) 

222 if cmp_version('4.7.5', version) == 1 \ 

223 or cmp_version(version, '5.1') == 1: 

224 

225 logger.warning( 

226 'Module pyrocko.rdseed has not been tested with ' 

227 'version %s of rdseed.' % version) 

228 

229 Programs.avail = True 

230 return Programs.avail 

231 

232 

233class SeedVolumeNotFound(Exception): 

234 pass 

235 

236 

237class SeedVolumeAccess(eventdata.EventDataAccess): 

238 

239 def __init__(self, seedvolume, datapile=None): 

240 

241 ''' 

242 Create new SEED Volume access object. 

243 

244 :param seedvolume: filename of seed volume 

245 :param datapile: if not ``None``, this should be a 

246 :py:class:`pile.Pile` object with data traces which are then used 

247 instead of the data provided by the SEED volume. 

248 (This is useful for dataless SEED volumes.) 

249 ''' 

250 

251 eventdata.EventDataAccess.__init__(self, datapile=datapile) 

252 self.tempdir = None 

253 Programs.check() 

254 

255 self.tempdir = None 

256 self.seedvolume = seedvolume 

257 if not os.path.isfile(self.seedvolume): 

258 raise SeedVolumeNotFound() 

259 

260 self.tempdir = tempfile.mkdtemp('', 'SeedVolumeAccess-') 

261 self.station_headers_file = os.path.join( 

262 self.tempdir, 'station_header_infos') 

263 self._unpack() 

264 self.shutil = shutil 

265 

266 def __del__(self): 

267 if self.tempdir: 

268 self.shutil.rmtree(self.tempdir) 

269 

270 def get_pile(self): 

271 if self._pile is None: 

272 fns = util.select_files([self.tempdir], include=r'\.SAC$') 

273 self._pile = pile.Pile() 

274 self._pile.load_files(fns, fileformat='sac') 

275 

276 return self._pile 

277 

278 def get_resp_file(self, tr): 

279 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % tr.nslc_id) 

280 if not os.path.exists(respfile): 

281 raise eventdata.NoRestitution( 

282 'no response information available for trace %s.%s.%s.%s' 

283 % tr.nslc_id) 

284 

285 return respfile 

286 

287 def get_stationxml(self): 

288 stations = self.get_pyrocko_stations() 

289 respfiles = [] 

290 for station in stations: 

291 for channel in station.get_channels(): 

292 nslc = station.nsl() + (channel.name,) 

293 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % nslc) 

294 respfiles.append(respfile) 

295 

296 from pyrocko.client.fdsn import resp 

297 sxml = resp.make_stationxml(stations, resp.iload(respfiles)) 

298 return sxml 

299 

300 def get_restitution(self, tr, allowed_methods): 

301 if 'evalresp' in allowed_methods: 

302 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % tr.nslc_id) 

303 if not os.path.exists(respfile): 

304 raise eventdata.NoRestitution( 

305 'no response information available for trace %s.%s.%s.%s' 

306 % tr.nslc_id) 

307 

308 trans = response.InverseEvalresp(respfile, tr) 

309 return trans 

310 else: 

311 raise eventdata.NoRestitution( 

312 'no response information available for trace %s.%s.%s.%s' 

313 % tr.nslc_id) 

314 

315 def get_pyrocko_response(self, tr, target): 

316 ''' 

317 Extract the frequency response as :py:class:`response.EvalResp` 

318 instance for *tr*. 

319 

320 :param tr: :py:class:`trace.Trace` instance 

321 ''' 

322 respfile = self.get_resp_file(tr) 

323 return response.Evalresp(respfile, tr, target) 

324 

325 def _unpack(self): 

326 input_fn = self.seedvolume 

327 output_dir = self.tempdir 

328 

329 def strerr(s): 

330 return '\n'.join(['rdseed: '+line.decode() 

331 for line in s.splitlines()]) 

332 try: 

333 

334 # seismograms: 

335 if self._pile is None: 

336 rdseed_proc = subprocess.Popen( 

337 [Programs.rdseed, '-f', input_fn, '-d', '-z', '3', '-o', 

338 '1', '-p', '-R', '-q', output_dir], 

339 stdin=subprocess.PIPE, 

340 stdout=subprocess.PIPE, 

341 stderr=subprocess.PIPE) 

342 

343 (out, err) = rdseed_proc.communicate() 

344 logging.info(strerr(err)) 

345 else: 

346 rdseed_proc = subprocess.Popen( 

347 [Programs.rdseed, '-f', input_fn, '-z', '3', '-p', '-R', 

348 '-q', output_dir], 

349 stdin=subprocess.PIPE, 

350 stdout=subprocess.PIPE, 

351 stderr=subprocess.PIPE) 

352 

353 (out, err) = rdseed_proc.communicate() 

354 logging.info(strerr(err)) 

355 

356 # event data: 

357 rdseed_proc = subprocess.Popen( 

358 [Programs.rdseed, '-f', input_fn, '-e', '-q', output_dir], 

359 stdin=subprocess.PIPE, 

360 stdout=subprocess.PIPE, 

361 stderr=subprocess.PIPE) 

362 

363 (out, err) = rdseed_proc.communicate() 

364 logging.info(strerr(err)) 

365 

366 # station summary information: 

367 rdseed_proc = subprocess.Popen( 

368 [Programs.rdseed, '-f', input_fn, '-S', '-q', output_dir], 

369 stdin=subprocess.PIPE, 

370 stdout=subprocess.PIPE, 

371 stderr=subprocess.PIPE) 

372 

373 (out, err) = rdseed_proc.communicate() 

374 logging.info(strerr(err)) 

375 

376 # station headers: 

377 rdseed_proc = subprocess.Popen( 

378 [Programs.rdseed, '-f', input_fn, '-s', '-q', output_dir], 

379 stdin=subprocess.PIPE, 

380 stdout=subprocess.PIPE, 

381 stderr=subprocess.PIPE) 

382 

383 (out, err) = rdseed_proc.communicate() 

384 with open(self.station_headers_file, 'w') as fout: 

385 fout.write(out.decode()) 

386 logging.info(strerr(err)) 

387 

388 except OSError as e: 

389 if e.errno == 2: 

390 reason = "Could not find executable: '%s'." % Programs.rdseed 

391 else: 

392 reason = str(e) 

393 

394 logging.fatal('Failed to unpack SEED volume. %s' % reason) 

395 sys.exit(1) 

396 

397 def _get_events_from_file(self): 

398 rdseed_event_file = os.path.join(self.tempdir, 'rdseed.events') 

399 if not os.path.isfile(rdseed_event_file): 

400 return [] 

401 

402 with open(rdseed_event_file, 'r') as f: 

403 events = [] 

404 for line in f: 

405 toks = line.split(', ') 

406 if len(toks) == 9: 

407 datetime = toks[1].split('.')[0] 

408 format = '%Y/%m/%d %H:%M:%S' 

409 secs = util.to_time_float(calendar.timegm( 

410 time.strptime(datetime, format))) 

411 e = model.Event( 

412 lat=float(toks[2]), 

413 lon=float(toks[3]), 

414 depth=float(toks[4])*1000., 

415 magnitude=float(toks[8]), 

416 time=secs) 

417 

418 events.append(e) 

419 else: 

420 raise Exception('Event description in unrecognized format') 

421 

422 return events 

423 

424 def _get_stations_from_file(self): 

425 stations = read_station_header_file(self.station_headers_file) 

426 return stations 

427 

428 def _insert_channel_descriptions(self, stations): 

429 # this is done beforehand in this class 

430 pass 

431 

432 

433def check_have_rdseed(): 

434 if not Programs.check(): 

435 raise ExternalProgramMissing( 

436 'rdseed is not installed or cannot be found') 

437 

438 

439if __name__ == '__main__': 

440 print(SeedVolumeAccess(sys.argv[1]).get_stationxml().dump_xml())