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

253 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Reader for RESP files. 

8''' 

9 

10import time 

11import re 

12import logging 

13 

14from pyrocko import util, guts 

15from pyrocko.io import io_common 

16from pyrocko.io import stationxml as sxml 

17 

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

19 

20 

21class RespError(io_common.FileLoadError): 

22 pass 

23 

24 

25def ppolezero(s): 

26 v = s.split() 

27 return sxml.PoleZero( 

28 number=int(v[0]), 

29 real=sxml.FloatNoUnit( 

30 value=float(v[1]), 

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

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

33 imaginary=sxml.FloatNoUnit( 

34 value=float(v[2]), 

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

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

37 

38 

39def pcfu(s): 

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

41 return sxml.FloatWithUnit( 

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

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

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

45 

46 

47def pnc(s): 

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

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

50 

51 

52def punit(s): 

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

54 

55 

56def psymmetry(s): 

57 return { 

58 b'A': 'NONE', 

59 b'B': 'ODD', 

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

61 

62 

63def ptftype(s): 

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

65 return 'LAPLACE (RADIANS/SECOND)' 

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

67 return 'LAPLACE (HERTZ)' 

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

69 return 'DIGITAL (Z-TRANSFORM)' 

70 else: 

71 raise RespError('Unknown PZ transfer function type.') 

72 

73 

74def pcftype(s): 

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

76 return 'ANALOG (RADIANS/SECOND)' 

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

78 return 'ANALOG (HERTZ)' 

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

80 return 'DIGITAL' 

81 else: 

82 raise RespError('Unknown cf transfer function type.') 

83 

84 

85def pblock_060(content): 

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

87 return stage_number, None 

88 

89 

90def pblock_053(content): 

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

92 

93 pzs = sxml.PolesZeros( 

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

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

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

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

98 normalization_frequency=sxml.Frequency( 

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

100 

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

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

103 

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

105 x.number = i 

106 

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

108 x.number = i 

109 

110 return stage_number, pzs 

111 

112 

113def pblock_043(content): 

114 stage_number = -1 

115 

116 pzs = sxml.PolesZeros( 

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

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

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

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

121 normalization_frequency=sxml.Frequency( 

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

123 

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

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

126 

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

128 x.number = i 

129 

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

131 x.number = i 

132 

133 return stage_number, pzs 

134 

135 

136def pblock_058(content): 

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

138 

139 gain = sxml.Gain( 

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

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

142 

143 return stage_number, gain 

144 

145 

146def pblock_048(content): 

147 stage_number = -1 

148 

149 gain = sxml.Gain( 

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

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

152 

153 return stage_number, gain 

154 

155 

156def pblock_054(content): 

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

158 

159 cfs = sxml.Coefficients( 

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

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

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

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

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

165 

166 return stage_number, cfs 

167 

168 

169def pblock_044(content): 

170 stage_number = -1 

171 

172 cfs = sxml.Coefficients( 

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

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

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

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

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

178 

179 return stage_number, cfs 

180 

181 

182def pblock_057(content): 

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

184 

185 deci = sxml.Decimation( 

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

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

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

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

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

191 

192 return stage_number, deci 

193 

194 

195def pblock_047(content): 

196 stage_number = -1 

197 

198 deci = sxml.Decimation( 

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

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

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

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

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

204 

205 return stage_number, deci 

206 

207 

208def pblock_061(content): 

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

210 

211 fir = sxml.FIR( 

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

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

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

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

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

217 

218 return stage_number, fir 

219 

220 

221def pblock_041(content): 

222 stage_number = -1 

223 

224 fir = sxml.FIR( 

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

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

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

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

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

230 

231 return stage_number, fir 

232 

233 

234bdefs = { 

235 b'050': { 

236 'name': 'Station Identifier Blockette', 

237 }, 

238 b'052': { 

239 'name': 'Channel Identifier Blockette', 

240 }, 

241 b'060': { 

242 'name': 'Response Reference Information', 

243 'parse': pblock_060, 

244 }, 

245 b'053': { 

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

247 'parse': pblock_053, 

248 }, 

249 b'043': { 

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

251 'parse': pblock_043, 

252 }, 

253 b'054': { 

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

255 'parse': pblock_054, 

256 }, 

257 b'044': { 

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

259 'parse': pblock_044, 

260 }, 

261 b'057': { 

262 'name': 'Decimation Blockette', 

263 'parse': pblock_057, 

264 }, 

265 b'047': { 

266 'name': 'Decimation Dictionary Blockette', 

267 'parse': pblock_047, 

268 }, 

269 b'058': { 

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

271 'parse': pblock_058, 

272 }, 

273 b'048': { 

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

275 'parse': pblock_048, 

276 }, 

277 b'061': { 

278 'name': 'FIR Response Blockette', 

279 'parse': pblock_061, 

280 }, 

281 b'041': { 

282 'name': 'FIR Dictionary Blockette', 

283 'parse': pblock_041, 

284 }, 

285} 

286 

287 

288def parse1(f): 

289 for line in f: 

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

291 m = re.match( 

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

293 if m: 

294 if m.group(2): 

295 pass 

296 

297 elif m.group(3): 

298 block = m.group(3) 

299 field = m.group(4) 

300 key = m.group(7) 

301 value = m.group(8) 

302 yield block, field, key, value 

303 

304 

305def parse2(f): 

306 current_b = None 

307 content = [] 

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

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

310 if current_b is not None: 

311 yield current_b, content 

312 

313 current_b = block 

314 content = [] 

315 

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

317 

318 if current_b is not None: 

319 yield current_b, content 

320 

321 

322def parse3(f): 

323 state = [None, None, []] 

324 for block, content in parse2(f): 

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

326 yield state 

327 state = [None, None, []] 

328 

329 if block == b'050': 

330 state[0] = content 

331 elif block == b'052': 

332 state[1] = content 

333 else: 

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

335 

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

337 yield state 

338 

339 

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

341 for field_, _, value in content: 

342 if field_ == field: 

343 return value 

344 else: 

345 if optional: 

346 return None 

347 elif default is not None: 

348 return default 

349 else: 

350 raise RespError('Key not found: %s' % field) 

351 

352 

353def getn(content, field): 

354 lst = [] 

355 for field_, _, value in content: 

356 if field_ == field: 

357 lst.append(value) 

358 return lst 

359 

360 

361def pdate(s): 

362 if len(s) < 17: 

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

364 

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

366 return None 

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

368 return None 

369 else: 

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

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

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

373 

374 return util.str_to_time( 

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

376 

377 

378def ploc(s): 

379 if s == b'??': 

380 return '' 

381 else: 

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

383 

384 

385def pcode(s): 

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

387 

388 

389def gett(lst, t): 

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

391 

392 

393def gett1o(lst, t): 

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

395 if len(lst) == 0: 

396 return None 

397 elif len(lst) == 1: 

398 return lst[0] 

399 else: 

400 raise RespError('Duplicate entry.') 

401 

402 

403def gett1(lst, t): 

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

405 if len(lst) == 0: 

406 raise RespError('Entry not found.') 

407 elif len(lst) == 1: 

408 return lst[0] 

409 else: 

410 raise RespError('Duplicate entry.') 

411 

412 

413class ChannelResponse(guts.Object): 

414 ''' 

415 Response information + channel codes and time span. 

416 ''' 

417 

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

419 start_date = guts.Timestamp.T() 

420 end_date = guts.Timestamp.T() 

421 response = sxml.Response.T() 

422 

423 

424def iload_fh(f): 

425 ''' 

426 Read RESP information from open file handle. 

427 ''' 

428 

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

430 nslc = ( 

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

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

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

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

435 

436 try: 

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

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

439 except util.TimeStrError as e: 

440 raise RespError('Invalid date in RESP information (%s).' % str(e)) 

441 

442 stage_elements = {} 

443 

444 istage = -1 

445 for block, content in rcs: 

446 if block not in bdefs: 

447 raise RespError('Unknown block type found: %s' % block) 

448 

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

450 if istage_temp != -1: 

451 istage = istage_temp 

452 

453 if x is None: 

454 continue 

455 

456 x.validate() 

457 if istage not in stage_elements: 

458 stage_elements[istage] = [] 

459 

460 stage_elements[istage].append(x) 

461 

462 istages = sorted(stage_elements.keys()) 

463 stages = [] 

464 totalgain = None 

465 for istage in istages: 

466 elements = stage_elements[istage] 

467 if istage == 0: 

468 totalgain = gett1(elements, sxml.Gain) 

469 else: 

470 stage = sxml.ResponseStage( 

471 number=istage, 

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

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

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

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

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

477 

478 stages.append(stage) 

479 

480 if stages: 

481 resp = sxml.Response( 

482 stage_list=stages) 

483 

484 if totalgain: 

485 totalgain_value = totalgain.value 

486 totalgain_frequency = totalgain.frequency 

487 

488 else: 

489 totalgain_value = 1. 

490 gain_frequencies = [] 

491 for stage in stages: 

492 totalgain_value *= stage.stage_gain.value 

493 gain_frequencies.append(stage.stage_gain.frequency) 

494 

495 totalgain_frequency = gain_frequencies[0] 

496 

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

498 logger.warning( 

499 'No total gain reported and inconsistent gain ' 

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

501 'omitting total gain and frequency from created ' 

502 'instrument sensitivity object.' % nslc) 

503 

504 totalgain_value = None 

505 totalgain_frequency = None 

506 

507 resp.instrument_sensitivity = sxml.Sensitivity( 

508 value=totalgain_value, 

509 frequency=totalgain_frequency, 

510 input_units=stages[0].input_units, 

511 output_units=stages[-1].output_units) 

512 

513 yield ChannelResponse( 

514 codes=nslc, 

515 start_date=tmin, 

516 end_date=tmax, 

517 response=resp) 

518 

519 else: 

520 logger.warning( 

521 'Incomplete response information for %s (%s - %s).', 

522 '.'.join(nslc), 

523 util.time_to_str(tmin), 

524 util.time_to_str(tmax)) 

525 

526 yield ChannelResponse( 

527 codes=nslc, 

528 start_date=tmin, 

529 end_date=tmax, 

530 response=None) 

531 

532 

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

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

535 

536 

537def make_stationxml(pyrocko_stations, channel_responses): 

538 ''' 

539 Create stationxml from pyrocko station list and RESP information. 

540 

541 :param pyrocko_stations: 

542 Station information. 

543 :type pyrocko_stations: 

544 :py:class:`list` of :py:class:`pyrocko.model.station.Station` 

545 

546 :param channel_responses: 

547 iterable yielding :py:class:`ChannelResponse` objects 

548 

549 :returns: :py:class:`pyrocko.io.stationxml.FDSNStationXML` object with 

550 merged information 

551 

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

553 is skipped and a warning is emitted. 

554 ''' 

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

556 networks = {} 

557 stations = {} 

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

559 pstation = pstations[net, sta, loc] 

560 if net not in networks: 

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

562 

563 if (net, sta) not in stations: 

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

565 code=sta, 

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

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

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

569 

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

571 

572 for cr in channel_responses: 

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

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

575 pstation = pstations[net, sta, loc] 

576 pchannel = pstation.get_channel(cha) 

577 extra = {} 

578 if pchannel is not None: 

579 if pchannel.azimuth is not None: 

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

581 

582 if pchannel.dip is not None: 

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

584 

585 channel = sxml.Channel( 

586 code=cha, 

587 location_code=loc, 

588 start_date=cr.start_date, 

589 end_date=cr.end_date, 

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

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

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

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

594 response=cr.response, 

595 **extra) 

596 

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

598 else: 

599 logger.warning('No station information for %s.%s.%s.' % 

600 (net, sta, loc)) 

601 

602 for station in stations.values(): 

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

604 

605 return sxml.FDSNStationXML( 

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

607 created=time.time(), 

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

609 

610 

611if __name__ == '__main__': 

612 import sys 

613 from pyrocko.model.station import load_stations 

614 

615 util.setup_logging(__name__) 

616 

617 if len(sys.argv) < 2: 

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

619 

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

621 

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

623 

624 print(sxml.dump_xml())