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

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

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

6'''

7Reader for RESP files.

8'''

10import time

11import re

12import logging

14from pyrocko import util, guts

15from pyrocko.io import io_common

16from pyrocko.io import stationxml as sxml

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

21class RespError(io_common.FileLoadError):

22 pass

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))

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)

47def pnc(s):

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

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

52def punit(s):

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

56def psymmetry(s):

57 return {

58 b'A': 'NONE',

59 b'B': 'ODD',

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

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.')

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.')

85def pblock_060(content):

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

87 return stage_number, None

90def pblock_053(content):

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

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())