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

# http://pyrocko.org - GPLv3 

# 

# The Pyrocko Developers, 21st Century 

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

from __future__ import absolute_import, division, print_function 

from builtins import map 

import time 

import re 

import logging 

from pyrocko import util, guts 

from pyrocko.io import io_common 

from pyrocko.io import stationxml as sxml 

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

class RespError(io_common.FileLoadError): 

pass 

def ppolezero(s): 

v = s.split() 

return sxml.PoleZero( 

number=int(v[0]), 

real=sxml.FloatNoUnit( 

value=float(v[1]), 

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

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

imaginary=sxml.FloatNoUnit( 

value=float(v[2]), 

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

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

def pcfu(s): 

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

return sxml.FloatWithUnit( 

value=float(v[-2]), 

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

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

def pnc(s): 

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

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

def punit(s): 

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

def psymmetry(s): 

return { 

b'A': 'NONE', 

b'B': 'ODD', 

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

def ptftype(s): 

if s.startswith(b'A'): 

return 'LAPLACE (RADIANS/SECOND)' 

elif s.startswith(b'B'): 

return 'LAPLACE (HERTZ)' 

elif s.startswith(b'D'): 

return 'DIGITAL (Z-TRANSFORM)' 

else: 

raise RespError('unknown pz transfer function type') 

def pcftype(s): 

if s.startswith(b'A'): 

return 'ANALOG (RADIANS/SECOND)' 

elif s.startswith(b'B'): 

return 'ANALOG (HERTZ)' 

elif s.startswith(b'D'): 

return 'DIGITAL' 

else: 

raise RespError('unknown cf transfer function type') 

def pblock_060(content): 

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

return stage_number, None 

def pblock_053(content): 

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

pzs = sxml.PolesZeros( 

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

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

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

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

normalization_frequency=sxml.Frequency( 

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

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

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

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

x.number = i 

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

x.number = i 

return stage_number, pzs 

def pblock_043(content): 

stage_number = -1 

pzs = sxml.PolesZeros( 

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

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

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

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

normalization_frequency=sxml.Frequency( 

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

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

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

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

x.number = i 

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

x.number = i 

return stage_number, pzs 

def pblock_058(content): 

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

gain = sxml.Gain( 

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

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

return stage_number, gain 

def pblock_048(content): 

stage_number = -1 

gain = sxml.Gain( 

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

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

return stage_number, gain 

def pblock_054(content): 

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

cfs = sxml.Coefficients( 

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

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

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

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

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

return stage_number, cfs 

def pblock_044(content): 

stage_number = -1 

cfs = sxml.Coefficients( 

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

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

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

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

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

return stage_number, cfs 

def pblock_057(content): 

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

deci = sxml.Decimation( 

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

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

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

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

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

return stage_number, deci 

def pblock_047(content): 

stage_number = -1 

deci = sxml.Decimation( 

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

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

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

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

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

return stage_number, deci 

def pblock_061(content): 

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

fir = sxml.FIR( 

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

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

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

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

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

return stage_number, fir 

def pblock_041(content): 

stage_number = -1 

fir = sxml.FIR( 

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

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

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

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

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

return stage_number, fir 

bdefs = { 

b'050': { 

'name': 'Station Identifier Blockette', 

}, 

b'052': { 

'name': 'Channel Identifier Blockette', 

}, 

b'060': { 

'name': 'Response Reference Information', 

'parse': pblock_060, 

}, 

b'053': { 

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

'parse': pblock_053, 

}, 

b'043': { 

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

'parse': pblock_043, 

}, 

b'054': { 

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

'parse': pblock_054, 

}, 

b'044': { 

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

'parse': pblock_044, 

}, 

b'057': { 

'name': 'Decimation Blockette', 

'parse': pblock_057, 

}, 

b'047': { 

'name': 'Decimation Dictionary Blockette', 

'parse': pblock_047, 

}, 

b'058': { 

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

'parse': pblock_058, 

}, 

b'048': { 

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

'parse': pblock_048, 

}, 

b'061': { 

'name': 'FIR Response Blockette', 

'parse': pblock_061, 

}, 

b'041': { 

'name': 'FIR Dictionary Blockette', 

'parse': pblock_041, 

}, 

} 

def parse1(f): 

for line in f: 

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

m = re.match( 

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

if m: 

if m.group(2): 

pass 

elif m.group(3): 

block = m.group(3) 

field = m.group(4) 

key = m.group(7) 

value = m.group(8) 

yield block, field, key, value 

def parse2(f): 

current_b = None 

content = [] 

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

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

if current_b is not None: 

yield current_b, content 

current_b = block 

content = [] 

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

if current_b is not None: 

yield current_b, content 

def parse3(f): 

state = [None, None, []] 

for block, content in parse2(f): 

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

yield state 

state = [None, None, []] 

if block == b'050': 

state[0] = content 

elif block == b'052': 

state[1] = content 

else: 

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

if state[0] and state[1]: 

yield state 

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

for field_, _, value in content: 

if field_ == field: 

return value 

else: 

if optional: 

return None 

elif default is not None: 

return default 

else: 

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

def getn(content, field): 

lst = [] 

for field_, _, value in content: 

if field_ == field: 

lst.append(value) 

return lst 

def pdate(s): 

if len(s) < 17: 

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

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

return None 

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

return None 

else: 

t = s.split(b',') 

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

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

return util.str_to_time( 

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

def ploc(s): 

if s == b'??': 

return '' 

else: 

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

def pcode(s): 

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

def gett(lst, t): 

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

def gett1o(lst, t): 

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

if len(lst) == 0: 

return None 

elif len(lst) == 1: 

return lst[0] 

else: 

raise RespError('duplicate entry') 

def gett1(lst, t): 

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

if len(lst) == 0: 

raise RespError('entry not found') 

elif len(lst) == 1: 

return lst[0] 

else: 

raise RespError('duplicate entry') 

class ChannelResponse(guts.Object): 

'''Response information + channel codes and time span.''' 

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

start_date = guts.Timestamp.T() 

end_date = guts.Timestamp.T() 

response = sxml.Response.T() 

def iload_fh(f): 

'''Read RESP information from open file handle.''' 

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

nslc = ( 

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

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

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

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

try: 

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

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

except util.TimeStrError as e: 

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

stage_elements = {} 

istage = -1 

for block, content in rcs: 

if block not in bdefs: 

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

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

if istage_temp != -1: 

istage = istage_temp 

if x is None: 

continue 

x.validate() 

if istage not in stage_elements: 

stage_elements[istage] = [] 

stage_elements[istage].append(x) 

istages = sorted(stage_elements.keys()) 

stages = [] 

totalgain = None 

for istage in istages: 

elements = stage_elements[istage] 

if istage == 0: 

totalgain = gett1(elements, sxml.Gain) 

else: 

stage = sxml.ResponseStage( 

number=istage, 

poles_zeros_list=gett(elements, sxml.PolesZeros), 

coefficients_list=gett(elements, sxml.Coefficients), 

fir=gett1o(elements, sxml.FIR), 

decimation=gett1o(elements, sxml.Decimation), 

stage_gain=gett1o(elements, sxml.Gain)) 

stages.append(stage) 

if stages: 

resp = sxml.Response( 

stage_list=stages) 

if totalgain: 

totalgain_value = totalgain.value 

totalgain_frequency = totalgain.frequency 

else: 

totalgain_value = 1. 

gain_frequencies = [] 

for stage in stages: 

totalgain_value *= stage.stage_gain.value 

gain_frequencies.append(stage.stage_gain.frequency) 

totalgain_frequency = gain_frequencies[0] 

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

logger.warn( 

'no total gain reported and inconsistent gain ' 

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

'omitting total gain and frequency from created ' 

'instrument sensitivity object' % nslc) 

totalgain_value = None 

totalgain_frequency = None 

resp.instrument_sensitivity = sxml.Sensitivity( 

value=totalgain_value, 

frequency=totalgain_frequency, 

input_units=stages[0].input_units, 

output_units=stages[-1].output_units) 

yield ChannelResponse( 

codes=nslc, 

start_date=tmin, 

end_date=tmax, 

response=resp) 

else: 

raise RespError('incomplete response information') 

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

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

def make_stationxml(pyrocko_stations, channel_responses): 

'''Create stationxml from pyrocko station list and RESP information. 

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

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

objects 

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

merged information 

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

is skipped and a warning is emitted. 

''' 

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

networks = {} 

stations = {} 

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

pstation = pstations[net, sta, loc] 

if net not in networks: 

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

if (net, sta) not in stations: 

stations[net, sta] = sxml.Station( 

code=sta, 

latitude=sxml.Latitude(pstation.lat), 

longitude=sxml.Longitude(pstation.lon), 

elevation=sxml.Distance(pstation.elevation)) 

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

for cr in channel_responses: 

net, sta, loc, cha = cr.codes 

if (net, sta, loc) in pstations: 

pstation = pstations[net, sta, loc] 

channel = sxml.Channel( 

code=cha, 

location_code=loc, 

start_date=cr.start_date, 

end_date=cr.end_date, 

latitude=sxml.Latitude(pstation.lat), 

longitude=sxml.Longitude(pstation.lon), 

elevation=sxml.Distance(pstation.elevation), 

depth=sxml.Distance(pstation.depth), 

response=cr.response) 

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

else: 

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

(net, sta, loc)) 

for station in stations.values(): 

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

return sxml.FDSNStationXML( 

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

created=time.time(), 

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

if __name__ == '__main__': 

import sys 

from pyrocko.model.station import load_stations 

util.setup_logging(__name__) 

if len(sys.argv) < 2: 

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

stations = load_stations(sys.argv[1]) 

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

print(sxml.dump_xml())