Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/plot.py: 14%
508 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 19:53 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 19:53 +0000
1import math
2import logging
3import random
5import numpy as num
6from matplotlib import cm, patches, colors as mcolors
8from pyrocko.guts import Tuple, Float, Int, String, List, Bool, StringChoice
10from pyrocko.plot import mpl_margins, mpl_color, mpl_init, mpl_graph_color
11from pyrocko.plot import beachball, hudson
13from grond.plot.config import PlotConfig
14from grond.plot.section import SectionPlotConfig, SectionPlot
15from grond.plot.collection import PlotItem
16from grond import meta, core, stats
17from matplotlib import pyplot as plt
19logger = logging.getLogger('grond.problem.plot')
21guts_prefix = 'grond'
24def cluster_label(icluster, perc):
25 if icluster == -1:
26 return 'Unclust. (%.0f%%)' % perc
27 else:
28 return 'Cluster %i (%.0f%%)' % (icluster, perc)
31def cluster_color(icluster):
32 if icluster == -1:
33 return mpl_color('aluminium3')
34 else:
35 return mpl_graph_color(icluster)
38def fixlim(lo, hi):
39 if lo == hi:
40 return lo - 1.0, hi + 1.0
41 else:
42 return lo, hi
45def eigh_sorted(mat):
46 evals, evecs = num.linalg.eigh(mat)
47 iorder = num.argsort(evals)
48 return evals[iorder], evecs[:, iorder]
51class JointparPlot(PlotConfig):
52 '''
53 Source problem parameter's tradeoff plots.
54 '''
56 name = 'jointpar'
57 size_cm = Tuple.T(2, Float.T(), default=(20., 20.))
58 misfit_cutoff = Float.T(optional=True)
59 ibootstrap = Int.T(optional=True)
60 color_parameter = String.T(default='misfit')
61 exclude = List.T(String.T())
62 include = List.T(String.T())
63 show_ellipses = Bool.T(default=False)
64 nsubplots = Int.T(default=6)
65 show_ticks = Bool.T(default=False)
66 show_reference = Bool.T(default=True)
67 show_axes = Bool.T(default=False)
69 def make(self, environ):
70 cm = environ.get_plot_collection_manager()
71 history = environ.get_history(subset='harvest')
72 optimiser = environ.get_optimiser()
73 sref = 'Dark gray boxes mark reference solution.' \
74 if self.show_reference else ''
76 mpl_init(fontsize=self.font_size)
77 cm.create_group_mpl(
78 self,
79 self.draw_figures(history, optimiser),
80 title=u'Jointpar Plot',
81 section='solution',
82 feather_icon='crosshair',
83 description=u'''
84Source problem parameter's scatter plots, to evaluate the resolution of source
85parameters and possible trade-offs between pairs of model parameters.
87A subset of model solutions (from harvest) is shown in two dimensions for all
88possible parameter pairs as points. The point color indicates the misfit for
89the model solution with cold colors (blue) for high misfit models and warm
90colors (red) for low misfit models. The plot ranges are defined by the given
91parameter bounds and shows the model space of the optimsation. %s''' % sref)
93 def draw_figures(self, history, optimiser):
95 color_parameter = self.color_parameter
96 exclude = self.exclude
97 include = self.include
98 nsubplots = self.nsubplots
99 figsize = self.size_inch
100 ibootstrap = 0 if self.ibootstrap is None else self.ibootstrap
101 misfit_cutoff = self.misfit_cutoff
102 show_ellipses = self.show_ellipses
103 msize = 1.5
104 cmap = 'coolwarm'
106 problem = history.problem
107 if not problem:
108 return []
110 models = history.models
112 exclude = list(exclude)
113 bounds = problem.get_combined_bounds()
114 for ipar in range(problem.ncombined):
115 par = problem.combined[ipar]
116 lo, hi = bounds[ipar]
117 if lo == hi:
118 exclude.append(par.name)
120 xref = problem.get_reference_model()
122 isort = history.get_sorted_misfits_idx(chain=ibootstrap)[::-1]
123 models = history.get_sorted_models(chain=ibootstrap)[::-1]
124 nmodels = history.nmodels
126 gms = history.get_sorted_misfits(chain=ibootstrap)[::-1]
127 if misfit_cutoff is not None:
128 ibest = gms < misfit_cutoff
129 gms = gms[ibest]
130 models = models[ibest]
132 kwargs = {}
134 if color_parameter == 'dist':
135 mx = num.mean(models, axis=0)
136 cov = num.cov(models.T)
137 mdists = core.mahalanobis_distance(models, mx, cov)
138 icolor = meta.ordersort(mdists)
140 elif color_parameter == 'misfit':
141 iorder = num.arange(nmodels)
142 icolor = iorder
144 elif color_parameter in problem.parameter_names:
145 ind = problem.name_to_index(color_parameter)
146 icolor = problem.extract(models, ind)
148 elif color_parameter in history.attribute_names:
149 icolor = history.get_attribute(color_parameter)[isort]
150 icolor_need = num.unique(icolor)
152 colors = []
153 for i in range(icolor_need[-1]+1):
154 colors.append(mpl_graph_color(i))
156 cmap = mcolors.ListedColormap(colors)
157 cmap.set_under(mpl_color('aluminium3'))
158 kwargs.update(dict(vmin=0, vmax=icolor_need[-1]))
159 else:
160 raise meta.GrondError(
161 'Invalid color_parameter: %s' % color_parameter)
163 smap = {}
164 iselected = 0
165 for ipar in range(problem.ncombined):
166 par = problem.combined[ipar]
167 if exclude and par.name in exclude or \
168 include and par.name not in include:
169 continue
171 smap[iselected] = ipar
172 iselected += 1
174 nselected = iselected
176 if nselected < 2:
177 logger.warning(
178 'Cannot draw joinpar figures with less than two parameters '
179 'selected.')
180 return []
182 nfig = (nselected - 2) // nsubplots + 1
184 figs = []
185 for ifig in range(nfig):
186 figs_row = []
187 for jfig in range(nfig):
188 if ifig >= jfig:
189 item = PlotItem(name='fig_%i_%i' % (ifig, jfig))
190 item.attributes['parameters'] = []
191 figs_row.append((item, plt.figure(figsize=figsize)))
192 else:
193 figs_row.append(None)
195 figs.append(figs_row)
197 for iselected in range(nselected):
198 ipar = smap[iselected]
199 ypar = problem.combined[ipar]
200 for jselected in range(iselected):
201 jpar = smap[jselected]
202 xpar = problem.combined[jpar]
204 ixg = (iselected - 1)
205 iyg = jselected
207 ix = ixg % nsubplots
208 iy = iyg % nsubplots
210 ifig = ixg // nsubplots
211 jfig = iyg // nsubplots
213 aind = (nsubplots, nsubplots, (ix * nsubplots) + iy + 1)
215 item, fig = figs[ifig][jfig]
217 tlist = item.attributes['parameters']
218 if xpar.name not in tlist:
219 tlist.append(xpar.name)
220 if ypar.name not in tlist:
221 tlist.append(ypar.name)
223 axes = fig.add_subplot(*aind)
225 axes.axvline(0., color=mpl_color('aluminium3'), lw=0.5)
226 axes.axhline(0., color=mpl_color('aluminium3'), lw=0.5)
227 for spine in axes.spines.values():
228 spine.set_edgecolor(mpl_color('aluminium5'))
229 spine.set_linewidth(0.5)
231 xmin, xmax = fixlim(*xpar.scaled(bounds[jpar]))
232 ymin, ymax = fixlim(*ypar.scaled(bounds[ipar]))
234 if not self.show_axes:
235 if ix == 0 or jselected + 1 == iselected:
236 for (xpos, xoff, x) in [
237 (0.0, 10., xmin),
238 (1.0, -10., xmax)]:
240 axes.annotate(
241 '%.3g%s' % (x, xpar.get_unit_suffix()),
242 xy=(xpos, 1.05),
243 xycoords='axes fraction',
244 xytext=(xoff, 5.),
245 textcoords='offset points',
246 verticalalignment='bottom',
247 horizontalalignment='left',
248 rotation=45.)
250 if iy == nsubplots - 1 or jselected + 1 == iselected:
251 for (ypos, yoff, y) in [
252 (0., 10., ymin),
253 (1.0, -10., ymax)]:
255 axes.annotate(
256 '%.3g%s' % (y, ypar.get_unit_suffix()),
257 xy=(1.0, ypos),
258 xycoords='axes fraction',
259 xytext=(5., yoff),
260 textcoords='offset points',
261 verticalalignment='bottom',
262 horizontalalignment='left',
263 rotation=45.)
265 axes.set_xlim(xmin, xmax)
266 axes.set_ylim(ymin, ymax)
268 if not self.show_axes:
269 if not self.show_ticks:
270 axes.get_xaxis().set_ticks([])
271 axes.get_yaxis().set_ticks([])
272 else:
273 axes.tick_params(length=4, which='both')
274 axes.get_yaxis().set_ticklabels([])
275 axes.get_xaxis().set_ticklabels([])
277 if iselected == nselected - 1 or ix == nsubplots - 1:
278 axes.annotate(
279 xpar.get_label(with_unit=False),
280 xy=(0.5, -0.05),
281 xycoords='axes fraction',
282 verticalalignment='top',
283 horizontalalignment='right',
284 rotation=45.)
286 if iy == 0:
287 axes.annotate(
288 ypar.get_label(with_unit=False),
289 xy=(-0.05, 0.5),
290 xycoords='axes fraction',
291 verticalalignment='top',
292 horizontalalignment='right',
293 rotation=45.)
294 else:
295 axes.set_xlabel(xpar.get_label())
296 axes.set_ylabel(ypar.get_label())
298 fx = problem.extract(models, jpar)
299 fy = problem.extract(models, ipar)
301 axes.scatter(
302 xpar.scaled(fx),
303 ypar.scaled(fy),
304 c=icolor,
305 s=msize, alpha=0.5, cmap=cmap, edgecolors='none', **kwargs)
307 if show_ellipses:
308 cov = num.cov((xpar.scaled(fx), ypar.scaled(fy)))
309 evals, evecs = eigh_sorted(cov)
310 evals = num.sqrt(evals)
311 ell = patches.Ellipse(
312 xy=(
313 num.mean(xpar.scaled(fx)),
314 num.mean(ypar.scaled(fy))),
315 width=evals[0] * 2,
316 height=evals[1] * 2,
317 angle=num.rad2deg(
318 num.arctan2(evecs[1][0], evecs[0][0])))
320 ell.set_facecolor('none')
321 ell.set_edgecolor(mpl_color('aluminium5'))
322 axes.add_artist(ell)
324 if self.show_reference:
325 fx = problem.extract(xref, jpar)
326 fy = problem.extract(xref, ipar)
328 ref_color = mpl_color('aluminium6')
329 ref_color_light = 'none'
330 axes.plot(
331 xpar.scaled(fx), ypar.scaled(fy), 's',
332 mew=1.5, ms=5, mfc=ref_color_light, mec=ref_color)
334 figs_flat = []
335 for figs_row in figs:
336 figs_flat.extend(
337 item_fig for item_fig in figs_row if item_fig is not None)
339 return figs_flat
342class HistogramPlot(PlotConfig):
343 '''
344 Histograms or Gaussian kernel densities (default) of all parameters
345 (marginal distributions of model parameters).
347 The histograms (by default shown as Gaussian kernel densities) show (red
348 curved solid line) the distributions of the parameters (marginals) along
349 with some characteristics: The red solid vertical line gives the median of
350 the distribution and the dashed red vertical line the mean value. Dark gray
351 vertical lines show reference values (given in the event.txt file). The
352 overlapping red-shaded areas show the 68% confidence intervals (innermost
353 area), the 90% confidence intervals (middle area) and the minimum and
354 maximum values (widest area). The plot ranges are defined by the given
355 parameter bounds and show the model space. Well resolved model parameters
356 show peaked distributions.
357 '''
359 name = 'histogram'
360 size_cm = Tuple.T(2, Float.T(), default=(12.5, 7.5))
361 exclude = List.T(String.T())
362 include = List.T(String.T())
363 method = StringChoice.T(
364 choices=['gaussian_kde', 'histogram'],
365 default='gaussian_kde')
366 show_reference = Bool.T(default=True)
368 def make(self, environ):
369 cm = environ.get_plot_collection_manager()
370 history = environ.get_history(subset='harvest')
372 mpl_init(fontsize=self.font_size)
373 cm.create_group_mpl(
374 self,
375 self.draw_figures(history),
376 title=u'Histogram',
377 section='solution',
378 feather_icon='bar-chart-2',
379 description=u'''
380Distribution of the problem's parameters.
382The histograms are shown either as Gaussian kernel densities (red curved solid
383line) or as bar plots the distributions of the parameters (marginals) along
384with some characteristics:
386The red solid vertical line gives the median of the distribution and the dashed
387red vertical line the mean value. Dark gray vertical lines show reference
388parameter values if given in the event.txt file. The overlapping red-shaded
389areas show the 68% confidence intervals (innermost area), the 90% confidence
390intervals (middle area) and the minimum and maximum values (widest area). The
391plot ranges are defined by the given parameter bounds and show the model
392space.
393''')
395 def draw_figures(self, history):
397 import scipy.stats
398 from grond.core import make_stats
400 exclude = self.exclude
401 include = self.include
402 figsize = self.size_inch
403 fontsize = self.font_size
404 method = self.method
405 ref_color = mpl_color('aluminium6')
406 stats_color = mpl_color('scarletred2')
407 bar_color = mpl_color('scarletred1')
408 stats_color3 = mpl_color('scarletred3')
410 problem = history.problem
412 models = history.models
414 bounds = problem.get_combined_bounds()
415 exclude = list(exclude)
416 for ipar in range(problem.ncombined):
417 par = problem.combined[ipar]
418 vmin, vmax = bounds[ipar]
419 if vmin == vmax:
420 exclude.append(par.name)
422 xref = problem.get_reference_model()
424 smap = {}
425 iselected = 0
426 for ipar in range(problem.ncombined):
427 par = problem.combined[ipar]
428 if exclude and par.name in exclude or \
429 include and par.name not in include:
430 continue
432 smap[iselected] = ipar
433 iselected += 1
435 nselected = iselected
436 del iselected
438 pnames = [
439 problem.combined[smap[iselected]].name
440 for iselected in range(nselected)]
442 rstats = make_stats(problem, models,
443 history.get_primary_chain_misfits(),
444 pnames=pnames)
446 for iselected in range(nselected):
447 ipar = smap[iselected]
448 par = problem.combined[ipar]
449 vs = problem.extract(models, ipar)
450 vmin, vmax = bounds[ipar]
452 fig = plt.figure(figsize=figsize)
453 labelpos = mpl_margins(
454 fig, nw=1, nh=1, w=7., bottom=5., top=1, units=fontsize)
456 axes = fig.add_subplot(1, 1, 1)
457 labelpos(axes, 2.5, 2.0)
458 axes.set_xlabel(par.get_label())
459 axes.set_ylabel('PDF')
460 axes.set_xlim(*fixlim(*par.scaled((vmin, vmax))))
462 if method == 'gaussian_kde':
463 try:
464 kde = scipy.stats.gaussian_kde(vs)
465 except Exception:
466 logger.warning(
467 'Cannot create plot histogram with gaussian_kde: '
468 'possibly all samples have the same value.')
469 continue
471 vps = num.linspace(vmin, vmax, 600)
472 pps = kde(vps)
474 axes.plot(
475 par.scaled(vps), par.inv_scaled(pps), color=stats_color)
477 elif method == 'histogram':
478 pps, edges = num.histogram(
479 vs,
480 bins=num.linspace(vmin, vmax, num=40),
481 density=True)
482 vps = 0.5 * (edges[:-1] + edges[1:])
484 axes.bar(par.scaled(vps), par.inv_scaled(pps),
485 par.scaled(2.*(vps - edges[:-1])),
486 color=bar_color)
488 pstats = rstats.parameter_stats_list[iselected]
490 axes.axvspan(
491 par.scaled(pstats.minimum),
492 par.scaled(pstats.maximum),
493 color=stats_color, alpha=0.1)
494 axes.axvspan(
495 par.scaled(pstats.percentile16),
496 par.scaled(pstats.percentile84),
497 color=stats_color, alpha=0.1)
498 axes.axvspan(
499 par.scaled(pstats.percentile5),
500 par.scaled(pstats.percentile95),
501 color=stats_color, alpha=0.1)
503 axes.axvline(
504 par.scaled(pstats.median),
505 color=stats_color3, alpha=0.5)
506 axes.axvline(
507 par.scaled(pstats.mean),
508 color=stats_color3, ls=':', alpha=0.5)
510 if self.show_reference:
511 axes.axvline(
512 par.scaled(problem.extract(xref, ipar)),
513 color=ref_color)
515 item = PlotItem(name=par.name)
516 item.attributes['parameters'] = [par.name]
517 yield item, fig
520class MTDecompositionPlot(PlotConfig):
521 '''
522 Moment tensor decomposition plot.
523 '''
525 name = 'mt_decomposition'
526 size_cm = Tuple.T(2, Float.T(), default=(15., 5.))
527 cluster_attribute = meta.StringID.T(
528 optional=True,
529 help='name of attribute to use as cluster IDs')
530 show_reference = Bool.T(default=True)
532 def make(self, environ):
533 cm = environ.get_plot_collection_manager()
534 history = environ.get_history(subset='harvest')
535 mpl_init(fontsize=self.font_size)
536 cm.create_group_mpl(
537 self,
538 self.draw_figures(history),
539 title=u'MT Decomposition',
540 section='solution',
541 feather_icon='sun',
542 description=u'''
543Moment tensor decomposition of the best-fitting solution into isotropic,
544deviatoric and best double couple components.
546Shown are the ensemble best, the ensemble mean%s and, if available, a reference
547mechanism. The symbol size indicates the relative strength of the components.
548The inversion result is consistent and stable if ensemble mean and ensemble
549best have similar symbol size and patterns.
550''' % (', cluster results' if self.cluster_attribute else ''))
552 def draw_figures(self, history):
554 fontsize = self.font_size
556 fig = plt.figure(figsize=self.size_inch)
557 axes = fig.add_subplot(1, 1, 1, aspect=1.0)
558 fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
560 problem = history.problem
561 models = history.models
563 if models.size == 0:
564 logger.warning('Empty models vector.')
565 return []
567 # gms = problem.combine_misfits(history.misfits)
568 # isort = num.argsort(gms)
569 # iorder = num.empty_like(isort)
570 # iorder[isort] = num.arange(iorder.size)[::-1]
572 ref_source = problem.base_source
574 mean_source = stats.get_mean_source(
575 problem, history.models)
577 best_source = history.get_best_source()
579 nlines_max = int(round(self.size_cm[1] / 5. * 4. - 1.0))
581 if self.cluster_attribute:
582 cluster_sources = history.mean_sources_by_cluster(
583 self.cluster_attribute)
584 else:
585 cluster_sources = []
587 def get_deco(source):
588 mt = source.pyrocko_moment_tensor()
589 return mt.standard_decomposition()
591 lines = []
592 lines.append(
593 ('Ensemble best', get_deco(best_source), mpl_color('aluminium5')))
595 lines.append(
596 ('Ensemble mean', get_deco(mean_source), mpl_color('aluminium5')))
598 for (icluster, perc, source) in cluster_sources:
599 if len(lines) < nlines_max - int(self.show_reference):
600 lines.append(
601 (cluster_label(icluster, perc),
602 get_deco(source),
603 cluster_color(icluster)))
604 else:
605 logger.warning(
606 'Skipping display of cluster %i because figure height is '
607 'too small. Figure height should be at least %g cm.' % (
608 icluster, (3 + len(cluster_sources)
609 + int(self.show_reference)) * 5/4.))
611 if self.show_reference:
612 lines.append(
613 ('Reference', get_deco(ref_source), mpl_color('aluminium3')))
615 moment_full_max = max(deco[-1][0] for (_, deco, _) in lines)
617 for xpos, label in [
618 (0., 'Full'),
619 (2., 'Isotropic'),
620 (4., 'Deviatoric'),
621 (6., 'CLVD'),
622 (8., 'DC')]:
624 axes.annotate(
625 label,
626 xy=(1 + xpos, nlines_max),
627 xycoords='data',
628 xytext=(0., 0.),
629 textcoords='offset points',
630 ha='center',
631 va='center',
632 color='black',
633 fontsize=fontsize)
635 for i, (label, deco, color_t) in enumerate(lines):
636 ypos = nlines_max - i - 1.0
638 [(moment_iso, ratio_iso, m_iso),
639 (moment_dc, ratio_dc, m_dc),
640 (moment_clvd, ratio_clvd, m_clvd),
641 (moment_devi, ratio_devi, m_devi),
642 (moment_full, ratio_full, m_full)] = deco
644 size0 = moment_full / moment_full_max
646 axes.annotate(
647 label,
648 xy=(-2., ypos),
649 xycoords='data',
650 xytext=(0., 0.),
651 textcoords='offset points',
652 ha='left',
653 va='center',
654 color='black',
655 fontsize=fontsize)
657 for xpos, mt_part, ratio, ops in [
658 (0., m_full, ratio_full, '-'),
659 (2., m_iso, ratio_iso, '='),
660 (4., m_devi, ratio_devi, '='),
661 (6., m_clvd, ratio_clvd, '+'),
662 (8., m_dc, ratio_dc, None)]:
664 if ratio > 1e-4:
665 try:
666 beachball.plot_beachball_mpl(
667 mt_part, axes,
668 beachball_type='full',
669 position=(1. + xpos, ypos),
670 size=0.9 * size0 * math.sqrt(ratio),
671 size_units='data',
672 color_t=color_t,
673 linewidth=1.0)
675 except beachball.BeachballError as e:
676 logger.warning(str(e))
678 axes.annotate(
679 'ERROR',
680 xy=(1. + xpos, ypos),
681 ha='center',
682 va='center',
683 color='red',
684 fontsize=fontsize)
686 else:
687 axes.annotate(
688 'N/A',
689 xy=(1. + xpos, ypos),
690 ha='center',
691 va='center',
692 color='black',
693 fontsize=fontsize)
695 if ops is not None:
696 axes.annotate(
697 ops,
698 xy=(2. + xpos, ypos),
699 ha='center',
700 va='center',
701 color='black',
702 fontsize=fontsize)
704 axes.axison = False
705 axes.set_xlim(-2.25, 9.75)
706 axes.set_ylim(-0.5, nlines_max+0.5)
708 item = PlotItem(name='main')
709 return [[item, fig]]
712class MTLocationPlot(SectionPlotConfig):
713 ''' MT location plot of the best solutions in three cross-sections. '''
714 name = 'location_mt'
715 beachball_type = StringChoice.T(
716 choices=['full', 'deviatoric', 'dc'],
717 default='dc')
718 normalisation_gamma = Float.T(
719 default=3.,
720 help='Normalisation of colors and alpha as :math:`x^\\gamma`.'
721 'A linear colormap/alpha with :math:`\\gamma=1`.')
723 def make(self, environ):
724 environ.setup_modelling()
725 cm = environ.get_plot_collection_manager()
726 history = environ.get_history(subset='harvest')
727 mpl_init(fontsize=self.font_size)
728 self._to_be_closed = []
729 cm.create_group_mpl(
730 self,
731 self.draw_figures(history),
732 title=u'MT Location',
733 section='solution',
734 feather_icon='target',
735 description=u'''
736Location plot of the ensemble of best solutions in three cross-sections.
738The coordinate range is defined by the search space given in the config file.
739Symbols show best double-couple mechanisms, and colors indicate low (red) and
740high (blue) misfit.
741''')
742 for obj in self._to_be_closed:
743 obj.close()
745 def draw_figures(self, history, color_p_axis=False):
746 from matplotlib import colors
748 color = 'black'
749 fontsize = self.font_size
750 markersize = fontsize * 1.5
751 beachballsize_small = markersize * 0.5
752 beachball_type = self.beachball_type
754 problem = history.problem
755 sp = SectionPlot(config=self)
756 self._to_be_closed.append(sp)
758 fig = sp.fig
759 axes_en = sp.axes_xy
760 axes_dn = sp.axes_zy
761 axes_ed = sp.axes_xz
763 bounds = problem.get_combined_bounds()
765 models = history.get_sorted_primary_models()[::-1]
767 iorder = num.arange(history.nmodels)
769 for parname, set_label, set_lim in [
770 ['east_shift', sp.set_xlabel, sp.set_xlim],
771 ['north_shift', sp.set_ylabel, sp.set_ylim],
772 ['depth', sp.set_zlabel, sp.set_zlim]]:
774 ipar = problem.name_to_index(parname)
775 par = problem.combined[ipar]
776 set_label(par.get_label())
777 xmin, xmax = fixlim(*par.scaled(bounds[ipar]))
778 set_lim(xmin, xmax)
780 if 'volume_change' in problem.parameter_names:
781 volumes = models[:, problem.name_to_index('volume_change')]
782 volume_max = volumes.max()
783 volume_min = volumes.min()
785 def scale_size(source):
786 if not hasattr(source, 'volume_change'):
787 return beachballsize_small
789 volume_change = source.volume_change
790 fac = (volume_change - volume_min) / (volume_max - volume_min)
791 return markersize * .25 + markersize * .5 * fac
793 for axes, xparname, yparname in [
794 (axes_en, 'east_shift', 'north_shift'),
795 (axes_dn, 'depth', 'north_shift'),
796 (axes_ed, 'east_shift', 'depth')]:
798 ixpar = problem.name_to_index(xparname)
799 iypar = problem.name_to_index(yparname)
801 xpar = problem.combined[ixpar]
802 ypar = problem.combined[iypar]
804 xmin, xmax = fixlim(*xpar.scaled(bounds[ixpar]))
805 ymin, ymax = fixlim(*ypar.scaled(bounds[iypar]))
807 try:
808 axes.set_facecolor(mpl_color('aluminium1'))
809 except AttributeError:
810 axes.patch.set_facecolor(mpl_color('aluminium1'))
812 rect = patches.Rectangle(
813 (xmin, ymin), xmax-xmin, ymax-ymin,
814 facecolor=mpl_color('white'),
815 edgecolor=mpl_color('aluminium2'))
817 axes.add_patch(rect)
819 # fxs = xpar.scaled(problem.extract(models, ixpar))
820 # fys = ypar.scaled(problem.extract(models, iypar))
822 # axes.set_xlim(*fixlim(num.min(fxs), num.max(fxs)))
823 # axes.set_ylim(*fixlim(num.min(fys), num.max(fys)))
825 cmap = cm.ScalarMappable(
826 norm=colors.PowerNorm(
827 gamma=self.normalisation_gamma,
828 vmin=iorder.min(),
829 vmax=iorder.max()),
831 cmap=plt.get_cmap('coolwarm'))
833 for ix, x in enumerate(models):
835 source = problem.get_source(x)
836 mt = source.pyrocko_moment_tensor(
837 store=problem.get_gf_store(problem.targets[0]),
838 target=problem.targets[0])
839 fx = problem.extract(x, ixpar)
840 fy = problem.extract(x, iypar)
841 sx, sy = xpar.scaled(fx), ypar.scaled(fy)
843 # TODO: Add rotation in cross-sections
844 color = cmap.to_rgba(iorder[ix])
846 alpha = (iorder[ix] - iorder.min()) / \
847 float(iorder.max() - iorder.min())
848 alpha = alpha**self.normalisation_gamma
850 try:
851 beachball.plot_beachball_mpl(
852 mt, axes,
853 beachball_type=beachball_type,
854 position=(sx, sy),
855 size=scale_size(source),
856 color_t=color,
857 color_p=color if color_p_axis else 'white',
858 alpha=alpha,
859 zorder=1,
860 linewidth=0.25)
862 except beachball.BeachballError as e:
863 logger.warning(str(e))
865 item = PlotItem(name='main')
866 return [[item, fig]]
869class MTFuzzyPlot(PlotConfig):
870 '''Fuzzy, propabalistic moment tensor plot '''
872 name = 'mt_fuzzy'
873 size_cm = Tuple.T(2, Float.T(), default=(10., 10.))
874 cluster_attribute = meta.StringID.T(
875 optional=True,
876 help='name of attribute to use as cluster IDs')
878 def make(self, environ):
879 cm = environ.get_plot_collection_manager()
880 history = environ.get_history(subset='harvest')
881 mpl_init(fontsize=self.font_size)
882 cm.create_group_mpl(
883 self,
884 self.draw_figures(history),
885 title=u'Fuzzy MT',
886 section='solution',
887 feather_icon='wind',
888 description=u'''
889A fuzzy moment tensor, illustrating the solution's uncertainty.
891The P wave radiation pattern strength of every ensemble solution is stacked for
892all ray spokes. The projection shows the stacked radiation pattern. If the
893variability of the ensemble solutions is small, the fuzzy plot has clearly
894separated black and white fields, consistent with the nodal lines of the %s
895best solution (indicated in red).
896''' % ('cluster' if self.cluster_attribute is not None else 'global'))
898 def draw_figures(self, history):
899 problem = history.problem
901 by_cluster = history.imodels_by_cluster(
902 self.cluster_attribute)
904 for icluster, percentage, imodels in by_cluster:
905 misfits = history.misfits[imodels]
906 models = history.models[imodels]
908 mts = []
909 for ix, x in enumerate(models):
910 source = problem.get_source(x)
911 mts.append(source.pyrocko_moment_tensor())
913 best_mt = stats.get_best_source(
914 problem, models, misfits).pyrocko_moment_tensor()
916 fig = plt.figure(figsize=self.size_inch)
917 fig.subplots_adjust(left=0., right=1., bottom=0., top=1.)
918 axes = fig.add_subplot(1, 1, 1, aspect=1.0)
920 if self.cluster_attribute is not None:
921 color = cluster_color(icluster)
922 else:
923 color = 'black'
925 beachball.plot_fuzzy_beachball_mpl_pixmap(
926 mts, axes, best_mt,
927 beachball_type='full',
928 size=8.*math.sqrt(percentage/100.),
929 position=(5., 5.),
930 color_t=color,
931 edgecolor='black',
932 best_color=mpl_color('scarletred2'))
934 if self.cluster_attribute is not None:
935 axes.annotate(
936 cluster_label(icluster, percentage),
937 xy=(5., 0.),
938 xycoords='data',
939 xytext=(0., self.font_size/2.),
940 textcoords='offset points',
941 ha='center',
942 va='bottom',
943 color='black',
944 fontsize=self.font_size)
946 axes.set_xlim(0., 10.)
947 axes.set_ylim(0., 10.)
948 axes.set_axis_off()
950 item = PlotItem(
951 name=(
952 'cluster_%i' % icluster
953 if icluster >= 0
954 else 'unclustered'))
956 yield [item, fig]
959class HudsonPlot(PlotConfig):
961 '''
962 Illustration of the solution distribution of decomposed moment tensor.
963 '''
965 name = 'hudson'
966 size_cm = Tuple.T(2, Float.T(), default=(17.5, 17.5*(3./4.)))
967 beachball_type = StringChoice.T(
968 choices=['full', 'deviatoric', 'dc'],
969 default='dc')
970 show_reference = Bool.T(default=True)
972 def make(self, environ):
973 cm = environ.get_plot_collection_manager()
974 history = environ.get_history(subset='harvest')
975 mpl_init(fontsize=self.font_size)
976 cm.create_group_mpl(
977 self,
978 self.draw_figures(history),
979 title=u'Hudson Plot',
980 section='solution',
981 feather_icon='box',
982 description=u'''
983Hudson's source type plot with the ensemble of bootstrap solutions.
985For about 10% of the solutions (randomly chosen), the focal mechanism is
986depicted, others are represented as dots. The square marks the global best
987fitting solution.
988''')
990 def draw_figures(self, history):
992 color = 'black'
993 fontsize = self.font_size
994 markersize = fontsize * 1.5
995 markersize_small = markersize * 0.2
996 beachballsize = markersize
997 beachballsize_small = beachballsize * 0.5
998 beachball_type = self.beachball_type
1000 problem = history.problem
1001 best_source = history.get_best_source()
1002 mean_source = history.get_mean_source()
1004 fig = plt.figure(figsize=self.size_inch)
1005 axes = fig.add_subplot(1, 1, 1)
1007 data = []
1008 for ix, x in enumerate(history.models):
1009 source = problem.get_source(x)
1010 mt = source.pyrocko_moment_tensor()
1011 u, v = hudson.project(mt)
1013 if random.random() < 0.1:
1014 try:
1015 beachball.plot_beachball_mpl(
1016 mt, axes,
1017 beachball_type=beachball_type,
1018 position=(u, v),
1019 size=beachballsize_small,
1020 color_t=color,
1021 alpha=0.5,
1022 zorder=1,
1023 linewidth=0.25)
1024 except beachball.BeachballError as e:
1025 logger.warning(str(e))
1027 else:
1028 data.append((u, v))
1030 if data:
1031 u, v = num.array(data).T
1032 axes.plot(
1033 u, v, 'o',
1034 color=color,
1035 ms=markersize_small,
1036 mec='none',
1037 mew=0,
1038 alpha=0.25,
1039 zorder=0)
1041 hudson.draw_axes(axes)
1043 mt = mean_source.pyrocko_moment_tensor()
1044 u, v = hudson.project(mt)
1046 try:
1047 beachball.plot_beachball_mpl(
1048 mt, axes,
1049 beachball_type=beachball_type,
1050 position=(u, v),
1051 size=beachballsize,
1052 color_t=color,
1053 zorder=2,
1054 linewidth=0.5)
1055 except beachball.BeachballError as e:
1056 logger.warning(str(e))
1058 mt = best_source.pyrocko_moment_tensor()
1059 u, v = hudson.project(mt)
1061 axes.plot(
1062 u, v, 's',
1063 markersize=markersize,
1064 mew=1,
1065 mec='black',
1066 mfc='none',
1067 zorder=-2)
1069 if self.show_reference:
1070 mt = problem.base_source.pyrocko_moment_tensor()
1071 u, v = hudson.project(mt)
1073 try:
1074 beachball.plot_beachball_mpl(
1075 mt, axes,
1076 beachball_type=beachball_type,
1077 position=(u, v),
1078 size=beachballsize,
1079 color_t='red',
1080 zorder=2,
1081 linewidth=0.5)
1082 except beachball.BeachballError as e:
1083 logger.warning(str(e))
1085 item = PlotItem(
1086 name='main')
1087 return [[item, fig]]
1090def get_plot_classes():
1091 return [
1092 JointparPlot,
1093 HistogramPlot,
1094 ]