1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import absolute_import 

6 

7import os 

8import sys 

9import subprocess 

10import tempfile 

11import calendar 

12import time 

13import re 

14import logging 

15import shutil 

16 

17from pyrocko import ExternalProgramMissing 

18from pyrocko import pile, model, util, response 

19from pyrocko.io import eventdata 

20 

21pjoin = os.path.join 

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

23 

24 

25def cmp(a, b): 

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

27 

28 

29def read_station_header_file(fn): 

30 

31 def m(i, *args): 

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

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

34 return True 

35 return False 

36 

37 def f(i): 

38 return float(toks[i]) 

39 

40 def s(i): 

41 if len(toks) > i: 

42 return toks[i] 

43 else: 

44 return '' 

45 

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

47 

48 stations = [] 

49 atsec, station, channel = None, None, None 

50 

51 for line in fi: 

52 toks = line.split() 

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

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

55 atsec = 'station' 

56 station = {'channels': []} 

57 stations.append(station) 

58 continue 

59 

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

61 atsec = 'channel' 

62 channel = {} 

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

64 continue 

65 

66 if atsec == 'station': 

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

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

69 

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

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

72 

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

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

75 

76 if atsec == 'channel': 

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

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

79 

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

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

82 

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

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

85 

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

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

88 

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

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

91 

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

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

94 

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

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

97 

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

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

100 

101 nsl_stations = {} 

102 for station in stations: 

103 for channel in station['channels']: 

104 def cs(k, default=None): 

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

106 

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

108 if nsl not in nsl_stations: 

109 nsl_stations[nsl] = model.Station( 

110 network=station['network'], 

111 station=station['station'], 

112 location=channel['location'], 

113 lat=cs('lat'), 

114 lon=cs('lon'), 

115 elevation=cs('elevation'), 

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

117 name=station['name']) 

118 

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

120 channel['channel'], 

121 azimuth=channel['azimuth'], 

122 dip=channel['dip'])) 

123 

124 return list(nsl_stations.values()) 

125 

126 

127def cmp_version(a, b): 

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

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

130 return cmp(ai, bi) 

131 

132 

133def dumb_parser(data): 

134 

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

136 

137 state = in_ws 

138 

139 rows = [] 

140 cols = [] 

141 accu = '' 

142 for c in data: 

143 if state == in_ws: 

144 if c == '"': 

145 new_state = in_str 

146 

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

148 new_state = in_kw 

149 

150 if state == in_kw: 

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

152 cols.append(accu) 

153 accu = '' 

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

155 rows.append(cols) 

156 cols = [] 

157 new_state = in_ws 

158 

159 if state == in_str: 

160 if c == '"': 

161 accu += c 

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

163 accu = '' 

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

165 rows.append(cols) 

166 cols = [] 

167 new_state = in_ws 

168 

169 state = new_state 

170 

171 if state in (in_kw, in_str): 

172 accu += c 

173 

174 if len(cols) != 0: 

175 rows.append(cols) 

176 

177 return rows 

178 

179 

180class Programs(object): 

181 rdseed = 'rdseed' 

182 avail = None 

183 

184 @staticmethod 

185 def check(): 

186 if Programs.avail is not None: 

187 return Programs.avail 

188 

189 else: 

190 try: 

191 rdseed_proc = subprocess.Popen( 

192 [Programs.rdseed], 

193 stdin=subprocess.PIPE, 

194 stdout=subprocess.PIPE, 

195 stderr=subprocess.PIPE) 

196 

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

198 

199 except OSError as e: 

200 if e.errno == 2: 

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

202 % Programs.rdseed 

203 else: 

204 reason = str(e) 

205 

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

207 Programs.avail = False 

208 return Programs.avail 

209 

210 ms = [re.search( 

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

212 for s in (err, out)] 

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

214 if not ms: 

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

216 else: 

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

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

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

220 

221 logger.warning( 

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

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

224 

225 Programs.avail = True 

226 return Programs.avail 

227 

228 

229class SeedVolumeNotFound(Exception): 

230 pass 

231 

232 

233class SeedVolumeAccess(eventdata.EventDataAccess): 

234 

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

236 

237 ''' 

238 Create new SEED Volume access object. 

239 

240 :param seedvolume: filename of seed volume 

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

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

243 instead of the data provided by the SEED volume. 

244 (This is useful for dataless SEED volumes.) 

245 ''' 

246 

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

248 self.tempdir = None 

249 Programs.check() 

250 

251 self.tempdir = None 

252 self.seedvolume = seedvolume 

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

254 raise SeedVolumeNotFound() 

255 

256 self.tempdir = tempfile.mkdtemp("", "SeedVolumeAccess-") 

257 self.station_headers_file = os.path.join( 

258 self.tempdir, 'station_header_infos') 

259 self._unpack() 

260 self.shutil = shutil 

261 

262 def __del__(self): 

263 if self.tempdir: 

264 self.shutil.rmtree(self.tempdir) 

265 

266 def get_pile(self): 

267 if self._pile is None: 

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

269 self._pile = pile.Pile() 

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

271 

272 return self._pile 

273 

274 def get_resp_file(self, tr): 

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

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

277 raise eventdata.NoRestitution( 

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

279 % tr.nslc_id) 

280 

281 return respfile 

282 

283 def get_stationxml(self): 

284 stations = self.get_pyrocko_stations() 

285 respfiles = [] 

286 for station in stations: 

287 for channel in station.get_channels(): 

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

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

290 respfiles.append(respfile) 

291 

292 from pyrocko.fdsn import resp 

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

294 return sxml 

295 

296 def get_restitution(self, tr, allowed_methods): 

297 if 'evalresp' in allowed_methods: 

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

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

300 raise eventdata.NoRestitution( 

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

302 % tr.nslc_id) 

303 

304 trans = response.InverseEvalresp(respfile, tr) 

305 return trans 

306 else: 

307 raise eventdata.NoRestitution( 

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

309 % tr.nslc_id) 

310 

311 def get_pyrocko_response(self, tr, target): 

312 ''' 

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

314 instance for *tr*. 

315 

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

317 ''' 

318 respfile = self.get_resp_file(tr) 

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

320 

321 def _unpack(self): 

322 input_fn = self.seedvolume 

323 output_dir = self.tempdir 

324 

325 def strerr(s): 

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

327 for line in s.splitlines()]) 

328 try: 

329 

330 # seismograms: 

331 if self._pile is None: 

332 rdseed_proc = subprocess.Popen( 

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

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

335 stdin=subprocess.PIPE, 

336 stdout=subprocess.PIPE, 

337 stderr=subprocess.PIPE) 

338 

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

340 logging.info(strerr(err)) 

341 else: 

342 rdseed_proc = subprocess.Popen( 

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

344 '-q', output_dir], 

345 stdin=subprocess.PIPE, 

346 stdout=subprocess.PIPE, 

347 stderr=subprocess.PIPE) 

348 

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

350 logging.info(strerr(err)) 

351 

352 # event data: 

353 rdseed_proc = subprocess.Popen( 

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

355 stdin=subprocess.PIPE, 

356 stdout=subprocess.PIPE, 

357 stderr=subprocess.PIPE) 

358 

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

360 logging.info(strerr(err)) 

361 

362 # station summary information: 

363 rdseed_proc = subprocess.Popen( 

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

365 stdin=subprocess.PIPE, 

366 stdout=subprocess.PIPE, 

367 stderr=subprocess.PIPE) 

368 

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

370 logging.info(strerr(err)) 

371 

372 # station headers: 

373 rdseed_proc = subprocess.Popen( 

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

375 stdin=subprocess.PIPE, 

376 stdout=subprocess.PIPE, 

377 stderr=subprocess.PIPE) 

378 

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

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

381 fout.write(out.decode()) 

382 logging.info(strerr(err)) 

383 

384 except OSError as e: 

385 if e.errno == 2: 

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

387 else: 

388 reason = str(e) 

389 

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

391 sys.exit(1) 

392 

393 def _get_events_from_file(self): 

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

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

396 return [] 

397 

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

399 events = [] 

400 for line in f: 

401 toks = line.split(', ') 

402 if len(toks) == 9: 

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

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

405 secs = util.to_time_float(calendar.timegm( 

406 time.strptime(datetime, format))) 

407 e = model.Event( 

408 lat=float(toks[2]), 

409 lon=float(toks[3]), 

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

411 magnitude=float(toks[8]), 

412 time=secs) 

413 

414 events.append(e) 

415 else: 

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

417 

418 return events 

419 

420 def _get_stations_from_file(self): 

421 stations = read_station_header_file(self.station_headers_file) 

422 return stations 

423 

424 def _insert_channel_descriptions(self, stations): 

425 # this is done beforehand in this class 

426 pass 

427 

428 

429def check_have_rdseed(): 

430 if not Programs.check(): 

431 raise ExternalProgramMissing( 

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

433 

434 

435if __name__ == '__main__': 

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