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):
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 if axes_phase:
115 dta = num.diff(num.log(ta))
116 iflat = num.nanargmin(num.abs(num.diff(dta)) + num.abs(dta[:-1]))
117 tp = num.unwrap(num.angle(tf))
118 ioff = int(num.round(tp[iflat] / (2.0*num.pi)))
119 tp -= ioff * 2.0 * num.pi
120 axes_phase.plot(f, tp/num.pi, label=label, **style)
121 else:
122 tp = [0.]
124 return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)
127def setup_axes(axes_amplitude=None, axes_phase=None):
128 '''
129 Configure axes in Bode plot style.
130 '''
132 if axes_amplitude is not None:
133 axes_amplitude.set_ylabel('Amplitude ratio')
134 axes_amplitude.set_xscale('log')
135 axes_amplitude.set_yscale('log')
136 axes_amplitude.grid(True)
137 axes_amplitude.axhline(1.0, lw=0.5, color='black')
138 if axes_phase is None:
139 axes_amplitude.set_xlabel('Frequency [Hz]')
140 axes_amplitude.set_xscale('log')
141 else:
142 axes_amplitude.xaxis.set_ticklabels([])
144 if axes_phase is not None:
145 axes_phase.set_ylabel('Phase [$\\pi$]')
146 axes_phase.set_xscale('log')
147 axes_phase.set_xlabel('Frequency [Hz]')
148 axes_phase.grid(True)
149 axes_phase.axhline(0.0, lw=0.5, color='black')
152def plot(
153 responses,
154 filename=None,
155 dpi=100,
156 fmin=0.01, fmax=100., nf=100,
157 normalize=False,
158 fontsize=10.,
159 figsize=None,
160 styles=None,
161 labels=None):
163 '''
164 Draw instrument responses in Bode plot style.
166 :param responses: instrument responses as
167 :py:class:`pyrocko.response.FrequencyResponse` objects
168 :param fmin: minimum frequency [Hz]
169 :param fmax: maximum frequency [Hz]
170 :param nf: number of frequencies where to evaluate the response
171 :param normalize: if ``True`` normalize flat part of response to be ``1``
172 :param styles: :py:class:`list` of :py:class:`dict` objects with keyword
173 arguments to be passed to matplotlib's
174 :py:meth:`matplotlib.axes.Axes.plot` function when drawing the response
175 lines. Length must match number of responses.
176 :param filename: file name to pass to matplotlib's ``savefig`` function. If
177 ``None``, the plot is shown with :py:func:`matplotlib.pyplot.show`.
178 :param fontsize: font size in points used in axis labels and legend
179 :param figsize: :py:class:`tuple`, ``(width, height)`` in inches
180 :param labels: :py:class:`list` of names to show in legend. Length must
181 correspond to number of responses.
182 '''
184 from matplotlib import pyplot as plt
185 from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize
186 from pyrocko.plot import graph_colors, to01
188 mpl_init(fontsize=fontsize)
190 if figsize is None:
191 figsize = mpl_papersize('a4', 'portrait')
193 fig = plt.figure(figsize=figsize)
194 labelpos = mpl_margins(
195 fig, w=7., h=5., units=fontsize, nw=1, nh=2, hspace=2.)
196 axes_amplitude = fig.add_subplot(2, 1, 1)
197 labelpos(axes_amplitude, 2., 1.5)
198 axes_phase = fig.add_subplot(2, 1, 2)
199 labelpos(axes_phase, 2., 1.5)
201 setup_axes(axes_amplitude, axes_phase)
203 if styles is None:
204 styles = [
205 dict(color=to01(graph_colors[i % len(graph_colors)]))
206 for i in range(len(responses))]
207 else:
208 assert len(styles) == len(responses)
210 if labels is None:
211 labels = [None] * len(responses)
212 else:
213 assert len(labels) == len(responses)
215 a_ranges, p_ranges = [], []
216 have_labels = False
217 for style, resp, label in zip(styles, responses, labels):
218 a_range, p_range = draw(
219 response=resp,
220 axes_amplitude=axes_amplitude,
221 axes_phase=axes_phase,
222 fmin=fmin, fmax=fmax, nf=nf,
223 normalize=normalize,
224 style=style,
225 label=label)
227 if label is not None:
228 have_labels = True
230 a_ranges.append(a_range)
231 p_ranges.append(p_range)
233 if have_labels:
234 axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))
236 if a_ranges:
237 a_ranges = num.array(a_ranges)
238 p_ranges = num.array(p_ranges)
240 amin, amax = num.nanmin(a_ranges), num.nanmax(a_ranges)
241 pmin, pmax = num.nanmin(p_ranges), num.nanmax(p_ranges)
243 amin *= 0.5
244 amax *= 2.0
246 pmin -= 0.5
247 pmax += 0.5
249 pmin = math.floor(pmin)
250 pmax = math.ceil(pmax)
252 axes_amplitude.set_ylim(amin, amax)
253 axes_phase.set_ylim(pmin, pmax)
254 axes_amplitude.set_xlim(fmin, fmax)
255 axes_phase.set_xlim(fmin, fmax)
257 if filename is not None:
258 fig.savefig(filename, dpi=dpi)
259 else:
260 plt.show()
263def tts(t):
264 if t is None:
265 return '?'
266 else:
267 return util.tts(t, format='%Y-%m-%d')
270def load_response_information(
271 filename, format,
272 nslc_patterns=None, fake_input_units=None, stages=None):
274 from pyrocko import pz, trace
275 from pyrocko.io import resp as fresp
277 resps = []
278 labels = []
279 if format == 'sacpz':
280 if fake_input_units is not None:
281 raise Exception(
282 'cannot guess true input units from plain SAC PZ files')
284 zeros, poles, constant = pz.read_sac_zpk(filename)
285 resp = trace.PoleZeroResponse(
286 zeros=zeros, poles=poles, constant=constant)
288 resps.append(resp)
289 labels.append(filename)
291 elif format == 'pf':
292 if fake_input_units is not None:
293 raise Exception(
294 'Cannot guess input units from plain response files.')
296 resp = guts.load(filename=filename)
297 resps.append(resp)
298 labels.append(filename)
300 elif format in ('resp', 'evalresp'):
301 for resp in list(fresp.iload_filename(filename)):
302 if nslc_patterns is not None and not util.match_nslc(
303 nslc_patterns, resp.codes):
304 continue
306 units = ''
307 in_units = None
308 if resp.response.instrument_sensitivity:
309 s = resp.response.instrument_sensitivity
310 in_units = fake_input_units or s.input_units.name
311 if s.input_units and s.output_units:
312 units = ', %s -> %s' % (in_units, s.output_units.name)
314 if format == 'resp':
315 resps.append(resp.response.get_pyrocko_response(
316 '.'.join(resp.codes),
317 fake_input_units=fake_input_units,
318 stages=stages).expect_one())
319 else:
320 target = {
321 'M/S': 'vel',
322 'M': 'dis',
323 }[in_units]
325 if resp.end_date is not None:
326 time = (resp.start_date + resp.end_date)*0.5
327 else:
328 time = resp.start_date + 3600*24*10
330 r = trace.Evalresp(
331 respfile=filename,
332 nslc_id=resp.codes,
333 target=target,
334 time=time,
335 stages=stages)
337 resps.append(r)
339 labels.append('%s (%s.%s.%s.%s, %s - %s%s)' % (
340 (filename, ) + resp.codes +
341 (tts(resp.start_date), tts(resp.end_date), units)))
343 elif format == 'stationxml':
344 from pyrocko.fdsn import station as fs
346 sx = fs.load_xml(filename=filename)
347 for network in sx.network_list:
348 for station in network.station_list:
349 for channel in station.channel_list:
350 nslc = (
351 network.code,
352 station.code,
353 channel.location_code,
354 channel.code)
356 if nslc_patterns is not None and not util.match_nslc(
357 nslc_patterns, nslc):
358 continue
360 if not channel.response:
361 logger.warning(
362 'no response for channel %s.%s.%s.%s given.'
363 % nslc)
364 continue
366 units = ''
367 if channel.response.instrument_sensitivity:
368 s = channel.response.instrument_sensitivity
369 if s.input_units and s.output_units:
370 units = ', %s -> %s' % (
371 fake_input_units or s.input_units.name,
372 s.output_units.name)
374 resps.append(channel.response.get_pyrocko_response(
375 '.'.join(nslc),
376 fake_input_units=fake_input_units,
377 stages=stages).expect_one())
379 labels.append(
380 '%s (%s.%s.%s.%s, %s - %s%s)' % (
381 (filename, ) + nslc +
382 (tts(channel.start_date),
383 tts(channel.end_date),
384 units)))
386 return resps, labels
389if __name__ == '__main__':
390 import sys
391 from optparse import OptionParser
393 util.setup_logging('pyrocko.plot.response.__main__', 'warning')
395 usage = 'python -m pyrocko.plot.response <filename> ... [options]'
397 description = '''Plot instrument responses (transfer functions).'''
399 allowed_formats = ['sacpz', 'resp', 'evalresp', 'stationxml', 'pf']
401 parser = OptionParser(
402 usage=usage,
403 description=description,
404 formatter=util.BetterHelpFormatter())
406 parser.add_option(
407 '--format',
408 dest='format',
409 default='sacpz',
410 choices=allowed_formats,
411 help='assume input files are of given FORMAT. Choices: %s' % (
412 ', '.join(allowed_formats)))
414 parser.add_option(
415 '--fmin',
416 dest='fmin',
417 type='float',
418 default=0.01,
419 help='minimum frequency [Hz], default: %default')
421 parser.add_option(
422 '--fmax',
423 dest='fmax',
424 type='float',
425 default=100.,
426 help='maximum frequency [Hz], default: %default')
428 parser.add_option(
429 '--normalize',
430 dest='normalize',
431 action='store_true',
432 help='normalize response to be 1 on flat part')
434 parser.add_option(
435 '--save',
436 dest='filename',
437 help='save figure to file with name FILENAME')
439 parser.add_option(
440 '--dpi',
441 dest='dpi',
442 type='float',
443 default=100.,
444 help='DPI setting for pixel image output, default: %default')
446 parser.add_option(
447 '--patterns',
448 dest='nslc_patterns',
449 metavar='NET.STA.LOC.CHA,...',
450 help='show only responses of channels matching any of the given '
451 'patterns')
453 parser.add_option(
454 '--stages',
455 dest='stages',
456 metavar='START:STOP',
457 help='show only responses selected stages')
459 parser.add_option(
460 '--fake-input-units',
461 dest='fake_input_units',
462 choices=['M', 'M/S', 'M/S**2'],
463 metavar='UNITS',
464 help='show converted response for given input units, choices: '
465 '["M", "M/S", "M/S**2"]')
467 (options, args) = parser.parse_args(sys.argv[1:])
469 if len(args) == 0:
470 parser.print_help()
471 sys.exit(1)
473 fns = args
475 resps = []
476 labels = []
478 stages = None
479 if options.stages:
480 stages = tuple(int(x) for x in options.stages.split(':', 1))
481 stages = stages[0]-1, stages[1]
483 for fn in fns:
485 if options.nslc_patterns is not None:
486 nslc_patterns = options.nslc_patterns.split(',')
487 else:
488 nslc_patterns = None
490 resps_this, labels_this = load_response_information(
491 fn, options.format, nslc_patterns,
492 fake_input_units=options.fake_input_units,
493 stages=stages)
495 resps.extend(resps_this)
496 labels.extend(labels_this)
498 plot(
499 resps,
500 fmin=options.fmin, fmax=options.fmax, nf=200,
501 normalize=options.normalize,
502 labels=labels, filename=options.filename, dpi=options.dpi)