Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/response.py: 72%

1# http://pyrocko.org - GPLv3

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]

106

107 if normalize:

108 tf = normalize_on_flat(f, tf)

109

110 ta = num.abs(tf)

111

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

124

125 style = dict(style)

126 style['color'] = color_pool[0][c_key]

127

128 ikey = list(color_pool[0].keys()).index(c_key) + 1

129 slabel = '[%i]' % ikey

130

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)

136

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'))

143

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'))

152

153 if show_breakpoints:

154 for br_frequency, br_change in response.construction():

155 if not (fmin <= br_frequency <= fmax):

156 continue

157

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)

164

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'))

173

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.]

183

184 return (num.min(ta), num.max(ta)), (num.min(tp)/num.pi, num.max(tp)/num.pi)

185

186

187def setup_axes(axes_amplitude=None, axes_phase=None):

188 '''

189 Configure axes in Bode plot style.

190 '''

191

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([])

203

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')

210

211

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):

224

225 '''

226 Draw instrument responses in Bode plot style.

227

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 '''

246

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

250

251 mpl_init(fontsize=fontsize)

252

253 if figsize is None:

254 figsize = mpl_papersize('a4', 'portrait')

255

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)

263

264 setup_axes(axes_amplitude, axes_phase)

265

266 if styles is None:

267 styles = [{} for i in range(len(responses))]

268 color_pool = ({}, [to01(color) for color in graph_colors])

269

270 else:

271 color_pool = None

272 assert len(styles) == len(responses)

273

274 if separate_combined_labels is not None:

275 label_pool = []

276 else:

277 label_pool = None

278

279 if labels is None:

280 labels = [None] * len(responses)

281 else:

282 assert len(labels) == len(responses)

283

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)

299

300 if label is not None:

301 have_labels = True

302

303 a_ranges.append(a_range)

304 p_ranges.append(p_range)

305

306 except InvalidResponseError as e:

307 logger.error(

308 'Error occured for response labeled "%s": %s' % (label, e))

309

310 if have_labels:

311 axes_amplitude.legend(loc='lower right', prop=dict(size=fontsize))

312

313 if separate_combined_labels == 'print':

314 for label in label_pool:

315 print(label)

316

317 if a_ranges:

318 a_ranges = num.array(a_ranges)

319 p_ranges = num.array(p_ranges)

320

321 amin, amax = num.nanmin(a_ranges), num.nanmax(a_ranges)

322 pmin, pmax = num.nanmin(p_ranges), num.nanmax(p_ranges)

323

324 amin *= 0.5

325 amax *= 2.0

326

327 pmin -= 0.5

328 pmax += 0.5

329

330 pmin = math.floor(pmin)

331 pmax = math.ceil(pmax)

332

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)

337

338 if filename is not None:

339 fig.savefig(filename, dpi=dpi)

340 else:

341 plt.show()

342

343 if separate_combined_labels == 'return':

344 return label_pool

345

346

347def tts(t):

348 if t is None:

349 return '?'

350 else:

351 return util.tts(t, format='%Y-%m-%d')

352

353

354def load_response_information(

355 filename, format,

356 nslc_patterns=None,

357 fake_input_units=None,

358 stages=None):

359

360 from pyrocko import pz, trace

361 from pyrocko.io import resp as fresp

362

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')

369

370 zeros, poles, constant = pz.read_sac_zpk(filename)

371 resp = trace.PoleZeroResponse(

372 zeros=zeros, poles=poles, constant=constant)

373

374 resps.append(resp)

375 labels.append(filename)

376

377 elif format == 'pf':

378 if fake_input_units is not None:

379 raise Exception(

380 'Cannot guess input units from plain response files.')

381

382 resp = guts.load(filename=filename)

383 resps.append(resp)

384 labels.append(filename)

385

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

391

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)

399

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]

410

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

415

416 r = trace.Evalresp(

417 respfile=filename,

418 nslc_id=resp.codes,

419 target=target,

420 time=time,

421 stages=stages)

422

423 resps.append(r)

424

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)))

428

429 elif format == 'stationxml':

430 from pyrocko.io import stationxml as fs

431

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)

441

442 if nslc_patterns is not None and not util.match_nslc(

443 nslc_patterns, nslc):

444 continue

445

446 if not channel.response:

447 logger.warning(

448 'no response for channel %s.%s.%s.%s given.'

449 % nslc)

450 continue

451

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)

459

460 resps.append(channel.response.get_pyrocko_response(

461 '.'.join(nslc),

462 fake_input_units=fake_input_units,

463 stages=stages).expect_one())

464

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)))

471

472 return resps, labels

473

474

475if __name__ == '__main__':

476 import sys

477 from optparse import OptionParser

478

479 util.setup_logging('pyrocko.plot.response.__main__', 'warning')

480

481 usage = 'python -m pyrocko.plot.response <filename> ... [options]'

482

483 description = '''Plot instrument responses (transfer functions).'''

484

485 allowed_formats = ['sacpz', 'resp', 'evalresp', 'stationxml', 'pf']

486

487 parser = OptionParser(

488 usage=usage,

489 description=description,

490 formatter=util.BetterHelpFormatter())

491

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)))

499

500 parser.add_option(

501 '--fmin',

502 dest='fmin',

503 type='float',

504 default=0.01,

505 help='minimum frequency [Hz], default: %default')

506

507 parser.add_option(

508 '--fmax',

509 dest='fmax',

510 type='float',

511 default=100.,

512 help='maximum frequency [Hz], default: %default')

513

514 parser.add_option(

515 '--normalize',

516 dest='normalize',

517 action='store_true',

518 help='normalize response to be 1 on flat part')

519

520 parser.add_option(

521 '--save',

522 dest='filename',

523 help='save figure to file with name FILENAME')

524

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')

531

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')

538

539 parser.add_option(

540 '--stages',

541 dest='stages',

542 metavar='START:STOP',

543 help='show only responses selected stages')

544

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"]')

552

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')

559

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.')

568

569 (options, args) = parser.parse_args(sys.argv[1:])

570

571 if len(args) == 0:

572 parser.print_help()

573 sys.exit(1)

574

575 fns = args

576

577 resps = []

578 labels = []

579

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]

584

585 for fn in fns:

586

587 if options.nslc_patterns is not None:

588 nslc_patterns = options.nslc_patterns.split(',')

589 else:

590 nslc_patterns = None

591

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)

596

597 resps.extend(resps_this)

598 labels.extend(labels_this)

599

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)