Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/seisan_response.py: 16%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Read instrument responses from `SEISAN <http://seisan.info/>`_ files. 

8''' 

9 

10import calendar 

11import logging 

12import numpy as num 

13from scipy import signal 

14 

15from pyrocko import util, response 

16 

17unpack_fixed = util.unpack_fixed 

18 

19logger = logging.getLogger('pyrocko.io.seisan_response') 

20 

21d2r = num.pi/180. 

22 

23 

24class SeisanResponseFileError(Exception): 

25 pass 

26 

27 

28class SeisanResponseFile(object): 

29 

30 def __init__(self): 

31 pass 

32 

33 def read(self, filename): 

34 

35 f = open(filename, 'rb') 

36 line = f.readline() 

37 line = str(line.decode('ascii')) 

38 

39 station, component, century, deltayear, doy, month, day, hr, mi, sec \ 

40 = unpack_fixed( 

41 'a5,a4,@1,i2,x1,i3,x1,i2,x1,i2,x1,i2,x1,i2,x1,f6', 

42 line[0:35], 

43 lambda s: {' ': 1900, '0': 1900, '1': 2000}[s]) 

44 

45 # is_accelerometer = line[6] == 'A' 

46 

47 latitude, longitude, elevation, filetype, cf_flag = \ 

48 unpack_fixed( 

49 'f8?,x1,f9?,x1,f5?,x2,@1,a1', 

50 line[50:80], 

51 lambda s: { 

52 ' ': 'gains-and-filters', 

53 't': 'tabulated', 

54 'p': 'poles-and-zeros'}[s.lower()]) 

55 

56 line = f.readline() 

57 line = str(line.decode('ascii')) 

58 

59 comment = line.strip() 

60 tmin = util.to_time_float(calendar.timegm( 

61 (century+deltayear, 1, doy, hr, mi, int(sec)))) + sec-int(sec) 

62 

63 if filetype == 'gains-and-filters': 

64 

65 line = f.readline() 

66 line = str(line.decode('ascii')) 

67 

68 period, damping, sensor_sensitivity, amplifier_gain, \ 

69 digitizer_gain, gain_1hz, filter1_corner, filter1_order, \ 

70 filter2_corner, filter2_order = unpack_fixed( 

71 'f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', 

72 line) 

73 

74 filter_defs = [ 

75 filter1_corner, 

76 filter1_order, 

77 filter2_corner, 

78 filter2_order] 

79 

80 line = f.readline() 

81 line = str(line.decode('ascii')) 

82 

83 filter_defs.extend( 

84 unpack_fixed('f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', line)) 

85 

86 filters = [] 

87 for order, corner in zip(filter_defs[1::2], filter_defs[0::2]): 

88 if order != 0.0: 

89 filters.append((order, corner)) 

90 

91 if filetype in ('gains-and-filters', 'tabulated'): 

92 data = ([], [], []) 

93 for iy in range(3): 

94 for ix in range(3): 

95 line = f.readline() 

96 line = str(line.decode('ascii')) 

97 

98 data[ix].extend(unpack_fixed( 

99 'f8,f8,f8,f8,f8,f8,f8,f8,f8,f8', line)) 

100 

101 response_table = num.array(data, dtype=float) 

102 

103 if filetype == 'poles-and-zeros': 

104 assert False, 'poles-and-zeros file type not implemented yet ' \ 

105 'for seisan response file format' 

106 

107 f.close() 

108 

109 if num.all(num.abs(response_table[2]) <= num.pi): 

110 logger.warning( 

111 'assuming tabulated phases are given in radians instead of ' 

112 'degrees') 

113 

114 cresp = response_table[1] * ( 

115 num.cos(response_table[2]) 

116 + 1.0j*num.sin(response_table[2])) 

117 else: 

118 cresp = response_table[1] * ( 

119 num.cos(response_table[2]*d2r) 

120 + 1.0j*num.sin(response_table[2]*d2r)) 

121 

122 self.station = station 

123 self.component = component 

124 self.tmin = tmin 

125 self.latitude = latitude 

126 self.longitude = longitude 

127 self.elevation = elevation 

128 self.filetype = filetype 

129 self.comment = comment 

130 self.period = period 

131 self.damping = damping 

132 self.sensor_sensitivity = sensor_sensitivity 

133 self.amplifier_gain = amplifier_gain 

134 self.digitizer_gain = digitizer_gain 

135 self.gain_1hz = gain_1hz 

136 self.filters = filters 

137 

138 self.sampled_response = response.SampledResponse( 

139 response_table[0], cresp) 

140 

141 self._check_tabulated_response(filename=filename) 

142 

143 def response(self, freqs, method='from-filetype', type='displacement'): 

144 

145 freqs = num.asarray(freqs) 

146 want_scalar = False 

147 if freqs.ndim == 0: 

148 freqs = num.array([freqs]) 

149 want_scalar = True 

150 

151 if method == 'from-filetype': 

152 method = self.filetype 

153 

154 if method == 'gains-and-filters': 

155 dresp = self._response_prototype_and_filters(freqs) \ 

156 / abs(self._response_prototype_and_filters([1.0])[0]) \ 

157 * self.gain_1hz 

158 

159 elif method == 'tabulated': 

160 dresp = self._response_tabulated(freqs) \ 

161 / abs(self._response_tabulated([1.0])[0]) * self.gain_1hz 

162 

163 elif method == 'poles-and-zeros': 

164 raise Exception( 

165 'fix me! poles-and-zeros response in seisan is broken: where ' 

166 'should the "normalization" come from?') 

167 

168 # dresp = self._response_from_poles_and_zeros(freqs) *normalization 

169 

170 elif method == 'gains': 

171 dresp = num.ones(freqs.size) * self.gain_total() * 1.0j * 2. \ 

172 * num.pi * freqs 

173 

174 else: 

175 assert False, 'invalid response method specified' 

176 

177 if type == 'velocity': 

178 dresp /= 1.0j * 2. * num.pi * freqs 

179 

180 if want_scalar: 

181 return dresp[0] 

182 else: 

183 return dresp 

184 

185 def gain_total(self): 

186 return self.sensor_sensitivity * 10.**(self.amplifier_gain/20.) \ 

187 * self.digitizer_gain 

188 

189 def _prototype_response_velocity(self, s): 

190 omega0 = 2. * num.pi / self.period 

191 return s**2/(omega0**2 + s**2 + 2.0*s*omega0*self.damping) 

192 

193 def _response_prototype_and_filters(self, freqs): 

194 freqs = num.asarray(freqs, dtype=float) 

195 iomega = 1.0j * 2. * num.pi * freqs 

196 

197 trans = iomega * self._prototype_response_velocity(iomega) 

198 

199 for (order, corner) in self.filters: 

200 if order < 0. or corner < 0.: 

201 b, a = signal.butter( 

202 abs(order), [abs(corner)], btype='high', analog=1) 

203 else: 

204 b, a = signal.butter( 

205 order, [corner], btype='low', analog=1) 

206 

207 trans *= signal.freqs(b, a, freqs)[1] 

208 

209 return trans 

210 

211 def _response_tabulated(self, freqs): 

212 freqs = num.asarray(freqs, dtype=float) 

213 return self.sampled_response.evaluate(freqs) 

214 

215 def _response_from_poles_and_zeros(self, freqs): 

216 assert False, 'poles-and-zeros file type not implemented yet for ' \ 

217 'seisan response file format' 

218 return None 

219 

220 def _check_tabulated_response(self, filename='?'): 

221 if self.filetype == 'gains-and-filters': 

222 freqs = self.sampled_response.frequencies() 

223 

224 trans_gaf = self.response(freqs, method='gains-and-filters') 

225 atrans_gaf = num.abs(trans_gaf) 

226 

227 trans_tab = self.response(freqs, method='tabulated') 

228 atrans_tab = num.abs(trans_tab) 

229 

230 if num.any(num.abs(atrans_gaf-atrans_tab) 

231 / num.abs(atrans_gaf+atrans_tab) > 1./100.): 

232 

233 logger.warning( 

234 'inconsistent amplitudes in tabulated response ' 

235 '(in file "%s") (max |a-b|/|a+b| is %g' % ( 

236 filename, 

237 num.amax(num.abs(atrans_gaf-atrans_tab) 

238 / abs(atrans_gaf+atrans_tab)))) 

239 

240 else: 

241 if num.any(num.abs(trans_gaf-trans_tab) 

242 > abs(trans_gaf+trans_tab)/100.): 

243 

244 logger.warning( 

245 'inconsistent phase values in tabulated response ' 

246 '(in file "%s"' % filename) 

247 

248 def __str__(self): 

249 

250 return '''--- Seisan Response File --- 

251station: %s 

252component: %s 

253start time: %s 

254latitude: %f 

255longitude: %f 

256elevation: %f 

257filetype: %s 

258comment: %s 

259sensor period: %g 

260sensor damping: %g 

261sensor sensitivity: %g 

262amplifier gain: %g 

263digitizer gain: %g 

264gain at 1 Hz: %g 

265filters: %s 

266''' % (self.station, self.component, util.time_to_str(self.tmin), 

267 self.latitude, self.longitude, self.elevation, self.filetype, 

268 self.comment, self.period, self.damping, self.sensor_sensitivity, 

269 self.amplifier_gain, self.digitizer_gain, self.gain_1hz, 

270 self.filters) 

271 

272 def plot_amplitudes(self, filename_pdf, type='displacement'): 

273 

274 from pyrocko.plot import gmtpy 

275 

276 p = gmtpy.LogLogPlot() 

277 f = self.sampled_response.frequencies() 

278 atab = num.abs(self.response(f, method='tabulated', type=type)) 

279 acom = num.abs(self.response(f, method='gains-and-filters', type=type)) 

280 aconst = num.abs(self.response(f, method='gains', type=type)) 

281 

282 p.plot((f, atab), '-W2p,red') 

283 p.plot((f, acom), '-W1p,blue') 

284 p.plot((f, aconst), '-W1p,green') 

285 p.save(filename_pdf) 

286 

287 def plot_phases(self, filename_pdf, type='displacement'): 

288 

289 from pyrocko.plot import gmtpy 

290 

291 p = gmtpy.LogLinPlot() 

292 f = self.sampled_response.frequencies() 

293 atab = num.unwrap(num.angle(self.response( 

294 f, method='tabulated', type=type))) / d2r 

295 acom = num.unwrap(num.angle(self.response( 

296 f, method='gains-and-filters', type=type))) / d2r 

297 aconst = num.unwrap(num.angle(self.response( 

298 f, method='gains', type=type))) / d2r 

299 

300 p.plot((f, atab), '-W2p,red') 

301 p.plot((f, acom), '-W1p,blue') 

302 p.plot((f, aconst), '-W1p,green') 

303 p.save(filename_pdf)