1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import time
7from collections import defaultdict
9from pyrocko import util, pz, response
10from pyrocko.io import io_common
11from pyrocko.io import stationxml as sxml
12from pyrocko.guts import Object, Tuple, String, Timestamp, Float
15guts_prefix = 'pf'
18class EnhancedSacPzError(io_common.FileLoadError):
19 pass
22class EnhancedSacPzResponse(Object):
23 codes = Tuple.T(4, String.T())
24 tmin = Timestamp.T(optional=True)
25 tmax = Timestamp.T(optional=True)
26 lat = Float.T()
27 lon = Float.T()
28 elevation = Float.T()
29 depth = Float.T()
30 dip = Float.T()
31 azimuth = Float.T()
32 input_unit = String.T()
33 output_unit = String.T()
34 response = response.PoleZeroResponse.T()
36 def spans(self, *args):
37 if len(args) == 0:
38 return True
39 elif len(args) == 1:
40 return ((self.tmin is None or
41 self.tmin <= args[0]) and
42 (self.tmax is None or
43 args[0] <= self.tmax))
45 elif len(args) == 2:
46 return ((self.tmin is None or
47 args[1] >= self.tmin) and
48 (self.tmax is None or
49 self.tmax >= args[0]))
52def make_stationxml_response(presponse, input_unit, output_unit):
53 return sxml.Response.from_pyrocko_pz_response(
54 presponse, input_unit, output_unit, 1.0)
57this_year = time.gmtime()[0]
60def dummy_aware_str_to_time(s, time_format='%Y-%m-%dT%H:%M:%S'):
61 try:
62 util.str_to_time(s, format=time_format)
63 except util.TimeStrError:
64 year = int(s[:4])
65 if year > this_year + 100:
66 return None # StationXML contained a dummy end date
68 raise
71def iload_fh(f, time_format='%Y-%m-%dT%H:%M:%S'):
72 zeros, poles, constant, comments = pz.read_sac_zpk(file=f,
73 get_comments=True)
74 d = {}
75 for line in comments:
76 toks = line.split(':', 1)
77 if len(toks) == 2:
78 temp = toks[0].strip('* \t')
79 for k in ('network', 'station', 'location', 'channel', 'start',
80 'end', 'latitude', 'longitude', 'depth', 'elevation',
81 'dip', 'azimuth', 'input unit', 'output unit'):
83 if temp.lower().startswith(k):
84 d[k] = toks[1].strip()
86 resp = response.PoleZeroResponse(zeros, poles, constant)
88 try:
89 yield EnhancedSacPzResponse(
90 codes=(d['network'], d['station'], d['location'], d['channel']),
91 tmin=util.str_to_time(d['start'], format=time_format),
92 tmax=dummy_aware_str_to_time(d['end']),
93 lat=float(d['latitude']),
94 lon=float(d['longitude']),
95 elevation=float(d['elevation']),
96 depth=float(d['depth']),
97 dip=float(d['dip']),
98 azimuth=float(d['azimuth']),
99 input_unit=d['input unit'],
100 output_unit=d['output unit'],
101 response=resp)
102 except KeyError as e:
103 raise EnhancedSacPzError(
104 'cannot get all required information "%s"' % e.args[0])
107iload_filename, iload_dirname, iload_glob, iload = util.make_iload_family(
108 iload_fh, 'SACPZ', ':py:class:`EnhancedSacPzResponse`')
111def make_stationxml(responses, inconsistencies='warn'):
112 '''
113 Create stationxml from "enhanced" SACPZ information.
115 :param responses: iterable yielding
116 :py:class:`EnhancedSacPzResponse` objects
117 :returns: :py:class:`pyrocko.fdsn.station.FDSNStationXML` object
118 '''
120 networks = {}
121 stations = {}
123 station_coords = defaultdict(list)
124 station_channels = defaultdict(list)
125 for cr in responses:
126 net, sta, loc, cha = cr.codes
127 station_coords[net, sta].append(
128 (cr.codes, cr.lat, cr.lon, cr.elevation))
130 station_channels[net, sta].append(sxml.Channel(
131 code=cha,
132 location_code=loc,
133 start_date=cr.tmin,
134 end_date=cr.tmax,
135 latitude=sxml.Latitude(cr.lat),
136 longitude=sxml.Longitude(cr.lon),
137 elevation=sxml.Distance(cr.elevation),
138 depth=sxml.Distance(cr.depth),
139 response=make_stationxml_response(
140 cr.response, cr.input_unit, cr.output_unit)))
142 for (net, sta), v in sorted(station_coords.items()):
143 lat, lon, elevation = util.consistency_merge(
144 v,
145 message='channel lat/lon/elevation values differ',
146 error=inconsistencies)
148 if net not in networks:
149 networks[net] = sxml.Network(code=net)
151 stations[net, sta] = sxml.Station(
152 code=sta,
153 latitude=sxml.Latitude(lat),
154 longitude=sxml.Longitude(lon),
155 elevation=sxml.Distance(elevation),
156 channel_list=sorted(
157 station_channels[net, sta],
158 key=lambda c: (c.location_code, c.code)))
160 networks[net].station_list.append(stations[net, sta])
162 return sxml.FDSNStationXML(
163 source='Converted from "enhanced" SACPZ information',
164 created=time.time(),
165 network_list=[networks[net_] for net_ in sorted(networks.keys())])
168if __name__ == '__main__':
170 import sys
172 util.setup_logging(__name__)
174 if len(sys.argv) < 2:
175 sys.exit('usage: python -m pyrocko.station.enhanced_sacpz <sacpz> ...')
177 sxml_in = make_stationxml(iload(sys.argv[1:]))
178 print(sxml_in.dump_xml())