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