1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import time
8import re
9import logging
11from pyrocko import util, guts
12from pyrocko.io import io_common
13from pyrocko.io import stationxml as sxml
15logger = logging.getLogger('pyrocko.io.resp')
18class RespError(io_common.FileLoadError):
19 pass
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))
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)
44def pnc(s):
45 v = list(map(float, s.split()))
46 return sxml.NumeratorCoefficient(i=int(v[0]), value=float(v[1]))
49def punit(s):
50 return str(s.split()[0].decode('ascii'))
53def psymmetry(s):
54 return {
55 b'A': 'NONE',
56 b'B': 'ODD',
57 b'C': 'EVEN'}[s.upper()]
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')
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')
82def pblock_060(content):
83 stage_number = int(get1(content, b'04'))
84 return stage_number, None
87def pblock_053(content):
88 stage_number = int(get1(content, b'04'))
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'))),
98 zero_list=list(map(ppolezero, getn(content, b'10-13'))),
99 pole_list=list(map(ppolezero, getn(content, b'15-18'))))
101 for i, x in enumerate(pzs.zero_list):
102 x.number = i
104 for i, x in enumerate(pzs.pole_list):
105 x.number = i
107 return stage_number, pzs
110def pblock_043(content):
111 stage_number = -1
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'))),
121 zero_list=list(map(ppolezero, getn(content, b'11-14'))),
122 pole_list=list(map(ppolezero, getn(content, b'16-19'))))
124 for i, x in enumerate(pzs.zero_list):
125 x.number = i
127 for i, x in enumerate(pzs.pole_list):
128 x.number = i
130 return stage_number, pzs
133def pblock_058(content):
134 stage_number = int(get1(content, b'03'))
136 gain = sxml.Gain(
137 value=float(get1(content, b'04')),
138 frequency=float(get1(content, b'05').split()[0]))
140 return stage_number, gain
143def pblock_048(content):
144 stage_number = -1
146 gain = sxml.Gain(
147 value=float(get1(content, b'05')),
148 frequency=float(get1(content, b'06').split()[0]))
150 return stage_number, gain
153def pblock_054(content):
154 stage_number = int(get1(content, b'04'))
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'))))
163 return stage_number, cfs
166def pblock_044(content):
167 stage_number = -1
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'))))
176 return stage_number, cfs
179def pblock_057(content):
180 stage_number = int(get1(content, b'03'))
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'))))
189 return stage_number, deci
192def pblock_047(content):
193 stage_number = -1
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'))))
202 return stage_number, deci
205def pblock_061(content):
206 stage_number = int(get1(content, b'03'))
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'))))
215 return stage_number, fir
218def pblock_041(content):
219 stage_number = -1
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'))))
228 return stage_number, fir
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}
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
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
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
310 current_b = block
311 content = []
313 content.append((field, key, value))
315 if current_b is not None:
316 yield current_b, content
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, []]
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))
333 if state[0] and state[1]:
334 yield state
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)
350def getn(content, field):
351 lst = []
352 for field_, _, value in content:
353 if field_ == field:
354 lst.append(value)
355 return lst
358def pdate(s):
359 if len(s) < 17:
360 s += b'0000,001,00:00:00'[len(s):]
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:])
371 return util.str_to_time(
372 str(s.decode('ascii')), format='%Y,%j,%H:%M:%S.OPTFRAC')
375def ploc(s):
376 if s == b'??':
377 return ''
378 else:
379 return str(s.decode('ascii'))
382def pcode(s):
383 return str(s.decode('ascii'))
386def gett(lst, t):
387 return [x for x in lst if isinstance(x, t)]
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')
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')
410class ChannelResponse(guts.Object):
411 '''
412 Response information + channel codes and time span.
413 '''
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()
421def iload_fh(f):
422 '''
423 Read RESP information from open file handle.
424 '''
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')))
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))
439 stage_elements = {}
441 istage = -1
442 for block, content in rcs:
443 if block not in bdefs:
444 raise RespError('unknown block type found: %s' % block)
446 istage_temp, x = bdefs[block]['parse'](content)
447 if istage_temp != -1:
448 istage = istage_temp
450 if x is None:
451 continue
453 x.validate()
454 if istage not in stage_elements:
455 stage_elements[istage] = []
457 stage_elements[istage].append(x)
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))
475 stages.append(stage)
477 if stages:
478 resp = sxml.Response(
479 stage_list=stages)
481 if totalgain:
482 totalgain_value = totalgain.value
483 totalgain_frequency = totalgain.frequency
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)
492 totalgain_frequency = gain_frequencies[0]
494 if not all(f == totalgain_frequency for f in gain_frequencies):
495 logger.warn(
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)
501 totalgain_value = None
502 totalgain_frequency = None
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)
510 yield ChannelResponse(
511 codes=nslc,
512 start_date=tmin,
513 end_date=tmax,
514 response=resp)
516 else:
517 raise RespError('incomplete response information')
520iload_filename, iload_dirname, iload_glob, iload = util.make_iload_family(
521 iload_fh, 'RESP', ':py:class:`ChannelResponse`')
524def make_stationxml(pyrocko_stations, channel_responses):
525 '''
526 Create stationxml from pyrocko station list and RESP information.
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
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)
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))
552 networks[net].station_list.append(stations[net, sta])
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)
564 if pchannel.dip is not None:
565 extra['dip'] = sxml.Dip(pchannel.dip)
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)
579 stations[net, sta].channel_list.append(channel)
580 else:
581 logger.warning('no station information for %s.%s.%s' %
582 (net, sta, loc))
584 for station in stations.values():
585 station.channel_list.sort(key=lambda c: (c.location_code, c.code))
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())])
593if __name__ == '__main__':
594 import sys
595 from pyrocko.model.station import load_stations
597 util.setup_logging(__name__)
599 if len(sys.argv) < 2:
600 sys.exit('usage: python -m pyrocko.fdsn.resp <stations> <resp> ...')
602 stations = load_stations(sys.argv[1])
604 sxml = make_stationxml(stations, iload(sys.argv[2:]))
606 print(sxml.dump_xml())