1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import numpy as num
8try:
9 from StringIO import StringIO as BytesIO
10except ImportError:
11 from io import BytesIO
13from pyrocko import trace
15d2r = math.pi/180.
18class SacPoleZeroError(Exception):
19 pass
22def read_sac_zpk(filename=None, file=None, string=None, get_comments=False):
23 '''
24 Read SAC Pole-Zero file.
26 :returns: ``(zeros, poles, constant)`` or
27 ``(zeros, poles, constant, comments)`` if ``get_comments`` is True.
28 '''
30 if filename is not None:
31 f = open(filename, 'rb')
33 elif file is not None:
34 f = file
36 elif string is not None:
37 f = BytesIO(string)
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
52 if toks[0][0] in '*#':
53 comments.append(line)
54 continue
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))
62 if toks[0].startswith('*'):
63 continue
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])))
80 if f != file:
81 f.close()
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))
96 if npoles_ < npoles:
97 poles.extend([complex(0.)]*(npoles-npoles_))
99 if nzeros_ < npoles:
100 zeros.extend([complex(0.)]*(nzeros-nzeros_))
102 if len(poles) == 0 and len(zeros) == 0:
103 raise SacPoleZeroError(
104 'No poles and zeros found in file "%s"' % (filename))
106 if not num.all(num.isfinite(poles)):
107 raise SacPoleZeroError(
108 'Not finite pole(s) found in pole-zero file "%s"'
109 % filename)
111 if not num.all(num.isfinite(zeros)):
112 raise SacPoleZeroError(
113 'Not finite zero(s) found in pole-zero file "%s"'
114 % filename)
116 if not num.isfinite(constant):
117 raise SacPoleZeroError(
118 'Ivalid constant (%g) found in pole-zero file "%s"'
119 % (constant, filename))
121 if get_comments:
122 return zeros, poles, constant, comments
123 else:
124 return zeros, poles, constant
127def write_sac_zpk(zeros, poles, constant, filename):
128 if hasattr(filename, 'write'):
129 f = filename
130 else:
131 f = open('w', filename)
133 def write_complex(x):
134 f.write('%12.8g %12.8g\n' % (complex(x).real, complex(x).imag))
136 f.write('POLES %i\n' % len(poles))
137 for p in poles:
138 if p != 0.0:
139 write_complex(p)
141 f.write('ZEROS %i\n' % len(zeros))
142 for z in zeros:
143 if z != 0.0:
144 write_complex(z)
146 f.write('CONSTANT %12.8g\n' % constant)
147 if not hasattr(filename, 'write'):
148 f.close()
151def evaluate(zeros, poles, constant, fmin=0.001, fmax=100., nf=100):
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)
161def evaluate_at(zeros, poles, constant, f):
162 jomeg = 1.0j * 2. * math.pi * f
164 a = constant
165 for z in zeros:
166 a *= jomeg-z
167 for p in poles:
168 a /= jomeg-p
170 return a
173def read_to_pyrocko_response(filename=None, file=None, string=None):
174 '''
175 Read SAC pole-zero file into Pyrocko response object.
177 :returns: Response as a :py:class:~pyrocko.trace.PoleZeroResponse` object.
178 '''
180 from pyrocko import trace
182 zeros, poles, constant = read_sac_zpk(
183 filename=filename, file=file, string=string)
184 return trace.PoleZeroResponse(zeros, poles, constant)
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.
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
201 :returns: Response as a :py:class:~pyrocko.io.stationxml.Response` object
202 with a single pole-zero response stage.
203 '''
205 from pyrocko.io import stationxml
207 presponse = read_to_pyrocko_response(
208 filename=filename, file=file, string=string)
210 return stationxml.Response.from_pyrocko_pz_response(
211 presponse, input_unit, output_unit, normalization_frequency)
214def plot_amplitudes_zpk(
215 zpks, filename_pdf,
216 fmin=0.001,
217 fmax=100.,
218 nf=100,
219 fnorm=None):
221 from pyrocko.plot import gmtpy
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)
229 amp = num.abs(h)
230 p.plot((f, amp), '-W2p,%s' % gmtpy.color(i))
232 p.save(filename_pdf)
235def plot_phases_zpk(zpks, filename_pdf, fmin=0.001, fmax=100., nf=100):
237 from pyrocko.plot import gmtpy
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))
245 p.save(filename_pdf)