1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import time
7import re
8import logging
10from pyrocko import util, guts
11from pyrocko.io import io_common
12from pyrocko.io import stationxml as sxml
14logger = logging.getLogger('pyrocko.io.resp')
17class RespError(io_common.FileLoadError):
18 pass
21def ppolezero(s):
22 v = s.split()
23 return sxml.PoleZero(
24 number=int(v[0]),
25 real=sxml.FloatNoUnit(
26 value=float(v[1]),
27 plus_error=float(v[3]) or None,
28 minus_error=float(v[3]) or None),
29 imaginary=sxml.FloatNoUnit(
30 value=float(v[2]),
31 plus_error=float(v[4]) or None,
32 minus_error=float(v[4]) or None))
35def pcfu(s):
36 v = list(map(float, s.split()))
37 return sxml.FloatWithUnit(
38 value=float(v[-2]),
39 plus_error=float(v[-1]) or None,
40 minus_error=float(v[-1]) or None)
43def pnc(s):
44 v = list(map(float, s.split()))
45 return sxml.NumeratorCoefficient(i=int(v[0]), value=float(v[1]))
48def punit(s):
49 return str(s.split()[0].decode('ascii'))
52def psymmetry(s):
53 return {
54 b'A': 'NONE',
55 b'B': 'ODD',
56 b'C': 'EVEN'}[s.upper()]
59def ptftype(s):
60 if s.startswith(b'A'):
61 return 'LAPLACE (RADIANS/SECOND)'
62 elif s.startswith(b'B'):
63 return 'LAPLACE (HERTZ)'
64 elif s.startswith(b'D'):
65 return 'DIGITAL (Z-TRANSFORM)'
66 else:
67 raise RespError('Unknown PZ transfer function type.')
70def pcftype(s):
71 if s.startswith(b'A'):
72 return 'ANALOG (RADIANS/SECOND)'
73 elif s.startswith(b'B'):
74 return 'ANALOG (HERTZ)'
75 elif s.startswith(b'D'):
76 return 'DIGITAL'
77 else:
78 raise RespError('Unknown cf transfer function type.')
81def pblock_060(content):
82 stage_number = int(get1(content, b'04'))
83 return stage_number, None
86def pblock_053(content):
87 stage_number = int(get1(content, b'04'))
89 pzs = sxml.PolesZeros(
90 pz_transfer_function_type=ptftype(get1(content, b'03')),
91 input_units=sxml.Units(name=punit(get1(content, b'05'))),
92 output_units=sxml.Units(name=punit(get1(content, b'06'))),
93 normalization_factor=float(get1(content, b'07')),
94 normalization_frequency=sxml.Frequency(
95 value=float(get1(content, b'08'))),
97 zero_list=list(map(ppolezero, getn(content, b'10-13'))),
98 pole_list=list(map(ppolezero, getn(content, b'15-18'))))
100 for i, x in enumerate(pzs.zero_list):
101 x.number = i
103 for i, x in enumerate(pzs.pole_list):
104 x.number = i
106 return stage_number, pzs
109def pblock_043(content):
110 stage_number = -1
112 pzs = sxml.PolesZeros(
113 pz_transfer_function_type=ptftype(get1(content, b'05')),
114 input_units=sxml.Units(name=punit(get1(content, b'06'))),
115 output_units=sxml.Units(name=punit(get1(content, b'07'))),
116 normalization_factor=float(get1(content, b'08')),
117 normalization_frequency=sxml.Frequency(
118 value=float(get1(content, b'09'))),
120 zero_list=list(map(ppolezero, getn(content, b'11-14'))),
121 pole_list=list(map(ppolezero, getn(content, b'16-19'))))
123 for i, x in enumerate(pzs.zero_list):
124 x.number = i
126 for i, x in enumerate(pzs.pole_list):
127 x.number = i
129 return stage_number, pzs
132def pblock_058(content):
133 stage_number = int(get1(content, b'03'))
135 gain = sxml.Gain(
136 value=float(get1(content, b'04')),
137 frequency=float(get1(content, b'05').split()[0]))
139 return stage_number, gain
142def pblock_048(content):
143 stage_number = -1
145 gain = sxml.Gain(
146 value=float(get1(content, b'05')),
147 frequency=float(get1(content, b'06').split()[0]))
149 return stage_number, gain
152def pblock_054(content):
153 stage_number = int(get1(content, b'04'))
155 cfs = sxml.Coefficients(
156 cf_transfer_function_type=pcftype(get1(content, b'03')),
157 input_units=sxml.Units(name=punit(get1(content, b'05'))),
158 output_units=sxml.Units(name=punit(get1(content, b'06'))),
159 numerator_list=list(map(pcfu, getn(content, b'08-09'))),
160 denominator_list=list(map(pcfu, getn(content, b'11-12'))))
162 return stage_number, cfs
165def pblock_044(content):
166 stage_number = -1
168 cfs = sxml.Coefficients(
169 cf_transfer_function_type=pcftype(get1(content, b'05')),
170 input_units=sxml.Units(name=punit(get1(content, b'06'))),
171 output_units=sxml.Units(name=punit(get1(content, b'07'))),
172 numerator_list=list(map(pcfu, getn(content, b'09-10'))),
173 denominator_list=list(map(pcfu, getn(content, b'12-13'))))
175 return stage_number, cfs
178def pblock_057(content):
179 stage_number = int(get1(content, b'03'))
181 deci = sxml.Decimation(
182 input_sample_rate=sxml.Frequency(value=float(get1(content, b'04'))),
183 factor=int(get1(content, b'05')),
184 offset=int(get1(content, b'06')),
185 delay=sxml.FloatWithUnit(value=float(get1(content, b'07'))),
186 correction=sxml.FloatWithUnit(value=float(get1(content, b'08'))))
188 return stage_number, deci
191def pblock_047(content):
192 stage_number = -1
194 deci = sxml.Decimation(
195 input_sample_rate=sxml.Frequency(value=float(get1(content, b'05'))),
196 factor=int(get1(content, b'06')),
197 offset=int(get1(content, b'07')),
198 delay=sxml.FloatWithUnit(value=float(get1(content, b'08'))),
199 correction=sxml.FloatWithUnit(value=float(get1(content, b'09'))))
201 return stage_number, deci
204def pblock_061(content):
205 stage_number = int(get1(content, b'03'))
207 fir = sxml.FIR(
208 name=get1(content, b'04', optional=True),
209 input_units=sxml.Units(name=punit(get1(content, b'06'))),
210 output_units=sxml.Units(name=punit(get1(content, b'07'))),
211 symmetry=psymmetry(get1(content, b'05')),
212 numerator_coefficient_list=list(map(pnc, getn(content, b'09'))))
214 return stage_number, fir
217def pblock_041(content):
218 stage_number = -1
220 fir = sxml.FIR(
221 name=get1(content, b'04', optional=True),
222 input_units=sxml.Units(name=punit(get1(content, b'06'))),
223 output_units=sxml.Units(name=punit(get1(content, b'07'))),
224 symmetry=psymmetry(get1(content, b'05')),
225 numerator_coefficient_list=list(map(pnc, getn(content, b'09'))))
227 return stage_number, fir
230bdefs = {
231 b'050': {
232 'name': 'Station Identifier Blockette',
233 },
234 b'052': {
235 'name': 'Channel Identifier Blockette',
236 },
237 b'060': {
238 'name': 'Response Reference Information',
239 'parse': pblock_060,
240 },
241 b'053': {
242 'name': 'Response (Poles & Zeros) Blockette',
243 'parse': pblock_053,
244 },
245 b'043': {
246 'name': 'Response (Poles & Zeros) Dictionary Blockette',
247 'parse': pblock_043,
248 },
249 b'054': {
250 'name': 'Response (Coefficients) Blockette',
251 'parse': pblock_054,
252 },
253 b'044': {
254 'name': 'Response (Coefficients) Dictionary Blockette',
255 'parse': pblock_044,
256 },
257 b'057': {
258 'name': 'Decimation Blockette',
259 'parse': pblock_057,
260 },
261 b'047': {
262 'name': 'Decimation Dictionary Blockette',
263 'parse': pblock_047,
264 },
265 b'058': {
266 'name': 'Channel Sensitivity/Gain Blockette',
267 'parse': pblock_058,
268 },
269 b'048': {
270 'name': 'Channel Sensitivity/Gain Dictionary Blockette',
271 'parse': pblock_048,
272 },
273 b'061': {
274 'name': 'FIR Response Blockette',
275 'parse': pblock_061,
276 },
277 b'041': {
278 'name': 'FIR Dictionary Blockette',
279 'parse': pblock_041,
280 },
281}
284def parse1(f):
285 for line in f:
286 line = line.rstrip(b'\r\n')
287 m = re.match(
288 br'\s*(#(.+)|B(\d\d\d)F(\d\d(-\d\d)?)\s+(([^:]+):\s*)?(.*))', line)
289 if m:
290 if m.group(2):
291 pass
293 elif m.group(3):
294 block = m.group(3)
295 field = m.group(4)
296 key = m.group(7)
297 value = m.group(8)
298 yield block, field, key, value
301def parse2(f):
302 current_b = None
303 content = []
304 for block, field, key, value in parse1(f):
305 if current_b != block or field == b'03':
306 if current_b is not None:
307 yield current_b, content
309 current_b = block
310 content = []
312 content.append((field, key, value))
314 if current_b is not None:
315 yield current_b, content
318def parse3(f):
319 state = [None, None, []]
320 for block, content in parse2(f):
321 if block == b'050' and state[0] and state[1]:
322 yield state
323 state = [None, None, []]
325 if block == b'050':
326 state[0] = content
327 elif block == b'052':
328 state[1] = content
329 else:
330 state[2].append((block, content))
332 if state[0] and state[1]:
333 yield state
336def get1(content, field, default=None, optional=False):
337 for field_, _, value in content:
338 if field_ == field:
339 return value
340 else:
341 if optional:
342 return None
343 elif default is not None:
344 return default
345 else:
346 raise RespError('Key not found: %s' % field)
349def getn(content, field):
350 lst = []
351 for field_, _, value in content:
352 if field_ == field:
353 lst.append(value)
354 return lst
357def pdate(s):
358 if len(s) < 17:
359 s += b'0000,001,00:00:00'[len(s):]
361 if s.startswith(b'2599') or s.startswith(b'2999'):
362 return None
363 elif s.lower().startswith(b'no'):
364 return None
365 else:
366 t = s.split(b',')
367 if len(t) > 2 and t[1] == b'000':
368 s = b','.join([t[0], b'001'] + t[2:])
370 return util.str_to_time(
371 str(s.decode('ascii')), format='%Y,%j,%H:%M:%S.OPTFRAC')
374def ploc(s):
375 if s == b'??':
376 return ''
377 else:
378 return str(s.decode('ascii'))
381def pcode(s):
382 return str(s.decode('ascii'))
385def gett(lst, t):
386 return [x for x in lst if isinstance(x, t)]
389def gett1o(lst, t):
390 lst = [x for x in lst if isinstance(x, t)]
391 if len(lst) == 0:
392 return None
393 elif len(lst) == 1:
394 return lst[0]
395 else:
396 raise RespError('Duplicate entry.')
399def gett1(lst, t):
400 lst = [x for x in lst if isinstance(x, t)]
401 if len(lst) == 0:
402 raise RespError('Entry not found.')
403 elif len(lst) == 1:
404 return lst[0]
405 else:
406 raise RespError('Duplicate entry.')
409class ChannelResponse(guts.Object):
410 '''
411 Response information + channel codes and time span.
412 '''
414 codes = guts.Tuple.T(4, guts.String.T(default=''))
415 start_date = guts.Timestamp.T()
416 end_date = guts.Timestamp.T()
417 response = sxml.Response.T()
420def iload_fh(f):
421 '''
422 Read RESP information from open file handle.
423 '''
425 for sc, cc, rcs in parse3(f):
426 nslc = (
427 pcode(get1(sc, b'16')),
428 pcode(get1(sc, b'03')),
429 ploc(get1(cc, b'03', b'')),
430 pcode(get1(cc, b'04')))
432 try:
433 tmin = pdate(get1(cc, b'22'))
434 tmax = pdate(get1(cc, b'23'))
435 except util.TimeStrError as e:
436 raise RespError('Invalid date in RESP information (%s).' % str(e))
438 stage_elements = {}
440 istage = -1
441 for block, content in rcs:
442 if block not in bdefs:
443 raise RespError('Unknown block type found: %s' % block)
445 istage_temp, x = bdefs[block]['parse'](content)
446 if istage_temp != -1:
447 istage = istage_temp
449 if x is None:
450 continue
452 x.validate()
453 if istage not in stage_elements:
454 stage_elements[istage] = []
456 stage_elements[istage].append(x)
458 istages = sorted(stage_elements.keys())
459 stages = []
460 totalgain = None
461 for istage in istages:
462 elements = stage_elements[istage]
463 if istage == 0:
464 totalgain = gett1(elements, sxml.Gain)
465 else:
466 stage = sxml.ResponseStage(
467 number=istage,
468 poles_zeros_list=gett(elements, sxml.PolesZeros),
469 coefficients_list=gett(elements, sxml.Coefficients),
470 fir=gett1o(elements, sxml.FIR),
471 decimation=gett1o(elements, sxml.Decimation),
472 stage_gain=gett1o(elements, sxml.Gain))
474 stages.append(stage)
476 if stages:
477 resp = sxml.Response(
478 stage_list=stages)
480 if totalgain:
481 totalgain_value = totalgain.value
482 totalgain_frequency = totalgain.frequency
484 else:
485 totalgain_value = 1.
486 gain_frequencies = []
487 for stage in stages:
488 totalgain_value *= stage.stage_gain.value
489 gain_frequencies.append(stage.stage_gain.frequency)
491 totalgain_frequency = gain_frequencies[0]
493 if not all(f == totalgain_frequency for f in gain_frequencies):
494 logger.warning(
495 'No total gain reported and inconsistent gain '
496 'frequency values found in resp file for %s.%s.%s.%s: '
497 'omitting total gain and frequency from created '
498 'instrument sensitivity object.' % nslc)
500 totalgain_value = None
501 totalgain_frequency = None
503 resp.instrument_sensitivity = sxml.Sensitivity(
504 value=totalgain_value,
505 frequency=totalgain_frequency,
506 input_units=stages[0].input_units,
507 output_units=stages[-1].output_units)
509 yield ChannelResponse(
510 codes=nslc,
511 start_date=tmin,
512 end_date=tmax,
513 response=resp)
515 else:
516 logger.warning(
517 'Incomplete response information for %s (%s - %s).',
518 '.'.join(nslc),
519 util.time_to_str(tmin),
520 util.time_to_str(tmax))
522 yield ChannelResponse(
523 codes=nslc,
524 start_date=tmin,
525 end_date=tmax,
526 response=None)
529iload_filename, iload_dirname, iload_glob, iload = util.make_iload_family(
530 iload_fh, 'RESP', ':py:class:`ChannelResponse`')
533def make_stationxml(pyrocko_stations, channel_responses):
534 '''
535 Create stationxml from pyrocko station list and RESP information.
537 :param pyrocko_stations: list of :py:class:`pyrocko.model.Station` objects
538 :param channel_responses: iterable yielding :py:class:`ChannelResponse`
539 objects
540 :returns: :py:class:`pyrocko.fdsn.station.FDSNStationXML` object with
541 merged information
543 If no station information is available for any response information, it
544 is skipped and a warning is emitted.
545 '''
546 pstations = dict((s.nsl(), s) for s in pyrocko_stations)
547 networks = {}
548 stations = {}
549 for (net, sta, loc) in sorted(pstations.keys()):
550 pstation = pstations[net, sta, loc]
551 if net not in networks:
552 networks[net] = sxml.Network(code=net)
554 if (net, sta) not in stations:
555 stations[net, sta] = sxml.Station(
556 code=sta,
557 latitude=sxml.Latitude(pstation.lat),
558 longitude=sxml.Longitude(pstation.lon),
559 elevation=sxml.Distance(pstation.elevation))
561 networks[net].station_list.append(stations[net, sta])
563 for cr in channel_responses:
564 net, sta, loc, cha = cr.codes
565 if (net, sta, loc) in pstations:
566 pstation = pstations[net, sta, loc]
567 pchannel = pstation.get_channel(cha)
568 extra = {}
569 if pchannel is not None:
570 if pchannel.azimuth is not None:
571 extra['azimuth'] = sxml.Azimuth(pchannel.azimuth)
573 if pchannel.dip is not None:
574 extra['dip'] = sxml.Dip(pchannel.dip)
576 channel = sxml.Channel(
577 code=cha,
578 location_code=loc,
579 start_date=cr.start_date,
580 end_date=cr.end_date,
581 latitude=sxml.Latitude(pstation.lat),
582 longitude=sxml.Longitude(pstation.lon),
583 elevation=sxml.Distance(pstation.elevation),
584 depth=sxml.Distance(pstation.depth),
585 response=cr.response,
586 **extra)
588 stations[net, sta].channel_list.append(channel)
589 else:
590 logger.warning('No station information for %s.%s.%s.' %
591 (net, sta, loc))
593 for station in stations.values():
594 station.channel_list.sort(key=lambda c: (c.location_code, c.code))
596 return sxml.FDSNStationXML(
597 source='Converted from Pyrocko stations file and RESP information',
598 created=time.time(),
599 network_list=[networks[net_] for net_ in sorted(networks.keys())])
602if __name__ == '__main__':
603 import sys
604 from pyrocko.model.station import load_stations
606 util.setup_logging(__name__)
608 if len(sys.argv) < 2:
609 sys.exit('usage: python -m pyrocko.fdsn.resp <stations> <resp> ...')
611 stations = load_stations(sys.argv[1])
613 sxml = make_stationxml(stations, iload(sys.argv[2:]))
615 print(sxml.dump_xml())