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. 

8 

9Example 

10 

11* :download:`test_response_plot.py </../../examples/test_response_plot.py>` 

12 

13:: 

14 

15 from pyrocko.plot import response 

16 from pyrocko.example import get_example_data 

17 

18 get_example_data('test_response.resp') 

19 

20 resps, labels = response.load_response_information( 

21 'test_response.resp', 'resp') 

22 

23 response.plot( 

24 responses=resps, labels=labels, filename='test_response.png', 

25 fmin=0.001, fmax=400., dpi=75.) 

26 

27 

28.. figure :: /static/test_response.png 

29 :align: center 

30 

31 Example response plot 

32''' 

33 

34import logging 

35import math 

36import hashlib 

37from io import BytesIO 

38 

39import numpy as num 

40 

41from pyrocko import util 

42from pyrocko import guts 

43from pyrocko.response import InvalidResponseError 

44 

45 

46logger = logging.getLogger('plot.response') 

47 

48 

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

53 

54 

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

65 

66 ''' 

67 Draw instrument response in Bode plot style to given Matplotlib axes 

68 

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

83 

84 f = num.exp(num.linspace(num.log(fmin), num.log(fmax), nf)) 

85 

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

93 

94 f = f[f <= resp_fmax] 

95 

96 if f.size == 0: 

97 return 

98 

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.fdsn import station 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)