1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import
7import calendar
8import logging
9import numpy as num
10from scipy import signal
12from pyrocko import util, trace
14unpack_fixed = util.unpack_fixed
16logger = logging.getLogger('pyrocko.io.seisan_response')
18d2r = num.pi/180.
21class SeisanResponseFileError(Exception):
22 pass
25class SeisanResponseFile(object):
27 def __init__(self):
28 pass
30 def read(self, filename):
32 f = open(filename, 'rb')
33 line = f.readline()
34 line = str(line.decode('ascii'))
36 station, component, century, deltayear, doy, month, day, hr, mi, sec \
37 = unpack_fixed(
38 'a5,a4,@1,i2,x1,i3,x1,i2,x1,i2,x1,i2,x1,i2,x1,f6',
39 line[0:35],
40 lambda s: {' ': 1900, '0': 1900, '1': 2000}[s])
42 # is_accelerometer = line[6] == 'A'
44 latitude, longitude, elevation, filetype, cf_flag = \
45 unpack_fixed(
46 'f8?,x1,f9?,x1,f5?,x2,@1,a1',
47 line[50:80],
48 lambda s: {
49 ' ': 'gains-and-filters',
50 't': 'tabulated',
51 'p': 'poles-and-zeros'}[s.lower()])
53 line = f.readline()
54 line = str(line.decode('ascii'))
56 comment = line.strip()
57 tmin = util.to_time_float(calendar.timegm(
58 (century+deltayear, 1, doy, hr, mi, int(sec)))) + sec-int(sec)
60 if filetype == 'gains-and-filters':
62 line = f.readline()
63 line = str(line.decode('ascii'))
65 period, damping, sensor_sensitivity, amplifier_gain, \
66 digitizer_gain, gain_1hz, filter1_corner, filter1_order, \
67 filter2_corner, filter2_order = unpack_fixed(
68 'f8,f8,f8,f8,f8,f8,f8,f8,f8,f8',
69 line)
71 filter_defs = [
72 filter1_corner,
73 filter1_order,
74 filter2_corner,
75 filter2_order]
77 line = f.readline()
78 line = str(line.decode('ascii'))
80 filter_defs.extend(
81 unpack_fixed('f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', line))
83 filters = []
84 for order, corner in zip(filter_defs[1::2], filter_defs[0::2]):
85 if order != 0.0:
86 filters.append((order, corner))
88 if filetype in ('gains-and-filters', 'tabulated'):
89 data = ([], [], [])
90 for iy in range(3):
91 for ix in range(3):
92 line = f.readline()
93 line = str(line.decode('ascii'))
95 data[ix].extend(unpack_fixed(
96 'f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', line))
98 response_table = num.array(data, dtype=float)
100 if filetype == 'poles-and-zeros':
101 assert False, 'poles-and-zeros file type not implemented yet ' \
102 'for seisan response file format'
104 f.close()
106 if num.all(num.abs(response_table[2]) <= num.pi):
107 logger.warning(
108 'assuming tabulated phases are given in radians instead of '
109 'degrees')
111 cresp = response_table[1] * (
112 num.cos(response_table[2])
113 + 1.0j*num.sin(response_table[2]))
114 else:
115 cresp = response_table[1] * (
116 num.cos(response_table[2]*d2r)
117 + 1.0j*num.sin(response_table[2]*d2r))
119 self.station = station
120 self.component = component
121 self.tmin = tmin
122 self.latitude = latitude
123 self.longitude = longitude
124 self.elevation = elevation
125 self.filetype = filetype
126 self.comment = comment
127 self.period = period
128 self.damping = damping
129 self.sensor_sensitivity = sensor_sensitivity
130 self.amplifier_gain = amplifier_gain
131 self.digitizer_gain = digitizer_gain
132 self.gain_1hz = gain_1hz
133 self.filters = filters
135 self.sampled_response = trace.SampledResponse(
136 response_table[0], cresp)
138 self._check_tabulated_response(filename=filename)
140 def response(self, freqs, method='from-filetype', type='displacement'):
142 freqs = num.asarray(freqs)
143 want_scalar = False
144 if freqs.ndim == 0:
145 freqs = num.array([freqs])
146 want_scalar = True
148 if method == 'from-filetype':
149 method = self.filetype
151 if method == 'gains-and-filters':
152 dresp = self._response_prototype_and_filters(freqs) \
153 / abs(self._response_prototype_and_filters([1.0])[0]) \
154 * self.gain_1hz
156 elif method == 'tabulated':
157 dresp = self._response_tabulated(freqs) \
158 / abs(self._response_tabulated([1.0])[0]) * self.gain_1hz
160 elif method == 'poles-and-zeros':
161 raise Exception(
162 'fix me! poles-and-zeros response in seisan is broken: where '
163 'should the "normalization" come from?')
165 # dresp = self._response_from_poles_and_zeros(freqs) *normalization
167 elif method == 'gains':
168 dresp = num.ones(freqs.size) * self.gain_total() * 1.0j * 2. \
169 * num.pi * freqs
171 else:
172 assert False, 'invalid response method specified'
174 if type == 'velocity':
175 dresp /= 1.0j * 2. * num.pi * freqs
177 if want_scalar:
178 return dresp[0]
179 else:
180 return dresp
182 def gain_total(self):
183 return self.sensor_sensitivity * 10.**(self.amplifier_gain/20.) \
184 * self.digitizer_gain
186 def _prototype_response_velocity(self, s):
187 omega0 = 2. * num.pi / self.period
188 return s**2/(omega0**2 + s**2 + 2.0*s*omega0*self.damping)
190 def _response_prototype_and_filters(self, freqs):
191 freqs = num.asarray(freqs, dtype=float)
192 iomega = 1.0j * 2. * num.pi * freqs
194 trans = iomega * self._prototype_response_velocity(iomega)
196 for (order, corner) in self.filters:
197 if order < 0. or corner < 0.:
198 b, a = signal.butter(
199 abs(order), [abs(corner)], btype='high', analog=1)
200 else:
201 b, a = signal.butter(
202 order, [corner], btype='low', analog=1)
204 trans *= signal.freqs(b, a, freqs)[1]
206 return trans
208 def _response_tabulated(self, freqs):
209 freqs = num.asarray(freqs, dtype=float)
210 return self.sampled_response.evaluate(freqs)
212 def _response_from_poles_and_zeros(self, freqs):
213 assert False, 'poles-and-zeros file type not implemented yet for ' \
214 'seisan response file format'
215 return None
217 def _check_tabulated_response(self, filename='?'):
218 if self.filetype == 'gains-and-filters':
219 freqs = self.sampled_response.frequencies()
221 trans_gaf = self.response(freqs, method='gains-and-filters')
222 atrans_gaf = num.abs(trans_gaf)
224 trans_tab = self.response(freqs, method='tabulated')
225 atrans_tab = num.abs(trans_tab)
227 if num.any(num.abs(atrans_gaf-atrans_tab)
228 / num.abs(atrans_gaf+atrans_tab) > 1./100.):
230 logger.warning(
231 'inconsistent amplitudes in tabulated response '
232 '(in file "%s") (max |a-b|/|a+b| is %g' % (
233 filename,
234 num.amax(num.abs(atrans_gaf-atrans_tab)
235 / abs(atrans_gaf+atrans_tab))))
237 else:
238 if num.any(num.abs(trans_gaf-trans_tab)
239 > abs(trans_gaf+trans_tab)/100.):
241 logger.warning(
242 'inconsistent phase values in tabulated response '
243 '(in file "%s"' % filename)
245 def __str__(self):
247 return '''--- Seisan Response File ---
248station: %s
249component: %s
250start time: %s
251latitude: %f
252longitude: %f
253elevation: %f
254filetype: %s
255comment: %s
256sensor period: %g
257sensor damping: %g
258sensor sensitivity: %g
259amplifier gain: %g
260digitizer gain: %g
261gain at 1 Hz: %g
262filters: %s
263''' % (self.station, self.component, util.time_to_str(self.tmin),
264 self.latitude, self.longitude, self.elevation, self.filetype,
265 self.comment, self.period, self.damping, self.sensor_sensitivity,
266 self.amplifier_gain, self.digitizer_gain, self.gain_1hz,
267 self.filters)
269 def plot_amplitudes(self, filename_pdf, type='displacement'):
271 from pyrocko.plot import gmtpy
273 p = gmtpy.LogLogPlot()
274 f = self.sampled_response.frequencies()
275 atab = num.abs(self.response(f, method='tabulated', type=type))
276 acom = num.abs(self.response(f, method='gains-and-filters', type=type))
277 aconst = num.abs(self.response(f, method='gains', type=type))
279 p.plot((f, atab), '-W2p,red')
280 p.plot((f, acom), '-W1p,blue')
281 p.plot((f, aconst), '-W1p,green')
282 p.save(filename_pdf)
284 def plot_phases(self, filename_pdf, type='displacement'):
286 from pyrocko.plot import gmtpy
288 p = gmtpy.LogLinPlot()
289 f = self.sampled_response.frequencies()
290 atab = num.unwrap(num.angle(self.response(
291 f, method='tabulated', type=type))) / d2r
292 acom = num.unwrap(num.angle(self.response(
293 f, method='gains-and-filters', type=type))) / d2r
294 aconst = num.unwrap(num.angle(self.response(
295 f, method='gains', type=type))) / d2r
297 p.plot((f, atab), '-W2p,red')
298 p.plot((f, acom), '-W1p,blue')
299 p.plot((f, aconst), '-W1p,green')
300 p.save(filename_pdf)