1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import calendar
7import logging
8import numpy as num
9from scipy import signal
11from pyrocko import util, response
13unpack_fixed = util.unpack_fixed
15logger = logging.getLogger('pyrocko.io.seisan_response')
17d2r = num.pi/180.
20class SeisanResponseFileError(Exception):
21 pass
24class SeisanResponseFile(object):
26 def __init__(self):
27 pass
29 def read(self, filename):
31 f = open(filename, 'rb')
32 line = f.readline()
33 line = str(line.decode('ascii'))
35 station, component, century, deltayear, doy, month, day, hr, mi, sec \
36 = unpack_fixed(
37 'a5,a4,@1,i2,x1,i3,x1,i2,x1,i2,x1,i2,x1,i2,x1,f6',
38 line[0:35],
39 lambda s: {' ': 1900, '0': 1900, '1': 2000}[s])
41 # is_accelerometer = line[6] == 'A'
43 latitude, longitude, elevation, filetype, cf_flag = \
44 unpack_fixed(
45 'f8?,x1,f9?,x1,f5?,x2,@1,a1',
46 line[50:80],
47 lambda s: {
48 ' ': 'gains-and-filters',
49 't': 'tabulated',
50 'p': 'poles-and-zeros'}[s.lower()])
52 line = f.readline()
53 line = str(line.decode('ascii'))
55 comment = line.strip()
56 tmin = util.to_time_float(calendar.timegm(
57 (century+deltayear, 1, doy, hr, mi, int(sec)))) + sec-int(sec)
59 if filetype == 'gains-and-filters':
61 line = f.readline()
62 line = str(line.decode('ascii'))
64 period, damping, sensor_sensitivity, amplifier_gain, \
65 digitizer_gain, gain_1hz, filter1_corner, filter1_order, \
66 filter2_corner, filter2_order = unpack_fixed(
67 'f8,f8,f8,f8,f8,f8,f8,f8,f8,f8',
68 line)
70 filter_defs = [
71 filter1_corner,
72 filter1_order,
73 filter2_corner,
74 filter2_order]
76 line = f.readline()
77 line = str(line.decode('ascii'))
79 filter_defs.extend(
80 unpack_fixed('f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', line))
82 filters = []
83 for order, corner in zip(filter_defs[1::2], filter_defs[0::2]):
84 if order != 0.0:
85 filters.append((order, corner))
87 if filetype in ('gains-and-filters', 'tabulated'):
88 data = ([], [], [])
89 for iy in range(3):
90 for ix in range(3):
91 line = f.readline()
92 line = str(line.decode('ascii'))
94 data[ix].extend(unpack_fixed(
95 'f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', line))
97 response_table = num.array(data, dtype=float)
99 if filetype == 'poles-and-zeros':
100 assert False, 'poles-and-zeros file type not implemented yet ' \
101 'for seisan response file format'
103 f.close()
105 if num.all(num.abs(response_table[2]) <= num.pi):
106 logger.warning(
107 'assuming tabulated phases are given in radians instead of '
108 'degrees')
110 cresp = response_table[1] * (
111 num.cos(response_table[2])
112 + 1.0j*num.sin(response_table[2]))
113 else:
114 cresp = response_table[1] * (
115 num.cos(response_table[2]*d2r)
116 + 1.0j*num.sin(response_table[2]*d2r))
118 self.station = station
119 self.component = component
120 self.tmin = tmin
121 self.latitude = latitude
122 self.longitude = longitude
123 self.elevation = elevation
124 self.filetype = filetype
125 self.comment = comment
126 self.period = period
127 self.damping = damping
128 self.sensor_sensitivity = sensor_sensitivity
129 self.amplifier_gain = amplifier_gain
130 self.digitizer_gain = digitizer_gain
131 self.gain_1hz = gain_1hz
132 self.filters = filters
134 self.sampled_response = response.SampledResponse(
135 response_table[0], cresp)
137 self._check_tabulated_response(filename=filename)
139 def response(self, freqs, method='from-filetype', type='displacement'):
141 freqs = num.asarray(freqs)
142 want_scalar = False
143 if freqs.ndim == 0:
144 freqs = num.array([freqs])
145 want_scalar = True
147 if method == 'from-filetype':
148 method = self.filetype
150 if method == 'gains-and-filters':
151 dresp = self._response_prototype_and_filters(freqs) \
152 / abs(self._response_prototype_and_filters([1.0])[0]) \
153 * self.gain_1hz
155 elif method == 'tabulated':
156 dresp = self._response_tabulated(freqs) \
157 / abs(self._response_tabulated([1.0])[0]) * self.gain_1hz
159 elif method == 'poles-and-zeros':
160 raise Exception(
161 'fix me! poles-and-zeros response in seisan is broken: where '
162 'should the "normalization" come from?')
164 # dresp = self._response_from_poles_and_zeros(freqs) *normalization
166 elif method == 'gains':
167 dresp = num.ones(freqs.size) * self.gain_total() * 1.0j * 2. \
168 * num.pi * freqs
170 else:
171 assert False, 'invalid response method specified'
173 if type == 'velocity':
174 dresp /= 1.0j * 2. * num.pi * freqs
176 if want_scalar:
177 return dresp[0]
178 else:
179 return dresp
181 def gain_total(self):
182 return self.sensor_sensitivity * 10.**(self.amplifier_gain/20.) \
183 * self.digitizer_gain
185 def _prototype_response_velocity(self, s):
186 omega0 = 2. * num.pi / self.period
187 return s**2/(omega0**2 + s**2 + 2.0*s*omega0*self.damping)
189 def _response_prototype_and_filters(self, freqs):
190 freqs = num.asarray(freqs, dtype=float)
191 iomega = 1.0j * 2. * num.pi * freqs
193 trans = iomega * self._prototype_response_velocity(iomega)
195 for (order, corner) in self.filters:
196 if order < 0. or corner < 0.:
197 b, a = signal.butter(
198 abs(order), [abs(corner)], btype='high', analog=1)
199 else:
200 b, a = signal.butter(
201 order, [corner], btype='low', analog=1)
203 trans *= signal.freqs(b, a, freqs)[1]
205 return trans
207 def _response_tabulated(self, freqs):
208 freqs = num.asarray(freqs, dtype=float)
209 return self.sampled_response.evaluate(freqs)
211 def _response_from_poles_and_zeros(self, freqs):
212 assert False, 'poles-and-zeros file type not implemented yet for ' \
213 'seisan response file format'
214 return None
216 def _check_tabulated_response(self, filename='?'):
217 if self.filetype == 'gains-and-filters':
218 freqs = self.sampled_response.frequencies()
220 trans_gaf = self.response(freqs, method='gains-and-filters')
221 atrans_gaf = num.abs(trans_gaf)
223 trans_tab = self.response(freqs, method='tabulated')
224 atrans_tab = num.abs(trans_tab)
226 if num.any(num.abs(atrans_gaf-atrans_tab)
227 / num.abs(atrans_gaf+atrans_tab) > 1./100.):
229 logger.warning(
230 'inconsistent amplitudes in tabulated response '
231 '(in file "%s") (max |a-b|/|a+b| is %g' % (
232 filename,
233 num.amax(num.abs(atrans_gaf-atrans_tab)
234 / abs(atrans_gaf+atrans_tab))))
236 else:
237 if num.any(num.abs(trans_gaf-trans_tab)
238 > abs(trans_gaf+trans_tab)/100.):
240 logger.warning(
241 'inconsistent phase values in tabulated response '
242 '(in file "%s"' % filename)
244 def __str__(self):
246 return '''--- Seisan Response File ---
247station: %s
248component: %s
249start time: %s
250latitude: %f
251longitude: %f
252elevation: %f
253filetype: %s
254comment: %s
255sensor period: %g
256sensor damping: %g
257sensor sensitivity: %g
258amplifier gain: %g
259digitizer gain: %g
260gain at 1 Hz: %g
261filters: %s
262''' % (self.station, self.component, util.time_to_str(self.tmin),
263 self.latitude, self.longitude, self.elevation, self.filetype,
264 self.comment, self.period, self.damping, self.sensor_sensitivity,
265 self.amplifier_gain, self.digitizer_gain, self.gain_1hz,
266 self.filters)
268 def plot_amplitudes(self, filename_pdf, type='displacement'):
270 from pyrocko.plot import gmtpy
272 p = gmtpy.LogLogPlot()
273 f = self.sampled_response.frequencies()
274 atab = num.abs(self.response(f, method='tabulated', type=type))
275 acom = num.abs(self.response(f, method='gains-and-filters', type=type))
276 aconst = num.abs(self.response(f, method='gains', type=type))
278 p.plot((f, atab), '-W2p,red')
279 p.plot((f, acom), '-W1p,blue')
280 p.plot((f, aconst), '-W1p,green')
281 p.save(filename_pdf)
283 def plot_phases(self, filename_pdf, type='displacement'):
285 from pyrocko.plot import gmtpy
287 p = gmtpy.LogLinPlot()
288 f = self.sampled_response.frequencies()
289 atab = num.unwrap(num.angle(self.response(
290 f, method='tabulated', type=type))) / d2r
291 acom = num.unwrap(num.angle(self.response(
292 f, method='gains-and-filters', type=type))) / d2r
293 aconst = num.unwrap(num.angle(self.response(
294 f, method='gains', type=type))) / d2r
296 p.plot((f, atab), '-W2p,red')
297 p.plot((f, acom), '-W1p,blue')
298 p.plot((f, aconst), '-W1p,green')
299 p.save(filename_pdf)