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
36import math
38import numpy as num
40from pyrocko import util
41from pyrocko import guts
44logger = logging.getLogger('plot.response')
47def normalize_on_flat(f, tf, factor=10000.):
48 df = num.diff(num.log(f))
49 tap = 1.0 / (1.0 + factor * (num.diff(num.log(num.abs(tf)))/df)**2)
50 return tf / (num.sum(num.abs(tf[1:]) * tap) / num.sum(tap))
53def draw(
54 response,
55 axes_amplitude=None, axes_phase=None,
56 fmin=0.01, fmax=100., nf=100,
57 normalize=False,
58 style={},
59 label=None,
60 show_breakpoints=False):
62 '''
63 Draw instrument response in Bode plot style to given Matplotlib axes
65 :param response: instrument response as a
66 :py:class:`pyrocko.response.FrequencyResponse` object
67 :param axes_amplitude: :py:class:`matplotlib.axes.Axes` object to use when
68 drawing the amplitude response
69 :param axes_phase: :py:class:`matplotlib.axes.Axes` object to use when
70 drawing the phase response
71 :param fmin: minimum frequency [Hz]
72 :param fmax: maximum frequency [Hz]
73 :param nf: number of frequencies where to evaluate the response
74 :param style: :py:class:`dict` with keyword arguments to tune the line
75 style
76 :param label: string to be passed to the ``label`` argument of
77 :py:meth:`matplotlib.axes.Axes.plot`
78 '''
80 f = num.exp(num.linspace(num.log(fmin), num.log(fmax), nf))
82 resp_fmax = response.get_fmax()
83 if resp_fmax is not None:
84 if fmax > resp_fmax:
85 logger.warning(
86 'Maximum frequency above range supported by response. '
87 'Clipping to supported%s.' % (
88 ' (%s)' % label if label else ''))
90 f = f[f <= resp_fmax]
92 if f.size == 0:
93 return
95 tf = response.evaluate(f)
96 ok = num.isfinite(tf)
97 if not num.all(ok):
98 logger.warning('NaN values present in evaluated response%s.' % (
99 ' (%s)' % label if label else ''))
100 f = f[ok]
101 tf = tf[ok]
103 if normalize:
104 tf = normalize_on_flat(f, tf)
106 ta = num.abs(tf)
108 if axes_amplitude:
109 axes_amplitude.plot(f, ta, label=label, **style)
110 for checkpoint in response.checkpoints:
111 axes_amplitude.plot(
112 checkpoint.frequency, checkpoint.value, 'o',
113 color=style.get('color', 'black'))
115 axes_amplitude.annotate(
116 '%.3g s' % (1.0/checkpoint.frequency)
117 if checkpoint.frequency < 1.0 else
118 '%.3g Hz' % checkpoint.frequency,
119 xy=(checkpoint.frequency, checkpoint.value),
120 xytext=(10, 10),
121 textcoords='offset points',
122 color=style.get('color', 'black'))
124 if show_breakpoints:
125 for br_frequency, br_change in response.construction():
126 if not (fmin <= br_frequency <= fmax):
127 continue
129 br_value = abs(response.evaluate1(br_frequency))
130 axes_amplitude.plot(
131 br_frequency, br_value, 'v' if br_change < 0 else '^',
132 mec=style.get('color', 'black'),
133 color='none',
134 ms=10)
136 axes_amplitude.annotate(
137 '%.3g s (%i)' % (1.0/br_frequency, br_change)
138 if br_frequency < 1.0 else
139 '%.3g Hz' % br_frequency,
140 xy=(br_frequency, br_value),
141 xytext=(10, 10),
142 textcoords='offset points',
143 color=style.get('color', 'black'))
145 if axes_phase:
146 dta = num.diff(num.log(ta))
147 iflat = num.nanargmin(num.abs(num.diff(dta)) + num.abs(dta[:-1]))
148 tp = num.unwrap(num.angle(tf))
149 ioff = int(num.round(tp[iflat] / (2.0*num.pi)))
150 tp -= ioff * 2.0 * num.pi
151 axes_phase.plot(f, tp/num.pi, label=label, **style)
152 else:
153 tp = [0.]
155 return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)
158def setup_axes(axes_amplitude=None, axes_phase=None):
159 '''
160 Configure axes in Bode plot style.
161 '''
163 if axes_amplitude is not None:
164 axes_amplitude.set_ylabel('Amplitude ratio')
165 axes_amplitude.set_xscale('log')
166 axes_amplitude.set_yscale('log')
167 axes_amplitude.grid(True)
168 axes_amplitude.axhline(1.0, lw=0.5, color='black')
169 if axes_phase is None:
170 axes_amplitude.set_xlabel('Frequency [Hz]')
171 axes_amplitude.set_xscale('log')
172 else:
173 axes_amplitude.xaxis.set_ticklabels([])
175 if axes_phase is not None:
176 axes_phase.set_ylabel('Phase [$\\pi$]')
177 axes_phase.set_xscale('log')
178 axes_phase.set_xlabel('Frequency [Hz]')
179 axes_phase.grid(True)
180 axes_phase.axhline(0.0, lw=0.5, color='black')
183def plot(
184 responses,
185 filename=None,
186 dpi=100,
187 fmin=0.01, fmax=100., nf=100,
188 normalize=False,
189 fontsize=10.,
190 figsize=None,
191 styles=None,
192 labels=None,
193 show_breakpoints=False):
195 '''
196 Draw instrument responses in Bode plot style.
198 :param responses: instrument responses as
199 :py:class:`pyrocko.response.FrequencyResponse` objects
200 :param fmin: minimum frequency [Hz]
201 :param fmax: maximum frequency [Hz]
202 :param nf: number of frequencies where to evaluate the response
203 :param normalize: if ``True`` normalize flat part of response to be ``1``
204 :param styles: :py:class:`list` of :py:class:`dict` objects with keyword
205 arguments to be passed to matplotlib's
206 :py:meth:`matplotlib.axes.Axes.plot` function when drawing the response
207 lines. Length must match number of responses.
208 :param filename: file name to pass to matplotlib's ``savefig`` function. If
209 ``None``, the plot is shown with :py:func:`matplotlib.pyplot.show`.
210 :param fontsize: font size in points used in axis labels and legend
211 :param figsize: :py:class:`tuple`, ``(width, height)`` in inches
212 :param labels: :py:class:`list` of names to show in legend. Length must
213 correspond to number of responses.
214 :param show_breakpoints: show breakpoints of pole-zero responses.
215 '''
217 from matplotlib import pyplot as plt
218 from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
219 from pyrocko.plot import graph_colors, to01
221 mpl_init(fontsize=fontsize)
223 if figsize is None:
224 figsize = mpl_papersize('a4', 'portrait')
226 fig = plt.figure(figsize=figsize)
227 labelpos = mpl_margins(
228 fig, w=7., h=5., units=fontsize, nw=1, nh=2, hspace=2.)
229 axes_amplitude = fig.add_subplot(2, 1, 1)
230 labelpos(axes_amplitude, 2., 1.5)
231 axes_phase = fig.add_subplot(2, 1, 2)
232 labelpos(axes_phase, 2., 1.5)
234 setup_axes(axes_amplitude, axes_phase)
236 if styles is None:
237 styles = [
238 dict(color=to01(graph_colors[i % len(graph_colors)]))
239 for i in range(len(responses))]
240 else:
241 assert len(styles) == len(responses)
243 if labels is None:
244 labels = [None] * len(responses)
245 else:
246 assert len(labels) == len(responses)
248 a_ranges, p_ranges = [], []
249 have_labels = False
250 for style, resp, label in zip(styles, responses, labels):
251 a_range, p_range = draw(
252 response=resp,
253 axes_amplitude=axes_amplitude,
254 axes_phase=axes_phase,
255 fmin=fmin, fmax=fmax, nf=nf,
256 normalize=normalize,
257 style=style,
258 label=label,
259 show_breakpoints=show_breakpoints)
261 if label is not None:
262 have_labels = True
264 a_ranges.append(a_range)
265 p_ranges.append(p_range)
267 if have_labels:
268 axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))
270 if a_ranges:
271 a_ranges = num.array(a_ranges)
272 p_ranges = num.array(p_ranges)
274 amin, amax = num.nanmin(a_ranges), num.nanmax(a_ranges)
275 pmin, pmax = num.nanmin(p_ranges), num.nanmax(p_ranges)
277 amin *= 0.5
278 amax *= 2.0
280 pmin -= 0.5
281 pmax += 0.5
283 pmin = math.floor(pmin)
284 pmax = math.ceil(pmax)
286 axes_amplitude.set_ylim(amin, amax)
287 axes_phase.set_ylim(pmin, pmax)
288 axes_amplitude.set_xlim(fmin, fmax)
289 axes_phase.set_xlim(fmin, fmax)
291 if filename is not None:
292 fig.savefig(filename, dpi=dpi)
293 else:
294 plt.show()
297def tts(t):
298 if t is None:
299 return '?'
300 else:
301 return util.tts(t, format='%Y-%m-%d')
304def load_response_information(
305 filename, format,
306 nslc_patterns=None, fake_input_units=None, stages=None):
308 from pyrocko import pz, trace
309 from pyrocko.io import resp as fresp
311 resps = []
312 labels = []
313 if format == 'sacpz':
314 if fake_input_units is not None:
315 raise Exception(
316 'cannot guess true input units from plain SAC PZ files')
318 zeros, poles, constant = pz.read_sac_zpk(filename)
319 resp = trace.PoleZeroResponse(
320 zeros=zeros, poles=poles, constant=constant)
322 resps.append(resp)
323 labels.append(filename)
325 elif format == 'pf':
326 if fake_input_units is not None:
327 raise Exception(
328 'Cannot guess input units from plain response files.')
330 resp = guts.load(filename=filename)
331 resps.append(resp)
332 labels.append(filename)
334 elif format in ('resp', 'evalresp'):
335 for resp in list(fresp.iload_filename(filename)):
336 if nslc_patterns is not None and not util.match_nslc(
337 nslc_patterns, resp.codes):
338 continue
340 units = ''
341 in_units = None
342 if resp.response.instrument_sensitivity:
343 s = resp.response.instrument_sensitivity
344 in_units = fake_input_units or s.input_units.name
345 if s.input_units and s.output_units:
346 units = ', %s -> %s' % (in_units, s.output_units.name)
348 if format == 'resp':
349 resps.append(resp.response.get_pyrocko_response(
350 '.'.join(resp.codes),
351 fake_input_units=fake_input_units,
352 stages=stages).expect_one())
353 else:
354 target = {
355 'M/S': 'vel',
356 'M': 'dis',
357 }[in_units]
359 if resp.end_date is not None:
360 time = (resp.start_date + resp.end_date)*0.5
361 else:
362 time = resp.start_date + 3600*24*10
364 r = trace.Evalresp(
365 respfile=filename,
366 nslc_id=resp.codes,
367 target=target,
368 time=time,
369 stages=stages)
371 resps.append(r)
373 labels.append('%s (%s.%s.%s.%s, %s - %s%s)' % (
374 (filename, ) + resp.codes +
375 (tts(resp.start_date), tts(resp.end_date), units)))
377 elif format == 'stationxml':
378 from pyrocko.fdsn import station as fs
380 sx = fs.load_xml(filename=filename)
381 for network in sx.network_list:
382 for station in network.station_list:
383 for channel in station.channel_list:
384 nslc = (
385 network.code,
386 station.code,
387 channel.location_code,
388 channel.code)
390 if nslc_patterns is not None and not util.match_nslc(
391 nslc_patterns, nslc):
392 continue
394 if not channel.response:
395 logger.warning(
396 'no response for channel %s.%s.%s.%s given.'
397 % nslc)
398 continue
400 units = ''
401 if channel.response.instrument_sensitivity:
402 s = channel.response.instrument_sensitivity
403 if s.input_units and s.output_units:
404 units = ', %s -> %s' % (
405 fake_input_units or s.input_units.name,
406 s.output_units.name)
408 resps.append(channel.response.get_pyrocko_response(
409 '.'.join(nslc),
410 fake_input_units=fake_input_units,
411 stages=stages).expect_one())
413 labels.append(
414 '%s (%s.%s.%s.%s, %s - %s%s)' % (
415 (filename, ) + nslc +
416 (tts(channel.start_date),
417 tts(channel.end_date),
418 units)))
420 return resps, labels
423if __name__ == '__main__':
424 import sys
425 from optparse import OptionParser
427 util.setup_logging('pyrocko.plot.response.__main__', 'warning')
429 usage = 'python -m pyrocko.plot.response <filename> ... [options]'
431 description = '''Plot instrument responses (transfer functions).'''
433 allowed_formats = ['sacpz', 'resp', 'evalresp', 'stationxml', 'pf']
435 parser = OptionParser(
436 usage=usage,
437 description=description,
438 formatter=util.BetterHelpFormatter())
440 parser.add_option(
441 '--format',
442 dest='format',
443 default='sacpz',
444 choices=allowed_formats,
445 help='assume input files are of given FORMAT. Choices: %s' % (
446 ', '.join(allowed_formats)))
448 parser.add_option(
449 '--fmin',
450 dest='fmin',
451 type='float',
452 default=0.01,
453 help='minimum frequency [Hz], default: %default')
455 parser.add_option(
456 '--fmax',
457 dest='fmax',
458 type='float',
459 default=100.,
460 help='maximum frequency [Hz], default: %default')
462 parser.add_option(
463 '--normalize',
464 dest='normalize',
465 action='store_true',
466 help='normalize response to be 1 on flat part')
468 parser.add_option(
469 '--save',
470 dest='filename',
471 help='save figure to file with name FILENAME')
473 parser.add_option(
474 '--dpi',
475 dest='dpi',
476 type='float',
477 default=100.,
478 help='DPI setting for pixel image output, default: %default')
480 parser.add_option(
481 '--patterns',
482 dest='nslc_patterns',
483 metavar='NET.STA.LOC.CHA,...',
484 help='show only responses of channels matching any of the given '
485 'patterns')
487 parser.add_option(
488 '--stages',
489 dest='stages',
490 metavar='START:STOP',
491 help='show only responses selected stages')
493 parser.add_option(
494 '--fake-input-units',
495 dest='fake_input_units',
496 choices=['M', 'M/S', 'M/S**2'],
497 metavar='UNITS',
498 help='show converted response for given input units, choices: '
499 '["M", "M/S", "M/S**2"]')
501 parser.add_option(
502 '--show-breakpoints',
503 dest='show_breakpoints',
504 action='store_true',
505 default=False,
506 help='show breakpoints of pole-zero responses')
508 (options, args) = parser.parse_args(sys.argv[1:])
510 if len(args) == 0:
511 parser.print_help()
512 sys.exit(1)
514 fns = args
516 resps = []
517 labels = []
519 stages = None
520 if options.stages:
521 stages = tuple(int(x) for x in options.stages.split(':', 1))
522 stages = stages[0]-1, stages[1]
524 for fn in fns:
526 if options.nslc_patterns is not None:
527 nslc_patterns = options.nslc_patterns.split(',')
528 else:
529 nslc_patterns = None
531 resps_this, labels_this = load_response_information(
532 fn, options.format, nslc_patterns,
533 fake_input_units=options.fake_input_units,
534 stages=stages)
536 resps.extend(resps_this)
537 labels.extend(labels_this)
539 plot(
540 resps,
541 fmin=options.fmin, fmax=options.fmax, nf=200,
542 normalize=options.normalize,
543 labels=labels, filename=options.filename, dpi=options.dpi,
544 show_breakpoints=options.show_breakpoints)