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