|
# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg---------- ''' This module contains functions to plot instrument response transfer functions in Bode plot style using Matplotlib.
Example
* :download:`test_response_plot.py </../../examples/test_response_plot.py>`
::
from pyrocko.plot import response from pyrocko.example import get_example_data
get_example_data('test_response.resp')
resps, labels = response.load_response_information( 'test_response.resp', 'resp')
response.plot( responses=resps, labels=labels, filename='test_response.png', fmin=0.001, fmax=400., dpi=75.)
.. figure :: /static/test_response.png :align: center
Example response plot ''' from __future__ import absolute_import from builtins import range
import logging
import numpy as num
from pyrocko import util from pyrocko import guts
logger = logging.getLogger('plot.response')
def normalize_on_flat(f, tf, factor=10000.): df = num.diff(num.log(f)) tap = 1.0 / (1.0 + factor * (num.diff(num.log(num.abs(tf)))/df)**2) return tf / (num.sum(num.abs(tf[1:]) * tap) / num.sum(tap))
def draw( response, axes_amplitude=None, axes_phase=None, fmin=0.01, fmax=100., nf=100, normalize=False, style={}, label=None):
'''Draw instrument response in Bode plot style to given Matplotlib axes
:param response: instrument response as a :py:class:`pyrocko.trace.FrequencyResponse` object :param axes_amplitude: :py:class:`matplotlib.axes.Axes` object to use when drawing the amplitude response :param axes_phase: :py:class:`matplotlib.axes.Axes` object to use when drawing the phase response :param fmin: minimum frequency [Hz] :param fmax: maximum frequency [Hz] :param nf: number of frequencies where to evaluate the response :param style: :py:class:`dict` with keyword arguments to tune the line style :param label: string to be passed to the ``label`` argument of :py:meth:`matplotlib.axes.Axes.plot` '''
f = num.exp(num.linspace(num.log(fmin), num.log(fmax), nf)) tf = response.evaluate(f)
if normalize: tf = normalize_on_flat(f, tf)
ta = num.abs(tf)
if axes_amplitude: axes_amplitude.plot(f, ta, label=label, **style)
if axes_phase: dta = num.diff(num.log(ta)) iflat = num.argmin(num.abs(num.diff(dta)) + num.abs(dta[:-1])) tp = num.unwrap(num.angle(tf)) ioff = int(num.round(tp[iflat] / (2.0*num.pi))) tp -= ioff * 2.0 * num.pi axes_phase.plot(f, tp/num.pi, label=label, **style) else: tp = [0.]
return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)
def setup_axes(axes_amplitude=None, axes_phase=None): '''Configure axes in Bode plot style.'''
if axes_amplitude is not None: axes_amplitude.set_ylabel('Amplitude ratio') axes_amplitude.set_xscale('log') axes_amplitude.set_yscale('log') axes_amplitude.grid(True) axes_amplitude.axhline(1.0, lw=0.5, color='black') if axes_phase is None: axes_amplitude.set_xlabel('Frequency [Hz]') axes_amplitude.set_xscale('log') else: axes_amplitude.xaxis.set_ticklabels([])
if axes_phase is not None: axes_phase.set_ylabel('Phase [$\\pi$]') axes_phase.set_xscale('log') axes_phase.set_xlabel('Frequency [Hz]') axes_phase.grid(True) axes_phase.axhline(0.0, lw=0.5, color='black')
def plot( responses, filename=None, dpi=100, fmin=0.01, fmax=100., nf=100, normalize=False, fontsize=10., figsize=None, styles=None, labels=None):
''' Draw instrument responses in Bode plot style.
:param responses: instrument responses as :py:class:`pyrocko.trace.FrequencyResponse` objects :param fmin: minimum frequency [Hz] :param fmax: maximum frequency [Hz] :param nf: number of frequencies where to evaluate the response :param normalize: if ``True`` normalize flat part of response to be ``1`` :param styles: :py:class:`list` of :py:class:`dict` objects with keyword arguments to be passed to matplotlib's :py:meth:`matplotlib.axes.Axes.plot` function when drawing the response lines. Length must match number of responses. :param filename: file name to pass to matplotlib's ``savefig`` function. If ``None``, the plot is shown with :py:func:`matplotlib.pyplot.show`. :param fontsize: font size in points used in axis labels and legend :param figsize: :py:class:`tuple`, ``(width, height)`` in inches :param labels: :py:class:`list` of names to show in legend. Length must correspond to number of responses. '''
from matplotlib import pyplot as plt from pyrocko.plot import mpl_init, mpl_margins, mpl_papersize from pyrocko.plot import graph_colors, to01
mpl_init(fontsize=fontsize)
if figsize is None: figsize = mpl_papersize('a4', 'portrait')
fig = plt.figure(figsize=figsize) labelpos = mpl_margins( fig, w=7., h=5., units=fontsize, nw=1, nh=2, hspace=2.) axes_amplitude = fig.add_subplot(2, 1, 1) labelpos(axes_amplitude, 2., 1.5) axes_phase = fig.add_subplot(2, 1, 2) labelpos(axes_phase, 2., 1.5)
setup_axes(axes_amplitude, axes_phase)
if styles is None: styles = [ dict(color=to01(graph_colors[i % len(graph_colors)])) for i in range(len(responses))] else: assert len(styles) == len(responses)
if labels is None: labels = [None] * len(responses) else: assert len(labels) == len(responses)
a_ranges, p_ranges = [], [] have_labels = False for style, resp, label in zip(styles, responses, labels): a_range, p_range = draw( response=resp, axes_amplitude=axes_amplitude, axes_phase=axes_phase, fmin=fmin, fmax=fmax, nf=nf, normalize=normalize, style=style, label=label)
if label is not None: have_labels = True
a_ranges.append(a_range) p_ranges.append(p_range)
if have_labels: axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))
a_ranges = num.array(a_ranges) p_ranges = num.array(p_ranges)
amin, amax = num.min(a_ranges), num.max(a_ranges) pmin, pmax = num.min(p_ranges), num.max(p_ranges)
amin *= 0.5 amax *= 2.0
pmin -= 0.5 pmax += 0.5
axes_amplitude.set_ylim(amin, amax) axes_phase.set_ylim(pmin, pmax) axes_amplitude.set_xlim(fmin, fmax) axes_phase.set_xlim(fmin, fmax)
if filename is not None: fig.savefig(filename, dpi=dpi) else: plt.show()
def tts(t): if t is None: return '?' else: return util.tts(t, format='%Y-%m-%d')
def load_response_information( filename, format, nslc_patterns=None, fake_input_units=None):
from pyrocko import pz, trace from pyrocko.io import resp as fresp
resps = [] labels = [] if format == 'sacpz': if fake_input_units is not None: raise Exception( 'cannot guess true input units from plain SAC PZ files')
zeros, poles, constant = pz.read_sac_zpk(filename) resp = trace.PoleZeroResponse( zeros=zeros, poles=poles, constant=constant)
resps.append(resp) labels.append(filename)
elif format == 'pf': if fake_input_units is not None: raise Exception( 'cannot guess true input units from plain response files')
resp = guts.load(filename=filename) resps.append(resp) labels.append(filename)
elif format == 'resp': for resp in list(fresp.iload_filename(filename)): if nslc_patterns is not None and not util.match_nslc( nslc_patterns, resp.codes): continue
units = '' if resp.response.instrument_sensitivity: s = resp.response.instrument_sensitivity if s.input_units and s.output_units: units = ', %s -> %s' % ( fake_input_units or s.input_units.name, s.output_units.name)
resps.append(resp.response.get_pyrocko_response( resp.codes, fake_input_units=fake_input_units))
labels.append('%s (%s.%s.%s.%s, %s - %s%s)' % ( (filename, ) + resp.codes + (tts(resp.start_date), tts(resp.end_date), units)))
elif format == 'stationxml': from pyrocko.fdsn import station as fs
sx = fs.load_xml(filename=filename) for network in sx.network_list: for station in network.station_list: for channel in station.channel_list: nslc = ( network.code, station.code, channel.location_code, channel.code)
if nslc_patterns is not None and not util.match_nslc( nslc_patterns, nslc): continue
if not channel.response: logger.warn( 'no response for channel %s.%s.%s.%s given.' % nslc) continue
units = '' if channel.response.instrument_sensitivity: s = channel.response.instrument_sensitivity if s.input_units and s.output_units: units = ', %s -> %s' % ( fake_input_units or s.input_units.name, s.output_units.name)
resps.append(channel.response.get_pyrocko_response( nslc, fake_input_units=fake_input_units))
labels.append( '%s (%s.%s.%s.%s, %s - %s%s)' % ( (filename, ) + nslc + (tts(channel.start_date), tts(channel.end_date), units)))
return resps, labels
if __name__ == '__main__': import sys from optparse import OptionParser
util.setup_logging('pyrocko.plot.response.__main__', 'warning')
usage = 'python -m pyrocko.plot.response <filename> ... [options]'
description = '''Plot instrument responses (transfer functions).'''
allowed_formats = ['sacpz', 'resp', 'stationxml', 'pf']
parser = OptionParser( usage=usage, description=description, formatter=util.BetterHelpFormatter())
parser.add_option( '--format', dest='format', default='sacpz', choices=allowed_formats, help='assume input files are of given FORMAT. Choices: %s' % ( ', '.join(allowed_formats)))
parser.add_option( '--fmin', dest='fmin', type='float', default=0.01, help='minimum frequency [Hz], default: %default')
parser.add_option( '--fmax', dest='fmax', type='float', default=100., help='maximum frequency [Hz], default: %default')
parser.add_option( '--normalize', dest='normalize', action='store_true', help='normalize response to be 1 on flat part')
parser.add_option( '--save', dest='filename', help='save figure to file with name FILENAME')
parser.add_option( '--dpi', dest='dpi', type='float', default=100., help='DPI setting for pixel image output, default: %default')
parser.add_option( '--patterns', dest='nslc_patterns', metavar='NET.STA.LOC.CHA,...', help='show only responses of channels matching any of the given ' 'patterns')
parser.add_option( '--fake-input-units', dest='fake_input_units', choices=['M', 'M/S', 'M/S**2'], metavar='UNITS', help='show converted response for given input units, choices: ' '["M", "M/S", "M/S**2"]')
(options, args) = parser.parse_args(sys.argv[1:])
if len(args) == 0: parser.print_help() sys.exit(1)
fns = args
resps = [] labels = []
for fn in fns:
if options.nslc_patterns is not None: nslc_patterns = options.nslc_patterns.split(',') else: nslc_patterns = None
resps_this, labels_this = load_response_information( fn, options.format, nslc_patterns, fake_input_units=options.fake_input_units)
resps.extend(resps_this) labels.extend(labels_this)
plot( resps, fmin=options.fmin, fmax=options.fmax, nf=200, normalize=options.normalize, labels=labels, filename=options.filename, dpi=options.dpi) |