1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5from __future__ import absolute_import, division, print_function 

6 

7import time 

8import re 

9import logging 

10 

11from pyrocko import util, guts 

12from pyrocko.io import io_common 

13from pyrocko.io import stationxml as sxml 

14 

15logger = logging.getLogger('pyrocko.io.resp') 

16 

17 

18class RespError(io_common.FileLoadError): 

19 pass 

20 

21 

22def ppolezero(s): 

23 v = s.split() 

24 return sxml.PoleZero( 

25 number=int(v[0]), 

26 real=sxml.FloatNoUnit( 

27 value=float(v[1]), 

28 plus_error=float(v[3]) or None, 

29 minus_error=float(v[3]) or None), 

30 imaginary=sxml.FloatNoUnit( 

31 value=float(v[2]), 

32 plus_error=float(v[4]) or None, 

33 minus_error=float(v[4]) or None)) 

34 

35 

36def pcfu(s): 

37 v = list(map(float, s.split())) 

38 return sxml.FloatWithUnit( 

39 value=float(v[-2]), 

40 plus_error=float(v[-1]) or None, 

41 minus_error=float(v[-1]) or None) 

42 

43 

44def pnc(s): 

45 v = list(map(float, s.split())) 

46 return sxml.NumeratorCoefficient(i=int(v[0]), value=float(v[1])) 

47 

48 

49def punit(s): 

50 return str(s.split()[0].decode('ascii')) 

51 

52 

53def psymmetry(s): 

54 return { 

55 b'A': 'NONE', 

56 b'B': 'ODD', 

57 b'C': 'EVEN'}[s.upper()] 

58 

59 

60def ptftype(s): 

61 if s.startswith(b'A'): 

62 return 'LAPLACE (RADIANS/SECOND)' 

63 elif s.startswith(b'B'): 

64 return 'LAPLACE (HERTZ)' 

65 elif s.startswith(b'D'): 

66 return 'DIGITAL (Z-TRANSFORM)' 

67 else: 

68 raise RespError('unknown pz transfer function type') 

69 

70 

71def pcftype(s): 

72 if s.startswith(b'A'): 

73 return 'ANALOG (RADIANS/SECOND)' 

74 elif s.startswith(b'B'): 

75 return 'ANALOG (HERTZ)' 

76 elif s.startswith(b'D'): 

77 return 'DIGITAL' 

78 else: 

79 raise RespError('unknown cf transfer function type') 

80 

81 

82def pblock_060(content): 

83 stage_number = int(get1(content, b'04')) 

84 return stage_number, None 

85 

86 

87def pblock_053(content): 

88 stage_number = int(get1(content, b'04')) 

89 

90 pzs = sxml.PolesZeros( 

91 pz_transfer_function_type=ptftype(get1(content, b'03')), 

92 input_units=sxml.Units(name=punit(get1(content, b'05'))), 

93 output_units=sxml.Units(name=punit(get1(content, b'06'))), 

94 normalization_factor=float(get1(content, b'07')), 

95 normalization_frequency=sxml.Frequency( 

96 value=float(get1(content, b'08'))), 

97 

98 zero_list=list(map(ppolezero, getn(content, b'10-13'))), 

99 pole_list=list(map(ppolezero, getn(content, b'15-18')))) 

100 

101 for i, x in enumerate(pzs.zero_list): 

102 x.number = i 

103 

104 for i, x in enumerate(pzs.pole_list): 

105 x.number = i 

106 

107 return stage_number, pzs 

108 

109 

110def pblock_043(content): 

111 stage_number = -1 

112 

113 pzs = sxml.PolesZeros( 

114 pz_transfer_function_type=ptftype(get1(content, b'05')), 

115 input_units=sxml.Units(name=punit(get1(content, b'06'))), 

116 output_units=sxml.Units(name=punit(get1(content, b'07'))), 

117 normalization_factor=float(get1(content, b'08')), 

118 normalization_frequency=sxml.Frequency( 

119 value=float(get1(content, b'09'))), 

120 

121 zero_list=list(map(ppolezero, getn(content, b'11-14'))), 

122 pole_list=list(map(ppolezero, getn(content, b'16-19')))) 

123 

124 for i, x in enumerate(pzs.zero_list): 

125 x.number = i 

126 

127 for i, x in enumerate(pzs.pole_list): 

128 x.number = i 

129 

130 return stage_number, pzs 

131 

132 

133def pblock_058(content): 

134 stage_number = int(get1(content, b'03')) 

135 

136 gain = sxml.Gain( 

137 value=float(get1(content, b'04')), 

138 frequency=float(get1(content, b'05').split()[0])) 

139 

140 return stage_number, gain 

141 

142 

143def pblock_048(content): 

144 stage_number = -1 

145 

146 gain = sxml.Gain( 

147 value=float(get1(content, b'05')), 

148 frequency=float(get1(content, b'06').split()[0])) 

149 

150 return stage_number, gain 

151 

152 

153def pblock_054(content): 

154 stage_number = int(get1(content, b'04')) 

155 

156 cfs = sxml.Coefficients( 

157 cf_transfer_function_type=pcftype(get1(content, b'03')), 

158 input_units=sxml.Units(name=punit(get1(content, b'05'))), 

159 output_units=sxml.Units(name=punit(get1(content, b'06'))), 

160 numerator_list=list(map(pcfu, getn(content, b'08-09'))), 

161 denominator_list=list(map(pcfu, getn(content, b'11-12')))) 

162 

163 return stage_number, cfs 

164 

165 

166def pblock_044(content): 

167 stage_number = -1 

168 

169 cfs = sxml.Coefficients( 

170 cf_transfer_function_type=pcftype(get1(content, b'05')), 

171 input_units=sxml.Units(name=punit(get1(content, b'06'))), 

172 output_units=sxml.Units(name=punit(get1(content, b'07'))), 

173 numerator_list=list(map(pcfu, getn(content, b'09-10'))), 

174 denominator_list=list(map(pcfu, getn(content, b'12-13')))) 

175 

176 return stage_number, cfs 

177 

178 

179def pblock_057(content): 

180 stage_number = int(get1(content, b'03')) 

181 

182 deci = sxml.Decimation( 

183 input_sample_rate=sxml.Frequency(value=float(get1(content, b'04'))), 

184 factor=int(get1(content, b'05')), 

185 offset=int(get1(content, b'06')), 

186 delay=sxml.FloatWithUnit(value=float(get1(content, b'07'))), 

187 correction=sxml.FloatWithUnit(value=float(get1(content, b'08')))) 

188 

189 return stage_number, deci 

190 

191 

192def pblock_047(content): 

193 stage_number = -1 

194 

195 deci = sxml.Decimation( 

196 input_sample_rate=sxml.Frequency(value=float(get1(content, b'05'))), 

197 factor=int(get1(content, b'06')), 

198 offset=int(get1(content, b'07')), 

199 delay=sxml.FloatWithUnit(value=float(get1(content, b'08'))), 

200 correction=sxml.FloatWithUnit(value=float(get1(content, b'09')))) 

201 

202 return stage_number, deci 

203 

204 

205def pblock_061(content): 

206 stage_number = int(get1(content, b'03')) 

207 

208 fir = sxml.FIR( 

209 name=get1(content, b'04', optional=True), 

210 input_units=sxml.Units(name=punit(get1(content, b'06'))), 

211 output_units=sxml.Units(name=punit(get1(content, b'07'))), 

212 symmetry=psymmetry(get1(content, b'05')), 

213 numerator_coefficient_list=list(map(pnc, getn(content, b'09')))) 

214 

215 return stage_number, fir 

216 

217 

218def pblock_041(content): 

219 stage_number = -1 

220 

221 fir = sxml.FIR( 

222 name=get1(content, b'04', optional=True), 

223 input_units=sxml.Units(name=punit(get1(content, b'06'))), 

224 output_units=sxml.Units(name=punit(get1(content, b'07'))), 

225 symmetry=psymmetry(get1(content, b'05')), 

226 numerator_coefficient_list=list(map(pnc, getn(content, b'09')))) 

227 

228 return stage_number, fir 

229 

230 

231bdefs = { 

232 b'050': { 

233 'name': 'Station Identifier Blockette', 

234 }, 

235 b'052': { 

236 'name': 'Channel Identifier Blockette', 

237 }, 

238 b'060': { 

239 'name': 'Response Reference Information', 

240 'parse': pblock_060, 

241 }, 

242 b'053': { 

243 'name': 'Response (Poles & Zeros) Blockette', 

244 'parse': pblock_053, 

245 }, 

246 b'043': { 

247 'name': 'Response (Poles & Zeros) Dictionary Blockette', 

248 'parse': pblock_043, 

249 }, 

250 b'054': { 

251 'name': 'Response (Coefficients) Blockette', 

252 'parse': pblock_054, 

253 }, 

254 b'044': { 

255 'name': 'Response (Coefficients) Dictionary Blockette', 

256 'parse': pblock_044, 

257 }, 

258 b'057': { 

259 'name': 'Decimation Blockette', 

260 'parse': pblock_057, 

261 }, 

262 b'047': { 

263 'name': 'Decimation Dictionary Blockette', 

264 'parse': pblock_047, 

265 }, 

266 b'058': { 

267 'name': 'Channel Sensitivity/Gain Blockette', 

268 'parse': pblock_058, 

269 }, 

270 b'048': { 

271 'name': 'Channel Sensitivity/Gain Dictionary Blockette', 

272 'parse': pblock_048, 

273 }, 

274 b'061': { 

275 'name': 'FIR Response Blockette', 

276 'parse': pblock_061, 

277 }, 

278 b'041': { 

279 'name': 'FIR Dictionary Blockette', 

280 'parse': pblock_041, 

281 }, 

282} 

283 

284 

285def parse1(f): 

286 for line in f: 

287 line = line.rstrip(b'\r\n') 

288 m = re.match( 

289 br'\s*(#(.+)|B(\d\d\d)F(\d\d(-\d\d)?)\s+(([^:]+):\s*)?(.*))', line) 

290 if m: 

291 if m.group(2): 

292 pass 

293 

294 elif m.group(3): 

295 block = m.group(3) 

296 field = m.group(4) 

297 key = m.group(7) 

298 value = m.group(8) 

299 yield block, field, key, value 

300 

301 

302def parse2(f): 

303 current_b = None 

304 content = [] 

305 for block, field, key, value in parse1(f): 

306 if current_b != block or field == b'03': 

307 if current_b is not None: 

308 yield current_b, content 

309 

310 current_b = block 

311 content = [] 

312 

313 content.append((field, key, value)) 

314 

315 if current_b is not None: 

316 yield current_b, content 

317 

318 

319def parse3(f): 

320 state = [None, None, []] 

321 for block, content in parse2(f): 

322 if block == b'050' and state[0] and state[1]: 

323 yield state 

324 state = [None, None, []] 

325 

326 if block == b'050': 

327 state[0] = content 

328 elif block == b'052': 

329 state[1] = content 

330 else: 

331 state[2].append((block, content)) 

332 

333 if state[0] and state[1]: 

334 yield state 

335 

336 

337def get1(content, field, default=None, optional=False): 

338 for field_, _, value in content: 

339 if field_ == field: 

340 return value 

341 else: 

342 if optional: 

343 return None 

344 elif default is not None: 

345 return default 

346 else: 

347 raise RespError('key not found: %s' % field) 

348 

349 

350def getn(content, field): 

351 lst = [] 

352 for field_, _, value in content: 

353 if field_ == field: 

354 lst.append(value) 

355 return lst 

356 

357 

358def pdate(s): 

359 if len(s) < 17: 

360 s += b'0000,001,00:00:00'[len(s):] 

361 

362 if s.startswith(b'2599') or s.startswith(b'2999'): 

363 return None 

364 elif s.lower().startswith(b'no'): 

365 return None 

366 else: 

367 t = s.split(b',') 

368 if len(t) > 2 and t[1] == b'000': 

369 s = b','.join([t[0], b'001'] + t[2:]) 

370 

371 return util.str_to_time( 

372 str(s.decode('ascii')), format='%Y,%j,%H:%M:%S.OPTFRAC') 

373 

374 

375def ploc(s): 

376 if s == b'??': 

377 return '' 

378 else: 

379 return str(s.decode('ascii')) 

380 

381 

382def pcode(s): 

383 return str(s.decode('ascii')) 

384 

385 

386def gett(lst, t): 

387 return [x for x in lst if isinstance(x, t)] 

388 

389 

390def gett1o(lst, t): 

391 lst = [x for x in lst if isinstance(x, t)] 

392 if len(lst) == 0: 

393 return None 

394 elif len(lst) == 1: 

395 return lst[0] 

396 else: 

397 raise RespError('duplicate entry') 

398 

399 

400def gett1(lst, t): 

401 lst = [x for x in lst if isinstance(x, t)] 

402 if len(lst) == 0: 

403 raise RespError('entry not found') 

404 elif len(lst) == 1: 

405 return lst[0] 

406 else: 

407 raise RespError('duplicate entry') 

408 

409 

410class ChannelResponse(guts.Object): 

411 ''' 

412 Response information + channel codes and time span. 

413 ''' 

414 

415 codes = guts.Tuple.T(4, guts.String.T(default='')) 

416 start_date = guts.Timestamp.T() 

417 end_date = guts.Timestamp.T() 

418 response = sxml.Response.T() 

419 

420 

421def iload_fh(f): 

422 ''' 

423 Read RESP information from open file handle. 

424 ''' 

425 

426 for sc, cc, rcs in parse3(f): 

427 nslc = ( 

428 pcode(get1(sc, b'16')), 

429 pcode(get1(sc, b'03')), 

430 ploc(get1(cc, b'03', b'')), 

431 pcode(get1(cc, b'04'))) 

432 

433 try: 

434 tmin = pdate(get1(cc, b'22')) 

435 tmax = pdate(get1(cc, b'23')) 

436 except util.TimeStrError as e: 

437 raise RespError('invalid date in RESP information. (%s)' % str(e)) 

438 

439 stage_elements = {} 

440 

441 istage = -1 

442 for block, content in rcs: 

443 if block not in bdefs: 

444 raise RespError('unknown block type found: %s' % block) 

445 

446 istage_temp, x = bdefs[block]['parse'](content) 

447 if istage_temp != -1: 

448 istage = istage_temp 

449 

450 if x is None: 

451 continue 

452 

453 x.validate() 

454 if istage not in stage_elements: 

455 stage_elements[istage] = [] 

456 

457 stage_elements[istage].append(x) 

458 

459 istages = sorted(stage_elements.keys()) 

460 stages = [] 

461 totalgain = None 

462 for istage in istages: 

463 elements = stage_elements[istage] 

464 if istage == 0: 

465 totalgain = gett1(elements, sxml.Gain) 

466 else: 

467 stage = sxml.ResponseStage( 

468 number=istage, 

469 poles_zeros_list=gett(elements, sxml.PolesZeros), 

470 coefficients_list=gett(elements, sxml.Coefficients), 

471 fir=gett1o(elements, sxml.FIR), 

472 decimation=gett1o(elements, sxml.Decimation), 

473 stage_gain=gett1o(elements, sxml.Gain)) 

474 

475 stages.append(stage) 

476 

477 if stages: 

478 resp = sxml.Response( 

479 stage_list=stages) 

480 

481 if totalgain: 

482 totalgain_value = totalgain.value 

483 totalgain_frequency = totalgain.frequency 

484 

485 else: 

486 totalgain_value = 1. 

487 gain_frequencies = [] 

488 for stage in stages: 

489 totalgain_value *= stage.stage_gain.value 

490 gain_frequencies.append(stage.stage_gain.frequency) 

491 

492 totalgain_frequency = gain_frequencies[0] 

493 

494 if not all(f == totalgain_frequency for f in gain_frequencies): 

495 logger.warning( 

496 'no total gain reported and inconsistent gain ' 

497 'frequency values found in resp file for %s.%s.%s.%s: ' 

498 'omitting total gain and frequency from created ' 

499 'instrument sensitivity object' % nslc) 

500 

501 totalgain_value = None 

502 totalgain_frequency = None 

503 

504 resp.instrument_sensitivity = sxml.Sensitivity( 

505 value=totalgain_value, 

506 frequency=totalgain_frequency, 

507 input_units=stages[0].input_units, 

508 output_units=stages[-1].output_units) 

509 

510 yield ChannelResponse( 

511 codes=nslc, 

512 start_date=tmin, 

513 end_date=tmax, 

514 response=resp) 

515 

516 else: 

517 raise RespError('incomplete response information') 

518 

519 

520iload_filename, iload_dirname, iload_glob, iload = util.make_iload_family( 

521 iload_fh, 'RESP', ':py:class:`ChannelResponse`') 

522 

523 

524def make_stationxml(pyrocko_stations, channel_responses): 

525 ''' 

526 Create stationxml from pyrocko station list and RESP information. 

527 

528 :param pyrocko_stations: list of :py:class:`pyrocko.model.Station` objects 

529 :param channel_responses: iterable yielding :py:class:`ChannelResponse` 

530 objects 

531 :returns: :py:class:`pyrocko.fdsn.station.FDSNStationXML` object with 

532 merged information 

533 

534 If no station information is available for any response information, it 

535 is skipped and a warning is emitted. 

536 ''' 

537 pstations = dict((s.nsl(), s) for s in pyrocko_stations) 

538 networks = {} 

539 stations = {} 

540 for (net, sta, loc) in sorted(pstations.keys()): 

541 pstation = pstations[net, sta, loc] 

542 if net not in networks: 

543 networks[net] = sxml.Network(code=net) 

544 

545 if (net, sta) not in stations: 

546 stations[net, sta] = sxml.Station( 

547 code=sta, 

548 latitude=sxml.Latitude(pstation.lat), 

549 longitude=sxml.Longitude(pstation.lon), 

550 elevation=sxml.Distance(pstation.elevation)) 

551 

552 networks[net].station_list.append(stations[net, sta]) 

553 

554 for cr in channel_responses: 

555 net, sta, loc, cha = cr.codes 

556 if (net, sta, loc) in pstations: 

557 pstation = pstations[net, sta, loc] 

558 pchannel = pstation.get_channel(cha) 

559 extra = {} 

560 if pchannel is not None: 

561 if pchannel.azimuth is not None: 

562 extra['azimuth'] = sxml.Azimuth(pchannel.azimuth) 

563 

564 if pchannel.dip is not None: 

565 extra['dip'] = sxml.Dip(pchannel.dip) 

566 

567 channel = sxml.Channel( 

568 code=cha, 

569 location_code=loc, 

570 start_date=cr.start_date, 

571 end_date=cr.end_date, 

572 latitude=sxml.Latitude(pstation.lat), 

573 longitude=sxml.Longitude(pstation.lon), 

574 elevation=sxml.Distance(pstation.elevation), 

575 depth=sxml.Distance(pstation.depth), 

576 response=cr.response, 

577 **extra) 

578 

579 stations[net, sta].channel_list.append(channel) 

580 else: 

581 logger.warning('no station information for %s.%s.%s' % 

582 (net, sta, loc)) 

583 

584 for station in stations.values(): 

585 station.channel_list.sort(key=lambda c: (c.location_code, c.code)) 

586 

587 return sxml.FDSNStationXML( 

588 source='Converted from Pyrocko stations file and RESP information', 

589 created=time.time(), 

590 network_list=[networks[net_] for net_ in sorted(networks.keys())]) 

591 

592 

593if __name__ == '__main__': 

594 import sys 

595 from pyrocko.model.station import load_stations 

596 

597 util.setup_logging(__name__) 

598 

599 if len(sys.argv) < 2: 

600 sys.exit('usage: python -m pyrocko.fdsn.resp <stations> <resp> ...') 

601 

602 stations = load_stations(sys.argv[1]) 

603 

604 sxml = make_stationxml(stations, iload(sys.argv[2:])) 

605 

606 print(sxml.dump_xml())