1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import math 

7import numpy as num 

8try: 

9 from StringIO import StringIO as BytesIO 

10except ImportError: 

11 from io import BytesIO 

12 

13from pyrocko import trace 

14 

15d2r = math.pi/180. 

16 

17 

18class SacPoleZeroError(Exception): 

19 pass 

20 

21 

22def read_sac_zpk(filename=None, file=None, string=None, get_comments=False): 

23 ''' 

24 Read SAC Pole-Zero file. 

25 

26 :returns: ``(zeros, poles, constant)`` or 

27 ``(zeros, poles, constant, comments)`` if ``get_comments`` is True. 

28 ''' 

29 

30 if filename is not None: 

31 f = open(filename, 'rb') 

32 

33 elif file is not None: 

34 f = file 

35 

36 elif string is not None: 

37 f = BytesIO(string) 

38 

39 sects = ('ZEROS', 'POLES', 'CONSTANT') 

40 sectdata = {'ZEROS': [], 'POLES': []} 

41 npoles = 0 

42 nzeros = 0 

43 constant = 1.0 

44 atsect = None 

45 comments = [] 

46 for iline, line in enumerate(f): 

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

48 toks = line.split() 

49 if len(toks) == 0: 

50 continue 

51 

52 if toks[0][0] in '*#': 

53 comments.append(line) 

54 continue 

55 

56 if len(toks) != 2: 

57 f.close() 

58 raise SacPoleZeroError( 

59 'Expected 2 tokens in line %i of file %s' 

60 % (iline+1, filename)) 

61 

62 if toks[0].startswith('*'): 

63 continue 

64 

65 lsect = toks[0].upper() 

66 if lsect in sects: 

67 atsect = lsect 

68 sectdata[atsect] = [] 

69 if lsect.upper() == 'ZEROS': 

70 nzeros = int(toks[1]) 

71 elif toks[0].upper() == 'POLES': 

72 npoles = int(toks[1]) 

73 elif toks[0].upper() == 'CONSTANT': 

74 constant = float(toks[1]) 

75 else: 

76 if atsect: 

77 sectdata[atsect].append( 

78 complex(float(toks[0]), float(toks[1]))) 

79 

80 if f != file: 

81 f.close() 

82 

83 poles = sectdata['POLES'] 

84 zeros = sectdata['ZEROS'] 

85 npoles_ = len(poles) 

86 nzeros_ = len(zeros) 

87 if npoles_ > npoles: 

88 raise SacPoleZeroError( 

89 'Expected %i poles but found %i in pole-zero file "%s"' 

90 % (npoles, npoles_, filename)) 

91 if nzeros_ > nzeros: 

92 raise SacPoleZeroError( 

93 'Expected %i zeros but found %i in pole-zero file "%s"' 

94 % (nzeros, nzeros_, filename)) 

95 

96 if npoles_ < npoles: 

97 poles.extend([complex(0.)]*(npoles-npoles_)) 

98 

99 if nzeros_ < npoles: 

100 zeros.extend([complex(0.)]*(nzeros-nzeros_)) 

101 

102 if len(poles) == 0 and len(zeros) == 0: 

103 raise SacPoleZeroError( 

104 'No poles and zeros found in file "%s"' % (filename)) 

105 

106 if not num.all(num.isfinite(poles)): 

107 raise SacPoleZeroError( 

108 'Not finite pole(s) found in pole-zero file "%s"' 

109 % filename) 

110 

111 if not num.all(num.isfinite(zeros)): 

112 raise SacPoleZeroError( 

113 'Not finite zero(s) found in pole-zero file "%s"' 

114 % filename) 

115 

116 if not num.isfinite(constant): 

117 raise SacPoleZeroError( 

118 'Ivalid constant (%g) found in pole-zero file "%s"' 

119 % (constant, filename)) 

120 

121 if get_comments: 

122 return zeros, poles, constant, comments 

123 else: 

124 return zeros, poles, constant 

125 

126 

127def write_sac_zpk(zeros, poles, constant, filename): 

128 if hasattr(filename, 'write'): 

129 f = filename 

130 else: 

131 f = open('w', filename) 

132 

133 def write_complex(x): 

134 f.write('%12.8g %12.8g\n' % (complex(x).real, complex(x).imag)) 

135 

136 f.write('POLES %i\n' % len(poles)) 

137 for p in poles: 

138 if p != 0.0: 

139 write_complex(p) 

140 

141 f.write('ZEROS %i\n' % len(zeros)) 

142 for z in zeros: 

143 if z != 0.0: 

144 write_complex(z) 

145 

146 f.write('CONSTANT %12.8g\n' % constant) 

147 if not hasattr(filename, 'write'): 

148 f.close() 

149 

150 

151def evaluate(zeros, poles, constant, fmin=0.001, fmax=100., nf=100): 

152 

153 logfmin = math.log(fmin) 

154 logfmax = math.log(fmax) 

155 logf = num.linspace(logfmin, logfmax, nf) 

156 f = num.exp(logf) 

157 trans = trace.PoleZeroResponse(zeros, poles, constant) 

158 return f, trans.evaluate(f) 

159 

160 

161def evaluate_at(zeros, poles, constant, f): 

162 jomeg = 1.0j * 2. * math.pi * f 

163 

164 a = constant 

165 for z in zeros: 

166 a *= jomeg-z 

167 for p in poles: 

168 a /= jomeg-p 

169 

170 return a 

171 

172 

173def read_to_pyrocko_response(filename=None, file=None, string=None): 

174 ''' 

175 Read SAC pole-zero file into Pyrocko response object. 

176 

177 :returns: Response as a :py:class:~pyrocko.trace.PoleZeroResponse` object. 

178 ''' 

179 

180 from pyrocko import trace 

181 

182 zeros, poles, constant = read_sac_zpk( 

183 filename=filename, file=file, string=string) 

184 return trace.PoleZeroResponse(zeros, poles, constant) 

185 

186 

187def read_to_stationxml_response( 

188 input_unit, output_unit, normalization_frequency=1.0, 

189 filename=None, file=None, string=None): 

190 ''' 

191 Read SAC pole-zero file into StationXML response object. 

192 

193 :param input_unit: Input unit to be reported in the StationXML response. 

194 :type input_unit: str 

195 :param output_unit: Output unit to be reported in the StationXML response. 

196 :type output_unit: str 

197 :param normalization_frequency: Frequency where the normalization factor 

198 for the StationXML response should be computed. 

199 :type normalization_frequency: float 

200 

201 :returns: Response as a :py:class:~pyrocko.io.stationxml.Response` object 

202 with a single pole-zero response stage. 

203 ''' 

204 

205 from pyrocko.io import stationxml 

206 

207 presponse = read_to_pyrocko_response( 

208 filename=filename, file=file, string=string) 

209 

210 return stationxml.Response.from_pyrocko_pz_response( 

211 presponse, input_unit, output_unit, normalization_frequency) 

212 

213 

214def plot_amplitudes_zpk( 

215 zpks, filename_pdf, 

216 fmin=0.001, 

217 fmax=100., 

218 nf=100, 

219 fnorm=None): 

220 

221 from pyrocko.plot import gmtpy 

222 

223 p = gmtpy.LogLogPlot(width=30*gmtpy.cm, yexp=0) 

224 for i, (zeros, poles, constant) in enumerate(zpks): 

225 f, h = evaluate(zeros, poles, constant, fmin, fmax, nf) 

226 if fnorm is not None: 

227 h /= evaluate_at(zeros, poles, constant, fnorm) 

228 

229 amp = num.abs(h) 

230 p.plot((f, amp), '-W2p,%s' % gmtpy.color(i)) 

231 

232 p.save(filename_pdf) 

233 

234 

235def plot_phases_zpk(zpks, filename_pdf, fmin=0.001, fmax=100., nf=100): 

236 

237 from pyrocko.plot import gmtpy 

238 

239 p = gmtpy.LogLinPlot(width=30*gmtpy.cm) 

240 for i, (zeros, poles, constant) in enumerate(zpks): 

241 f, h = evaluate(zeros, poles, constant, fmin, fmax, nf) 

242 phase = num.unwrap(num.angle(h)) / d2r 

243 p.plot((f, phase), '-W1p,%s' % gmtpy.color(i)) 

244 

245 p.save(filename_pdf)