Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/dynamic_rupture.py: 82%
633 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-23 12:35 +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.path as op
9import re
10import logging
11import tempfile
12import shutil
13import inspect
15import numpy as num
16from scipy.interpolate import RegularGridInterpolator as scrgi
18from matplotlib import pyplot as plt, patheffects
19from matplotlib.ticker import FuncFormatter
21from pyrocko import orthodrome as pod
22from pyrocko.guts import Object
23from pyrocko.plot import mpl_init, mpl_papersize, mpl_color, AutoScaler, \
24 gmtpy, mpl_get_cmap
25from pyrocko.plot.automap import Map, NoTopo
26from pyrocko.gf import PseudoDynamicRupture
27from pyrocko.gf.seismosizer import map_anchor
28from pyrocko.dataset.topo.tile import Tile
30logger = logging.getLogger(__name__)
32km = 1e3
34d2m = 111180.
35m2d = 1. / d2m
37cm2inch = gmtpy.cm/gmtpy.inch
39d2r = num.pi / 180.
40r2d = 1. / d2r
43def km_formatter(v, p):
44 return '%g' % (v / km)
47def get_kwargs(cls):
48 sig = inspect.signature(cls.__init__)
49 kwargs = {}
51 if issubclass(cls.T._cls, RuptureMap):
52 for p in Map.T.properties:
53 kwargs[p.name] = p._default
55 for par in sig.parameters.values():
56 if par.default is inspect._empty:
57 continue
58 kwargs[par.name] = par.default
60 return kwargs
63def _save_grid(lats, lons, data, filename):
64 '''
65 Save lat-lon gridded data as gmt .grd file.
67 :param lats:
68 Grid latitude coordinates in [deg].
69 :type lats:
70 iterable
72 :param lons:
73 Grid longitude coordinates in [deg].
74 :type lons:
75 iterable
77 :param data:
78 Grid data of any kind.
79 :type data:
80 :py:class:`~numpy.ndarray`, ``(n_lons, n_lats)``
82 :param filename:
83 Filename of the written grid file.
84 :type filename:
85 str
86 '''
88 gmtpy.savegrd(lons, lats, data, filename=filename, naming='lonlat')
91def _mplcmap_to_gmtcpt_code(mplcmap, steps=256):
92 '''
93 Get gmt readable R/G/A code from a given matplotlib colormap.
95 :param mplcmap:
96 Name of the demanded matplotlib colormap.
97 :type mplcmap:
98 str
100 :returns:
101 Series of comma seperate R/G/B values for gmtpy usage.
102 :rtype:
103 str
104 '''
106 cmap = mpl_get_cmap(mplcmap)
108 rgbas = [cmap(i) for i in num.linspace(0, 255, steps).astype(num.int64)]
110 return ','.join(['%d/%d/%d' % (
111 c[0] * 255, c[1] * 255, c[2] * 255) for c in rgbas])
114def make_cpt_path(gmt, base_cmap):
115 return gmt.tempfilename('my_%s.cpt' % base_cmap)
118def make_colormap(
119 gmt,
120 vmin,
121 vmax,
122 C=None,
123 cmap=None,
124 space=False):
125 '''
126 Create gmt-readable colormap cpt file called my_<cmap>.cpt.
128 :param vmin:
129 Minimum value covered by the colormap.
130 :type vmin:
131 float
133 :param vmax:
134 Maximum value covered by the colormap.
135 :type vmax:
136 float
138 :param C:
139 Comma seperated R/G/B values for cmap definition.
140 :type C:
141 str
143 :param cmap:
144 Name of the colormap. Colormap is stored as "my_<cmap>.cpt".
145 If name is equivalent to a matplotlib colormap, R/G/B strings are
146 extracted from this colormap.
147 :type cmap:
148 str
150 :param space:
151 If ``True``, the range of the colormap is broadened below vmin
152 and above vmax.
153 :type space: bool
154 '''
156 scaler = AutoScaler(mode='min-max')
157 scale = scaler.make_scale((vmin, vmax))
159 incr = scale[2]
160 margin = 0.
162 if vmin == vmax:
163 space = True
165 if space:
166 margin = incr
168 msg = ('Please give either a valid color code or a'
169 ' valid matplotlib colormap name.')
171 if C is None and cmap is None:
172 raise ValueError(msg)
174 if C is None:
175 try:
176 C = _mplcmap_to_gmtcpt_code(cmap)
177 except ValueError:
178 raise ValueError(msg)
180 if cmap is None:
181 logger.warning(
182 'No colormap name given. Uses temporary filename instead')
183 cmap = 'temp_cmap'
185 return gmt.makecpt(
186 C=C,
187 D='o',
188 T='%g/%g/%g' % (
189 vmin - margin, vmax + margin, incr),
190 Z=True,
191 out_filename=make_cpt_path(gmt, cmap),
192 suppress_defaults=True)
195def xy_to_latlon(source, x, y):
196 '''
197 Convert x and y relative coordinates on extended ruptures into latlon.
199 :param source:
200 Extended source class, on which the given point is located.
201 :type source:
202 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or
203 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
205 :param x:
206 Relative point coordinate along strike (range: -1:1).
207 :type x:
208 :py:class:`float` or :py:class:`~numpy.ndarray`
210 :param y:
211 Relative downdip point coordinate (range: -1:1).
212 :type y:
213 :py:class:`float` or :py:class:`~numpy.ndarray`
215 :returns:
216 Latitude and Longitude of the given point in [deg].
217 :rtype:
218 :py:class:`tuple` of :py:class:`float`
219 '''
221 s = source
222 ax, ay = map_anchor[s.anchor]
224 length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
225 strike, dip = s.strike * d2r, s.dip * d2r
227 northing = num.cos(strike)*length - num.cos(dip) * num.sin(strike)*width
228 easting = num.sin(strike)*length + num.cos(dip) * num.cos(strike)*width
230 northing += source.north_shift
231 easting += source.east_shift
233 return pod.ne_to_latlon(s.lat, s.lon, northing, easting)
236def xy_to_lw(source, x, y):
237 '''
238 Convert relative coordinates on extended ruptures into length and width.
240 :param source:
241 Extended source, on which the given points are located.
242 :type source:
243 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or
244 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
246 :param x:
247 Relative point coordinates along strike (range: -1:1).
248 :type x:
249 float or :py:class:`~numpy.ndarray`
251 :param y:
252 Relative downdip point coordinates (range: -1:1).
253 :type y:
254 float or :py:class:`~numpy.ndarray`
256 :returns:
257 Length and downdip width of the given points from the anchor in [m].
258 :rtype:
259 :py:class:`tuple` of :py:class:`float`
260 '''
262 length, width = source.length, source.width
264 ax, ay = map_anchor[source.anchor]
266 lengths = (x - ax) / 2. * length
267 widths = (y - ay) / 2. * width
269 return lengths, widths
272cbar_anchor = {
273 'center': 'MC',
274 'center_left': 'ML',
275 'center_right': 'MR',
276 'top': 'TC',
277 'top_left': 'TL',
278 'top_right': 'TR',
279 'bottom': 'BC',
280 'bottom_left': 'BL',
281 'bottom_right': 'BR'}
284cbar_helper = {
285 'traction': {
286 'unit': 'MPa',
287 'factor': 1e-6},
288 'tx': {
289 'unit': 'MPa',
290 'factor': 1e-6},
291 'ty': {
292 'unit': 'MPa',
293 'factor': 1e-6},
294 'tz': {
295 'unit': 'MPa',
296 'factor': 1e-6},
297 'time': {
298 'unit': 's',
299 'factor': 1.},
300 'strike': {
301 'unit': '°',
302 'factor': 1.},
303 'dip': {
304 'unit': '°',
305 'factor': 1.},
306 'vr': {
307 'unit': 'km/s',
308 'factor': 1e-3},
309 'length': {
310 'unit': 'km',
311 'factor': 1e-3},
312 'width': {
313 'unit': 'km',
314 'factor': 1e-3}
315}
318fonttype = 'Helvetica'
320c2disl = dict([('x', 0), ('y', 1), ('z', 2)])
323def _make_gmt_conf(fontcolor, size):
324 '''
325 Update different gmt parameters depending on figure size and fontcolor.
327 :param fontcolor:
328 GMT readable colorcode or colorstring for the font.
329 :type fontcolor:
330 str
332 :param size:
333 Tuple of the figure size (width, height) in [cm].
334 :type size:
335 :py:class:`tuple` of :py:class:`float`
337 :returns:
338 Best fitting fontsize, GMT configuration parameter
339 :rtype:
340 float, dict
341 '''
343 color = fontcolor
344 fontsize = num.max(size)
346 font = '%gp,%s' % (fontsize, fonttype)
348 pen_size = fontsize / 16.
349 tick_size = num.min(size) / 200.
351 return fontsize, dict(
352 MAP_FRAME_PEN='%s' % color,
353 MAP_TICK_PEN_PRIMARY='%gp,%s' % (pen_size, color),
354 MAP_TICK_PEN_SECONDARY='%gp,%s' % (pen_size, color),
355 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size,
356 MAP_TICK_LENGTH_SECONDARY='%gc' % (tick_size * 3),
357 FONT_ANNOT_PRIMARY='%s-Bold,%s' % (font, color),
358 FONT_LABEL='%s-Bold,%s' % (font, color),
359 FONT_TITLE='%s-Bold,%s' % (font, color),
360 PS_CHAR_ENCODING='ISOLatin1+',
361 MAP_FRAME_TYPE='fancy',
362 FORMAT_GEO_MAP='D',
363 PS_PAGE_ORIENTATION='portrait',
364 MAP_GRID_PEN_PRIMARY='thinnest,%s' % color,
365 MAP_ANNOT_OBLIQUE='6')
368class SourceError(Exception):
369 pass
372class RuptureMap(Map):
373 '''
374 Map plotting of attributes and results of the
375 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
376 '''
378 def __init__(
379 self,
380 source=None,
381 fontcolor='darkslategrey',
382 width=20.,
383 height=14.,
384 margins=None,
385 color_wet=(216, 242, 254),
386 color_dry=(238, 236, 230),
387 topo_cpt_wet='light_sea_uniform',
388 topo_cpt_dry='light_land_uniform',
389 show_cities=False,
390 *args,
391 **kwargs):
393 size = (width, height)
394 fontsize, gmt_config = _make_gmt_conf(fontcolor, size)
396 if margins is None:
397 margins = [
398 fontsize * 0.15, num.min(size) / 200.,
399 num.min(size) / 200., fontsize * 0.05]
401 Map.__init__(self, *args, margins=margins, width=width, height=height,
402 gmt_config=gmt_config,
403 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet,
404 color_wet=color_wet, color_dry=color_dry,
405 **kwargs)
407 if show_cities:
408 self.draw_cities()
410 self._source = source
411 self._fontcolor = fontcolor
412 self._fontsize = fontsize
413 self.axes_layout = 'WeSn'
415 @property
416 def size(self):
417 '''
418 Figure size in [cm].
419 '''
420 return (self.width, self.height)
422 @property
423 def font(self):
424 '''
425 Font style (size and type).
426 '''
428 return '%sp,%s' % (self._fontsize, fonttype)
430 @property
431 def source(self):
432 '''
433 PseudoDynamicRupture whose attributes are plotted.
435 Note, that source.patches attribute needs to be calculated in advance.
436 '''
438 if self._source is None:
439 raise SourceError('No source given. Please define it!')
441 if not isinstance(self._source, PseudoDynamicRupture):
442 raise SourceError('This function is only capable for a source of'
443 ' type: %s' % type(PseudoDynamicRupture()))
445 if not self._source.patches:
446 raise TypeError('No source patches are defined. Please run'
447 '"source.discretize_patches()" on your source')
449 return self._source
451 @source.setter
452 def source(self, source):
453 self._source = source
455 def _get_topotile(self):
456 if self._dems is None:
457 self._setup()
459 try:
460 t, _ = self._get_topo_tile('land')
461 except NoTopo:
462 wesn = self.wesn
464 nx = int(num.ceil(
465 self.width * cm2inch * self.topo_resolution_max))
466 ny = int(num.ceil(
467 self.height * cm2inch * self.topo_resolution_max))
469 data = num.zeros((nx, ny))
471 t = Tile(wesn[0], wesn[2],
472 (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny,
473 data)
475 return t
477 def _patches_to_lw(self):
478 '''
479 Generate regular rect. length-width grid based on the patch distance.
481 Prerequesite is a regular grid of patches (constant lengths and widths)
482 Both coordinates are given relative to the source anchor point in [m]
483 The grid is extended from the patch centres to the edges of the source.
485 :returns:
486 Lengths along strike, widths downdip in [m].
487 :rtype:
488 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray`
489 '''
491 source = self.source
492 patches = source.patches
494 patch_l, patch_w = patches[0].length, patches[0].width
496 patch_lengths = num.concatenate((
497 num.array([0.]),
498 num.array([il*patch_l+patch_l/2. for il in range(source.nx)]),
499 num.array([patch_l * source.nx])))
501 patch_widths = num.concatenate((
502 num.array([0.]),
503 num.array([iw*patch_w+patch_w/2. for iw in range(source.ny)]),
504 num.array([patch_w * source.ny])))
506 ax, ay = map_anchor[source.anchor]
508 patch_lengths -= source.length * (ax + 1.) / 2.
509 patch_widths -= source.width * (ay + 1.) / 2.
511 return patch_lengths, patch_widths
513 def _xy_to_lw(self, x, y):
514 '''
515 Generate regular rect. length-width grid based on the xy coordinates.
517 Prerequesite is a regular grid with constant dx and dy. x and y are
518 relative coordinates on the rupture plane (range -1:1) along strike and
519 downdip.
520 Length and width are obtained relative to the source anchor point
521 in [m].
523 :returns:
524 Lengths along strike, widths downdip in [m].
525 :rtype:
526 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray`
527 '''
529 x, y = num.unique(x), num.unique(y)
530 dx, dy = x[1] - x[0], y[1] - y[0]
532 if any(num.abs(num.diff(x) - dx) >= 1e-6):
533 raise ValueError('Regular grid with constant spacing needed.'
534 ' Please check the x coordinates.')
536 if any(num.abs(num.diff(y) - dy) >= 1e-6):
537 raise ValueError('Regular grid with constant spacing needed.'
538 ' Please check the y coordinates.')
540 return xy_to_lw(self.source, x, y)
542 def _tile_to_lw(self, ref_lat, ref_lon,
543 north_shift=0., east_shift=0., strike=0.):
545 '''
546 Coordinate transformation from the topo tile grid into length-width.
548 The topotile latlon grid is rotated into the length width grid. The
549 width is defined here as its horizontal component. The resulting grid
550 is used for interpolation of grid data.
552 :param ref_lat:
553 Reference latitude, from which length-width relative coordinates
554 grid are calculated.
555 :type ref_lat:
556 float
558 :param ref_lon:
559 Reference longitude, from which length-width relative coordinates
560 grid are calculated.
561 :type ref_lon:
562 float
564 :param north_shift:
565 North shift of the reference point in [m].
566 :type north_shift:
567 float
569 :param east_shift:
570 East shift of the reference point [m].
571 :type east_shift:
572 float
574 :param strike:
575 striking of the length axis compared to the North axis in [deg].
576 :type strike:
577 float
579 :returns:
580 Topotile grid nodes as array of length-width coordinates.
581 :rtype:
582 :py:class:`~numpy.ndarray`, ``(n_nodes, 2)``
583 '''
585 t = self._get_topotile()
586 grid_lats = t.y()
587 grid_lons = t.x()
589 meshgrid_lons, meshgrid_lats = num.meshgrid(grid_lons, grid_lats)
591 grid_northing, grid_easting = pod.latlon_to_ne_numpy(
592 ref_lat, ref_lon, meshgrid_lats.flatten(), meshgrid_lons.flatten())
594 grid_northing -= north_shift
595 grid_easting -= east_shift
597 strike *= d2r
598 sin, cos = num.sin(strike), num.cos(strike)
599 points_out = num.zeros((grid_northing.shape[0], 2))
600 points_out[:, 0] = -sin * grid_northing + cos * grid_easting
601 points_out[:, 1] = cos * grid_northing + sin * grid_easting
603 return points_out
605 def _prep_patch_grid_data(self, data):
606 '''
607 Extend patch data from patch centres to the outer source edges.
609 Patch data is always defined in the centre of the patches. For
610 interpolation the data is extended here to the edges of the rupture
611 plane.
613 :param data:
614 Patch wise data.
615 :type data:
616 :py:class:`~numpy.ndarray`
618 :returns:
619 Extended data array.
620 :rtype:
621 :py:class:`~numpy.ndarray`
622 '''
624 shape = (self.source.nx + 2, self.source.ny + 2)
625 data = data.reshape(self.source.nx, self.source.ny).copy()
627 data_new = num.zeros(shape)
628 data_new[1:-1, 1:-1] = data
629 data_new[1:-1, 0] = data[:, 0]
630 data_new[1:-1, -1] = data[:, -1]
631 data_new[0, 1:-1] = data[0, :]
632 data_new[-1, 1:-1] = data[-1, :]
634 for i, j in zip([-1, -1, 0, 0], [-1, 0, -1, 0]):
635 data_new[i, j] = data[i, j]
637 return data_new
639 def _regular_data_to_grid(self, lengths, widths, data, filename):
640 '''
641 Interpolate regular data onto topotile grid.
643 Regular gridded data is interpolated onto the latlon grid of the
644 topotile. It is then stored as a gmt-readable .grd-file.
646 :param lengths:
647 Grid coordinates along strike relative to anchor in [m].
648 :type lengths:
649 :py:class:`~numpy.ndarray`
651 :param widths:
652 Grid coordinates downdip relative to anchor in [m].
653 :type widths:
654 :py:class:`~numpy.ndarray`
656 :param data:
657 Data grid array.
658 :type data:
659 :py:class:`~numpy.ndarray`
661 :param filename:
662 Filename, where grid is stored.
663 :type filename:
664 str
665 '''
667 source = self.source
669 interpolator = scrgi(
670 (widths * num.cos(d2r * source.dip), lengths),
671 data.T,
672 bounds_error=False,
673 method='nearest')
675 points_out = self._tile_to_lw(
676 ref_lat=source.lat,
677 ref_lon=source.lon,
678 north_shift=source.north_shift,
679 east_shift=source.east_shift,
680 strike=source.strike)
682 t = self._get_topotile()
683 t.data = num.zeros_like(t.data, dtype=float)
684 t.data[:] = num.nan
686 t.data = interpolator(points_out).reshape(t.data.shape)
688 _save_grid(t.y(), t.x(), t.data, filename=filename)
690 def patch_data_to_grid(self, data, *args, **kwargs):
691 '''
692 Generate a grid file based on regular patch wise data.
694 :param data:
695 Patchwise grid data.
696 :type data:
697 :py:class:`~numpy.ndarray`
698 '''
700 lengths, widths = self._patches_to_lw()
702 data_new = self._prep_patch_grid_data(data)
704 self._regular_data_to_grid(lengths, widths, data_new, *args, **kwargs)
706 def xy_data_to_grid(self, x, y, data, *args, **kwargs):
707 '''
708 Generate a grid file based on gridded data using xy coordinates.
710 Convert a grid based on relative fault coordinates (range -1:1) along
711 strike (x) and downdip (y) into a .grd file.
713 :param x:
714 Relative point coordinate along strike (range: -1:1).
715 :type x:
716 float or :py:class:`~numpy.ndarray`
718 :param y:
719 Relative downdip point coordinate (range: -1:1).
720 :type y:
721 float or :py:class:`~numpy.ndarray`
723 :param data:
724 Patchwise grid data.
725 :type data:
726 :py:class:`~numpy.ndarray`
727 '''
729 lengths, widths = self._xy_to_lw(x, y)
731 self._regular_data_to_grid(
732 lengths, widths, data.reshape((lengths.shape[0], widths.shape[0])),
733 *args, **kwargs)
735 def draw_image(self, gridfile, cmap, cbar=True, **kwargs):
736 '''
737 Draw grid data as image and include, if whished, a colorbar.
739 :param gridfile:
740 File of the grid which shall be plotted.
741 :type gridfile:
742 str
744 :param cmap:
745 Name of the colormap, which shall be used. A .cpt-file
746 "my_<cmap>.cpt" must exist.
747 :type cmap:
748 str
750 :param cbar:
751 If ``True``, a colorbar corresponding to the grid data is
752 added. Keyword arguments are parsed to it.
753 :type cbar:
754 bool
755 '''
757 self.gmt.grdimage(
758 gridfile,
759 C=make_cpt_path(self.gmt, cmap),
760 E='200',
761 Q=True,
762 n='+t0.0',
763 *self.jxyr)
765 if cbar:
766 self.draw_colorbar(cmap=cmap, **kwargs)
768 def draw_contour(
769 self,
770 gridfile,
771 contour_int,
772 anot_int,
773 angle=None,
774 unit='',
775 color='',
776 style='',
777 **kwargs):
779 '''
780 Draw grid data as contour lines.
782 :param gridfile:
783 File of the grid which shall be plotted.
784 :type gridfile:
785 str
787 :param contour_int:
788 Interval of contour lines in units of the gridfile.
789 :type contour_int:
790 float
792 :param anot_int:
793 Interval of labelled contour lines in units of the gridfile. Must
794 be a integer multiple of contour_int.
795 :type anot_int:
796 float
798 :param angle:
799 Rotation angle of the labels in [deg].
800 :type angle:
801 float
803 :param unit:
804 Name of the unit in the grid file. It will be displayed behind the
805 label on labelled contour lines.
806 :type unit:
807 str
809 :param color:
810 GMT readable color code or string of the contour lines.
811 :type color:
812 str
814 :param style:
815 Line style of the contour lines. If not given, solid lines are
816 plotted.
817 :type style:
818 str
819 '''
821 pen_size = self._fontsize / 40.
823 if not color:
824 color = self._fontcolor
826 a_string = '%g+f%s,%s+r%gc+u%s' % (
827 anot_int, self.font, color, pen_size*4, unit)
828 if angle:
829 a_string += '+a%g' % angle
830 c_string = '%g' % contour_int
832 if kwargs:
833 kwargs['A'], kwargs['C'] = a_string, c_string
834 else:
835 kwargs = dict(A=a_string, C=c_string)
837 if style:
838 style = ',' + style
840 args = ['-Wc%gp,%s%s+s' % (pen_size, color, style)]
842 self.gmt.grdcontour(
843 gridfile,
844 S='10',
845 W='a%gp,%s%s+s' % (pen_size*4, color, style),
846 *self.jxyr + args,
847 **kwargs)
849 def draw_colorbar(self, cmap, label='', anchor='top_right', **kwargs):
850 '''
851 Draw a colorbar based on a existing colormap.
853 :param cmap:
854 Name of the colormap, which shall be used. A .cpt-file
855 "my_<cmap>.cpt" must exist.
856 :type cmap:
857 str
859 :param label:
860 Title of the colorbar.
861 :type label:
862 str
864 :param anchor:
865 Placement of the colorbar. Combine 'top', 'center' and 'bottom'
866 with 'left', None for middle and 'right'.
867 :type anchor:
868 str
869 '''
871 if not kwargs:
872 kwargs = {}
874 if label:
875 kwargs['B'] = 'af+l%s' % label
877 kwargs['C'] = make_cpt_path(self.gmt, cmap)
878 a_str = cbar_anchor[anchor]
880 w = self.width / 3.
881 h = w / 10.
883 lgap = rgap = w / 10.
884 bgap, tgap = h, h / 10.
886 dx, dy = 2.5 * lgap, 2. * tgap
888 if 'bottom' in anchor:
889 dy += 4 * h
891 self.gmt.psscale(
892 D='j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy),
893 F='+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap),
894 *self.jxyr,
895 **kwargs)
897 def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs):
898 '''
899 Draw vectors based on two grid files.
901 Two grid files for vector lengths in x and y need to be given. The
902 function calls gmt.grdvector. All arguments defined for this function
903 in gmt can be passed as keyword arguments. Different standard settings
904 are applied if not defined differently.
906 :param x_gridfile:
907 File of the grid defining vector lengths in x.
908 :type x_gridfile:
909 str
911 :param y_gridfile:
912 File of the grid defining vector lengths in y.
913 :type y_gridfile:
914 str
916 :param vcolor:
917 Vector face color as defined in "G" option.
918 :type vcolor:
919 str
920 '''
922 kwargs['S'] = kwargs.get('S', 'il1.')
923 kwargs['I'] = kwargs.get('I', 'x20')
924 kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black')
925 kwargs['Q'] = kwargs.get('Q', '4c+e+n5km+h1')
927 self.gmt.grdvector(
928 x_gridfile, y_gridfile,
929 G='%s' % 'lightgrey' if not vcolor else vcolor,
930 *self.jxyr,
931 **kwargs)
933 def draw_dynamic_data(self, data, **kwargs):
934 '''
935 Draw an image of any data gridded on the patches e.g dislocation.
937 :param data:
938 Patchwise grid data.
939 :type data:
940 :py:class:`~numpy.ndarray`
941 '''
943 plot_data = data
945 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
947 clim = kwargs.pop('clim', (plot_data.min(), plot_data.max()))
949 cpt_path = make_cpt_path(self.gmt, kwargs['cmap'])
950 if not op.exists(cpt_path):
951 make_colormap(self.gmt, clim[0], clim[1],
952 cmap=kwargs['cmap'], space=False)
954 tmp_grd_file = self.gmt.tempfilename('tmpdata.grd')
955 self.patch_data_to_grid(plot_data, tmp_grd_file)
956 self.draw_image(tmp_grd_file, **kwargs)
958 def draw_patch_parameter(self, attribute, **kwargs):
959 '''
960 Draw an image of a chosen patch attribute e.g traction.
962 :param attribute:
963 Patch attribute, which is plotted. All patch attributes can be
964 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`)
965 and also ``traction``, ``tx``, ``ty`` or ``tz`` to display the
966 length or the single components of the traction vector.
967 :type attribute:
968 str
969 '''
971 a = attribute
972 source = self.source
974 if a == 'traction':
975 data = num.linalg.norm(source.get_tractions(), axis=1)
976 elif a == 'tx':
977 data = source.get_tractions()[:, 0]
978 elif a == 'ty':
979 data = source.get_tractions()[:, 1]
980 elif a == 'tz':
981 data = source.get_tractions()[:, 2]
982 else:
983 data = source.get_patch_attribute(attribute)
985 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
986 data *= factor
988 kwargs['label'] = kwargs.get(
989 'label',
990 '%s [%s]' % (a, cbar_helper[a]['unit']))
992 self.draw_dynamic_data(data, **kwargs)
994 def draw_time_contour(self, store, clevel=[], **kwargs):
995 '''
996 Draw high contour lines of the rupture front propgation time.
998 :param store:
999 Greens function store, which is used for time calculation.
1000 :type store:
1001 :py:class:`~pyrocko.gf.store.Store`
1002 :param clevel:
1003 List of times, for which contour lines are drawn.
1004 :type clevel:
1005 :py:class:`list` of :py:class:`float`
1006 '''
1008 _, _, _, _, points_xy = self.source._discretize_points(store, cs='xyz')
1009 _, _, times, _, _, _ = self.source.get_vr_time_interpolators(store)
1011 scaler = AutoScaler(mode='0-max', approx_ticks=8)
1012 scale = scaler.make_scale([num.min(times), num.max(times)])
1014 if clevel:
1015 if len(clevel) > 1:
1016 kwargs['anot_int'] = num.min(num.diff(clevel))
1017 else:
1018 kwargs['anot_int'] = clevel[0]
1020 kwargs['contour_int'] = kwargs['anot_int']
1021 kwargs['L'] = '0/%g' % num.max(clevel)
1023 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1024 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1025 kwargs['unit'] = kwargs.get('unit', cbar_helper['time']['unit'])
1026 kwargs['L'] = kwargs.get('L', '0/%g' % (num.max(times) + 1.))
1027 kwargs['G'] = kwargs.get('G', 'n2/3c')
1029 tmp_grd_file = self.gmt.tempfilename('tmpdata.grd')
1030 self.xy_data_to_grid(points_xy[:, 0], points_xy[:, 1],
1031 times, tmp_grd_file)
1032 self.draw_contour(tmp_grd_file, **kwargs)
1034 def draw_points(self, lats, lons, symbol='point', size=None, **kwargs):
1035 '''
1036 Draw points at given locations.
1038 :param lats:
1039 Point latitude coordinates in [deg].
1040 :type lats:
1041 iterable of :py:class:`float`
1043 :param lons:
1044 Point longitude coordinates in [deg].
1045 :type lons:
1046 iterable of :py:class:`float`
1048 :param symbol:
1049 Define symbol of the points (``'star', 'circle', 'point',
1050 'triangle'``) - default is ``point``.
1051 :type symbol:
1052 str
1054 :param size:
1055 Size of the points in [points].
1056 :type size:
1057 float
1058 '''
1060 sym_to_gmt = dict(
1061 star='a',
1062 circle='c',
1063 point='p',
1064 triangle='t')
1066 lats = num.atleast_1d(lats)
1067 lons = num.atleast_1d(lons)
1069 if lats.shape[0] != lons.shape[0]:
1070 raise IndexError('lats and lons do not have the same shape!')
1072 if size is None:
1073 size = self._fontsize / 3.
1075 kwargs['S'] = kwargs.get('S', sym_to_gmt[symbol] + '%gp' % size)
1076 kwargs['G'] = kwargs.get('G', gmtpy.color('scarletred2'))
1077 kwargs['W'] = kwargs.get('W', '2p,%s' % self._fontcolor)
1079 self.gmt.psxy(
1080 in_columns=[lons, lats],
1081 *self.jxyr,
1082 **kwargs)
1084 def draw_nucleation_point(self, **kwargs):
1085 '''
1086 Plot the nucleation point onto the map.
1087 '''
1089 nlat, nlon = xy_to_latlon(
1090 self.source, self.source.nucleation_x, self.source.nucleation_y)
1092 self.draw_points(nlat, nlon, **kwargs)
1094 def draw_dislocation(self, time=None, component='', **kwargs):
1095 '''
1096 Draw dislocation onto map at any time.
1098 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given
1099 component the patchwise dislocation is plotted onto the map.
1101 :param time:
1102 Time after origin, for which dislocation is computed. If ``None``,
1103 ``tmax`` is taken.
1104 :type time:
1105 float
1107 :param component:
1108 Dislocation component, which shall be plotted: ``x`` along strike,
1109 ``y`` along updip, ``z`` normal. If ``None``, the
1110 length of the dislocation vector is plotted.
1111 :type component: str
1112 '''
1114 disl = self.source.get_slip(time=time)
1116 if component:
1117 data = disl[:, c2disl[component]]
1118 else:
1119 data = num.linalg.norm(disl, axis=1)
1121 kwargs['label'] = kwargs.get(
1122 'label', 'u%s [m]' % (component))
1124 self.draw_dynamic_data(data, **kwargs)
1126 def draw_dislocation_contour(
1127 self, time=None, component=None, clevel=[], **kwargs):
1128 '''
1129 Draw dislocation contour onto map at any time.
1131 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given
1132 component the patchwise dislocation is plotted as contour onto the map.
1134 :param time:
1135 Time after origin, for which dislocation is computed. If ``None``,
1136 ``tmax`` is taken.
1137 :type time:
1138 float
1140 :param component:
1141 Dislocation component, which shall be plotted: ``x`` along strike,
1142 ``y`` along updip, ``z`` normal``. If ``None``, the length of the
1143 dislocation vector is plotted.
1144 :type component: str
1146 :param clevel:
1147 Times, for which contour lines are drawn.
1148 :type clevel:
1149 :py:class:`list` of :py:class:`float`
1150 '''
1152 disl = self.source.get_slip(time=time)
1154 if component:
1155 data = disl[:, c2disl[component]]
1156 else:
1157 data = num.linalg.norm(disl, axis=1)
1159 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1160 scale = scaler.make_scale([num.min(data), num.max(data)])
1162 if clevel:
1163 if len(clevel) > 1:
1164 kwargs['anot_int'] = num.min(num.diff(clevel))
1165 else:
1166 kwargs['anot_int'] = clevel[0]
1168 kwargs['contour_int'] = kwargs['anot_int']
1169 kwargs['L'] = '%g/%g' % (
1170 num.min(clevel) - kwargs['contour_int'],
1171 num.max(clevel) + kwargs['contour_int'])
1172 else:
1173 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1174 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1175 kwargs['L'] = kwargs.get('L', '%g/%g' % (
1176 num.min(data) - 1., num.max(data) + 1.))
1178 kwargs['unit'] = kwargs.get('unit', ' m')
1179 kwargs['G'] = kwargs.get('G', 'n2/3c')
1181 tmp_grd_file = self.gmt.tempfilename('tmpdata.grd')
1182 self.patch_data_to_grid(data, tmp_grd_file)
1183 self.draw_contour(tmp_grd_file, **kwargs)
1185 def draw_dislocation_vector(self, time=None, **kwargs):
1186 '''
1187 Draw vector arrows onto map indicating direction of dislocation.
1189 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1190 and given component the dislocation is plotted as vectors onto the map.
1192 :param time:
1193 Time after origin [s], for which dislocation is computed. If
1194 ``None``, ``tmax`` is used.
1195 :type time:
1196 float
1197 '''
1199 disl = self.source.get_slip(time=time)
1201 p_strike = self.source.get_patch_attribute('strike') * d2r
1202 p_dip = self.source.get_patch_attribute('dip') * d2r
1204 disl_dh = num.cos(p_dip) * disl[:, 1]
1205 disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh
1206 disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh
1208 tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')]
1210 self.patch_data_to_grid(disl_n, tmp_grd_files[0])
1211 self.patch_data_to_grid(disl_e, tmp_grd_files[1])
1213 self.draw_vector(
1214 tmp_grd_files[1], tmp_grd_files[0],
1215 **kwargs)
1217 def draw_top_edge(self, **kwargs):
1218 '''
1219 Indicate rupture top edge on map.
1220 '''
1222 outline = self.source.outline(cs='latlondepth')
1223 top_edge = outline[:2, :]
1225 kwargs = kwargs or {}
1226 kwargs['W'] = kwargs.get(
1227 'W', '%gp,%s' % (self._fontsize / 10., gmtpy.color('orange3')))
1229 self.gmt.psxy(
1230 in_columns=[top_edge[:, 1], top_edge[:, 0]],
1231 *self.jxyr,
1232 **kwargs)
1235class RuptureView(Object):
1236 '''
1237 Plot of attributes and results of the
1238 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1239 '''
1241 _patches_to_lw = RuptureMap._patches_to_lw
1243 def __init__(self, source=None, figsize=None, fontsize=None):
1244 self._source = source
1245 self._axes = None
1247 self.figsize = figsize or mpl_papersize('halfletter', 'landscape')
1248 self.fontsize = fontsize or 10
1249 mpl_init(fontsize=self.fontsize)
1251 self._fig = None
1252 self._axes = None
1253 self._is_1d = False
1255 @property
1256 def source(self):
1257 '''
1258 PseudoDynamicRupture whose attributes are plotted.
1260 Note, that source.patches attribute needs to be calculated for
1261 :type source: :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1262 '''
1264 if self._source is None:
1265 raise SourceError('No source given. Please define it!')
1267 if not isinstance(self._source, PseudoDynamicRupture):
1268 raise SourceError('This function is only capable for a source of'
1269 ' type: %s' % type(PseudoDynamicRupture()))
1271 if not self._source.patches:
1272 raise TypeError('No source patches are defined. Please run'
1273 '\"discretize_patches\" on your source')
1275 return self._source
1277 @source.setter
1278 def source(self, source):
1279 self._source = source
1281 def _setup(
1282 self,
1283 title='',
1284 xlabel='',
1285 ylabel='',
1286 aspect=1.,
1287 spatial_plot=True,
1288 **kwargs):
1290 if self._fig is not None and self._axes is not None:
1291 return
1293 self._fig = plt.figure(figsize=self.figsize)
1295 self._axes = self._fig.add_subplot(1, 1, 1, aspect=aspect)
1296 ax = self._axes
1298 if ax is not None:
1299 ax.set_title(title)
1300 ax.grid(alpha=.3)
1301 ax.set_xlabel(xlabel)
1302 ax.set_ylabel(ylabel)
1304 if spatial_plot:
1305 ax.xaxis.set_major_formatter(FuncFormatter(km_formatter))
1306 ax.yaxis.set_major_formatter(FuncFormatter(km_formatter))
1308 def _clear_all(self):
1309 plt.cla()
1310 plt.clf()
1311 plt.close()
1313 self._fig, self._axes = None, None
1315 def _draw_scatter(self, x, y, *args, **kwargs):
1316 default_kwargs = dict(
1317 linewidth=0,
1318 marker='o',
1319 markerfacecolor=mpl_color('skyblue2'),
1320 markersize=6.,
1321 markeredgecolor=mpl_color('skyblue3'))
1322 default_kwargs.update(kwargs)
1324 if self._axes is not None:
1325 self._axes.plot(x, y, *args, **default_kwargs)
1327 def _draw_image(self, length, width, data, *args, **kwargs):
1328 if self._axes is not None:
1329 if 'extent' not in kwargs:
1330 kwargs['extent'] = [
1331 num.min(length), num.max(length),
1332 num.max(width), num.min(width)]
1334 im = self._axes.imshow(
1335 data,
1336 *args,
1337 interpolation='none',
1338 vmin=kwargs.get('clim', [None])[0],
1339 vmax=kwargs.get('clim', [None, None])[1],
1340 **kwargs)
1342 del kwargs['extent']
1343 if 'aspect' in kwargs:
1344 del kwargs['aspect']
1345 if 'clim' in kwargs:
1346 del kwargs['clim']
1347 if 'cmap' in kwargs:
1348 del kwargs['cmap']
1350 plt.colorbar(
1351 im, *args, shrink=0.9, pad=0.03, aspect=15., **kwargs)
1353 def _draw_contour(self, x, y, data, clevel=None, unit='', *args, **kwargs):
1354 setup_kwargs = dict(
1355 xlabel=kwargs.pop('xlabel', 'along strike [km]'),
1356 ylabel=kwargs.pop('ylabel', 'down dip [km]'),
1357 title=kwargs.pop('title', ''),
1358 aspect=kwargs.pop('aspect', 1))
1360 self._setup(**setup_kwargs)
1362 if self._axes is not None:
1363 if clevel is None:
1364 scaler = AutoScaler(mode='min-max')
1365 scale = scaler.make_scale([data.min(), data.max()])
1367 clevel = num.arange(scale[0], scale[1] + scale[2], scale[2])
1369 if not isinstance(clevel, num.ndarray):
1370 clevel = num.array([clevel])
1372 clevel = clevel[clevel < data.max()]
1374 cont = self._axes.contour(
1375 x, y, data, clevel, *args,
1376 linewidths=1.5, **kwargs)
1378 clabels = self._axes.clabel(
1379 cont, clevel[::2], *args,
1380 inline=1, fmt='%g' + '%s' % unit,
1381 inline_spacing=15, rightside_up=True, use_clabeltext=True,
1382 **kwargs)
1384 plt.setp(clabels, path_effects=[
1385 patheffects.withStroke(linewidth=1.25, foreground='beige'),
1386 patheffects.Normal()])
1388 def draw_points(self, length, width, *args, **kwargs):
1389 '''
1390 Draw a point onto the figure.
1392 Args and kwargs can be defined according to
1393 :py:func:`matplotlib.pyplot.scatter`.
1395 :param length:
1396 Point(s) coordinate on the rupture plane along strike relative to
1397 the anchor point in [m].
1398 :type length:
1399 float, :py:class:`~numpy.ndarray`
1401 :param width:
1402 Point(s) coordinate on the rupture plane along downdip relative to
1403 the anchor point in [m].
1404 :type width:
1405 float, :py:class:`~numpy.ndarray`
1406 '''
1408 if self._axes is not None:
1409 kwargs['color'] = kwargs.get('color', mpl_color('scarletred2'))
1410 kwargs['s'] = kwargs.get('s', 100.)
1412 self._axes.scatter(length, width, *args, **kwargs)
1414 def draw_dynamic_data(self, data, **kwargs):
1415 '''
1416 Draw an image of any data gridded on the patches e.g dislocation.
1418 :param data:
1419 Patchwise grid data.
1420 :type data:
1421 :py:class:`~numpy.ndarray`
1422 '''
1424 plot_data = data
1426 anchor_x, anchor_y = map_anchor[self.source.anchor]
1428 length, width = xy_to_lw(
1429 self.source, num.array([-1., 1.]), num.array([-1., 1.]))
1431 data = plot_data.reshape(self.source.nx, self.source.ny).T
1433 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
1435 setup_kwargs = dict(
1436 xlabel='along strike [km]',
1437 ylabel='down dip [km]',
1438 title='',
1439 aspect=1)
1440 setup_kwargs.update(kwargs)
1442 kwargs = {k: v for k, v in kwargs.items() if k not in
1443 ('xlabel', 'ylabel', 'title')}
1444 self._setup(**setup_kwargs)
1445 self._draw_image(length=length, width=width, data=data, **kwargs)
1447 def draw_patch_parameter(self, attribute, **kwargs):
1448 '''
1449 Draw an image of a chosen patch attribute e.g traction.
1451 :param attribute:
1452 Patch attribute, which is plotted. All patch attributes can be
1453 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`)
1454 and also ``'traction', 'tx', 'ty', 'tz'`` to display the length
1455 or the single components of the traction vector.
1456 :type attribute:
1457 str
1458 '''
1460 a = attribute
1461 source = self.source
1463 if a == 'traction':
1464 data = num.linalg.norm(source.get_tractions(), axis=1)
1465 elif a == 'tx':
1466 data = source.get_tractions()[:, 0]
1467 elif a == 'ty':
1468 data = source.get_tractions()[:, 1]
1469 elif a == 'tz':
1470 data = source.get_tractions()[:, 2]
1471 else:
1472 data = source.get_patch_attribute(attribute)
1474 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1475 data *= factor
1477 kwargs['label'] = kwargs.get(
1478 'label',
1479 '%s [%s]' % (a, cbar_helper[a]['unit']))
1481 return self.draw_dynamic_data(data=data, **kwargs)
1483 def draw_time_contour(self, store, clevel=[], **kwargs):
1484 '''
1485 Draw high resolution contours of the rupture front propgation time
1487 :param store:
1488 Greens function store, which is used for time calculation.
1489 :type store:
1490 :py:class:`~pyrocko.gf.store.Store`
1492 :param clevel:
1493 Levels of the contour lines. If no levels are given, they are
1494 automatically computed based on ``tmin`` and ``tmax``.
1495 :type clevel:
1496 list
1497 '''
1499 source = self.source
1500 default_kwargs = dict(
1501 colors='#474747'
1502 )
1503 default_kwargs.update(kwargs)
1505 _, _, _, _, points_xy = source._discretize_points(store, cs='xyz')
1506 _, _, _, _, time_interpolator, _ = source.get_vr_time_interpolators(
1507 store)
1509 times = time_interpolator.values
1511 scaler = AutoScaler(mode='min-max', approx_ticks=8)
1512 scale = scaler.make_scale([times.min(), times.max()])
1514 if len(clevel) == 0:
1515 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1517 points_x = time_interpolator.grid[0]
1518 points_y = time_interpolator.grid[1]
1520 self._draw_contour(points_x, points_y, data=times.T,
1521 clevel=clevel, unit='s', **default_kwargs)
1523 def draw_nucleation_point(self, **kwargs):
1524 '''
1525 Draw the nucleation point onto the map.
1526 '''
1528 nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y
1529 length, width = xy_to_lw(self.source, nuc_x, nuc_y)
1530 self.draw_points(length, width, marker='o', **kwargs)
1532 def draw_dislocation(self, time=None, component='', **kwargs):
1533 '''
1534 Draw dislocation onto map at any time.
1536 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1537 and given component the patchwise dislocation is plotted onto the map.
1539 :param time:
1540 Time after origin [s], for which dislocation is computed.
1541 If ``None``, ``tmax`` is taken.
1542 :type time:
1543 float
1545 :param component:
1546 Dislocation component, which shall be plotted: ``x``
1547 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
1548 length of the dislocation vector is plotted
1549 :type component: str
1550 '''
1552 disl = self.source.get_slip(time=time)
1554 if component:
1555 data = disl[:, c2disl[component]]
1556 else:
1557 data = num.linalg.norm(disl, axis=1)
1559 kwargs['label'] = kwargs.get(
1560 'label', 'u%s [m]' % (component))
1562 self.draw_dynamic_data(data, **kwargs)
1564 def draw_dislocation_contour(
1565 self, time=None, component=None, clevel=[], **kwargs):
1566 '''
1567 Draw dislocation contour onto map at any time.
1569 For a given time (if time is ``None``, ``tmax`` is used) and given
1570 component the patchwise dislocation is plotted as contour onto the map.
1572 :param time:
1573 Time after origin, for which dislocation is computed. If
1574 ``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``
1581 is given, the length of the dislocation vector is plotted.
1582 :type component:
1583 str
1584 '''
1586 disl = self.source.get_slip(time=time)
1588 if component:
1589 data = disl[:, c2disl[component]]
1590 else:
1591 data = num.linalg.norm(disl, axis=1)
1593 data = data.reshape(self.source.ny, self.source.nx, order='F')
1595 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1596 scale = scaler.make_scale([data.min(), data.max()])
1598 if len(clevel) == 0:
1599 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1601 anchor_x, anchor_y = map_anchor[self.source.anchor]
1603 length, width = self._patches_to_lw()
1604 length, width = length[1:-1], width[1:-1]
1606 kwargs['colors'] = kwargs.get('colors', '#474747')
1608 self._setup(**kwargs)
1609 self._draw_contour(
1610 length, width, data=data, clevel=clevel, unit='m', **kwargs)
1612 def draw_source_dynamics(
1613 self, variable, store, deltat=None, *args, **kwargs):
1614 '''
1615 Display dynamic source parameter.
1617 Fast inspection possibility for the cumulative moment and the source
1618 time function approximation (assuming equal paths between different
1619 patches and observation point - valid for an observation point in the
1620 far field perpendicular to the source strike), so the cumulative moment
1621 rate function.
1623 :param variable:
1624 Dynamic parameter, which shall be plotted. Choose
1625 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1626 :type variable:
1627 str
1629 :param store:
1630 Greens function store, whose store.config.deltat defines
1631 the time increment between two parameter snapshots. If ``store`` is
1632 not given, the time increment is defined is taken from ``deltat``.
1633 :type store:
1634 :py:class:`~pyrocko.gf.store.Store`
1636 :param deltat:
1637 Time increment between two parameter snapshots. If not
1638 given, store.config.deltat is used to define ``deltat``.
1639 :type deltat:
1640 float
1641 '''
1643 v = variable
1645 data, times = self.source.get_moment_rate(store=store, deltat=deltat)
1646 deltat = num.concatenate([(num.diff(times)[0], ), num.diff(times)])
1648 if v in ('moment_rate', 'stf'):
1649 name, unit = 'dM/dt', 'Nm/s'
1650 elif v in ('cumulative_moment', 'moment'):
1651 data = num.cumsum(data * deltat)
1652 name, unit = 'M', 'Nm'
1653 else:
1654 raise ValueError('No dynamic data for given variable %s found' % v)
1656 self._setup(xlabel='time [s]',
1657 ylabel='%s / %.2g %s' % (name, data.max(), unit),
1658 aspect='auto',
1659 spatial_plot=False)
1660 self._draw_scatter(times, data/num.max(data), *args, **kwargs)
1661 self._is_1d = True
1663 def draw_patch_dynamics(
1664 self, variable, nx, ny, store=None, deltat=None, *args, **kwargs):
1665 '''
1666 Display dynamic boundary element / patch parameter.
1668 Fast inspection possibility for different dynamic parameter for a
1669 single patch / boundary element. The chosen parameter is plotted for
1670 the chosen patch.
1672 :param variable:
1673 Dynamic parameter, which shall be plotted. Choose
1674 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment').
1675 :type variable:
1676 str
1678 :param nx:
1679 Patch index along strike (range: 0:source.nx - 1).
1680 :type nx:
1681 int
1683 :param nx:
1684 Patch index downdip (range: 0:source.ny - 1).
1685 :type nx:
1686 int
1688 :param store:
1689 Greens function store, whose store.config.deltat defines
1690 the time increment between two parameter snapshots. If ``store`` is
1691 not given, the time increment is defined is taken from ``deltat``.
1692 :type store:
1693 :py:class:`~pyrocko.gf.store.Store`
1695 :param deltat:
1696 Time increment between two parameter snapshots. If not given,
1697 store.config.deltat is used to define ``deltat``.
1698 :type deltat:
1699 float
1700 '''
1702 v = variable
1703 source = self.source
1704 idx = nx * source.ny + nx
1706 m = re.match(r'dislocation_([xyz])', v)
1708 if v in ('moment_rate', 'cumulative_moment', 'moment', 'stf'):
1709 data, times = source.get_moment_rate_patches(
1710 store=store, deltat=deltat)
1711 elif 'dislocation' in v or 'slip_rate' == v:
1712 ddisloc, times = source.get_delta_slip(store=store, deltat=deltat)
1714 if v in ('moment_rate', 'stf'):
1715 data, times = source.get_moment_rate_patches(
1716 store=store, deltat=deltat)
1717 name, unit = 'dM/dt', 'Nm/s'
1718 elif v in ('cumulative_moment', 'moment'):
1719 data, times = source.get_moment_rate_patches(
1720 store=store, deltat=deltat)
1721 data = num.cumsum(data, axis=1)
1722 name, unit = 'M', 'Nm'
1723 elif v == 'slip_rate':
1724 data, times = source.get_slip_rate(store=store, deltat=deltat)
1725 name, unit = 'du/dt', 'm/s'
1726 elif v == 'dislocation':
1727 data, times = source.get_delta_slip(store=store, deltat=deltat)
1728 data = num.linalg.norm(num.cumsum(data, axis=2), axis=1)
1729 name, unit = 'du', 'm'
1730 elif m:
1731 data, times = source.get_delta_slip(store=store, deltat=deltat)
1732 data = num.cumsum(data, axis=2)[:, c2disl[m.group(1)], :]
1733 name, unit = 'du%s' % m.group(1), 'm'
1734 else:
1735 raise ValueError('No dynamic data for given variable %s found' % v)
1737 self._setup(xlabel='time [s]',
1738 ylabel='%s / %.2g %s' % (name, num.max(data), unit),
1739 aspect='auto',
1740 spatial_plot=False)
1741 self._draw_scatter(times, data[idx, :]/num.max(data),
1742 *args, **kwargs)
1743 self._is_1d = True
1745 def finalize(self):
1746 if self._is_1d:
1747 return
1749 length, width = xy_to_lw(
1750 self.source, num.array([-1., 1.]), num.array([1., -1.]))
1752 self._axes.set_xlim(length)
1753 self._axes.set_ylim(width)
1755 def gcf(self):
1756 self.finalize()
1757 return self._fig
1759 def save(self, filename, dpi=None):
1760 '''
1761 Save plot to file.
1763 :param filename:
1764 Filename and path, where the plot is stored.
1765 :type filename:
1766 str
1768 :param dpi:
1769 Resolution of the output plot in [dpi].
1770 :type dpi:
1771 int
1772 '''
1774 self.finalize()
1775 try:
1776 self._fig.savefig(fname=filename, dpi=dpi, bbox_inches='tight')
1777 except TypeError:
1778 self._fig.savefig(filename=filename, dpi=dpi, bbox_inches='tight')
1780 self._clear_all()
1782 def show_plot(self):
1783 '''
1784 Show plot.
1785 '''
1787 self.finalize()
1788 plt.show()
1789 self._clear_all()
1792def render_movie(fn_path, output_path, framerate=20):
1793 '''
1794 Generate a mp4 movie based on given png files using
1795 `ffmpeg <https://ffmpeg.org>`_.
1797 Render a movie based on a set of given .png files in fn_path. All files
1798 must have a filename specified by ``fn_path`` (e.g. giving ``fn_path`` with
1799 ``/temp/f%04.png`` a valid png filename would be ``/temp/f0001.png``). The
1800 files must have a numbering, indicating their order in the movie.
1802 :param fn_path:
1803 Path and fileformat specification of the input .png files.
1804 :type fn_path:
1805 str
1807 :param output_path:
1808 Path and filename of the output ``.mp4`` movie file.
1809 :type output_path:
1810 str
1812 :param deltat:
1813 Time between individual frames (``1 / framerate``) in [s].
1814 :type deltat:
1815 float
1816 '''
1818 try:
1819 check_call(['ffmpeg', '-loglevel', 'panic'])
1820 except CalledProcessError:
1821 pass
1822 except (TypeError, IOError):
1823 logger.warning(
1824 'Package ffmpeg needed for movie rendering. Please install it '
1825 '(e.g. on linux distr. via sudo apt-get ffmpeg) and retry.')
1826 return False
1828 try:
1829 check_call([
1830 'ffmpeg', '-loglevel', 'info', '-y',
1831 '-framerate', '%g' % framerate,
1832 '-i', fn_path,
1833 '-vcodec', 'libx264',
1834 '-preset', 'medium',
1835 '-tune', 'animation',
1836 '-pix_fmt', 'yuv420p',
1837 '-movflags', '+faststart',
1838 '-filter:v', "crop='floor(in_w/2)*2:floor(in_h/2)*2'",
1839 '-crf', '15',
1840 output_path])
1842 return True
1843 except CalledProcessError as e:
1844 logger.warning(e)
1845 return False
1848def render_gif(fn, output_path, loops=-1):
1849 '''
1850 Generate a gif based on a given mp4 using ffmpeg.
1852 Render a gif based on a given .mp4 movie file in ``fn`` path.
1854 :param fn:
1855 Path and file name of the input .mp4 file.
1856 :type fn:
1857 str
1859 :param output_path:
1860 Path and filename of the output animated gif file.
1861 :type output_path:
1862 str
1863 :param loops:
1864 Number of gif repeats (loops). ``-1`` is not repetition, ``0``
1865 infinite.
1866 :type loops:
1867 int
1868 '''
1870 try:
1871 check_call(['ffmpeg', '-loglevel', 'panic'])
1872 except CalledProcessError:
1873 pass
1874 except (TypeError, IOError):
1875 logger.warning(
1876 'Package ffmpeg needed for movie rendering. Please install it '
1877 '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.')
1878 return False
1880 try:
1881 check_call([
1882 'ffmpeg', '-hide_banner', '-loglevel', 'panic', '-y', '-i',
1883 fn,
1884 '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse',
1885 '-loop', '%d' % loops,
1886 output_path])
1888 return True
1889 except CalledProcessError as e:
1890 logger.warning(e)
1891 return False
1894def rupture_movie(
1895 source,
1896 store,
1897 variable='dislocation',
1898 draw_time_contours=False,
1899 fn_path='.',
1900 prefix='',
1901 plot_type='map',
1902 deltat=None,
1903 framerate=None,
1904 store_images=False,
1905 render_as_gif=False,
1906 gif_loops=-1,
1907 **kwargs):
1908 '''
1909 Generate a movie based on a given source for dynamic parameter.
1911 Create a MPEG-4 movie or gif of one of the following dynamic parameters
1912 (``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1913 (along updip), ``dislocation_z`` (normal), ``slip_rate``, ``moment_rate``).
1914 If desired, the single snap shots can be stored as images as well.
1915 ``kwargs`` have to be given according to the chosen ``plot_type``.
1917 :param source:
1918 Pseudo dynamic rupture, for which the movie is produced.
1919 :type source:
1920 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
1922 :param store:
1923 Greens function store, which is used for time calculation. If
1924 ``deltat`` is not given, it is taken from the store.config.deltat
1925 :type store:
1926 :py:class:`~pyrocko.gf.store.Store`
1928 :param variable:
1929 Dynamic parameter, which shall be plotted. Choose between
1930 ``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1931 (along updip), ``dislocation_z`` (normal), ``slip_rate`` and
1932 ``moment_rate``, default ``dislocation``.
1933 :type variable:
1934 str
1936 :param draw_time_contours:
1937 If ``True``, corresponding isochrones are drawn on the each plots.
1938 :type draw_time_contours:
1939 bool
1941 :param fn_path:
1942 Absolut or relative path, where movie (and optional images) are stored.
1943 :type fn_path:
1944 str
1946 :param prefix:
1947 File prefix used for the movie (and optional image) files.
1948 :type prefix:
1949 str
1951 :param plot_type:
1952 Choice of plot type: ``map``, ``view`` (map plot using
1953 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap`
1954 or plane view using
1955 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`).
1956 :type plot_type:
1957 str
1959 :param deltat:
1960 Time between parameter snapshots. If not given, store.config.deltat is
1961 used to define ``deltat``.
1962 :type deltat:
1963 float
1965 :param store_images:
1966 Choice to store the single .png parameter snapshots in ``fn_path`` or
1967 not.
1968 :type store_images:
1969 bool
1971 :param render_as_gif:
1972 If ``True``, the movie is converted into a gif. If ``False``, the movie
1973 is returned as mp4.
1974 :type render_as_gif:
1975 bool
1977 :param gif_loops:
1978 If ``render_as_gif`` is ``True``, a gif with ``gif_loops`` number of
1979 loops (repetitions) is returned. ``-1`` is no repetition, ``0``
1980 infinite.
1981 :type gif_loops:
1982 int
1983 '''
1985 v = variable
1986 assert plot_type in ('map', 'view')
1988 if not source.patches:
1989 source.discretize_patches(store, interpolation='multilinear')
1991 if source.coef_mat is None:
1992 source.calc_coef_mat()
1994 prefix = prefix or v
1995 deltat = deltat or store.config.deltat
1996 framerate = max(framerate or int(1./deltat), 1)
1998 if v == 'moment_rate':
1999 data, times = source.get_moment_rate_patches(deltat=deltat)
2000 name, unit = 'dM/dt', 'Nm/s'
2001 elif 'dislocation' in v or 'slip_rate' == v:
2002 ddisloc, times = source.get_delta_slip(deltat=deltat)
2003 else:
2004 raise ValueError('No dynamic data for given variable %s found' % v)
2006 deltat = times[1] - times[0]
2008 m = re.match(r'dislocation_([xyz])', v)
2009 if m:
2010 data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]]
2011 name, unit = 'du%s' % m.group(1), 'm'
2012 elif v == 'dislocation':
2013 data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2)
2014 name, unit = 'du', 'm'
2015 elif v == 'slip_rate':
2016 data = num.linalg.norm(ddisloc, axis=2) / deltat
2017 name, unit = 'du/dt', 'm/s'
2019 if plot_type == 'map':
2020 plt_base = RuptureMap
2021 elif plot_type == 'view':
2022 plt_base = RuptureView
2023 else:
2024 raise AttributeError('invalid type: %s' % plot_type)
2026 attrs_base = get_kwargs(plt_base)
2027 kwargs_base = dict([k for k in kwargs.items() if k[0] in attrs_base])
2028 kwargs_plt = dict([k for k in kwargs.items() if k[0] not in kwargs_base])
2030 if 'clim' in kwargs_plt:
2031 data = num.clip(data, kwargs_plt['clim'][0], kwargs_plt['clim'][1])
2032 else:
2033 kwargs_plt['clim'] = [num.min(data), num.max(data)]
2035 if 'label' not in kwargs_plt:
2036 vmax = num.max(num.abs(kwargs_plt['clim']))
2037 data /= vmax
2039 kwargs_plt['label'] = '%s / %.2g %s' % (name, vmax, unit)
2040 kwargs_plt['clim'] = [i / vmax for i in kwargs_plt['clim']]
2042 with tempfile.TemporaryDirectory(suffix='pyrocko') as temp_path:
2043 fns_temp = [op.join(temp_path, 'f%09d.png' % (it + 1))
2044 for it, _ in enumerate(times)]
2045 fn_temp_path = op.join(temp_path, 'f%09d.png')
2047 for it, (t, ft) in enumerate(zip(times, fns_temp)):
2048 plot_data = data[:, it]
2050 plt = plt_base(source=source, **kwargs_base)
2051 plt.draw_dynamic_data(plot_data, **kwargs_plt)
2052 plt.draw_nucleation_point()
2054 if draw_time_contours:
2055 plt.draw_time_contour(store, clevel=[t])
2057 plt.save(ft)
2059 fn_mp4 = op.join(temp_path, 'movie.mp4')
2060 return_code = render_movie(
2061 fn_temp_path,
2062 output_path=fn_mp4,
2063 framerate=framerate)
2065 if render_as_gif and return_code:
2066 render_gif(fn=fn_mp4, output_path=op.join(
2067 fn_path, '%s_%s_gif.gif' % (prefix, plot_type)),
2068 loops=gif_loops)
2070 elif return_code:
2071 shutil.move(fn_mp4, op.join(
2072 fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
2074 else:
2075 logger.error('ffmpeg failed. Exit')
2077 if store_images:
2078 fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t))
2079 for t in times]
2081 for ft, f in zip(fns_temp, fns):
2082 shutil.move(ft, f)
2085__all__ = [
2086 'make_colormap',
2087 'xy_to_latlon',
2088 'xy_to_lw',
2089 'SourceError',
2090 'RuptureMap',
2091 'RuptureView',
2092 'rupture_movie',
2093 'render_movie']