1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5'''
6This module contains functions to plot instrument response transfer functions
7in Bode plot style using Matplotlib.
9Example
11* :download:`test_response_plot.py </../../examples/test_response_plot.py>`
13::
15 from pyrocko.plot import response
16 from pyrocko.example import get_example_data
18 get_example_data('test_response.resp')
20 resps, labels = response.load_response_information(
21 'test_response.resp', 'resp')
23 response.plot(
24 responses=resps, labels=labels, filename='test_response.png',
25 fmin=0.001, fmax=400., dpi=75.)
28.. figure :: /static/test_response.png
29 :align: center
31 Example response plot
32'''
33from __future__ import absolute_import
35import logging
37import numpy as num
39from pyrocko import util
40from pyrocko import guts
43logger = logging.getLogger('plot.response')
46def normalize_on_flat(f, tf, factor=10000.):
47 df = num.diff(num.log(f))
48 tap = 1.0 / (1.0 + factor * (num.diff(num.log(num.abs(tf)))/df)**2)
49 return tf / (num.sum(num.abs(tf[1:]) * tap) / num.sum(tap))
52def draw(
53 response,
54 axes_amplitude=None, axes_phase=None,
55 fmin=0.01, fmax=100., nf=100,
56 normalize=False,
57 style={},
58 label=None):
60 '''
61 Draw instrument response in Bode plot style to given Matplotlib axes
63 :param response: instrument response as a
64 :py:class:`pyrocko.trace.FrequencyResponse` object
65 :param axes_amplitude: :py:class:`matplotlib.axes.Axes` object to use when
66 drawing the amplitude response
67 :param axes_phase: :py:class:`matplotlib.axes.Axes` object to use when
68 drawing the phase response
69 :param fmin: minimum frequency [Hz]
70 :param fmax: maximum frequency [Hz]
71 :param nf: number of frequencies where to evaluate the response
72 :param style: :py:class:`dict` with keyword arguments to tune the line
73 style
74 :param label: string to be passed to the ``label`` argument of
75 :py:meth:`matplotlib.axes.Axes.plot`
76 '''
78 f = num.exp(num.linspace(num.log(fmin), num.log(fmax), nf))
79 tf = response.evaluate(f)
81 if normalize:
82 tf = normalize_on_flat(f, tf)
84 ta = num.abs(tf)
86 if axes_amplitude:
87 axes_amplitude.plot(f, ta, label=label, **style)
89 if axes_phase:
90 dta = num.diff(num.log(ta))
91 iflat = num.argmin(num.abs(num.diff(dta)) + num.abs(dta[:-1]))
92 tp = num.unwrap(num.angle(tf))
93 ioff = int(num.round(tp[iflat] / (2.0*num.pi)))
94 tp -= ioff * 2.0 * num.pi
95 axes_phase.plot(f, tp/num.pi, label=label, **style)
96 else:
97 tp = [0.]
99 return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)
102def setup_axes(axes_amplitude=None, axes_phase=None):
103 '''
104 Configure axes in Bode plot style.
105 '''
107 if axes_amplitude is not None:
108 axes_amplitude.set_ylabel('Amplitude ratio')
109 axes_amplitude.set_xscale('log')
110 axes_amplitude.set_yscale('log')
111 axes_amplitude.grid(True)
112 axes_amplitude.axhline(1.0, lw=0.5, color='black')
113 if axes_phase is None:
114 axes_amplitude.set_xlabel('Frequency [Hz]')
115 axes_amplitude.set_xscale('log')
116 else:
117 axes_amplitude.xaxis.set_ticklabels([])
119 if axes_phase is not None:
120 axes_phase.set_ylabel('Phase [$\\pi$]')
121 axes_phase.set_xscale('log')
122 axes_phase.set_xlabel('Frequency [Hz]')
123 axes_phase.grid(True)
124 axes_phase.axhline(0.0, lw=0.5, color='black')
127def plot(
128 responses,
129 filename=None,
130 dpi=100,
131 fmin=0.01, fmax=100., nf=100,
132 normalize=False,
133 fontsize=10.,
134 figsize=None,
135 styles=None,
136 labels=None):
138 '''
139 Draw instrument responses in Bode plot style.
141 :param responses: instrument responses as
142 :py:class:`pyrocko.trace.FrequencyResponse` objects
143 :param fmin: minimum frequency [Hz]
144 :param fmax: maximum frequency [Hz]
145 :param nf: number of frequencies where to evaluate the response
146 :param normalize: if ``True`` normalize flat part of response to be ``1``
147 :param styles: :py:class:`list` of :py:class:`dict` objects with keyword
148 arguments to be passed to matplotlib's
149 :py:meth:`matplotlib.axes.Axes.plot` function when drawing the response
150 lines. Length must match number of responses.
151 :param filename: file name to pass to matplotlib's ``savefig`` function. If
152 ``None``, the plot is shown with :py:func:`matplotlib.pyplot.show`.
153 :param fontsize: font size in points used in axis labels and legend
154 :param figsize: :py:class:`tuple`, ``(width, height)`` in inches
155 :param labels: :py:class:`list` of names to show in legend. Length must
156 correspond to number of responses.
157 '''
159 from matplotlib import pyplot as plt
160 from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
161 from pyrocko.plot import graph_colors, to01
163 mpl_init(fontsize=fontsize)
165 if figsize is None:
166 figsize = mpl_papersize('a4', 'portrait')
168 fig = plt.figure(figsize=figsize)
169 labelpos = mpl_margins(
170 fig, w=7., h=5., units=fontsize, nw=1, nh=2, hspace=2.)
171 axes_amplitude = fig.add_subplot(2, 1, 1)
172 labelpos(axes_amplitude, 2., 1.5)
173 axes_phase = fig.add_subplot(2, 1, 2)
174 labelpos(axes_phase, 2., 1.5)
176 setup_axes(axes_amplitude, axes_phase)
178 if styles is None:
179 styles = [
180 dict(color=to01(graph_colors[i % len(graph_colors)]))
181 for i in range(len(responses))]
182 else:
183 assert len(styles) == len(responses)
185 if labels is None:
186 labels = [None] * len(responses)
187 else:
188 assert len(labels) == len(responses)
190 a_ranges, p_ranges = [], []
191 have_labels = False
192 for style, resp, label in zip(styles, responses, labels):
193 a_range, p_range = draw(
194 response=resp,
195 axes_amplitude=axes_amplitude,
196 axes_phase=axes_phase,
197 fmin=fmin, fmax=fmax, nf=nf,
198 normalize=normalize,
199 style=style,
200 label=label)
202 if label is not None:
203 have_labels = True
205 a_ranges.append(a_range)
206 p_ranges.append(p_range)
208 if have_labels:
209 axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))
211 if a_ranges:
212 a_ranges = num.array(a_ranges)
213 p_ranges = num.array(p_ranges)
215 amin, amax = num.min(a_ranges), num.max(a_ranges)
216 pmin, pmax = num.min(p_ranges), num.max(p_ranges)
218 amin *= 0.5
219 amax *= 2.0
221 pmin -= 0.5
222 pmax += 0.5
224 axes_amplitude.set_ylim(amin, amax)
225 axes_phase.set_ylim(pmin, pmax)
226 axes_amplitude.set_xlim(fmin, fmax)
227 axes_phase.set_xlim(fmin, fmax)
229 if filename is not None:
230 fig.savefig(filename, dpi=dpi)
231 else:
232 plt.show()
235def tts(t):
236 if t is None:
237 return '?'
238 else:
239 return util.tts(t, format='%Y-%m-%d')
242def load_response_information(
243 filename, format, nslc_patterns=None, fake_input_units=None):
245 from pyrocko import pz, trace
246 from pyrocko.io import resp as fresp
248 resps = []
249 labels = []
250 if format == 'sacpz':
251 if fake_input_units is not None:
252 raise Exception(
253 'cannot guess true input units from plain SAC PZ files')
255 zeros, poles, constant = pz.read_sac_zpk(filename)
256 resp = trace.PoleZeroResponse(
257 zeros=zeros, poles=poles, constant=constant)
259 resps.append(resp)
260 labels.append(filename)
262 elif format == 'pf':
263 if fake_input_units is not None:
264 raise Exception(
265 'Cannot guess input units from plain response files.')
267 resp = guts.load(filename=filename)
268 resps.append(resp)
269 labels.append(filename)
271 elif format == 'resp':
272 for resp in list(fresp.iload_filename(filename)):
273 if nslc_patterns is not None and not util.match_nslc(
274 nslc_patterns, resp.codes):
275 continue
277 units = ''
278 if resp.response.instrument_sensitivity:
279 s = resp.response.instrument_sensitivity
280 if s.input_units and s.output_units:
281 units = ', %s -> %s' % (
282 fake_input_units or s.input_units.name,
283 s.output_units.name)
285 resps.append(resp.response.get_pyrocko_response(
286 resp.codes, fake_input_units=fake_input_units))
288 labels.append('%s (%s.%s.%s.%s, %s - %s%s)' % (
289 (filename, ) + resp.codes +
290 (tts(resp.start_date), tts(resp.end_date), units)))
292 elif format == 'stationxml':
293 from pyrocko.fdsn import station as fs
295 sx = fs.load_xml(filename=filename)
296 for network in sx.network_list:
297 for station in network.station_list:
298 for channel in station.channel_list:
299 nslc = (
300 network.code,
301 station.code,
302 channel.location_code,
303 channel.code)
305 if nslc_patterns is not None and not util.match_nslc(
306 nslc_patterns, nslc):
307 continue
309 if not channel.response:
310 logger.warn(
311 'no response for channel %s.%s.%s.%s given.'
312 % nslc)
313 continue
315 units = ''
316 if channel.response.instrument_sensitivity:
317 s = channel.response.instrument_sensitivity
318 if s.input_units and s.output_units:
319 units = ', %s -> %s' % (
320 fake_input_units or s.input_units.name,
321 s.output_units.name)
323 resps.append(channel.response.get_pyrocko_response(
324 nslc, fake_input_units=fake_input_units))
326 labels.append(
327 '%s (%s.%s.%s.%s, %s - %s%s)' % (
328 (filename, ) + nslc +
329 (tts(channel.start_date),
330 tts(channel.end_date),
331 units)))
333 return resps, labels
336if __name__ == '__main__':
337 import sys
338 from optparse import OptionParser
340 util.setup_logging('pyrocko.plot.response.__main__', 'warning')
342 usage = 'python -m pyrocko.plot.response <filename> ... [options]'
344 description = '''Plot instrument responses (transfer functions).'''
346 allowed_formats = ['sacpz', 'resp', 'stationxml', 'pf']
348 parser = OptionParser(
349 usage=usage,
350 description=description,
351 formatter=util.BetterHelpFormatter())
353 parser.add_option(
354 '--format',
355 dest='format',
356 default='sacpz',
357 choices=allowed_formats,
358 help='assume input files are of given FORMAT. Choices: %s' % (
359 ', '.join(allowed_formats)))
361 parser.add_option(
362 '--fmin',
363 dest='fmin',
364 type='float',
365 default=0.01,
366 help='minimum frequency [Hz], default: %default')
368 parser.add_option(
369 '--fmax',
370 dest='fmax',
371 type='float',
372 default=100.,
373 help='maximum frequency [Hz], default: %default')
375 parser.add_option(
376 '--normalize',
377 dest='normalize',
378 action='store_true',
379 help='normalize response to be 1 on flat part')
381 parser.add_option(
382 '--save',
383 dest='filename',
384 help='save figure to file with name FILENAME')
386 parser.add_option(
387 '--dpi',
388 dest='dpi',
389 type='float',
390 default=100.,
391 help='DPI setting for pixel image output, default: %default')
393 parser.add_option(
394 '--patterns',
395 dest='nslc_patterns',
396 metavar='NET.STA.LOC.CHA,...',
397 help='show only responses of channels matching any of the given '
398 'patterns')
400 parser.add_option(
401 '--fake-input-units',
402 dest='fake_input_units',
403 choices=['M', 'M/S', 'M/S**2'],
404 metavar='UNITS',
405 help='show converted response for given input units, choices: '
406 '["M", "M/S", "M/S**2"]')
408 (options, args) = parser.parse_args(sys.argv[1:])
410 if len(args) == 0:
411 parser.print_help()
412 sys.exit(1)
414 fns = args
416 resps = []
417 labels = []
419 for fn in fns:
421 if options.nslc_patterns is not None:
422 nslc_patterns = options.nslc_patterns.split(',')
423 else:
424 nslc_patterns = None
426 resps_this, labels_this = load_response_information(
427 fn, options.format, nslc_patterns,
428 fake_input_units=options.fake_input_units)
430 resps.extend(resps_this)
431 labels.extend(labels_this)
433 plot(
434 resps,
435 fmin=options.fmin, fmax=options.fmax, nf=200,
436 normalize=options.normalize,
437 labels=labels, filename=options.filename, dpi=options.dpi)