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'''
34import logging
35import math
36import hashlib
37from io import BytesIO
39import numpy as num
41from pyrocko import util
42from pyrocko import guts
43from pyrocko.response import InvalidResponseError
46logger = logging.getLogger('plot.response')
49def normalize_on_flat(f, tf, factor=10000.):
50 df = num.diff(num.log(f))
51 tap = 1.0 / (1.0 + factor * (num.diff(num.log(num.abs(tf)))/df)**2)
52 return tf / (num.sum(num.abs(tf[1:]) * tap) / num.sum(tap))
55def draw(
56 response,
57 axes_amplitude=None, axes_phase=None,
58 fmin=0.01, fmax=100., nf=100,
59 normalize=False,
60 style={},
61 label=None,
62 show_breakpoints=False,
63 color_pool=None,
64 label_pool=None):
66 '''
67 Draw instrument response in Bode plot style to given Matplotlib axes
69 :param response: instrument response as a
70 :py:class:`pyrocko.response.FrequencyResponse` object
71 :param axes_amplitude: :py:class:`matplotlib.axes.Axes` object to use when
72 drawing the amplitude response
73 :param axes_phase: :py:class:`matplotlib.axes.Axes` object to use when
74 drawing the phase response
75 :param fmin: minimum frequency [Hz]
76 :param fmax: maximum frequency [Hz]
77 :param nf: number of frequencies where to evaluate the response
78 :param style: :py:class:`dict` with keyword arguments to tune the line
79 style
80 :param label: string to be passed to the ``label`` argument of
81 :py:meth:`matplotlib.axes.Axes.plot`
82 '''
84 f = num.exp(num.linspace(num.log(fmin), num.log(fmax), nf))
86 resp_fmax = response.get_fmax()
87 if resp_fmax is not None:
88 if fmax > resp_fmax:
89 logger.warning(
90 'Maximum frequency above range supported by response. '
91 'Clipping to supported%s.' % (
92 ' (%s)' % label if label else ''))
94 f = f[f <= resp_fmax]
96 if f.size == 0:
97 return
99 tf = response.evaluate(f)
100 ok = num.isfinite(tf)
101 if not num.all(ok):
102 logger.warning('NaN values present in evaluated response%s.' % (
103 ' (%s)' % label if label else ''))
104 f = f[ok]
105 tf = tf[ok]
107 if normalize:
108 tf = normalize_on_flat(f, tf)
110 ta = num.abs(tf)
112 if color_pool is not None:
113 fh = BytesIO()
114 f.dump(fh)
115 tf.dump(fh)
116 c_key = hashlib.sha1(fh.getvalue()).hexdigest()
117 fh.close()
118 if c_key not in color_pool[0]:
119 color_pool[0][c_key] \
120 = color_pool[1][len(color_pool[0]) % len(color_pool[1])]
121 new = True
122 else:
123 new = False
125 style = dict(style)
126 style['color'] = color_pool[0][c_key]
128 ikey = list(color_pool[0].keys()).index(c_key) + 1
129 slabel = '[%i]' % ikey
131 if label_pool is not None:
132 label_pool.append('(%s) %s' % (slabel, label))
133 eff_label = slabel if new else None
134 else:
135 eff_label = '%s %s' % (slabel, label)
137 if axes_amplitude:
138 axes_amplitude.plot(f, ta, label=eff_label, **style)
139 for checkpoint in response.checkpoints:
140 axes_amplitude.plot(
141 checkpoint.frequency, checkpoint.value, 'o',
142 color=style.get('color', 'black'))
144 axes_amplitude.annotate(
145 '%.3g s' % (1.0/checkpoint.frequency)
146 if checkpoint.frequency < 1.0 else
147 '%.3g Hz' % checkpoint.frequency,
148 xy=(checkpoint.frequency, checkpoint.value),
149 xytext=(10, 10),
150 textcoords='offset points',
151 color=style.get('color', 'black'))
153 if show_breakpoints:
154 for br_frequency, br_change in response.construction():
155 if not (fmin <= br_frequency <= fmax):
156 continue
158 br_value = abs(response.evaluate1(br_frequency))
159 axes_amplitude.plot(
160 br_frequency, br_value, 'v' if br_change < 0 else '^',
161 mec=style.get('color', 'black'),
162 color='none',
163 ms=10)
165 axes_amplitude.annotate(
166 '%.3g s (%i)' % (1.0/br_frequency, br_change)
167 if br_frequency < 1.0 else
168 '%.3g Hz' % br_frequency,
169 xy=(br_frequency, br_value),
170 xytext=(10, 10),
171 textcoords='offset points',
172 color=style.get('color', 'black'))
174 if axes_phase:
175 dta = num.diff(num.log(ta))
176 iflat = num.nanargmin(num.abs(num.diff(dta)) + num.abs(dta[:-1]))
177 tp = num.unwrap(num.angle(tf))
178 ioff = int(num.round(tp[iflat] / (2.0*num.pi)))
179 tp -= ioff * 2.0 * num.pi
180 axes_phase.plot(f, tp/num.pi, label=eff_label, **style)
181 else:
182 tp = [0.]
184 return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)
187def setup_axes(axes_amplitude=None, axes_phase=None):
188 '''
189 Configure axes in Bode plot style.
190 '''
192 if axes_amplitude is not None:
193 axes_amplitude.set_ylabel('Amplitude ratio')
194 axes_amplitude.set_xscale('log')
195 axes_amplitude.set_yscale('log')
196 axes_amplitude.grid(True)
197 axes_amplitude.axhline(1.0, lw=0.5, color='black')
198 if axes_phase is None:
199 axes_amplitude.set_xlabel('Frequency [Hz]')
200 axes_amplitude.set_xscale('log')
201 else:
202 axes_amplitude.xaxis.set_ticklabels([])
204 if axes_phase is not None:
205 axes_phase.set_ylabel('Phase [$\\pi$]')
206 axes_phase.set_xscale('log')
207 axes_phase.set_xlabel('Frequency [Hz]')
208 axes_phase.grid(True)
209 axes_phase.axhline(0.0, lw=0.5, color='black')
212def plot(
213 responses,
214 filename=None,
215 dpi=100,
216 fmin=0.01, fmax=100., nf=100,
217 normalize=False,
218 fontsize=10.,
219 figsize=None,
220 styles=None,
221 labels=None,
222 show_breakpoints=False,
223 separate_combined_labels=None):
225 '''
226 Draw instrument responses in Bode plot style.
228 :param responses: instrument responses as
229 :py:class:`pyrocko.response.FrequencyResponse` objects
230 :param fmin: minimum frequency [Hz]
231 :param fmax: maximum frequency [Hz]
232 :param nf: number of frequencies where to evaluate the response
233 :param normalize: if ``True`` normalize flat part of response to be ``1``
234 :param styles: :py:class:`list` of :py:class:`dict` objects with keyword
235 arguments to be passed to matplotlib's
236 :py:meth:`matplotlib.axes.Axes.plot` function when drawing the response
237 lines. Length must match number of responses.
238 :param filename: file name to pass to matplotlib's ``savefig`` function. If
239 ``None``, the plot is shown with :py:func:`matplotlib.pyplot.show`.
240 :param fontsize: font size in points used in axis labels and legend
241 :param figsize: :py:class:`tuple`, ``(width, height)`` in inches
242 :param labels: :py:class:`list` of names to show in legend. Length must
243 correspond to number of responses.
244 :param show_breakpoints: show breakpoints of pole-zero responses.
245 '''
247 from matplotlib import pyplot as plt
248 from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
249 from pyrocko.plot import graph_colors, to01
251 mpl_init(fontsize=fontsize)
253 if figsize is None:
254 figsize = mpl_papersize('a4', 'portrait')
256 fig = plt.figure(figsize=figsize)
257 labelpos = mpl_margins(
258 fig, w=7., h=5., units=fontsize, nw=1, nh=2, hspace=2.)
259 axes_amplitude = fig.add_subplot(2, 1, 1)
260 labelpos(axes_amplitude, 2., 1.5)
261 axes_phase = fig.add_subplot(2, 1, 2)
262 labelpos(axes_phase, 2., 1.5)
264 setup_axes(axes_amplitude, axes_phase)
266 if styles is None:
267 styles = [{} for i in range(len(responses))]
268 color_pool = ({}, [to01(color) for color in graph_colors])
270 else:
271 color_pool = None
272 assert len(styles) == len(responses)
274 if separate_combined_labels is not None:
275 label_pool = []
276 else:
277 label_pool = None
279 if labels is None:
280 labels = [None] * len(responses)
281 else:
282 assert len(labels) == len(responses)
284 a_ranges, p_ranges = [], []
285 have_labels = False
286 for style, resp, label in zip(styles, responses, labels):
287 try:
288 a_range, p_range = draw(
289 response=resp,
290 axes_amplitude=axes_amplitude,
291 axes_phase=axes_phase,
292 fmin=fmin, fmax=fmax, nf=nf,
293 normalize=normalize,
294 style=style,
295 label=label,
296 show_breakpoints=show_breakpoints,
297 color_pool=color_pool,
298 label_pool=label_pool)
300 if label is not None:
301 have_labels = True
303 a_ranges.append(a_range)
304 p_ranges.append(p_range)
306 except InvalidResponseError as e:
307 logger.error(
308 'Error occured for response labeled "%s": %s' % (label, e))
310 if have_labels:
311 axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))
313 if separate_combined_labels == 'print':
314 for label in label_pool:
315 print(label)
317 if a_ranges:
318 a_ranges = num.array(a_ranges)
319 p_ranges = num.array(p_ranges)
321 amin, amax = num.nanmin(a_ranges), num.nanmax(a_ranges)
322 pmin, pmax = num.nanmin(p_ranges), num.nanmax(p_ranges)
324 amin *= 0.5
325 amax *= 2.0
327 pmin -= 0.5
328 pmax += 0.5
330 pmin = math.floor(pmin)
331 pmax = math.ceil(pmax)
333 axes_amplitude.set_ylim(amin, amax)
334 axes_phase.set_ylim(pmin, pmax)
335 axes_amplitude.set_xlim(fmin, fmax)
336 axes_phase.set_xlim(fmin, fmax)
338 if filename is not None:
339 fig.savefig(filename, dpi=dpi)
340 else:
341 plt.show()
343 if separate_combined_labels == 'return':
344 return label_pool
347def tts(t):
348 if t is None:
349 return '?'
350 else:
351 return util.tts(t, format='%Y-%m-%d')
354def load_response_information(
355 filename, format,
356 nslc_patterns=None,
357 fake_input_units=None,
358 stages=None):
360 from pyrocko import pz, trace
361 from pyrocko.io import resp as fresp
363 resps = []
364 labels = []
365 if format == 'sacpz':
366 if fake_input_units is not None:
367 raise Exception(
368 'cannot guess true input units from plain SAC PZ files')
370 zeros, poles, constant = pz.read_sac_zpk(filename)
371 resp = trace.PoleZeroResponse(
372 zeros=zeros, poles=poles, constant=constant)
374 resps.append(resp)
375 labels.append(filename)
377 elif format == 'pf':
378 if fake_input_units is not None:
379 raise Exception(
380 'Cannot guess input units from plain response files.')
382 resp = guts.load(filename=filename)
383 resps.append(resp)
384 labels.append(filename)
386 elif format in ('resp', 'evalresp'):
387 for resp in list(fresp.iload_filename(filename)):
388 if nslc_patterns is not None and not util.match_nslc(
389 nslc_patterns, resp.codes):
390 continue
392 units = ''
393 in_units = None
394 if resp.response.instrument_sensitivity:
395 s = resp.response.instrument_sensitivity
396 in_units = fake_input_units or s.input_units.name
397 if s.input_units and s.output_units:
398 units = ', %s -> %s' % (in_units, s.output_units.name)
400 if format == 'resp':
401 resps.append(resp.response.get_pyrocko_response(
402 '.'.join(resp.codes),
403 fake_input_units=fake_input_units,
404 stages=stages).expect_one())
405 else:
406 target = {
407 'M/S': 'vel',
408 'M': 'dis',
409 }[in_units]
411 if resp.end_date is not None:
412 time = (resp.start_date + resp.end_date)*0.5
413 else:
414 time = resp.start_date + 3600*24*10
416 r = trace.Evalresp(
417 respfile=filename,
418 nslc_id=resp.codes,
419 target=target,
420 time=time,
421 stages=stages)
423 resps.append(r)
425 labels.append('%s (%s.%s.%s.%s, %s - %s%s)' % (
426 (filename, ) + resp.codes +
427 (tts(resp.start_date), tts(resp.end_date), units)))
429 elif format == 'stationxml':
430 from pyrocko.fdsn import station as fs
432 sx = fs.load_xml(filename=filename)
433 for network in sx.network_list:
434 for station in network.station_list:
435 for channel in station.channel_list:
436 nslc = (
437 network.code,
438 station.code,
439 channel.location_code,
440 channel.code)
442 if nslc_patterns is not None and not util.match_nslc(
443 nslc_patterns, nslc):
444 continue
446 if not channel.response:
447 logger.warning(
448 'no response for channel %s.%s.%s.%s given.'
449 % nslc)
450 continue
452 units = ''
453 if channel.response.instrument_sensitivity:
454 s = channel.response.instrument_sensitivity
455 if s.input_units and s.output_units:
456 units = ', %s -> %s' % (
457 fake_input_units or s.input_units.name,
458 s.output_units.name)
460 resps.append(channel.response.get_pyrocko_response(
461 '.'.join(nslc),
462 fake_input_units=fake_input_units,
463 stages=stages).expect_one())
465 labels.append(
466 '%s (%s.%s.%s.%s, %s - %s%s)' % (
467 (filename, ) + nslc +
468 (tts(channel.start_date),
469 tts(channel.end_date),
470 units)))
472 return resps, labels
475if __name__ == '__main__':
476 import sys
477 from optparse import OptionParser
479 util.setup_logging('pyrocko.plot.response.__main__', 'warning')
481 usage = 'python -m pyrocko.plot.response <filename> ... [options]'
483 description = '''Plot instrument responses (transfer functions).'''
485 allowed_formats = ['sacpz', 'resp', 'evalresp', 'stationxml', 'pf']
487 parser = OptionParser(
488 usage=usage,
489 description=description,
490 formatter=util.BetterHelpFormatter())
492 parser.add_option(
493 '--format',
494 dest='format',
495 default='sacpz',
496 choices=allowed_formats,
497 help='assume input files are of given FORMAT. Choices: %s' % (
498 ', '.join(allowed_formats)))
500 parser.add_option(
501 '--fmin',
502 dest='fmin',
503 type='float',
504 default=0.01,
505 help='minimum frequency [Hz], default: %default')
507 parser.add_option(
508 '--fmax',
509 dest='fmax',
510 type='float',
511 default=100.,
512 help='maximum frequency [Hz], default: %default')
514 parser.add_option(
515 '--normalize',
516 dest='normalize',
517 action='store_true',
518 help='normalize response to be 1 on flat part')
520 parser.add_option(
521 '--save',
522 dest='filename',
523 help='save figure to file with name FILENAME')
525 parser.add_option(
526 '--dpi',
527 dest='dpi',
528 type='float',
529 default=100.,
530 help='DPI setting for pixel image output, default: %default')
532 parser.add_option(
533 '--patterns',
534 dest='nslc_patterns',
535 metavar='NET.STA.LOC.CHA,...',
536 help='show only responses of channels matching any of the given '
537 'patterns')
539 parser.add_option(
540 '--stages',
541 dest='stages',
542 metavar='START:STOP',
543 help='show only responses selected stages')
545 parser.add_option(
546 '--fake-input-units',
547 dest='fake_input_units',
548 choices=['M', 'M/S', 'M/S**2'],
549 metavar='UNITS',
550 help='show converted response for given input units, choices: '
551 '["M", "M/S", "M/S**2"]')
553 parser.add_option(
554 '--show-breakpoints',
555 dest='show_breakpoints',
556 action='store_true',
557 default=False,
558 help='show breakpoints of pole-zero responses')
560 parser.add_option(
561 '--index-labels',
562 dest='index_labels',
563 action='store_true',
564 default=False,
565 help='label graphs only by index and print details to terminal '
566 'to save space when many labels would be shown. Aggregate '
567 'identical responses under a common index.')
569 (options, args) = parser.parse_args(sys.argv[1:])
571 if len(args) == 0:
572 parser.print_help()
573 sys.exit(1)
575 fns = args
577 resps = []
578 labels = []
580 stages = None
581 if options.stages:
582 stages = tuple(int(x) for x in options.stages.split(':', 1))
583 stages = stages[0]-1, stages[1]
585 for fn in fns:
587 if options.nslc_patterns is not None:
588 nslc_patterns = options.nslc_patterns.split(',')
589 else:
590 nslc_patterns = None
592 resps_this, labels_this = load_response_information(
593 fn, options.format, nslc_patterns,
594 fake_input_units=options.fake_input_units,
595 stages=stages)
597 resps.extend(resps_this)
598 labels.extend(labels_this)
600 plot(
601 resps,
602 fmin=options.fmin, fmax=options.fmax, nf=200,
603 normalize=options.normalize,
604 labels=labels, filename=options.filename, dpi=options.dpi,
605 show_breakpoints=options.show_breakpoints,
606 separate_combined_labels='print' if options.index_labels else None)