Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/waveform/plot.py: 9%
633 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:25 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:25 +0000
1import logging
2import math
3from collections import defaultdict
5import numpy as num
7from matplotlib import pyplot as plt
8from matplotlib import cm, patches
10from pyrocko import plot, gf, trace
11from pyrocko.plot import mpl_init, mpl_color
12from pyrocko.guts import Tuple, Float, Int, String
14from grond import core, meta
15from .target import WaveformMisfitResult, WaveformMisfitTarget
16from ..plot import StationDistributionPlot
18from grond.plot.config import PlotConfig
19from grond.plot.common import light
20from grond.plot.collection import PlotItem
22guts_prefix = 'grond'
23km = 1e3
24d2r = math.pi / 180.
26logger = logging.getLogger('targets.waveform.plot')
29def make_norm_trace(a, b, exponent):
30 tmin = max(a.tmin, b.tmin)
31 tmax = min(a.tmax, b.tmax)
32 c = a.chop(tmin, tmax, inplace=False)
33 bc = b.chop(tmin, tmax, inplace=False)
34 c.set_ydata(num.abs(c.get_ydata() - bc.get_ydata())**exponent)
35 return c
38def amp_spec_max(spec_trs, key):
39 amaxs = {}
40 for spec_tr in spec_trs:
41 amax = num.max(num.abs(spec_tr.ydata))
42 k = key(spec_tr)
43 if k not in amaxs:
44 amaxs[k] = amax
45 else:
46 amaxs[k] = max(amaxs[k], amax)
48 return amaxs
51def plot_trace(axes, tr, **kwargs):
52 return axes.plot(tr.get_xdata(), tr.get_ydata(), **kwargs)
55def plot_taper(axes, t, taper, **kwargs):
56 y = num.ones(t.size) * 0.9
57 taper(y, t[0], t[1] - t[0])
58 y2 = num.concatenate((y, -y[::-1]))
59 t2 = num.concatenate((t, t[::-1]))
60 axes.fill(t2, y2, **kwargs)
63def plot_dtrace(axes, tr, space, mi, ma, **kwargs):
64 t = tr.get_xdata()
65 y = tr.get_ydata()
66 y2 = (num.concatenate((y, num.zeros(y.size))) - mi) / \
67 (ma - mi) * space - (1.0 + space)
68 t2 = num.concatenate((t, t[::-1]))
69 axes.fill(
70 t2, y2,
71 clip_on=False,
72 **kwargs)
75def plot_spectrum(
76 axes, spec_syn, spec_obs, fmin, fmax, space, mi, ma,
77 syn_color='red', obs_color='black',
78 syn_lw=1.5, obs_lw=1.0, color_vline='gray', fontsize=9.):
80 fpad = (fmax - fmin) / 6.
82 for spec, color, lw in [
83 (spec_syn, syn_color, syn_lw),
84 (spec_obs, obs_color, obs_lw)]:
86 f = spec.get_xdata()
87 mask = num.logical_and(fmin - fpad <= f, f <= fmax + fpad)
89 f = f[mask]
90 y = num.abs(spec.get_ydata())[mask]
92 y2 = (num.concatenate((y, num.zeros(y.size))) - mi) / \
93 (ma - mi) * space - (1.0 + space)
94 f2 = num.concatenate((f, f[::-1]))
95 axes2 = axes.twiny()
96 axes2.set_axis_off()
98 axes2.set_xlim(fmin - fpad * 5, fmax + fpad * 5)
100 axes2.plot(f2, y2, clip_on=False, color=color, lw=lw)
101 axes2.fill(f2, y2, alpha=0.1, clip_on=False, color=color)
103 axes2.plot([fmin, fmin], [-1.0 - space, -1.0], color=color_vline)
104 axes2.plot([fmax, fmax], [-1.0 - space, -1.0], color=color_vline)
106 for (text, fx, ha) in [
107 ('%.3g Hz' % fmin, fmin, 'right'),
108 ('%.3g Hz' % fmax, fmax, 'left')]:
110 axes2.annotate(
111 text,
112 xy=(fx, -1.0),
113 xycoords='data',
114 xytext=(
115 fontsize * 0.4 * [-1, 1][ha == 'left'],
116 -fontsize * 0.2),
117 textcoords='offset points',
118 ha=ha,
119 va='top',
120 color=color_vline,
121 fontsize=fontsize)
124def plot_dtrace_vline(axes, t, space, **kwargs):
125 axes.plot([t, t], [-1.0 - space, -1.0], **kwargs)
128def layout(source, targets, nxmax, nymax):
129 if nxmax == 1 and nymax == 1:
130 nx = len(targets)
131 ny = 1
132 nxx = 1
133 nyy = 1
134 frame_to_target = {}
135 for ix, target in enumerate(sorted(targets, key=lambda t: t.path)):
136 frame_to_target[0, ix] = target
138 return frame_to_target, nx, ny, nxx, nyy
140 nframes = len(targets)
142 nx = int(math.ceil(math.sqrt(nframes)))
143 ny = (nframes - 1) // nx + 1
145 nxx = (nx - 1) // nxmax + 1
146 nyy = (ny - 1) // nymax + 1
148 # nz = nxx * nyy
150 xs = num.arange(nx) // ((max(2, nx) - 1.0) / 2.)
151 ys = num.arange(ny) // ((max(2, ny) - 1.0) / 2.)
153 xs -= num.mean(xs)
154 ys -= num.mean(ys)
156 fxs = num.tile(xs, ny)
157 fys = num.repeat(ys, nx)
159 data = []
161 for target in targets:
162 azi = source.azibazi_to(target)[0]
163 dist = source.distance_to(target)
164 x = dist * num.sin(num.deg2rad(azi))
165 y = dist * num.cos(num.deg2rad(azi))
166 data.append((x, y, dist))
168 gxs, gys, dists = num.array(data, dtype=float).T
170 iorder = num.argsort(dists)
172 gxs = gxs[iorder]
173 gys = gys[iorder]
174 targets_sorted = [targets[ii] for ii in iorder]
176 gxs -= num.mean(gxs)
177 gys -= num.mean(gys)
179 gmax = max(num.max(num.abs(gys)), num.max(num.abs(gxs)))
180 if gmax == 0.:
181 gmax = 1.
183 gxs /= gmax
184 gys /= gmax
186 dists = num.sqrt(
187 (fxs[num.newaxis, :] - gxs[:, num.newaxis])**2 +
188 (fys[num.newaxis, :] - gys[:, num.newaxis])**2)
190 distmax = num.max(dists)
192 availmask = num.ones(dists.shape[1], dtype=bool)
193 frame_to_target = {}
194 for itarget, target in enumerate(targets_sorted):
195 iframe = num.argmin(
196 num.where(availmask, dists[itarget], distmax + 1.))
197 availmask[iframe] = False
198 iy, ix = num.unravel_index(iframe, (ny, nx))
199 frame_to_target[iy, ix] = target
201 return frame_to_target, nx, ny, nxx, nyy
204class CheckWaveformsPlot(PlotConfig):
205 ''' Plot for checking the waveforms fit with a number of synthetics '''
206 name = 'check_waveform'
208 size_cm = Tuple.T(
209 2, Float.T(),
210 default=(9., 7.5),
211 help='width and length of the figure in cm')
212 n_random_synthetics = Int.T(
213 default=10,
214 help='Number of Synthetics to generate')
216 def make(self, environ):
217 cm = environ.get_plot_collection_manager()
218 mpl_init(fontsize=self.font_size)
220 environ.setup_modelling()
222 problem = environ.get_problem()
223 results_list = []
224 sources = []
225 if self.n_random_synthetics == 0:
226 x = problem.preconstrain(problem.get_reference_model())
227 sources.append(problem.base_source)
228 results = problem.evaluate(x)
229 results_list.append(results)
231 else:
232 for _ in range(self.n_random_synthetics):
233 x = problem.get_random_model()
234 sources.append(problem.get_source(x))
235 results = problem.evaluate(x)
236 results_list.append(results)
238 cm.create_group_mpl(self, self.draw_figures(
239 sources, problem.targets, results_list),
240 title=u'Waveform Check',
241 section='checks',
242 feather_icon='activity',
243 description=u'''
244Plot to judge waveform time window settings and source model parameter ranges.
246For each waveform target, observed and synthetic waveforms are shown. For the
247latter, models are randomly drawn from the configured parameter search space.
249The top panel shows the observed waveform; filtered (faint gray), and filtered
250and tapered (black). The colored outline around the observed trace shows the
251taper position for each drawn model in a different color. The middle panel
252shows the filtered synthetic waveforms of the drawn models and the bottom plot
253shows the corresponding filtered and tapered synthetic waveforms. The colors of
254taper and synthetic traces are consistent for each random model. The given time
255is relative to the reference event origin time.
256''')
258 def draw_figures(self, sources, targets, results_list):
259 results_list = list(zip(*results_list))
260 for itarget, target, results in zip(
261 range(len(targets)), targets, results_list):
263 if isinstance(target, WaveformMisfitTarget) and results:
264 item = PlotItem(name='t%i' % itarget)
265 item.attributes['targets'] = [target.string_id()]
266 fig = self.draw_figure(sources, target, results)
267 if fig is not None:
268 yield item, fig
270 def draw_figure(self, sources, target, results):
271 t0_mean = num.mean([s.time for s in sources])
273 # distances = [
274 # s.distance_to(target) for s in sources]
276 # distance_min = num.min(distances)
277 # distance_max = num.max(distances)
279 yabsmaxs = []
281 for result in results:
282 if isinstance(result, WaveformMisfitResult):
283 yabsmaxs.append(
284 num.max(num.abs(
285 result.filtered_obs.get_ydata())))
287 if yabsmaxs:
288 yabsmax = max(yabsmaxs) or 1.0
289 else:
290 yabsmax = None
292 fontsize = self.font_size
294 fig = plt.figure(figsize=self.size_inch)
296 labelpos = plot.mpl_margins(
297 fig, nw=1, nh=1, w=1., h=5.,
298 units=fontsize)
300 axes = fig.add_subplot(1, 1, 1)
302 labelpos(axes, 2.5, 2.0)
303 axes.set_frame_on(False)
304 axes.set_ylim(1., 4.)
305 axes.get_yaxis().set_visible(False)
306 axes.set_title('%s' % target.string_id())
307 axes.set_xlabel('Time [s]')
308 ii = 0
310 for source, result in zip(sources, results):
311 if not isinstance(result, WaveformMisfitResult):
312 continue
314 if result.tobs_shift != 0.0:
315 t0 = result.tsyn_pick
316 else:
317 t0 = t0_mean
319 t = result.filtered_obs.get_xdata()
320 ydata = result.filtered_obs.get_ydata() / yabsmax
321 axes.plot(
322 t-t0, ydata*0.40 + 3.5, color='black', lw=1.0)
324 color = plot.mpl_graph_color(ii)
326 t = result.filtered_syn.get_xdata()
327 ydata = result.filtered_syn.get_ydata()
328 ydata = ydata / (num.max(num.abs(ydata)) or 1.0)
330 axes.plot(t-t0, ydata*0.47 + 2.5, color=color, alpha=0.5, lw=1.0)
332 t = result.processed_syn.get_xdata()
333 ydata = result.processed_syn.get_ydata()
334 ydata = ydata / (num.max(num.abs(ydata)) or 1.0)
336 axes.plot(t-t0, ydata*0.47 + 1.5, color=color, alpha=0.5, lw=1.0)
337 if result.tobs_shift != 0.0:
338 axes.axvline(
339 result.tsyn_pick - t0,
340 color=(0.7, 0.7, 0.7),
341 zorder=2)
343 t = result.processed_syn.get_xdata()
344 taper = result.taper
346 y = num.ones(t.size) * 0.9
347 taper(y, t[0], t[1] - t[0])
348 y2 = num.concatenate((y, -y[::-1]))
349 t2 = num.concatenate((t, t[::-1]))
350 axes.plot(t2-t0, y2 * 0.47 + 3.5, color=color, alpha=0.2, lw=1.0)
351 ii += 1
353 return fig
356class FitsWaveformEnsemblePlot(PlotConfig):
357 ''' Plot showing all waveform fits for the ensemble of solutions'''
359 name = 'fits_waveform_ensemble'
360 size_cm = Tuple.T(
361 2, Float.T(),
362 default=(9., 5.),
363 help='width and length of the figure in cm')
364 nx = Int.T(
365 default=1,
366 help='horizontal number of subplots on every page')
367 ny = Int.T(
368 default=1,
369 help='vertical number of subplots on every page')
370 misfit_cutoff = Float.T(
371 optional=True,
372 help='Plot fits for models up to this misfit value')
373 color_parameter = String.T(
374 default='misfit',
375 help='Choice of value to color, options: dist and misfit')
376 font_size = Float.T(
377 default=8,
378 help='Font Size of all fonts, except title')
379 font_size_title = Float.T(
380 default=10,
381 help='Font size of title')
383 def make(self, environ):
384 cm = environ.get_plot_collection_manager()
385 mpl_init(fontsize=self.font_size)
386 environ.setup_modelling()
387 ds = environ.get_dataset()
388 history = environ.get_history(subset='harvest')
389 optimiser = environ.get_optimiser()
391 cm.create_group_mpl(
392 self,
393 self.draw_figures(ds, history, optimiser),
394 title=u'Waveform fits for the ensemble',
395 section='fits',
396 feather_icon='activity',
397 description=u'''
398Plot showing waveform (attribute) fits for the ensemble of solutions.
400Waveform fits for every nth model in the ensemble of bootstrap solutions.
401Depending on the target configuration different types of comparisons are
402possible: (i) time domain waveform differences, (ii) amplitude spectra, (iii)
403envelopes, (iv) cross correlation functions. Each waveform plot gives a number
404of details:
4061) Target information (left side, from top to bottom) gives station name with
407component, distance to source, azimuth of station with respect to source,
408target weight, target misfit and starting time of the waveform relative to the
409origin time.
4112) The background gray area shows the applied taper function.
4133) The waveforms shown are: the restituted and filtered observed trace without
414tapering (light grey) and the same trace with tapering and processing (dark
415gray), the synthetic trace (light red) and the filtered, tapered and (if
416enabled) shifted and processed synthetic trace (colored). The colors of the
417synthetic traces indicate how well the corresponding models fit in the global
418weighting scheme (when all bootstrap weights are equal), from better fit (red)
419to worse fit (blue). The amplitudes of the traces are scaled according to the
420target weight (small weight, small amplitude) and normed relative to the
421maximum amplitude of the targets of the corresponding normalisation family.
4234) The bottom panel shows, depending on the type of comparison, sample-wise
424residuals for time domain comparisons (red filled), spectra of observed and
425synthetic traces for amplitude spectrum comparisons, or cross correlation
426traces.''')
428 def draw_figures(self, ds, history, optimiser):
430 color_parameter = self.color_parameter
431 misfit_cutoff = self.misfit_cutoff
432 fontsize = self.font_size
433 fontsize_title = self.font_size_title
435 nxmax = self.nx
436 nymax = self.ny
438 problem = history.problem
440 for target in problem.targets:
441 target.set_dataset(ds)
443 target_index = {}
444 i = 0
445 for target in problem.targets:
446 target_index[target] = i, i+target.nmisfits
447 i += target.nmisfits
449 gms = history.get_sorted_primary_misfits()[::-1]
450 models = history.get_sorted_primary_models()[::-1]
452 if misfit_cutoff is not None:
453 ibest = gms < misfit_cutoff
454 gms = gms[ibest]
455 models = models[ibest]
457 gms = gms[::10]
458 models = models[::10]
460 nmodels = models.shape[0]
461 if color_parameter == 'dist':
462 mx = num.mean(models, axis=0)
463 cov = num.cov(models.T)
464 mdists = core.mahalanobis_distance(models, mx, cov)
465 icolor = meta.ordersort(mdists)
467 elif color_parameter == 'misfit':
468 iorder = num.arange(nmodels)
469 icolor = iorder
471 elif color_parameter in problem.parameter_names:
472 ind = problem.name_to_index(color_parameter)
473 icolor = problem.extract(models, ind)
475 target_to_results = defaultdict(list)
476 all_syn_trs = []
478 dtraces = []
479 for imodel in range(nmodels):
480 model = models[imodel, :]
482 source = problem.get_source(model)
483 results = problem.evaluate(model)
485 dtraces.append([])
487 for target, result in zip(problem.targets, results):
488 w = target.get_combined_weight()
490 if isinstance(result, gf.SeismosizerError) or \
491 not isinstance(target, WaveformMisfitTarget) or \
492 not num.all(num.isfinite(w)):
494 dtraces[-1].extend([None] * target.nmisfits)
495 continue
497 itarget, itarget_end = target_index[target]
498 assert itarget_end == itarget + 1
500 if target.misfit_config.domain == 'cc_max_norm':
501 tref = (
502 result.filtered_obs.tmin + result.filtered_obs.tmax) \
503 * 0.5
505 for tr_filt, tr_proc, tshift in (
506 (result.filtered_obs,
507 result.processed_obs,
508 0.),
509 (result.filtered_syn,
510 result.processed_syn,
511 result.tshift)):
513 norm = num.sum(num.abs(tr_proc.ydata)) \
514 / tr_proc.data_len()
515 tr_filt.ydata /= norm
516 tr_proc.ydata /= norm
518 tr_filt.shift(tshift)
519 tr_proc.shift(tshift)
521 ctr = result.cc
522 ctr.shift(tref)
524 dtrace = ctr
526 else:
527 for tr in (
528 result.filtered_obs,
529 result.filtered_syn,
530 result.processed_obs,
531 result.processed_syn):
533 tr.ydata *= w
535 if result.tshift is not None and result.tshift != 0.0:
536 # result.filtered_syn.shift(result.tshift)
537 result.processed_syn.shift(result.tshift)
539 dtrace = make_norm_trace(
540 result.processed_syn, result.processed_obs,
541 problem.norm_exponent)
543 target_to_results[target].append(result)
545 dtrace.meta = dict(
546 normalisation_family=target.normalisation_family,
547 path=target.path)
549 dtraces[-1].append(dtrace)
551 result.processed_syn.meta = dict(
552 normalisation_family=target.normalisation_family,
553 path=target.path)
555 all_syn_trs.append(result.processed_syn)
557 if not all_syn_trs:
558 logger.warning('No traces to show!')
559 return
561 def skey(tr):
562 return tr.meta['normalisation_family'], tr.meta['path']
564 trace_minmaxs = trace.minmax(all_syn_trs, skey)
566 dtraces_all = []
567 for dtraces_group in dtraces:
568 dtraces_all.extend(dtraces_group)
570 dminmaxs = trace.minmax([
571 dtrace_ for dtrace_ in dtraces_all if dtrace_ is not None], skey)
573 for tr in dtraces_all:
574 if tr:
575 dmin, dmax = dminmaxs[skey(tr)]
576 tr.ydata /= max(abs(dmin), abs(dmax))
578 cg_to_targets = meta.gather(
579 problem.waveform_targets,
580 lambda t: (t.path, t.codes[3]),
581 filter=lambda t: t in target_to_results)
583 cgs = sorted(cg_to_targets.keys())
585 from matplotlib import colors
586 cmap = cm.ScalarMappable(
587 norm=colors.Normalize(vmin=num.min(icolor), vmax=num.max(icolor)),
588 cmap=plt.get_cmap('coolwarm'))
590 imodel_to_color = []
591 for imodel in range(nmodels):
592 imodel_to_color.append(cmap.to_rgba(icolor[imodel]))
594 for cg in cgs:
595 targets = cg_to_targets[cg]
597 frame_to_target, nx, ny, nxx, nyy = layout(
598 source, targets, nxmax, nymax)
600 figures = {}
601 for iy in range(ny):
602 for ix in range(nx):
603 if (iy, ix) not in frame_to_target:
604 continue
606 ixx = ix // nxmax
607 iyy = iy // nymax
608 if (iyy, ixx) not in figures:
609 title = '_'.join(x for x in cg if x)
610 item = PlotItem(
611 name='fig_%s_%i_%i' % (title, ixx, iyy))
612 item.attributes['targets'] = []
613 figures[iyy, ixx] = (
614 item, plt.figure(figsize=self.size_inch))
616 figures[iyy, ixx][1].subplots_adjust(
617 left=0.03,
618 right=1.0 - 0.03,
619 bottom=0.03,
620 top=1.0 - 0.06,
621 wspace=0.2,
622 hspace=0.2)
624 item, fig = figures[iyy, ixx]
626 target = frame_to_target[iy, ix]
628 item.attributes['targets'].append(target.string_id())
630 amin, amax = trace_minmaxs[
631 target.normalisation_family, target.path]
632 absmax = max(abs(amin), abs(amax))
634 ny_this = nymax # min(ny, nymax)
635 nx_this = nxmax # min(nx, nxmax)
636 i_this = (iy % ny_this) * nx_this + (ix % nx_this) + 1
638 axes2 = fig.add_subplot(ny_this, nx_this, i_this)
640 space = 0.5
641 space_factor = 1.0 + space
642 axes2.set_axis_off()
643 axes2.set_ylim(-1.05 * space_factor, 1.05)
645 axes = axes2.twinx()
646 axes.set_axis_off()
648 if target.misfit_config.domain == 'cc_max_norm':
649 axes.set_ylim(-10. * space_factor, 10.)
650 else:
651 axes.set_ylim(-absmax*1.33 * space_factor, absmax*1.33)
653 itarget, itarget_end = target_index[target]
654 assert itarget_end == itarget + 1
656 for imodel, result in enumerate(target_to_results[target]):
658 syn_color = imodel_to_color[imodel]
660 dtrace = dtraces[imodel][itarget]
662 tap_color_annot = (0.35, 0.35, 0.25)
663 tap_color_edge = (0.85, 0.85, 0.80)
664 tap_color_fill = (0.95, 0.95, 0.90)
666 plot_taper(
667 axes2,
668 result.processed_obs.get_xdata(),
669 result.taper,
670 fc=tap_color_fill, ec=tap_color_edge, alpha=0.2)
672 obs_color = mpl_color('aluminium5')
673 obs_color_light = light(obs_color, 0.5)
675 plot_dtrace(
676 axes2, dtrace, space, 0., 1.,
677 fc='none',
678 ec=syn_color)
680 # plot_trace(
681 # axes, result.filtered_syn,
682 # color=syn_color_light, lw=1.0)
684 if imodel == 0:
685 plot_trace(
686 axes, result.filtered_obs,
687 color=obs_color_light, lw=0.75)
689 plot_trace(
690 axes, result.processed_syn,
691 color=syn_color, lw=1.0, alpha=0.3)
693 plot_trace(
694 axes, result.processed_obs,
695 color=obs_color, lw=0.75, alpha=0.3)
697 if imodel != 0:
698 continue
699 xdata = result.filtered_obs.get_xdata()
700 axes.set_xlim(xdata[0], xdata[-1])
702 tmarks = [
703 result.processed_obs.tmin,
704 result.processed_obs.tmax]
706 for tmark in tmarks:
707 axes2.plot(
708 [tmark, tmark], [-0.9, 0.1],
709 color=tap_color_annot)
711 dur = tmarks[1] - tmarks[0]
712 for tmark, text, ha in [
713 (tmarks[0],
714 '$\\,$ ' + meta.str_duration(
715 tmarks[0] - source.time),
716 'left'),
717 (tmarks[1],
718 '$\\Delta$ ' + meta.str_duration(
719 dur),
720 'right')]:
722 axes2.annotate(
723 text,
724 xy=(tmark, -0.9),
725 xycoords='data',
726 xytext=(
727 fontsize*0.4 * [-1, 1][ha == 'left'],
728 fontsize*0.2),
729 textcoords='offset points',
730 ha=ha,
731 va='bottom',
732 color=tap_color_annot,
733 fontsize=fontsize)
735 axes2.set_xlim(
736 tmarks[0] - dur*0.1, tmarks[1] + dur*0.1)
738 scale_string = None
740 if target.misfit_config.domain == 'cc_max_norm':
741 scale_string = 'Syn/obs scales differ!'
743 infos = []
744 if scale_string:
745 infos.append(scale_string)
747 if self.nx == 1 and self.ny == 1:
748 infos.append(target.string_id())
749 else:
750 infos.append('.'.join(x for x in target.codes if x))
751 dist = source.distance_to(target)
752 azi = source.azibazi_to(target)[0]
753 infos.append(meta.str_dist(dist))
754 infos.append(u'%.0f\u00B0' % azi)
755 axes2.annotate(
756 '\n'.join(infos),
757 xy=(0., 1.),
758 xycoords='axes fraction',
759 xytext=(2., 2.),
760 textcoords='offset points',
761 ha='left',
762 va='top',
763 fontsize=fontsize,
764 fontstyle='normal')
766 if (self.nx == 1 and self.ny == 1):
767 yield item, fig
768 del figures[iyy, ixx]
770 if not (self.nx == 1 and self.ny == 1):
771 for (iyy, ixx), (_, fig) in figures.items():
772 title = '.'.join(x for x in cg if x)
773 if len(figures) > 1:
774 title += ' (%i/%i, %i/%i)' % (iyy+1, nyy, ixx+1, nxx)
776 fig.suptitle(title, fontsize=fontsize_title)
778 for item, fig in figures.values():
779 yield item, fig
782class FitsWaveformPlot(PlotConfig):
783 ''' Plot showing the waveform fits for the best model '''
784 name = 'fits_waveform'
785 size_cm = Tuple.T(
786 2, Float.T(),
787 default=(9., 5.),
788 help='width and length of the figure in cm')
789 nx = Int.T(
790 default=1,
791 help='horizontal number of subplots on every page')
792 ny = Int.T(
793 default=1,
794 help='vertical number of subplots on every page')
795 font_size = Float.T(
796 default=8,
797 help='Font Size of all fonts, except title')
798 font_size_title = Float.T(
799 default=10,
800 help='Font Size of title')
802 def make(self, environ):
803 cm = environ.get_plot_collection_manager()
804 mpl_init(fontsize=self.font_size)
805 environ.setup_modelling()
806 ds = environ.get_dataset()
807 optimiser = environ.get_optimiser()
809 environ.setup_modelling()
811 history = environ.get_history(subset='harvest')
812 cm.create_group_mpl(
813 self,
814 self.draw_figures(ds, history, optimiser),
815 title=u'Waveform fits for best model',
816 section='fits',
817 feather_icon='activity',
818 description=u'''
819Plot showing observed and synthetic waveform (attributes) for the best fitting
820model.
822Best model's waveform fits for all targets. Depending on the target
823configurations different types of comparisons are possible: (i) time domain
824waveform differences, (ii) amplitude spectra, (iii) envelopes, (iv) cross
825correlation functions. Each waveform plot gives a number of details:
8271) Target information (left side, from top to bottom) gives station name with
828component, distance to source, azimuth of station with respect to source,
829target weight, target misfit and starting time of the waveform relative to the
830origin time.
8322) The background gray area shows the applied taper function.
8343) The waveforms shown are: the restituted and filtered observed trace without
835tapering (light grey) and the same trace with tapering and processing (dark
836gray), the synthetic trace (light red) and the filtered, tapered and (if
837enabled) shifted and processed synthetic target trace (red). The traces are
838scaled according to the target weight (small weight, small amplitude) and
839normed relative to the maximum amplitude of the targets of the corresponding
840normalisation family.
8424) The bottom panel shows, depending on the type of comparison, sample-wise
843residuals for time domain comparisons (red filled), spectra of observed and
844synthetic traces for amplitude spectrum comparisons, or cross correlation
845traces.
8475) Colored boxes on the upper right show the relative weight of the target
848within the entire dataset of the optimisation (top box, orange) and the
849relative misfit contribution to the global misfit of the optimisation (bottom
850box, red).
851''')
853 def draw_figures(self, ds, history, optimiser):
855 fontsize = self.font_size
856 fontsize_title = self.font_size_title
858 nxmax = self.nx
859 nymax = self.ny
861 problem = history.problem
863 for target in problem.targets:
864 target.set_dataset(ds)
866 target_index = {}
867 i = 0
868 for target in problem.targets:
869 target_index[target] = i, i+target.nmisfits
870 i += target.nmisfits
872 xbest = history.get_best_model()
873 misfits = history.misfits[history.get_sorted_misfits_idx(chain=0), ...]
875 ws = problem.get_target_weights()
877 gcms = problem.combine_misfits(
878 misfits[:1, :, :],
879 extra_correlated_weights=optimiser.get_correlated_weights(problem),
880 get_contributions=True)[0, :]
882 w_max = num.nanmax(ws)
883 gcm_max = num.nanmax(gcms)
885 source = problem.get_source(xbest)
887 target_to_result = {}
888 all_syn_trs = []
889 all_syn_specs = []
890 results = problem.evaluate(xbest)
892 dtraces = []
893 for target, result in zip(problem.targets, results):
894 if not isinstance(result, WaveformMisfitResult):
895 dtraces.extend([None] * target.nmisfits)
896 continue
898 itarget, itarget_end = target_index[target]
899 assert itarget_end == itarget + 1
901 w = target.get_combined_weight()
903 if target.misfit_config.domain == 'cc_max_norm':
904 tref = (
905 result.filtered_obs.tmin + result.filtered_obs.tmax) * 0.5
906 for tr_filt, tr_proc, tshift in (
907 (result.filtered_obs,
908 result.processed_obs,
909 0.),
910 (result.filtered_syn,
911 result.processed_syn,
912 result.tshift)):
914 norm = num.sum(num.abs(tr_proc.ydata)) / tr_proc.data_len()
915 tr_filt.ydata /= norm
916 tr_proc.ydata /= norm
918 tr_filt.shift(tshift)
919 tr_proc.shift(tshift)
921 ctr = result.cc
922 ctr.shift(tref)
924 dtrace = ctr
926 else:
927 for tr in (
928 result.filtered_obs,
929 result.filtered_syn,
930 result.processed_obs,
931 result.processed_syn):
933 tr.ydata *= w
935 for spec in (
936 result.spectrum_obs,
937 result.spectrum_syn):
939 if spec is not None:
940 spec.ydata *= w
942 if result.tshift is not None and result.tshift != 0.0:
943 # result.filtered_syn.shift(result.tshift)
944 result.processed_syn.shift(result.tshift)
946 dtrace = make_norm_trace(
947 result.processed_syn, result.processed_obs,
948 problem.norm_exponent)
950 target_to_result[target] = result
952 dtrace.meta = dict(
953 normalisation_family=target.normalisation_family,
954 path=target.path)
955 dtraces.append(dtrace)
957 result.processed_syn.meta = dict(
958 normalisation_family=target.normalisation_family,
959 path=target.path)
961 all_syn_trs.append(result.processed_syn)
963 if result.spectrum_syn:
964 result.spectrum_syn.meta = dict(
965 normalisation_family=target.normalisation_family,
966 path=target.path)
968 all_syn_specs.append(result.spectrum_syn)
970 if not all_syn_trs:
971 logger.warning('No traces to show!')
972 return
974 def skey(tr):
975 return tr.meta['normalisation_family'], tr.meta['path']
977 trace_minmaxs = trace.minmax(all_syn_trs, skey)
979 amp_spec_maxs = amp_spec_max(all_syn_specs, skey)
981 dminmaxs = trace.minmax([x for x in dtraces if x is not None], skey)
983 for tr in dtraces:
984 if tr:
985 dmin, dmax = dminmaxs[skey(tr)]
986 tr.ydata /= max(abs(dmin), abs(dmax))
988 cg_to_targets = meta.gather(
989 problem.waveform_targets,
990 lambda t: (t.path, t.codes[3]),
991 filter=lambda t: t in target_to_result)
993 cgs = sorted(cg_to_targets.keys())
995 for cg in cgs:
996 targets = cg_to_targets[cg]
998 frame_to_target, nx, ny, nxx, nyy = layout(
999 source, targets, nxmax, nymax)
1001 figures = {}
1002 for iy in range(ny):
1003 for ix in range(nx):
1004 if (iy, ix) not in frame_to_target:
1005 continue
1007 ixx = ix // nxmax
1008 iyy = iy // nymax
1009 if (iyy, ixx) not in figures:
1010 title = '_'.join(x for x in cg if x)
1011 item = PlotItem(
1012 name='fig_%s_%i_%i' % (title, ixx, iyy))
1013 item.attributes['targets'] = []
1014 figures[iyy, ixx] = (
1015 item, plt.figure(figsize=self.size_inch))
1017 figures[iyy, ixx][1].subplots_adjust(
1018 left=0.03,
1019 right=1.0 - 0.03,
1020 bottom=0.03,
1021 top=1.0 - 0.06,
1022 wspace=0.2,
1023 hspace=0.2)
1025 item, fig = figures[iyy, ixx]
1027 target = frame_to_target[iy, ix]
1029 item.attributes['targets'].append(target.string_id())
1031 amin, amax = trace_minmaxs[
1032 target.normalisation_family, target.path]
1033 absmax = max(abs(amin), abs(amax))
1035 ny_this = nymax # min(ny, nymax)
1036 nx_this = nxmax # min(nx, nxmax)
1037 i_this = (iy % ny_this) * nx_this + (ix % nx_this) + 1
1039 axes2 = fig.add_subplot(ny_this, nx_this, i_this)
1041 space = 0.5
1042 space_factor = 1.0 + space
1043 axes2.set_axis_off()
1044 axes2.set_ylim(-1.05 * space_factor, 1.05)
1046 axes = axes2.twinx()
1047 axes.set_axis_off()
1049 if target.misfit_config.domain == 'cc_max_norm':
1050 axes.set_ylim(-10. * space_factor, 10.)
1051 else:
1052 axes.set_ylim(
1053 -absmax * 1.33 * space_factor, absmax * 1.33)
1055 itarget, itarget_end = target_index[target]
1056 assert itarget_end == itarget + 1
1058 result = target_to_result[target]
1060 dtrace = dtraces[itarget]
1062 tap_color_annot = (0.35, 0.35, 0.25)
1063 tap_color_edge = (0.85, 0.85, 0.80)
1064 tap_color_fill = (0.95, 0.95, 0.90)
1066 plot_taper(
1067 axes2, result.processed_obs.get_xdata(), result.taper,
1068 fc=tap_color_fill, ec=tap_color_edge)
1070 obs_color = mpl_color('aluminium5')
1071 obs_color_light = light(obs_color, 0.5)
1073 syn_color = mpl_color('scarletred2')
1074 syn_color_light = light(syn_color, 0.5)
1076 misfit_color = mpl_color('scarletred2')
1077 weight_color = mpl_color('chocolate2')
1079 cc_color = mpl_color('aluminium5')
1081 if target.misfit_config.domain == 'cc_max_norm':
1082 tref = (result.filtered_obs.tmin +
1083 result.filtered_obs.tmax) * 0.5
1085 plot_dtrace(
1086 axes2, dtrace, space, -1., 1.,
1087 fc=light(cc_color, 0.5),
1088 ec=cc_color)
1090 plot_dtrace_vline(
1091 axes2, tref, space, color=tap_color_annot)
1093 elif target.misfit_config.domain == 'frequency_domain':
1095 asmax = amp_spec_maxs[
1096 target.normalisation_family, target.path]
1097 fmin, fmax = \
1098 target.misfit_config.get_full_frequency_range()
1100 plot_spectrum(
1101 axes2,
1102 result.spectrum_syn,
1103 result.spectrum_obs,
1104 fmin, fmax,
1105 space, 0., asmax,
1106 syn_color=syn_color,
1107 obs_color=obs_color,
1108 syn_lw=1.0,
1109 obs_lw=0.75,
1110 color_vline=tap_color_annot,
1111 fontsize=fontsize)
1113 else:
1114 plot_dtrace(
1115 axes2, dtrace, space, 0., 1.,
1116 fc=light(misfit_color, 0.3),
1117 ec=misfit_color)
1119 plot_trace(
1120 axes, result.filtered_syn,
1121 color=syn_color_light, lw=1.0)
1123 plot_trace(
1124 axes, result.filtered_obs,
1125 color=obs_color_light, lw=0.75)
1127 plot_trace(
1128 axes, result.processed_syn,
1129 color=syn_color, lw=1.0)
1131 plot_trace(
1132 axes, result.processed_obs,
1133 color=obs_color, lw=0.75)
1135 # xdata = result.filtered_obs.get_xdata()
1137 tmarks = [
1138 result.processed_obs.tmin,
1139 result.processed_obs.tmax]
1141 for tmark in tmarks:
1142 axes2.plot(
1143 [tmark, tmark], [-0.9, 0.1], color=tap_color_annot)
1145 dur = tmarks[1] - tmarks[0]
1146 for tmark, text, ha in [
1147 (tmarks[0],
1148 '$\\,$ ' + meta.str_duration(
1149 tmarks[0] - source.time),
1150 'left'),
1151 (tmarks[1],
1152 '$\\Delta$ ' + meta.str_duration(dur),
1153 'right')]:
1155 axes2.annotate(
1156 text,
1157 xy=(tmark, -0.9),
1158 xycoords='data',
1159 xytext=(
1160 fontsize * 0.4 * [-1, 1][ha == 'left'],
1161 fontsize * 0.2),
1162 textcoords='offset points',
1163 ha=ha,
1164 va='bottom',
1165 color=tap_color_annot,
1166 fontsize=fontsize)
1168 axes2.set_xlim(tmarks[0] - dur*0.1, tmarks[1] + dur*0.1)
1170 rel_w = ws[itarget] / w_max
1171 rel_c = gcms[itarget] / gcm_max
1173 sw = 0.25
1174 sh = 0.1
1175 ph = 0.01
1177 for (ih, rw, facecolor, edgecolor) in [
1178 (0, rel_w, light(weight_color, 0.5),
1179 weight_color),
1180 (1, rel_c, light(misfit_color, 0.5),
1181 misfit_color)]:
1183 bar = patches.Rectangle(
1184 (1.0 - rw * sw, 1.0 - (ih + 1) * sh + ph),
1185 rw * sw,
1186 sh - 2 * ph,
1187 facecolor=facecolor, edgecolor=edgecolor,
1188 zorder=10,
1189 transform=axes.transAxes, clip_on=False)
1191 axes.add_patch(bar)
1193 scale_string = None
1195 if target.misfit_config.domain == 'cc_max_norm':
1196 scale_string = 'Syn/obs scales differ!'
1198 infos = []
1199 if scale_string:
1200 infos.append(scale_string)
1202 if self.nx == 1 and self.ny == 1:
1203 infos.append(target.string_id())
1204 else:
1205 infos.append('.'.join(x for x in target.codes if x))
1207 dist = source.distance_to(target)
1208 azi = source.azibazi_to(target)[0]
1209 infos.append(meta.str_dist(dist))
1210 infos.append('%.0f\u00B0' % azi)
1211 infos.append('%.3g' % ws[itarget])
1212 infos.append('%.3g' % gcms[itarget])
1213 axes2.annotate(
1214 '\n'.join(infos),
1215 xy=(0., 1.),
1216 xycoords='axes fraction',
1217 xytext=(2., 2.),
1218 textcoords='offset points',
1219 ha='left',
1220 va='top',
1221 fontsize=fontsize,
1222 fontstyle='normal')
1224 if (self.nx == 1 and self.ny == 1):
1225 yield item, fig
1226 del figures[iyy, ixx]
1228 if not (self.nx == 1 and self.ny == 1):
1229 for (iyy, ixx), (_, fig) in figures.items():
1230 title = '.'.join(x for x in cg if x)
1231 if len(figures) > 1:
1232 title += ' (%i/%i, %i/%i)' % (
1233 iyy + 1, nyy, ixx + 1, nxx)
1235 fig.suptitle(title, fontsize=fontsize_title)
1237 for item, fig in figures.values():
1238 yield item, fig
1241class WaveformStationDistribution(StationDistributionPlot):
1242 ''' Plot showing all waveform fits for the ensemble of solutions'''
1244 name = 'seismic_stations'
1246 def make(self, environ):
1247 from grond.problems.base import ProblemDataNotAvailable
1248 from grond.environment import NoRundirAvailable
1250 cm = environ.get_plot_collection_manager()
1251 mpl_init(fontsize=self.font_size)
1253 problem = environ.get_problem()
1254 dataset = environ.get_dataset()
1255 try:
1256 history = environ.get_history(subset='harvest')
1257 except (NoRundirAvailable, ProblemDataNotAvailable):
1258 history = None
1260 cm.create_group_mpl(
1261 self,
1262 self.draw_figures(problem, dataset, history),
1263 title=u'Seismic station locations',
1264 section='checks',
1265 feather_icon='target',
1266 description=u'''
1267Plot showing seismic station locations and attributes.
1269Station locations in dependence of distance and azimuth are shown. The center
1270of the plot corresponds to the origin of the search space, not to the optimised
1271location of the source.
1272''')
1274 def draw_figures(self, problem, dataset, history):
1276 target_index = {}
1277 i = 0
1278 for target in problem.targets:
1279 target_index[target] = i, i+target.nmisfits
1280 i += target.nmisfits
1282 ws = problem.get_target_weights()
1284 if history:
1285 misfits = history.misfits[history.get_sorted_misfits_idx(), ...]
1286 gcms = problem.combine_misfits(
1287 misfits[:1, :, :], get_contributions=True)[0, :]
1289 event = problem.base_source
1291 cg_to_targets = meta.gather(
1292 problem.waveform_targets,
1293 lambda t: (t.path, t.codes[3]))
1295 cgs = sorted(cg_to_targets.keys())
1297 for cg in cgs:
1298 cg_str = '.'.join(cg)
1300 targets = cg_to_targets[cg]
1301 if len(targets) == 0:
1302 continue
1304 assert all(target_index[target][0] == target_index[target][1] - 1
1305 for target in targets)
1307 itargets = num.array(
1308 [target_index[target][0] for target in targets])
1310 labels = ['.'.join(x for x in t.codes[:3] if x) for t in targets]
1312 azimuths = num.array([event.azibazi_to(t)[0] for t in targets])
1313 distances = num.array([t.distance_to(event) for t in targets])
1315 item = PlotItem(
1316 name='seismic_stations_weights_%s' % cg_str,
1317 title=u'Station weights (%s)' % cg_str,
1318 description=u'\n\nMarkers are scaled according to the '
1319 u'weighting factor of the corresponding target\'s '
1320 u'contribution in the misfit function.')
1321 fig, ax, legend = self.plot_station_distribution(
1322 azimuths, distances, ws[itargets], labels,
1323 legend_title='Weight')
1325 yield (item, fig)
1327 if history:
1328 item = PlotItem(
1329 name='seismic_stations_contributions_%s' % cg_str,
1330 title=u'Station misfit contributions (%s)' % cg_str,
1331 description=u'\n\nMarkers are scaled according to their '
1332 u'misfit contribution for the globally best '
1333 u'source model.')
1334 fig, ax, legend = self.plot_station_distribution(
1335 azimuths, distances, gcms[itargets], labels,
1336 legend_title='Contribution')
1338 yield (item, fig)
1341def get_plot_classes():
1342 return [
1343 WaveformStationDistribution,
1344 CheckWaveformsPlot,
1345 FitsWaveformPlot,
1346 FitsWaveformEnsemblePlot
1347 ]