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
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,
59 show_breakpoints=False):
61 '''
62 Draw instrument response in Bode plot style to given Matplotlib axes
64 :param response: instrument response as a
65 :py:class:`pyrocko.response.FrequencyResponse` object
66 :param axes_amplitude: :py:class:`matplotlib.axes.Axes` object to use when
67 drawing the amplitude response
68 :param axes_phase: :py:class:`matplotlib.axes.Axes` object to use when
69 drawing the phase response
70 :param fmin: minimum frequency [Hz]
71 :param fmax: maximum frequency [Hz]
72 :param nf: number of frequencies where to evaluate the response
73 :param style: :py:class:`dict` with keyword arguments to tune the line
74 style
75 :param label: string to be passed to the ``label`` argument of
76 :py:meth:`matplotlib.axes.Axes.plot`
77 '''
79 f = num.exp(num.linspace(num.log(fmin), num.log(fmax), nf))
81 resp_fmax = response.get_fmax()
82 if resp_fmax is not None:
83 if fmax > resp_fmax:
84 logger.warning(
85 'Maximum frequency above range supported by response. '
86 'Clipping to supported%s.' % (
87 ' (%s)' % label if label else ''))
89 f = f[f <= resp_fmax]
91 if f.size == 0:
92 return
94 tf = response.evaluate(f)
95 ok = num.isfinite(tf)
96 if not num.all(ok):
97 logger.warning('NaN values present in evaluated response%s.' % (
98 ' (%s)' % label if label else ''))
99 f = f[ok]
100 tf = tf[ok]
102 if normalize:
103 tf = normalize_on_flat(f, tf)
105 ta = num.abs(tf)
107 if axes_amplitude:
108 axes_amplitude.plot(f, ta, label=label, **style)
109 for checkpoint in response.checkpoints:
110 axes_amplitude.plot(
111 checkpoint.frequency, checkpoint.value, 'o',
112 color=style.get('color', 'black'))
114 axes_amplitude.annotate(
115 '%.3g s' % (1.0/checkpoint.frequency)
116 if checkpoint.frequency < 1.0 else
117 '%.3g Hz' % checkpoint.frequency,
118 xy=(checkpoint.frequency, checkpoint.value),
119 xytext=(10, 10),
120 textcoords='offset points',
121 color=style.get('color', 'black'))
123 if show_breakpoints:
124 for br_frequency, br_change in response.construction():
125 if not (fmin <= br_frequency <= fmax):
126 continue
128 br_value = abs(response.evaluate1(br_frequency))
129 axes_amplitude.plot(
130 br_frequency, br_value, 'v' if br_change < 0 else '^',
131 mec=style.get('color', 'black'),
132 color='none',
133 ms=10)
135 axes_amplitude.annotate(
136 '%.3g s (%i)' % (1.0/br_frequency, br_change)
137 if br_frequency < 1.0 else
138 '%.3g Hz' % br_frequency,
139 xy=(br_frequency, br_value),
140 xytext=(10, 10),
141 textcoords='offset points',
142 color=style.get('color', 'black'))
144 if axes_phase:
145 dta = num.diff(num.log(ta))
146 iflat = num.nanargmin(num.abs(num.diff(dta)) + num.abs(dta[:-1]))
147 tp = num.unwrap(num.angle(tf))
148 ioff = int(num.round(tp[iflat] / (2.0*num.pi)))
149 tp -= ioff * 2.0 * num.pi
150 axes_phase.plot(f, tp/num.pi, label=label, **style)
151 else:
152 tp = [0.]
154 return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)
157def setup_axes(axes_amplitude=None, axes_phase=None):
158 '''
159 Configure axes in Bode plot style.
160 '''
162 if axes_amplitude is not None:
163 axes_amplitude.set_ylabel('Amplitude ratio')
164 axes_amplitude.set_xscale('log')
165 axes_amplitude.set_yscale('log')
166 axes_amplitude.grid(True)
167 axes_amplitude.axhline(1.0, lw=0.5, color='black')
168 if axes_phase is None:
169 axes_amplitude.set_xlabel('Frequency [Hz]')
170 axes_amplitude.set_xscale('log')
171 else:
172 axes_amplitude.xaxis.set_ticklabels([])
174 if axes_phase is not None:
175 axes_phase.set_ylabel('Phase [$\\pi$]')
176 axes_phase.set_xscale('log')
177 axes_phase.set_xlabel('Frequency [Hz]')
178 axes_phase.grid(True)
179 axes_phase.axhline(0.0, lw=0.5, color='black')
182def plot(
183 responses,
184 filename=None,
185 dpi=100,
186 fmin=0.01, fmax=100., nf=100,
187 normalize=False,
188 fontsize=10.,
189 figsize=None,
190 styles=None,
191 labels=None,
192 show_breakpoints=False):
194 '''
195 Draw instrument responses in Bode plot style.
197 :param responses: instrument responses as
198 :py:class:`pyrocko.response.FrequencyResponse` objects
199 :param fmin: minimum frequency [Hz]
200 :param fmax: maximum frequency [Hz]
201 :param nf: number of frequencies where to evaluate the response
202 :param normalize: if ``True`` normalize flat part of response to be ``1``
203 :param styles: :py:class:`list` of :py:class:`dict` objects with keyword
204 arguments to be passed to matplotlib's
205 :py:meth:`matplotlib.axes.Axes.plot` function when drawing the response
206 lines. Length must match number of responses.
207 :param filename: file name to pass to matplotlib's ``savefig`` function. If
208 ``None``, the plot is shown with :py:func:`matplotlib.pyplot.show`.
209 :param fontsize: font size in points used in axis labels and legend
210 :param figsize: :py:class:`tuple`, ``(width, height)`` in inches
211 :param labels: :py:class:`list` of names to show in legend. Length must
212 correspond to number of responses.
213 :param show_breakpoints: show breakpoints of pole-zero responses.
214 '''
216 from matplotlib import pyplot as plt
217 from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
218 from pyrocko.plot import graph_colors, to01
220 mpl_init(fontsize=fontsize)
222 if figsize is None:
223 figsize = mpl_papersize('a4', 'portrait')
225 fig = plt.figure(figsize=figsize)
226 labelpos = mpl_margins(
227 fig, w=7., h=5., units=fontsize, nw=1, nh=2, hspace=2.)
228 axes_amplitude = fig.add_subplot(2, 1, 1)
229 labelpos(axes_amplitude, 2., 1.5)
230 axes_phase = fig.add_subplot(2, 1, 2)
231 labelpos(axes_phase, 2., 1.5)
233 setup_axes(axes_amplitude, axes_phase)
235 if styles is None:
236 styles = [
237 dict(color=to01(graph_colors[i % len(graph_colors)]))
238 for i in range(len(responses))]
239 else:
240 assert len(styles) == len(responses)
242 if labels is None:
243 labels = [None] * len(responses)
244 else:
245 assert len(labels) == len(responses)
247 a_ranges, p_ranges = [], []
248 have_labels = False
249 for style, resp, label in zip(styles, responses, labels):
250 a_range, p_range = draw(
251 response=resp,
252 axes_amplitude=axes_amplitude,
253 axes_phase=axes_phase,
254 fmin=fmin, fmax=fmax, nf=nf,
255 normalize=normalize,
256 style=style,
257 label=label,
258 show_breakpoints=show_breakpoints)
260 if label is not None:
261 have_labels = True
263 a_ranges.append(a_range)
264 p_ranges.append(p_range)
266 if have_labels:
267 axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))
269 if a_ranges:
270 a_ranges = num.array(a_ranges)
271 p_ranges = num.array(p_ranges)
273 amin, amax = num.nanmin(a_ranges), num.nanmax(a_ranges)
274 pmin, pmax = num.nanmin(p_ranges), num.nanmax(p_ranges)
276 amin *= 0.5
277 amax *= 2.0
279 pmin -= 0.5
280 pmax += 0.5
282 pmin = math.floor(pmin)
283 pmax = math.ceil(pmax)
285 axes_amplitude.set_ylim(amin, amax)
286 axes_phase.set_ylim(pmin, pmax)
287 axes_amplitude.set_xlim(fmin, fmax)
288 axes_phase.set_xlim(fmin, fmax)
290 if filename is not None:
291 fig.savefig(filename, dpi=dpi)
292 else:
293 plt.show()
296def tts(t):
297 if t is None:
298 return '?'
299 else:
300 return util.tts(t, format='%Y-%m-%d')
303def load_response_information(
304 filename, format,
305 nslc_patterns=None, fake_input_units=None, stages=None):
307 from pyrocko import pz, trace
308 from pyrocko.io import resp as fresp
310 resps = []
311 labels = []
312 if format == 'sacpz':
313 if fake_input_units is not None:
314 raise Exception(
315 'cannot guess true input units from plain SAC PZ files')
317 zeros, poles, constant = pz.read_sac_zpk(filename)
318 resp = trace.PoleZeroResponse(
319 zeros=zeros, poles=poles, constant=constant)
321 resps.append(resp)
322 labels.append(filename)
324 elif format == 'pf':
325 if fake_input_units is not None:
326 raise Exception(
327 'Cannot guess input units from plain response files.')
329 resp = guts.load(filename=filename)
330 resps.append(resp)
331 labels.append(filename)
333 elif format in ('resp', 'evalresp'):
334 for resp in list(fresp.iload_filename(filename)):
335 if nslc_patterns is not None and not util.match_nslc(
336 nslc_patterns, resp.codes):
337 continue
339 units = ''
340 in_units = None
341 if resp.response.instrument_sensitivity:
342 s = resp.response.instrument_sensitivity
343 in_units = fake_input_units or s.input_units.name
344 if s.input_units and s.output_units:
345 units = ', %s -> %s' % (in_units, s.output_units.name)
347 if format == 'resp':
348 resps.append(resp.response.get_pyrocko_response(
349 '.'.join(resp.codes),
350 fake_input_units=fake_input_units,
351 stages=stages).expect_one())
352 else:
353 target = {
354 'M/S': 'vel',
355 'M': 'dis',
356 }[in_units]
358 if resp.end_date is not None:
359 time = (resp.start_date + resp.end_date)*0.5
360 else:
361 time = resp.start_date + 3600*24*10
363 r = trace.Evalresp(
364 respfile=filename,
365 nslc_id=resp.codes,
366 target=target,
367 time=time,
368 stages=stages)
370 resps.append(r)
372 labels.append('%s (%s.%s.%s.%s, %s - %s%s)' % (
373 (filename, ) + resp.codes +
374 (tts(resp.start_date), tts(resp.end_date), units)))
376 elif format == 'stationxml':
377 from pyrocko.fdsn import station as fs
379 sx = fs.load_xml(filename=filename)
380 for network in sx.network_list:
381 for station in network.station_list:
382 for channel in station.channel_list:
383 nslc = (
384 network.code,
385 station.code,
386 channel.location_code,
387 channel.code)
389 if nslc_patterns is not None and not util.match_nslc(
390 nslc_patterns, nslc):
391 continue
393 if not channel.response:
394 logger.warning(
395 'no response for channel %s.%s.%s.%s given.'
396 % nslc)
397 continue
399 units = ''
400 if channel.response.instrument_sensitivity:
401 s = channel.response.instrument_sensitivity
402 if s.input_units and s.output_units:
403 units = ', %s -> %s' % (
404 fake_input_units or s.input_units.name,
405 s.output_units.name)
407 resps.append(channel.response.get_pyrocko_response(
408 '.'.join(nslc),
409 fake_input_units=fake_input_units,
410 stages=stages).expect_one())
412 labels.append(
413 '%s (%s.%s.%s.%s, %s - %s%s)' % (
414 (filename, ) + nslc +
415 (tts(channel.start_date),
416 tts(channel.end_date),
417 units)))
419 return resps, labels
422if __name__ == '__main__':
423 import sys
424 from optparse import OptionParser
426 util.setup_logging('pyrocko.plot.response.__main__', 'warning')
428 usage = 'python -m pyrocko.plot.response <filename> ... [options]'
430 description = '''Plot instrument responses (transfer functions).'''
432 allowed_formats = ['sacpz', 'resp', 'evalresp', 'stationxml', 'pf']
434 parser = OptionParser(
435 usage=usage,
436 description=description,
437 formatter=util.BetterHelpFormatter())
439 parser.add_option(
440 '--format',
441 dest='format',
442 default='sacpz',
443 choices=allowed_formats,
444 help='assume input files are of given FORMAT. Choices: %s' % (
445 ', '.join(allowed_formats)))
447 parser.add_option(
448 '--fmin',
449 dest='fmin',
450 type='float',
451 default=0.01,
452 help='minimum frequency [Hz], default: %default')
454 parser.add_option(
455 '--fmax',
456 dest='fmax',
457 type='float',
458 default=100.,
459 help='maximum frequency [Hz], default: %default')
461 parser.add_option(
462 '--normalize',
463 dest='normalize',
464 action='store_true',
465 help='normalize response to be 1 on flat part')
467 parser.add_option(
468 '--save',
469 dest='filename',
470 help='save figure to file with name FILENAME')
472 parser.add_option(
473 '--dpi',
474 dest='dpi',
475 type='float',
476 default=100.,
477 help='DPI setting for pixel image output, default: %default')
479 parser.add_option(
480 '--patterns',
481 dest='nslc_patterns',
482 metavar='NET.STA.LOC.CHA,...',
483 help='show only responses of channels matching any of the given '
484 'patterns')
486 parser.add_option(
487 '--stages',
488 dest='stages',
489 metavar='START:STOP',
490 help='show only responses selected stages')
492 parser.add_option(
493 '--fake-input-units',
494 dest='fake_input_units',
495 choices=['M', 'M/S', 'M/S**2'],
496 metavar='UNITS',
497 help='show converted response for given input units, choices: '
498 '["M", "M/S", "M/S**2"]')
500 parser.add_option(
501 '--show-breakpoints',
502 dest='show_breakpoints',
503 action='store_true',
504 default=False,
505 help='show breakpoints of pole-zero responses')
507 (options, args) = parser.parse_args(sys.argv[1:])
509 if len(args) == 0:
510 parser.print_help()
511 sys.exit(1)
513 fns = args
515 resps = []
516 labels = []
518 stages = None
519 if options.stages:
520 stages = tuple(int(x) for x in options.stages.split(':', 1))
521 stages = stages[0]-1, stages[1]
523 for fn in fns:
525 if options.nslc_patterns is not None:
526 nslc_patterns = options.nslc_patterns.split(',')
527 else:
528 nslc_patterns = None
530 resps_this, labels_this = load_response_information(
531 fn, options.format, nslc_patterns,
532 fake_input_units=options.fake_input_units,
533 stages=stages)
535 resps.extend(resps_this)
536 labels.extend(labels_this)
538 plot(
539 resps,
540 fmin=options.fmin, fmax=options.fmax, nf=200,
541 normalize=options.normalize,
542 labels=labels, filename=options.filename, dpi=options.dpi,
543 show_breakpoints=options.show_breakpoints)