Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/dynamic_rupture.py: 82%
648 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-04 09:52 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from subprocess import check_call, CalledProcessError
8import os
9import os.path as op
10import re
11import logging
12import tempfile
13import shutil
14import inspect
16import numpy as num
17from scipy.interpolate import RegularGridInterpolator as scrgi
19from matplotlib import pyplot as plt, patheffects
20from matplotlib.ticker import FuncFormatter
22from pyrocko import orthodrome as pod
23from pyrocko.guts import Object
24from pyrocko.plot import mpl_init, mpl_papersize, mpl_color, AutoScaler, \
25 gmtpy, mpl_get_cmap
26from pyrocko.plot.automap import Map, NoTopo
27from pyrocko.gf import PseudoDynamicRupture
28from pyrocko.gf.seismosizer import map_anchor
29from pyrocko.dataset.topo.tile import Tile
31logger = logging.getLogger(__name__)
33km = 1e3
35d2m = 111180.
36m2d = 1. / d2m
38cm2inch = gmtpy.cm/gmtpy.inch
40d2r = num.pi / 180.
41r2d = 1. / d2r
44def km_formatter(v, p):
45 return '%g' % (v / km)
48def get_kwargs(cls):
49 sig = inspect.signature(cls.__init__)
50 kwargs = {}
52 if issubclass(cls.T._cls, RuptureMap):
53 for p in Map.T.properties:
54 kwargs[p.name] = p._default
56 for par in sig.parameters.values():
57 if par.default is inspect._empty:
58 continue
59 kwargs[par.name] = par.default
61 return kwargs
64def _save_grid(lats, lons, data, filename):
65 '''
66 Save lat-lon gridded data as gmt .grd file.
68 :param lats:
69 Grid latitude coordinates in [deg].
70 :type lats:
71 iterable
73 :param lons:
74 Grid longitude coordinates in [deg].
75 :type lons:
76 iterable
78 :param data:
79 Grid data of any kind.
80 :type data:
81 :py:class:`~numpy.ndarray`, ``(n_lons, n_lats)``
83 :param filename:
84 Filename of the written grid file.
85 :type filename:
86 str
87 '''
89 gmtpy.savegrd(lons, lats, data, filename=filename, naming='lonlat')
92def _mplcmap_to_gmtcpt_code(mplcmap, steps=256):
93 '''
94 Get gmt readable R/G/A code from a given matplotlib colormap.
96 :param mplcmap:
97 Name of the demanded matplotlib colormap.
98 :type mplcmap:
99 str
101 :returns:
102 Series of comma seperate R/G/B values for gmtpy usage.
103 :rtype:
104 str
105 '''
107 cmap = mpl_get_cmap(mplcmap)
109 rgbas = [cmap(i) for i in num.linspace(0, 255, steps).astype(num.int64)]
111 return ','.join(['%d/%d/%d' % (
112 c[0] * 255, c[1] * 255, c[2] * 255) for c in rgbas])
115def make_colormap(
116 gmt,
117 vmin,
118 vmax,
119 C=None,
120 cmap=None,
121 space=False):
122 '''
123 Create gmt-readable colormap cpt file called my_<cmap>.cpt.
125 :param vmin:
126 Minimum value covered by the colormap.
127 :type vmin:
128 float
130 :param vmax:
131 Maximum value covered by the colormap.
132 :type vmax:
133 float
135 :param C:
136 Comma seperated R/G/B values for cmap definition.
137 :type C:
138 str
140 :param cmap:
141 Name of the colormap. Colormap is stored as "my_<cmap>.cpt".
142 If name is equivalent to a matplotlib colormap, R/G/B strings are
143 extracted from this colormap.
144 :type cmap:
145 str
147 :param space:
148 If ``True``, the range of the colormap is broadened below vmin
149 and above vmax.
150 :type space: bool
151 '''
153 scaler = AutoScaler(mode='min-max')
154 scale = scaler.make_scale((vmin, vmax))
156 incr = scale[2]
157 margin = 0.
159 if vmin == vmax:
160 space = True
162 if space:
163 margin = incr
165 msg = ('Please give either a valid color code or a'
166 ' valid matplotlib colormap name.')
168 if C is None and cmap is None:
169 raise ValueError(msg)
171 if C is None:
172 try:
173 C = _mplcmap_to_gmtcpt_code(cmap)
174 except ValueError:
175 raise ValueError(msg)
177 if cmap is None:
178 logger.warning(
179 'No colormap name given. Uses temporary filename instead')
180 cmap = 'temp_cmap'
182 return gmt.makecpt(
183 C=C,
184 D='o',
185 T='%g/%g/%g' % (
186 vmin - margin, vmax + margin, incr),
187 Z=True,
188 out_filename='my_%s.cpt' % cmap,
189 suppress_defaults=True)
192def clear_temp(gridfiles=[], cpts=[]):
193 '''
194 Clear all temporary needed grid and colormap cpt files.
196 :param gridfiles:
197 List of all "...grd" files, which shall be deleted.
198 :type gridfiles:
199 list
201 :param cpts:
202 Cmaps, whose corresponding "my_<cmap>.cpt" file shall be deleted.
203 :type cpts:
204 list
205 '''
207 for fil in gridfiles:
208 try:
209 os.remove(fil)
210 except OSError:
211 continue
212 for fil in cpts:
213 try:
214 os.remove('my_%s.cpt' % fil)
215 except OSError:
216 continue
219def xy_to_latlon(source, x, y):
220 '''
221 Convert x and y relative coordinates on extended ruptures into latlon.
223 :param source:
224 Extended source class, on which the given point is located.
225 :type source:
226 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or
227 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
229 :param x:
230 Relative point coordinate along strike (range: -1:1).
231 :type x:
232 :py:class:`float` or :py:class:`~numpy.ndarray`
234 :param y:
235 Relative downdip point coordinate (range: -1:1).
236 :type y:
237 :py:class:`float` or :py:class:`~numpy.ndarray`
239 :returns:
240 Latitude and Longitude of the given point in [deg].
241 :rtype:
242 :py:class:`tuple` of :py:class:`float`
243 '''
245 s = source
246 ax, ay = map_anchor[s.anchor]
248 length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
249 strike, dip = s.strike * d2r, s.dip * d2r
251 northing = num.cos(strike)*length - num.cos(dip) * num.sin(strike)*width
252 easting = num.sin(strike)*length + num.cos(dip) * num.cos(strike)*width
254 northing += source.north_shift
255 easting += source.east_shift
257 return pod.ne_to_latlon(s.lat, s.lon, northing, easting)
260def xy_to_lw(source, x, y):
261 '''
262 Convert relative coordinates on extended ruptures into length and width.
264 :param source:
265 Extended source, on which the given points are located.
266 :type source:
267 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or
268 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
270 :param x:
271 Relative point coordinates along strike (range: -1:1).
272 :type x:
273 float or :py:class:`~numpy.ndarray`
275 :param y:
276 Relative downdip point coordinates (range: -1:1).
277 :type y:
278 float or :py:class:`~numpy.ndarray`
280 :returns:
281 Length and downdip width of the given points from the anchor in [m].
282 :rtype:
283 :py:class:`tuple` of :py:class:`float`
284 '''
286 length, width = source.length, source.width
288 ax, ay = map_anchor[source.anchor]
290 lengths = (x - ax) / 2. * length
291 widths = (y - ay) / 2. * width
293 return lengths, widths
296cbar_anchor = {
297 'center': 'MC',
298 'center_left': 'ML',
299 'center_right': 'MR',
300 'top': 'TC',
301 'top_left': 'TL',
302 'top_right': 'TR',
303 'bottom': 'BC',
304 'bottom_left': 'BL',
305 'bottom_right': 'BR'}
308cbar_helper = {
309 'traction': {
310 'unit': 'MPa',
311 'factor': 1e-6},
312 'tx': {
313 'unit': 'MPa',
314 'factor': 1e-6},
315 'ty': {
316 'unit': 'MPa',
317 'factor': 1e-6},
318 'tz': {
319 'unit': 'MPa',
320 'factor': 1e-6},
321 'time': {
322 'unit': 's',
323 'factor': 1.},
324 'strike': {
325 'unit': '°',
326 'factor': 1.},
327 'dip': {
328 'unit': '°',
329 'factor': 1.},
330 'vr': {
331 'unit': 'km/s',
332 'factor': 1e-3},
333 'length': {
334 'unit': 'km',
335 'factor': 1e-3},
336 'width': {
337 'unit': 'km',
338 'factor': 1e-3}
339}
342fonttype = 'Helvetica'
344c2disl = dict([('x', 0), ('y', 1), ('z', 2)])
347def _make_gmt_conf(fontcolor, size):
348 '''
349 Update different gmt parameters depending on figure size and fontcolor.
351 :param fontcolor:
352 GMT readable colorcode or colorstring for the font.
353 :type fontcolor:
354 str
356 :param size:
357 Tuple of the figure size (width, height) in [cm].
358 :type size:
359 :py:class:`tuple` of :py:class:`float`
361 :returns:
362 Best fitting fontsize, GMT configuration parameter
363 :rtype:
364 float, dict
365 '''
367 color = fontcolor
368 fontsize = num.max(size)
370 font = '%gp,%s' % (fontsize, fonttype)
372 pen_size = fontsize / 16.
373 tick_size = num.min(size) / 200.
375 return fontsize, dict(
376 MAP_FRAME_PEN='%s' % color,
377 MAP_TICK_PEN_PRIMARY='%gp,%s' % (pen_size, color),
378 MAP_TICK_PEN_SECONDARY='%gp,%s' % (pen_size, color),
379 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size,
380 MAP_TICK_LENGTH_SECONDARY='%gc' % (tick_size * 3),
381 FONT_ANNOT_PRIMARY='%s-Bold,%s' % (font, color),
382 FONT_LABEL='%s-Bold,%s' % (font, color),
383 FONT_TITLE='%s-Bold,%s' % (font, color),
384 PS_CHAR_ENCODING='ISOLatin1+',
385 MAP_FRAME_TYPE='fancy',
386 FORMAT_GEO_MAP='D',
387 PS_PAGE_ORIENTATION='portrait',
388 MAP_GRID_PEN_PRIMARY='thinnest,%s' % color,
389 MAP_ANNOT_OBLIQUE='6')
392class SourceError(Exception):
393 pass
396class RuptureMap(Map):
397 '''
398 Map plotting of attributes and results of the
399 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
400 '''
402 def __init__(
403 self,
404 source=None,
405 fontcolor='darkslategrey',
406 width=20.,
407 height=14.,
408 margins=None,
409 color_wet=(216, 242, 254),
410 color_dry=(238, 236, 230),
411 topo_cpt_wet='light_sea_uniform',
412 topo_cpt_dry='light_land_uniform',
413 show_cities=False,
414 *args,
415 **kwargs):
417 size = (width, height)
418 fontsize, gmt_config = _make_gmt_conf(fontcolor, size)
420 if margins is None:
421 margins = [
422 fontsize * 0.15, num.min(size) / 200.,
423 num.min(size) / 200., fontsize * 0.05]
425 Map.__init__(self, *args, margins=margins, width=width, height=height,
426 gmt_config=gmt_config,
427 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet,
428 color_wet=color_wet, color_dry=color_dry,
429 **kwargs)
431 if show_cities:
432 self.draw_cities()
434 self._source = source
435 self._fontcolor = fontcolor
436 self._fontsize = fontsize
437 self.axes_layout = 'WeSn'
439 @property
440 def size(self):
441 '''
442 Figure size in [cm].
443 '''
444 return (self.width, self.height)
446 @property
447 def font(self):
448 '''
449 Font style (size and type).
450 '''
452 return '%sp,%s' % (self._fontsize, fonttype)
454 @property
455 def source(self):
456 '''
457 PseudoDynamicRupture whose attributes are plotted.
459 Note, that source.patches attribute needs to be calculated in advance.
460 '''
462 if self._source is None:
463 raise SourceError('No source given. Please define it!')
465 if not isinstance(self._source, PseudoDynamicRupture):
466 raise SourceError('This function is only capable for a source of'
467 ' type: %s' % type(PseudoDynamicRupture()))
469 if not self._source.patches:
470 raise TypeError('No source patches are defined. Please run'
471 '"source.discretize_patches()" on your source')
473 return self._source
475 @source.setter
476 def source(self, source):
477 self._source = source
479 def _get_topotile(self):
480 if self._dems is None:
481 self._setup()
483 try:
484 t, _ = self._get_topo_tile('land')
485 except NoTopo:
486 wesn = self.wesn
488 nx = int(num.ceil(
489 self.width * cm2inch * self.topo_resolution_max))
490 ny = int(num.ceil(
491 self.height * cm2inch * self.topo_resolution_max))
493 data = num.zeros((nx, ny))
495 t = Tile(wesn[0], wesn[2],
496 (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny,
497 data)
499 return t
501 def _patches_to_lw(self):
502 '''
503 Generate regular rect. length-width grid based on the patch distance.
505 Prerequesite is a regular grid of patches (constant lengths and widths)
506 Both coordinates are given relative to the source anchor point in [m]
507 The grid is extended from the patch centres to the edges of the source.
509 :returns:
510 Lengths along strike, widths downdip in [m].
511 :rtype:
512 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray`
513 '''
515 source = self.source
516 patches = source.patches
518 patch_l, patch_w = patches[0].length, patches[0].width
520 patch_lengths = num.concatenate((
521 num.array([0.]),
522 num.array([il*patch_l+patch_l/2. for il in range(source.nx)]),
523 num.array([patch_l * source.nx])))
525 patch_widths = num.concatenate((
526 num.array([0.]),
527 num.array([iw*patch_w+patch_w/2. for iw in range(source.ny)]),
528 num.array([patch_w * source.ny])))
530 ax, ay = map_anchor[source.anchor]
532 patch_lengths -= source.length * (ax + 1.) / 2.
533 patch_widths -= source.width * (ay + 1.) / 2.
535 return patch_lengths, patch_widths
537 def _xy_to_lw(self, x, y):
538 '''
539 Generate regular rect. length-width grid based on the xy coordinates.
541 Prerequesite is a regular grid with constant dx and dy. x and y are
542 relative coordinates on the rupture plane (range -1:1) along strike and
543 downdip.
544 Length and width are obtained relative to the source anchor point
545 in [m].
547 :returns:
548 Lengths along strike, widths downdip in [m].
549 :rtype:
550 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray`
551 '''
553 x, y = num.unique(x), num.unique(y)
554 dx, dy = x[1] - x[0], y[1] - y[0]
556 if any(num.abs(num.diff(x) - dx) >= 1e-6):
557 raise ValueError('Regular grid with constant spacing needed.'
558 ' Please check the x coordinates.')
560 if any(num.abs(num.diff(y) - dy) >= 1e-6):
561 raise ValueError('Regular grid with constant spacing needed.'
562 ' Please check the y coordinates.')
564 return xy_to_lw(self.source, x, y)
566 def _tile_to_lw(self, ref_lat, ref_lon,
567 north_shift=0., east_shift=0., strike=0.):
569 '''
570 Coordinate transformation from the topo tile grid into length-width.
572 The topotile latlon grid is rotated into the length width grid. The
573 width is defined here as its horizontal component. The resulting grid
574 is used for interpolation of grid data.
576 :param ref_lat:
577 Reference latitude, from which length-width relative coordinates
578 grid are calculated.
579 :type ref_lat:
580 float
582 :param ref_lon:
583 Reference longitude, from which length-width relative coordinates
584 grid are calculated.
585 :type ref_lon:
586 float
588 :param north_shift:
589 North shift of the reference point in [m].
590 :type north_shift:
591 float
593 :param east_shift:
594 East shift of the reference point [m].
595 :type east_shift:
596 float
598 :param strike:
599 striking of the length axis compared to the North axis in [deg].
600 :type strike:
601 float
603 :returns:
604 Topotile grid nodes as array of length-width coordinates.
605 :rtype:
606 :py:class:`~numpy.ndarray`, ``(n_nodes, 2)``
607 '''
609 t = self._get_topotile()
610 grid_lats = t.y()
611 grid_lons = t.x()
613 meshgrid_lons, meshgrid_lats = num.meshgrid(grid_lons, grid_lats)
615 grid_northing, grid_easting = pod.latlon_to_ne_numpy(
616 ref_lat, ref_lon, meshgrid_lats.flatten(), meshgrid_lons.flatten())
618 grid_northing -= north_shift
619 grid_easting -= east_shift
621 strike *= d2r
622 sin, cos = num.sin(strike), num.cos(strike)
623 points_out = num.zeros((grid_northing.shape[0], 2))
624 points_out[:, 0] = -sin * grid_northing + cos * grid_easting
625 points_out[:, 1] = cos * grid_northing + sin * grid_easting
627 return points_out
629 def _prep_patch_grid_data(self, data):
630 '''
631 Extend patch data from patch centres to the outer source edges.
633 Patch data is always defined in the centre of the patches. For
634 interpolation the data is extended here to the edges of the rupture
635 plane.
637 :param data:
638 Patch wise data.
639 :type data:
640 :py:class:`~numpy.ndarray`
642 :returns:
643 Extended data array.
644 :rtype:
645 :py:class:`~numpy.ndarray`
646 '''
648 shape = (self.source.nx + 2, self.source.ny + 2)
649 data = data.reshape(self.source.nx, self.source.ny).copy()
651 data_new = num.zeros(shape)
652 data_new[1:-1, 1:-1] = data
653 data_new[1:-1, 0] = data[:, 0]
654 data_new[1:-1, -1] = data[:, -1]
655 data_new[0, 1:-1] = data[0, :]
656 data_new[-1, 1:-1] = data[-1, :]
658 for i, j in zip([-1, -1, 0, 0], [-1, 0, -1, 0]):
659 data_new[i, j] = data[i, j]
661 return data_new
663 def _regular_data_to_grid(self, lengths, widths, data, filename):
664 '''
665 Interpolate regular data onto topotile grid.
667 Regular gridded data is interpolated onto the latlon grid of the
668 topotile. It is then stored as a gmt-readable .grd-file.
670 :param lengths:
671 Grid coordinates along strike relative to anchor in [m].
672 :type lengths:
673 :py:class:`~numpy.ndarray`
675 :param widths:
676 Grid coordinates downdip relative to anchor in [m].
677 :type widths:
678 :py:class:`~numpy.ndarray`
680 :param data:
681 Data grid array.
682 :type data:
683 :py:class:`~numpy.ndarray`
685 :param filename:
686 Filename, where grid is stored.
687 :type filename:
688 str
689 '''
691 source = self.source
693 interpolator = scrgi(
694 (widths * num.cos(d2r * source.dip), lengths),
695 data.T,
696 bounds_error=False,
697 method='nearest')
699 points_out = self._tile_to_lw(
700 ref_lat=source.lat,
701 ref_lon=source.lon,
702 north_shift=source.north_shift,
703 east_shift=source.east_shift,
704 strike=source.strike)
706 t = self._get_topotile()
707 t.data = num.zeros_like(t.data, dtype=float)
708 t.data[:] = num.nan
710 t.data = interpolator(points_out).reshape(t.data.shape)
712 _save_grid(t.y(), t.x(), t.data, filename=filename)
714 def patch_data_to_grid(self, data, *args, **kwargs):
715 '''
716 Generate a grid file based on regular patch wise data.
718 :param data:
719 Patchwise grid data.
720 :type data:
721 :py:class:`~numpy.ndarray`
722 '''
724 lengths, widths = self._patches_to_lw()
726 data_new = self._prep_patch_grid_data(data)
728 self._regular_data_to_grid(lengths, widths, data_new, *args, **kwargs)
730 def xy_data_to_grid(self, x, y, data, *args, **kwargs):
731 '''
732 Generate a grid file based on gridded data using xy coordinates.
734 Convert a grid based on relative fault coordinates (range -1:1) along
735 strike (x) and downdip (y) into a .grd file.
737 :param x:
738 Relative point coordinate along strike (range: -1:1).
739 :type x:
740 float or :py:class:`~numpy.ndarray`
742 :param y:
743 Relative downdip point coordinate (range: -1:1).
744 :type y:
745 float or :py:class:`~numpy.ndarray`
747 :param data:
748 Patchwise grid data.
749 :type data:
750 :py:class:`~numpy.ndarray`
751 '''
753 lengths, widths = self._xy_to_lw(x, y)
755 self._regular_data_to_grid(
756 lengths, widths, data.reshape((lengths.shape[0], widths.shape[0])),
757 *args, **kwargs)
759 def draw_image(self, gridfile, cmap, cbar=True, **kwargs):
760 '''
761 Draw grid data as image and include, if whished, a colorbar.
763 :param gridfile:
764 File of the grid which shall be plotted.
765 :type gridfile:
766 str
768 :param cmap:
769 Name of the colormap, which shall be used. A .cpt-file
770 "my_<cmap>.cpt" must exist.
771 :type cmap:
772 str
774 :param cbar:
775 If ``True``, a colorbar corresponding to the grid data is
776 added. Keyword arguments are parsed to it.
777 :type cbar:
778 bool
779 '''
781 self.gmt.grdimage(
782 gridfile,
783 C='my_%s.cpt' % cmap,
784 E='200',
785 Q=True,
786 n='+t0.0',
787 *self.jxyr)
789 if cbar:
790 self.draw_colorbar(cmap=cmap, **kwargs)
792 def draw_contour(
793 self,
794 gridfile,
795 contour_int,
796 anot_int,
797 angle=None,
798 unit='',
799 color='',
800 style='',
801 **kwargs):
803 '''
804 Draw grid data as contour lines.
806 :param gridfile:
807 File of the grid which shall be plotted.
808 :type gridfile:
809 str
811 :param contour_int:
812 Interval of contour lines in units of the gridfile.
813 :type contour_int:
814 float
816 :param anot_int:
817 Interval of labelled contour lines in units of the gridfile. Must
818 be a integer multiple of contour_int.
819 :type anot_int:
820 float
822 :param angle:
823 Rotation angle of the labels in [deg].
824 :type angle:
825 float
827 :param unit:
828 Name of the unit in the grid file. It will be displayed behind the
829 label on labelled contour lines.
830 :type unit:
831 str
833 :param color:
834 GMT readable color code or string of the contour lines.
835 :type color:
836 str
838 :param style:
839 Line style of the contour lines. If not given, solid lines are
840 plotted.
841 :type style:
842 str
843 '''
845 pen_size = self._fontsize / 40.
847 if not color:
848 color = self._fontcolor
850 a_string = '%g+f%s,%s+r%gc+u%s' % (
851 anot_int, self.font, color, pen_size*4, unit)
852 if angle:
853 a_string += '+a%g' % angle
854 c_string = '%g' % contour_int
856 if kwargs:
857 kwargs['A'], kwargs['C'] = a_string, c_string
858 else:
859 kwargs = dict(A=a_string, C=c_string)
861 if style:
862 style = ',' + style
864 args = ['-Wc%gp,%s%s+s' % (pen_size, color, style)]
866 self.gmt.grdcontour(
867 gridfile,
868 S='10',
869 W='a%gp,%s%s+s' % (pen_size*4, color, style),
870 *self.jxyr + args,
871 **kwargs)
873 def draw_colorbar(self, cmap, label='', anchor='top_right', **kwargs):
874 '''
875 Draw a colorbar based on a existing colormap.
877 :param cmap:
878 Name of the colormap, which shall be used. A .cpt-file
879 "my_<cmap>.cpt" must exist.
880 :type cmap:
881 str
883 :param label:
884 Title of the colorbar.
885 :type label:
886 str
888 :param anchor:
889 Placement of the colorbar. Combine 'top', 'center' and 'bottom'
890 with 'left', None for middle and 'right'.
891 :type anchor:
892 str
893 '''
895 if not kwargs:
896 kwargs = {}
898 if label:
899 kwargs['B'] = 'af+l%s' % label
901 kwargs['C'] = 'my_%s.cpt' % cmap
902 a_str = cbar_anchor[anchor]
904 w = self.width / 3.
905 h = w / 10.
907 lgap = rgap = w / 10.
908 bgap, tgap = h, h / 10.
910 dx, dy = 2.5 * lgap, 2. * tgap
912 if 'bottom' in anchor:
913 dy += 4 * h
915 self.gmt.psscale(
916 D='j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy),
917 F='+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap),
918 *self.jxyr,
919 **kwargs)
921 def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs):
922 '''
923 Draw vectors based on two grid files.
925 Two grid files for vector lengths in x and y need to be given. The
926 function calls gmt.grdvector. All arguments defined for this function
927 in gmt can be passed as keyword arguments. Different standard settings
928 are applied if not defined differently.
930 :param x_gridfile:
931 File of the grid defining vector lengths in x.
932 :type x_gridfile:
933 str
935 :param y_gridfile:
936 File of the grid defining vector lengths in y.
937 :type y_gridfile:
938 str
940 :param vcolor:
941 Vector face color as defined in "G" option.
942 :type vcolor:
943 str
944 '''
946 kwargs['S'] = kwargs.get('S', 'il1.')
947 kwargs['I'] = kwargs.get('I', 'x20')
948 kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black')
949 kwargs['Q'] = kwargs.get('Q', '4c+e+n5km+h1')
951 self.gmt.grdvector(
952 x_gridfile, y_gridfile,
953 G='%s' % 'lightgrey' if not vcolor else vcolor,
954 *self.jxyr,
955 **kwargs)
957 def draw_dynamic_data(self, data, **kwargs):
958 '''
959 Draw an image of any data gridded on the patches e.g dislocation.
961 :param data:
962 Patchwise grid data.
963 :type data:
964 :py:class:`~numpy.ndarray`
965 '''
967 plot_data = data
969 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
971 clim = kwargs.pop('clim', (plot_data.min(), plot_data.max()))
973 cpt = []
974 if not op.exists('my_%s.cpt' % kwargs['cmap']):
975 make_colormap(self.gmt, clim[0], clim[1],
976 cmap=kwargs['cmap'], space=False)
977 cpt = [kwargs['cmap']]
979 tmp_grd_file = 'tmpdata.grd'
980 self.patch_data_to_grid(plot_data, tmp_grd_file)
981 self.draw_image(tmp_grd_file, **kwargs)
983 clear_temp(gridfiles=[tmp_grd_file], cpts=cpt)
985 def draw_patch_parameter(self, attribute, **kwargs):
986 '''
987 Draw an image of a chosen patch attribute e.g traction.
989 :param attribute:
990 Patch attribute, which is plotted. All patch attributes can be
991 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`)
992 and also ``traction``, ``tx``, ``ty`` or ``tz`` to display the
993 length or the single components of the traction vector.
994 :type attribute:
995 str
996 '''
998 a = attribute
999 source = self.source
1001 if a == 'traction':
1002 data = num.linalg.norm(source.get_tractions(), axis=1)
1003 elif a == 'tx':
1004 data = source.get_tractions()[:, 0]
1005 elif a == 'ty':
1006 data = source.get_tractions()[:, 1]
1007 elif a == 'tz':
1008 data = source.get_tractions()[:, 2]
1009 else:
1010 data = source.get_patch_attribute(attribute)
1012 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1013 data *= factor
1015 kwargs['label'] = kwargs.get(
1016 'label',
1017 '%s [%s]' % (a, cbar_helper[a]['unit']))
1019 self.draw_dynamic_data(data, **kwargs)
1021 def draw_time_contour(self, store, clevel=[], **kwargs):
1022 '''
1023 Draw high contour lines of the rupture front propgation time.
1025 :param store:
1026 Greens function store, which is used for time calculation.
1027 :type store:
1028 :py:class:`~pyrocko.gf.store.Store`
1029 :param clevel:
1030 List of times, for which contour lines are drawn.
1031 :type clevel:
1032 :py:class:`list` of :py:class:`float`
1033 '''
1035 _, _, _, _, points_xy = self.source._discretize_points(store, cs='xyz')
1036 _, _, times, _, _, _ = self.source.get_vr_time_interpolators(store)
1038 scaler = AutoScaler(mode='0-max', approx_ticks=8)
1039 scale = scaler.make_scale([num.min(times), num.max(times)])
1041 if clevel:
1042 if len(clevel) > 1:
1043 kwargs['anot_int'] = num.min(num.diff(clevel))
1044 else:
1045 kwargs['anot_int'] = clevel[0]
1047 kwargs['contour_int'] = kwargs['anot_int']
1048 kwargs['L'] = '0/%g' % num.max(clevel)
1050 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1051 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1052 kwargs['unit'] = kwargs.get('unit', cbar_helper['time']['unit'])
1053 kwargs['L'] = kwargs.get('L', '0/%g' % (num.max(times) + 1.))
1054 kwargs['G'] = kwargs.get('G', 'n2/3c')
1056 tmp_grd_file = 'tmpdata.grd'
1057 self.xy_data_to_grid(points_xy[:, 0], points_xy[:, 1],
1058 times, tmp_grd_file)
1059 self.draw_contour(tmp_grd_file, **kwargs)
1061 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
1063 def draw_points(self, lats, lons, symbol='point', size=None, **kwargs):
1064 '''
1065 Draw points at given locations.
1067 :param lats:
1068 Point latitude coordinates in [deg].
1069 :type lats:
1070 iterable of :py:class:`float`
1072 :param lons:
1073 Point longitude coordinates in [deg].
1074 :type lons:
1075 iterable of :py:class:`float`
1077 :param symbol:
1078 Define symbol of the points (``'star', 'circle', 'point',
1079 'triangle'``) - default is ``point``.
1080 :type symbol:
1081 str
1083 :param size:
1084 Size of the points in [points].
1085 :type size:
1086 float
1087 '''
1089 sym_to_gmt = dict(
1090 star='a',
1091 circle='c',
1092 point='p',
1093 triangle='t')
1095 lats = num.atleast_1d(lats)
1096 lons = num.atleast_1d(lons)
1098 if lats.shape[0] != lons.shape[0]:
1099 raise IndexError('lats and lons do not have the same shape!')
1101 if size is None:
1102 size = self._fontsize / 3.
1104 kwargs['S'] = kwargs.get('S', sym_to_gmt[symbol] + '%gp' % size)
1105 kwargs['G'] = kwargs.get('G', gmtpy.color('scarletred2'))
1106 kwargs['W'] = kwargs.get('W', '2p,%s' % self._fontcolor)
1108 self.gmt.psxy(
1109 in_columns=[lons, lats],
1110 *self.jxyr,
1111 **kwargs)
1113 def draw_nucleation_point(self, **kwargs):
1114 '''
1115 Plot the nucleation point onto the map.
1116 '''
1118 nlat, nlon = xy_to_latlon(
1119 self.source, self.source.nucleation_x, self.source.nucleation_y)
1121 self.draw_points(nlat, nlon, **kwargs)
1123 def draw_dislocation(self, time=None, component='', **kwargs):
1124 '''
1125 Draw dislocation onto map at any time.
1127 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given
1128 component the patchwise dislocation is plotted onto the map.
1130 :param time:
1131 Time after origin, for which dislocation is computed. If ``None``,
1132 ``tmax`` is taken.
1133 :type time:
1134 float
1136 :param component:
1137 Dislocation component, which shall be plotted: ``x`` along strike,
1138 ``y`` along updip, ``z`` normal. If ``None``, the
1139 length of the dislocation vector is plotted.
1140 :type component: str
1141 '''
1143 disl = self.source.get_slip(time=time)
1145 if component:
1146 data = disl[:, c2disl[component]]
1147 else:
1148 data = num.linalg.norm(disl, axis=1)
1150 kwargs['label'] = kwargs.get(
1151 'label', 'u%s [m]' % (component))
1153 self.draw_dynamic_data(data, **kwargs)
1155 def draw_dislocation_contour(
1156 self, time=None, component=None, clevel=[], **kwargs):
1157 '''
1158 Draw dislocation contour onto map at any time.
1160 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given
1161 component the patchwise dislocation is plotted as contour onto the map.
1163 :param time:
1164 Time after origin, for which dislocation is computed. If ``None``,
1165 ``tmax`` is taken.
1166 :type time:
1167 float
1169 :param component:
1170 Dislocation component, which shall be plotted: ``x`` along strike,
1171 ``y`` along updip, ``z`` normal``. If ``None``, the length of the
1172 dislocation vector is plotted.
1173 :type component: str
1175 :param clevel:
1176 Times, for which contour lines are drawn.
1177 :type clevel:
1178 :py:class:`list` of :py:class:`float`
1179 '''
1181 disl = self.source.get_slip(time=time)
1183 if component:
1184 data = disl[:, c2disl[component]]
1185 else:
1186 data = num.linalg.norm(disl, axis=1)
1188 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1189 scale = scaler.make_scale([num.min(data), num.max(data)])
1191 if clevel:
1192 if len(clevel) > 1:
1193 kwargs['anot_int'] = num.min(num.diff(clevel))
1194 else:
1195 kwargs['anot_int'] = clevel[0]
1197 kwargs['contour_int'] = kwargs['anot_int']
1198 kwargs['L'] = '%g/%g' % (
1199 num.min(clevel) - kwargs['contour_int'],
1200 num.max(clevel) + kwargs['contour_int'])
1201 else:
1202 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1203 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1204 kwargs['L'] = kwargs.get('L', '%g/%g' % (
1205 num.min(data) - 1., num.max(data) + 1.))
1207 kwargs['unit'] = kwargs.get('unit', ' m')
1208 kwargs['G'] = kwargs.get('G', 'n2/3c')
1210 tmp_grd_file = 'tmpdata.grd'
1211 self.patch_data_to_grid(data, tmp_grd_file)
1212 self.draw_contour(tmp_grd_file, **kwargs)
1214 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
1216 def draw_dislocation_vector(self, time=None, **kwargs):
1217 '''
1218 Draw vector arrows onto map indicating direction of dislocation.
1220 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1221 and given component the dislocation is plotted as vectors onto the map.
1223 :param time:
1224 Time after origin [s], for which dislocation is computed. If
1225 ``None``, ``tmax`` is used.
1226 :type time:
1227 float
1228 '''
1230 disl = self.source.get_slip(time=time)
1232 p_strike = self.source.get_patch_attribute('strike') * d2r
1233 p_dip = self.source.get_patch_attribute('dip') * d2r
1235 disl_dh = num.cos(p_dip) * disl[:, 1]
1236 disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh
1237 disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh
1239 tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')]
1241 self.patch_data_to_grid(disl_n, tmp_grd_files[0])
1242 self.patch_data_to_grid(disl_e, tmp_grd_files[1])
1244 self.draw_vector(
1245 tmp_grd_files[1], tmp_grd_files[0],
1246 **kwargs)
1248 clear_temp(gridfiles=tmp_grd_files, cpts=[])
1250 def draw_top_edge(self, **kwargs):
1251 '''
1252 Indicate rupture top edge on map.
1253 '''
1255 outline = self.source.outline(cs='latlondepth')
1256 top_edge = outline[:2, :]
1258 kwargs = kwargs or {}
1259 kwargs['W'] = kwargs.get(
1260 'W', '%gp,%s' % (self._fontsize / 10., gmtpy.color('orange3')))
1262 self.gmt.psxy(
1263 in_columns=[top_edge[:, 1], top_edge[:, 0]],
1264 *self.jxyr,
1265 **kwargs)
1268class RuptureView(Object):
1269 '''
1270 Plot of attributes and results of the
1271 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1272 '''
1274 _patches_to_lw = RuptureMap._patches_to_lw
1276 def __init__(self, source=None, figsize=None, fontsize=None):
1277 self._source = source
1278 self._axes = None
1280 self.figsize = figsize or mpl_papersize('halfletter', 'landscape')
1281 self.fontsize = fontsize or 10
1282 mpl_init(fontsize=self.fontsize)
1284 self._fig = None
1285 self._axes = None
1286 self._is_1d = False
1288 @property
1289 def source(self):
1290 '''
1291 PseudoDynamicRupture whose attributes are plotted.
1293 Note, that source.patches attribute needs to be calculated for
1294 :type source: :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1295 '''
1297 if self._source is None:
1298 raise SourceError('No source given. Please define it!')
1300 if not isinstance(self._source, PseudoDynamicRupture):
1301 raise SourceError('This function is only capable for a source of'
1302 ' type: %s' % type(PseudoDynamicRupture()))
1304 if not self._source.patches:
1305 raise TypeError('No source patches are defined. Please run'
1306 '\"discretize_patches\" on your source')
1308 return self._source
1310 @source.setter
1311 def source(self, source):
1312 self._source = source
1314 def _setup(
1315 self,
1316 title='',
1317 xlabel='',
1318 ylabel='',
1319 aspect=1.,
1320 spatial_plot=True,
1321 **kwargs):
1323 if self._fig is not None and self._axes is not None:
1324 return
1326 self._fig = plt.figure(figsize=self.figsize)
1328 self._axes = self._fig.add_subplot(1, 1, 1, aspect=aspect)
1329 ax = self._axes
1331 if ax is not None:
1332 ax.set_title(title)
1333 ax.grid(alpha=.3)
1334 ax.set_xlabel(xlabel)
1335 ax.set_ylabel(ylabel)
1337 if spatial_plot:
1338 ax.xaxis.set_major_formatter(FuncFormatter(km_formatter))
1339 ax.yaxis.set_major_formatter(FuncFormatter(km_formatter))
1341 def _clear_all(self):
1342 plt.cla()
1343 plt.clf()
1344 plt.close()
1346 self._fig, self._axes = None, None
1348 def _draw_scatter(self, x, y, *args, **kwargs):
1349 default_kwargs = dict(
1350 linewidth=0,
1351 marker='o',
1352 markerfacecolor=mpl_color('skyblue2'),
1353 markersize=6.,
1354 markeredgecolor=mpl_color('skyblue3'))
1355 default_kwargs.update(kwargs)
1357 if self._axes is not None:
1358 self._axes.plot(x, y, *args, **default_kwargs)
1360 def _draw_image(self, length, width, data, *args, **kwargs):
1361 if self._axes is not None:
1362 if 'extent' not in kwargs:
1363 kwargs['extent'] = [
1364 num.min(length), num.max(length),
1365 num.max(width), num.min(width)]
1367 im = self._axes.imshow(
1368 data,
1369 *args,
1370 interpolation='none',
1371 vmin=kwargs.get('clim', [None])[0],
1372 vmax=kwargs.get('clim', [None, None])[1],
1373 **kwargs)
1375 del kwargs['extent']
1376 if 'aspect' in kwargs:
1377 del kwargs['aspect']
1378 if 'clim' in kwargs:
1379 del kwargs['clim']
1380 if 'cmap' in kwargs:
1381 del kwargs['cmap']
1383 plt.colorbar(
1384 im, *args, shrink=0.9, pad=0.03, aspect=15., **kwargs)
1386 def _draw_contour(self, x, y, data, clevel=None, unit='', *args, **kwargs):
1387 setup_kwargs = dict(
1388 xlabel=kwargs.pop('xlabel', 'along strike [km]'),
1389 ylabel=kwargs.pop('ylabel', 'down dip [km]'),
1390 title=kwargs.pop('title', ''),
1391 aspect=kwargs.pop('aspect', 1))
1393 self._setup(**setup_kwargs)
1395 if self._axes is not None:
1396 if clevel is None:
1397 scaler = AutoScaler(mode='min-max')
1398 scale = scaler.make_scale([data.min(), data.max()])
1400 clevel = num.arange(scale[0], scale[1] + scale[2], scale[2])
1402 if not isinstance(clevel, num.ndarray):
1403 clevel = num.array([clevel])
1405 clevel = clevel[clevel < data.max()]
1407 cont = self._axes.contour(
1408 x, y, data, clevel, *args,
1409 linewidths=1.5, **kwargs)
1411 clabels = self._axes.clabel(
1412 cont, clevel[::2], *args,
1413 inline=1, fmt='%g' + '%s' % unit,
1414 inline_spacing=15, rightside_up=True, use_clabeltext=True,
1415 **kwargs)
1417 plt.setp(clabels, path_effects=[
1418 patheffects.withStroke(linewidth=1.25, foreground='beige'),
1419 patheffects.Normal()])
1421 def draw_points(self, length, width, *args, **kwargs):
1422 '''
1423 Draw a point onto the figure.
1425 Args and kwargs can be defined according to
1426 :py:func:`matplotlib.pyplot.scatter`.
1428 :param length:
1429 Point(s) coordinate on the rupture plane along strike relative to
1430 the anchor point in [m].
1431 :type length:
1432 float, :py:class:`~numpy.ndarray`
1434 :param width:
1435 Point(s) coordinate on the rupture plane along downdip relative to
1436 the anchor point in [m].
1437 :type width:
1438 float, :py:class:`~numpy.ndarray`
1439 '''
1441 if self._axes is not None:
1442 kwargs['color'] = kwargs.get('color', mpl_color('scarletred2'))
1443 kwargs['s'] = kwargs.get('s', 100.)
1445 self._axes.scatter(length, width, *args, **kwargs)
1447 def draw_dynamic_data(self, data, **kwargs):
1448 '''
1449 Draw an image of any data gridded on the patches e.g dislocation.
1451 :param data:
1452 Patchwise grid data.
1453 :type data:
1454 :py:class:`~numpy.ndarray`
1455 '''
1457 plot_data = data
1459 anchor_x, anchor_y = map_anchor[self.source.anchor]
1461 length, width = xy_to_lw(
1462 self.source, num.array([-1., 1.]), num.array([-1., 1.]))
1464 data = plot_data.reshape(self.source.nx, self.source.ny).T
1466 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
1468 setup_kwargs = dict(
1469 xlabel='along strike [km]',
1470 ylabel='down dip [km]',
1471 title='',
1472 aspect=1)
1473 setup_kwargs.update(kwargs)
1475 kwargs = {k: v for k, v in kwargs.items() if k not in
1476 ('xlabel', 'ylabel', 'title')}
1477 self._setup(**setup_kwargs)
1478 self._draw_image(length=length, width=width, data=data, **kwargs)
1480 def draw_patch_parameter(self, attribute, **kwargs):
1481 '''
1482 Draw an image of a chosen patch attribute e.g traction.
1484 :param attribute:
1485 Patch attribute, which is plotted. All patch attributes can be
1486 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`)
1487 and also ``'traction', 'tx', 'ty', 'tz'`` to display the length
1488 or the single components of the traction vector.
1489 :type attribute:
1490 str
1491 '''
1493 a = attribute
1494 source = self.source
1496 if a == 'traction':
1497 data = num.linalg.norm(source.get_tractions(), axis=1)
1498 elif a == 'tx':
1499 data = source.get_tractions()[:, 0]
1500 elif a == 'ty':
1501 data = source.get_tractions()[:, 1]
1502 elif a == 'tz':
1503 data = source.get_tractions()[:, 2]
1504 else:
1505 data = source.get_patch_attribute(attribute)
1507 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1508 data *= factor
1510 kwargs['label'] = kwargs.get(
1511 'label',
1512 '%s [%s]' % (a, cbar_helper[a]['unit']))
1514 return self.draw_dynamic_data(data=data, **kwargs)
1516 def draw_time_contour(self, store, clevel=[], **kwargs):
1517 '''
1518 Draw high resolution contours of the rupture front propgation time
1520 :param store:
1521 Greens function store, which is used for time calculation.
1522 :type store:
1523 :py:class:`~pyrocko.gf.store.Store`
1525 :param clevel:
1526 Levels of the contour lines. If no levels are given, they are
1527 automatically computed based on ``tmin`` and ``tmax``.
1528 :type clevel:
1529 list
1530 '''
1532 source = self.source
1533 default_kwargs = dict(
1534 colors='#474747'
1535 )
1536 default_kwargs.update(kwargs)
1538 _, _, _, _, points_xy = source._discretize_points(store, cs='xyz')
1539 _, _, _, _, time_interpolator, _ = source.get_vr_time_interpolators(
1540 store)
1542 times = time_interpolator.values
1544 scaler = AutoScaler(mode='min-max', approx_ticks=8)
1545 scale = scaler.make_scale([times.min(), times.max()])
1547 if len(clevel) == 0:
1548 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1550 points_x = time_interpolator.grid[0]
1551 points_y = time_interpolator.grid[1]
1553 self._draw_contour(points_x, points_y, data=times.T,
1554 clevel=clevel, unit='s', **default_kwargs)
1556 def draw_nucleation_point(self, **kwargs):
1557 '''
1558 Draw the nucleation point onto the map.
1559 '''
1561 nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y
1562 length, width = xy_to_lw(self.source, nuc_x, nuc_y)
1563 self.draw_points(length, width, marker='o', **kwargs)
1565 def draw_dislocation(self, time=None, component='', **kwargs):
1566 '''
1567 Draw dislocation onto map at any time.
1569 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1570 and given component the patchwise dislocation is plotted onto the map.
1572 :param time:
1573 Time after origin [s], for which dislocation is computed.
1574 If ``None``, ``tmax`` is taken.
1575 :type time:
1576 float
1578 :param component:
1579 Dislocation component, which shall be plotted: ``x``
1580 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
1581 length of the dislocation vector is plotted
1582 :type component: str
1583 '''
1585 disl = self.source.get_slip(time=time)
1587 if component:
1588 data = disl[:, c2disl[component]]
1589 else:
1590 data = num.linalg.norm(disl, axis=1)
1592 kwargs['label'] = kwargs.get(
1593 'label', 'u%s [m]' % (component))
1595 self.draw_dynamic_data(data, **kwargs)
1597 def draw_dislocation_contour(
1598 self, time=None, component=None, clevel=[], **kwargs):
1599 '''
1600 Draw dislocation contour onto map at any time.
1602 For a given time (if time is ``None``, ``tmax`` is used) and given
1603 component the patchwise dislocation is plotted as contour onto the map.
1605 :param time:
1606 Time after origin, for which dislocation is computed. If
1607 ``None``, ``tmax`` is taken.
1608 :type time:
1609 float
1611 :param component:
1612 Dislocation component, which shall be plotted. ``x``
1613 along strike, ``y`` along updip, ``z`` - normal. If ``None``
1614 is given, the length of the dislocation vector is plotted.
1615 :type component:
1616 str
1617 '''
1619 disl = self.source.get_slip(time=time)
1621 if component:
1622 data = disl[:, c2disl[component]]
1623 else:
1624 data = num.linalg.norm(disl, axis=1)
1626 data = data.reshape(self.source.ny, self.source.nx, order='F')
1628 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1629 scale = scaler.make_scale([data.min(), data.max()])
1631 if len(clevel) == 0:
1632 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1634 anchor_x, anchor_y = map_anchor[self.source.anchor]
1636 length, width = self._patches_to_lw()
1637 length, width = length[1:-1], width[1:-1]
1639 kwargs['colors'] = kwargs.get('colors', '#474747')
1641 self._setup(**kwargs)
1642 self._draw_contour(
1643 length, width, data=data, clevel=clevel, unit='m', **kwargs)
1645 def draw_source_dynamics(
1646 self, variable, store, deltat=None, *args, **kwargs):
1647 '''
1648 Display dynamic source parameter.
1650 Fast inspection possibility for the cumulative moment and the source
1651 time function approximation (assuming equal paths between different
1652 patches and observation point - valid for an observation point in the
1653 far field perpendicular to the source strike), so the cumulative moment
1654 rate function.
1656 :param variable:
1657 Dynamic parameter, which shall be plotted. Choose
1658 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1659 :type variable:
1660 str
1662 :param store:
1663 Greens function store, whose store.config.deltat defines
1664 the time increment between two parameter snapshots. If ``store`` is
1665 not given, the time increment is defined is taken from ``deltat``.
1666 :type store:
1667 :py:class:`~pyrocko.gf.store.Store`
1669 :param deltat:
1670 Time increment between two parameter snapshots. If not
1671 given, store.config.deltat is used to define ``deltat``.
1672 :type deltat:
1673 float
1674 '''
1676 v = variable
1678 data, times = self.source.get_moment_rate(store=store, deltat=deltat)
1679 deltat = num.concatenate([(num.diff(times)[0], ), num.diff(times)])
1681 if v in ('moment_rate', 'stf'):
1682 name, unit = 'dM/dt', 'Nm/s'
1683 elif v in ('cumulative_moment', 'moment'):
1684 data = num.cumsum(data * deltat)
1685 name, unit = 'M', 'Nm'
1686 else:
1687 raise ValueError('No dynamic data for given variable %s found' % v)
1689 self._setup(xlabel='time [s]',
1690 ylabel='%s / %.2g %s' % (name, data.max(), unit),
1691 aspect='auto',
1692 spatial_plot=False)
1693 self._draw_scatter(times, data/num.max(data), *args, **kwargs)
1694 self._is_1d = True
1696 def draw_patch_dynamics(
1697 self, variable, nx, ny, store=None, deltat=None, *args, **kwargs):
1698 '''
1699 Display dynamic boundary element / patch parameter.
1701 Fast inspection possibility for different dynamic parameter for a
1702 single patch / boundary element. The chosen parameter is plotted for
1703 the chosen patch.
1705 :param variable:
1706 Dynamic parameter, which shall be plotted. Choose
1707 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment').
1708 :type variable:
1709 str
1711 :param nx:
1712 Patch index along strike (range: 0:source.nx - 1).
1713 :type nx:
1714 int
1716 :param nx:
1717 Patch index downdip (range: 0:source.ny - 1).
1718 :type nx:
1719 int
1721 :param store:
1722 Greens function store, whose store.config.deltat defines
1723 the time increment between two parameter snapshots. If ``store`` is
1724 not given, the time increment is defined is taken from ``deltat``.
1725 :type store:
1726 :py:class:`~pyrocko.gf.store.Store`
1728 :param deltat:
1729 Time increment between two parameter snapshots. If not given,
1730 store.config.deltat is used to define ``deltat``.
1731 :type deltat:
1732 float
1733 '''
1735 v = variable
1736 source = self.source
1737 idx = nx * source.ny + nx
1739 m = re.match(r'dislocation_([xyz])', v)
1741 if v in ('moment_rate', 'cumulative_moment', 'moment', 'stf'):
1742 data, times = source.get_moment_rate_patches(
1743 store=store, deltat=deltat)
1744 elif 'dislocation' in v or 'slip_rate' == v:
1745 ddisloc, times = source.get_delta_slip(store=store, deltat=deltat)
1747 if v in ('moment_rate', 'stf'):
1748 data, times = source.get_moment_rate_patches(
1749 store=store, deltat=deltat)
1750 name, unit = 'dM/dt', 'Nm/s'
1751 elif v in ('cumulative_moment', 'moment'):
1752 data, times = source.get_moment_rate_patches(
1753 store=store, deltat=deltat)
1754 data = num.cumsum(data, axis=1)
1755 name, unit = 'M', 'Nm'
1756 elif v == 'slip_rate':
1757 data, times = source.get_slip_rate(store=store, deltat=deltat)
1758 name, unit = 'du/dt', 'm/s'
1759 elif v == 'dislocation':
1760 data, times = source.get_delta_slip(store=store, deltat=deltat)
1761 data = num.linalg.norm(num.cumsum(data, axis=2), axis=1)
1762 name, unit = 'du', 'm'
1763 elif m:
1764 data, times = source.get_delta_slip(store=store, deltat=deltat)
1765 data = num.cumsum(data, axis=2)[:, c2disl[m.group(1)], :]
1766 name, unit = 'du%s' % m.group(1), 'm'
1767 else:
1768 raise ValueError('No dynamic data for given variable %s found' % v)
1770 self._setup(xlabel='time [s]',
1771 ylabel='%s / %.2g %s' % (name, num.max(data), unit),
1772 aspect='auto',
1773 spatial_plot=False)
1774 self._draw_scatter(times, data[idx, :]/num.max(data),
1775 *args, **kwargs)
1776 self._is_1d = True
1778 def finalize(self):
1779 if self._is_1d:
1780 return
1782 length, width = xy_to_lw(
1783 self.source, num.array([-1., 1.]), num.array([1., -1.]))
1785 self._axes.set_xlim(length)
1786 self._axes.set_ylim(width)
1788 def gcf(self):
1789 self.finalize()
1790 return self._fig
1792 def save(self, filename, dpi=None):
1793 '''
1794 Save plot to file.
1796 :param filename:
1797 Filename and path, where the plot is stored.
1798 :type filename:
1799 str
1801 :param dpi:
1802 Resolution of the output plot in [dpi].
1803 :type dpi:
1804 int
1805 '''
1807 self.finalize()
1808 try:
1809 self._fig.savefig(fname=filename, dpi=dpi, bbox_inches='tight')
1810 except TypeError:
1811 self._fig.savefig(filename=filename, dpi=dpi, bbox_inches='tight')
1813 self._clear_all()
1815 def show_plot(self):
1816 '''
1817 Show plot.
1818 '''
1820 self.finalize()
1821 plt.show()
1822 self._clear_all()
1825def render_movie(fn_path, output_path, framerate=20):
1826 '''
1827 Generate a mp4 movie based on given png files using
1828 `ffmpeg <https://ffmpeg.org>`_.
1830 Render a movie based on a set of given .png files in fn_path. All files
1831 must have a filename specified by ``fn_path`` (e.g. giving ``fn_path`` with
1832 ``/temp/f%04.png`` a valid png filename would be ``/temp/f0001.png``). The
1833 files must have a numbering, indicating their order in the movie.
1835 :param fn_path:
1836 Path and fileformat specification of the input .png files.
1837 :type fn_path:
1838 str
1840 :param output_path:
1841 Path and filename of the output ``.mp4`` movie file.
1842 :type output_path:
1843 str
1845 :param deltat:
1846 Time between individual frames (``1 / framerate``) in [s].
1847 :type deltat:
1848 float
1849 '''
1851 try:
1852 check_call(['ffmpeg', '-loglevel', 'panic'])
1853 except CalledProcessError:
1854 pass
1855 except (TypeError, IOError):
1856 logger.warning(
1857 'Package ffmpeg needed for movie rendering. Please install it '
1858 '(e.g. on linux distr. via sudo apt-get ffmpeg) and retry.')
1859 return False
1861 try:
1862 check_call([
1863 'ffmpeg', '-loglevel', 'info', '-y',
1864 '-framerate', '%g' % framerate,
1865 '-i', fn_path,
1866 '-vcodec', 'libx264',
1867 '-preset', 'medium',
1868 '-tune', 'animation',
1869 '-pix_fmt', 'yuv420p',
1870 '-movflags', '+faststart',
1871 '-filter:v', "crop='floor(in_w/2)*2:floor(in_h/2)*2'",
1872 '-crf', '15',
1873 output_path])
1875 return True
1876 except CalledProcessError as e:
1877 logger.warning(e)
1878 return False
1881def render_gif(fn, output_path, loops=-1):
1882 '''
1883 Generate a gif based on a given mp4 using ffmpeg.
1885 Render a gif based on a given .mp4 movie file in ``fn`` path.
1887 :param fn:
1888 Path and file name of the input .mp4 file.
1889 :type fn:
1890 str
1892 :param output_path:
1893 Path and filename of the output animated gif file.
1894 :type output_path:
1895 str
1896 :param loops:
1897 Number of gif repeats (loops). ``-1`` is not repetition, ``0``
1898 infinite.
1899 :type loops:
1900 int
1901 '''
1903 try:
1904 check_call(['ffmpeg', '-loglevel', 'panic'])
1905 except CalledProcessError:
1906 pass
1907 except (TypeError, IOError):
1908 logger.warning(
1909 'Package ffmpeg needed for movie rendering. Please install it '
1910 '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.')
1911 return False
1913 try:
1914 check_call([
1915 'ffmpeg', '-hide_banner', '-loglevel', 'panic', '-y', '-i',
1916 fn,
1917 '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse',
1918 '-loop', '%d' % loops,
1919 output_path])
1921 return True
1922 except CalledProcessError as e:
1923 logger.warning(e)
1924 return False
1927def rupture_movie(
1928 source,
1929 store,
1930 variable='dislocation',
1931 draw_time_contours=False,
1932 fn_path='.',
1933 prefix='',
1934 plot_type='map',
1935 deltat=None,
1936 framerate=None,
1937 store_images=False,
1938 render_as_gif=False,
1939 gif_loops=-1,
1940 **kwargs):
1941 '''
1942 Generate a movie based on a given source for dynamic parameter.
1944 Create a MPEG-4 movie or gif of one of the following dynamic parameters
1945 (``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1946 (along updip), ``dislocation_z`` (normal), ``slip_rate``, ``moment_rate``).
1947 If desired, the single snap shots can be stored as images as well.
1948 ``kwargs`` have to be given according to the chosen ``plot_type``.
1950 :param source:
1951 Pseudo dynamic rupture, for which the movie is produced.
1952 :type source:
1953 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
1955 :param store:
1956 Greens function store, which is used for time calculation. If
1957 ``deltat`` is not given, it is taken from the store.config.deltat
1958 :type store:
1959 :py:class:`~pyrocko.gf.store.Store`
1961 :param variable:
1962 Dynamic parameter, which shall be plotted. Choose between
1963 ``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1964 (along updip), ``dislocation_z`` (normal), ``slip_rate`` and
1965 ``moment_rate``, default ``dislocation``.
1966 :type variable:
1967 str
1969 :param draw_time_contours:
1970 If ``True``, corresponding isochrones are drawn on the each plots.
1971 :type draw_time_contours:
1972 bool
1974 :param fn_path:
1975 Absolut or relative path, where movie (and optional images) are stored.
1976 :type fn_path:
1977 str
1979 :param prefix:
1980 File prefix used for the movie (and optional image) files.
1981 :type prefix:
1982 str
1984 :param plot_type:
1985 Choice of plot type: ``map``, ``view`` (map plot using
1986 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap`
1987 or plane view using
1988 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`).
1989 :type plot_type:
1990 str
1992 :param deltat:
1993 Time between parameter snapshots. If not given, store.config.deltat is
1994 used to define ``deltat``.
1995 :type deltat:
1996 float
1998 :param store_images:
1999 Choice to store the single .png parameter snapshots in ``fn_path`` or
2000 not.
2001 :type store_images:
2002 bool
2004 :param render_as_gif:
2005 If ``True``, the movie is converted into a gif. If ``False``, the movie
2006 is returned as mp4.
2007 :type render_as_gif:
2008 bool
2010 :param gif_loops:
2011 If ``render_as_gif`` is ``True``, a gif with ``gif_loops`` number of
2012 loops (repetitions) is returned. ``-1`` is no repetition, ``0``
2013 infinite.
2014 :type gif_loops:
2015 int
2016 '''
2018 v = variable
2019 assert plot_type in ('map', 'view')
2021 if not source.patches:
2022 source.discretize_patches(store, interpolation='multilinear')
2024 if source.coef_mat is None:
2025 source.calc_coef_mat()
2027 prefix = prefix or v
2028 deltat = deltat or store.config.deltat
2029 framerate = max(framerate or int(1./deltat), 1)
2031 if v == 'moment_rate':
2032 data, times = source.get_moment_rate_patches(deltat=deltat)
2033 name, unit = 'dM/dt', 'Nm/s'
2034 elif 'dislocation' in v or 'slip_rate' == v:
2035 ddisloc, times = source.get_delta_slip(deltat=deltat)
2036 else:
2037 raise ValueError('No dynamic data for given variable %s found' % v)
2039 deltat = times[1] - times[0]
2041 m = re.match(r'dislocation_([xyz])', v)
2042 if m:
2043 data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]]
2044 name, unit = 'du%s' % m.group(1), 'm'
2045 elif v == 'dislocation':
2046 data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2)
2047 name, unit = 'du', 'm'
2048 elif v == 'slip_rate':
2049 data = num.linalg.norm(ddisloc, axis=2) / deltat
2050 name, unit = 'du/dt', 'm/s'
2052 if plot_type == 'map':
2053 plt_base = RuptureMap
2054 elif plot_type == 'view':
2055 plt_base = RuptureView
2056 else:
2057 raise AttributeError('invalid type: %s' % plot_type)
2059 attrs_base = get_kwargs(plt_base)
2060 kwargs_base = dict([k for k in kwargs.items() if k[0] in attrs_base])
2061 kwargs_plt = dict([k for k in kwargs.items() if k[0] not in kwargs_base])
2063 if 'clim' in kwargs_plt:
2064 data = num.clip(data, kwargs_plt['clim'][0], kwargs_plt['clim'][1])
2065 else:
2066 kwargs_plt['clim'] = [num.min(data), num.max(data)]
2068 if 'label' not in kwargs_plt:
2069 vmax = num.max(num.abs(kwargs_plt['clim']))
2070 data /= vmax
2072 kwargs_plt['label'] = '%s / %.2g %s' % (name, vmax, unit)
2073 kwargs_plt['clim'] = [i / vmax for i in kwargs_plt['clim']]
2075 with tempfile.TemporaryDirectory(suffix='pyrocko') as temp_path:
2076 fns_temp = [op.join(temp_path, 'f%09d.png' % (it + 1))
2077 for it, _ in enumerate(times)]
2078 fn_temp_path = op.join(temp_path, 'f%09d.png')
2080 for it, (t, ft) in enumerate(zip(times, fns_temp)):
2081 plot_data = data[:, it]
2083 plt = plt_base(source=source, **kwargs_base)
2084 plt.draw_dynamic_data(plot_data, **kwargs_plt)
2085 plt.draw_nucleation_point()
2087 if draw_time_contours:
2088 plt.draw_time_contour(store, clevel=[t])
2090 plt.save(ft)
2092 fn_mp4 = op.join(temp_path, 'movie.mp4')
2093 return_code = render_movie(
2094 fn_temp_path,
2095 output_path=fn_mp4,
2096 framerate=framerate)
2098 if render_as_gif and return_code:
2099 render_gif(fn=fn_mp4, output_path=op.join(
2100 fn_path, '%s_%s_gif.gif' % (prefix, plot_type)),
2101 loops=gif_loops)
2103 elif return_code:
2104 shutil.move(fn_mp4, op.join(
2105 fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
2107 else:
2108 logger.error('ffmpeg failed. Exit')
2110 if store_images:
2111 fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t))
2112 for t in times]
2114 for ft, f in zip(fns_temp, fns):
2115 shutil.move(ft, f)
2118__all__ = [
2119 'make_colormap',
2120 'clear_temp',
2121 'xy_to_latlon',
2122 'xy_to_lw',
2123 'SourceError',
2124 'RuptureMap',
2125 'RuptureView',
2126 'rupture_movie',
2127 'render_movie']