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 cm, 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, gmtpy
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 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 = cm.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_colormap(
115 gmt,
116 vmin,
117 vmax,
118 C=None,
119 cmap=None,
120 space=False):
121 '''
122 Create gmt-readable colormap cpt file called my_<cmap>.cpt.
124 :param vmin:
125 Minimum value covered by the colormap.
126 :type vmin:
127 float
129 :param vmax:
130 Maximum value covered by the colormap.
131 :type vmax:
132 float
134 :param C:
135 Comma seperated R/G/B values for cmap definition.
136 :type C:
137 optional, str
139 :param cmap:
140 Name of the colormap. Colormap is stored as "my_<cmap>.cpt".
141 If name is equivalent to a matplotlib colormap, R/G/B strings are
142 extracted from this colormap.
143 :type cmap:
144 optional, str
146 :param space:
147 If ``True``, the range of the colormap is broadened below vmin
148 and above vmax.
149 :type space: optional, bool
150 '''
152 scaler = AutoScaler(mode='min-max')
153 scale = scaler.make_scale((vmin, vmax))
155 incr = scale[2]
156 margin = 0.
158 if vmin == vmax:
159 space = True
161 if space:
162 margin = incr
164 msg = ('Please give either a valid color code or a'
165 ' valid matplotlib colormap name.')
167 if C is None and cmap is None:
168 raise ValueError(msg)
170 if C is None:
171 try:
172 C = _mplcmap_to_gmtcpt_code(cmap)
173 except ValueError:
174 raise ValueError(msg)
176 if cmap is None:
177 logger.warning(
178 'No colormap name given. Uses temporary filename instead')
179 cmap = 'temp_cmap'
181 return gmt.makecpt(
182 C=C,
183 D='o',
184 T='%g/%g/%g' % (
185 vmin - margin, vmax + margin, incr),
186 Z=True,
187 out_filename='my_%s.cpt' % cmap,
188 suppress_defaults=True)
191def clear_temp(gridfiles=[], cpts=[]):
192 '''
193 Clear all temporary needed grid and colormap cpt files.
195 :param gridfiles:
196 List of all "...grd" files, which shall be deleted.
197 :type gridfiles:
198 optional, list
200 :param cpts:
201 Cmaps, whose corresponding "my_<cmap>.cpt" file shall be deleted.
202 :type cpts:
203 optional, list
204 '''
206 for fil in gridfiles:
207 try:
208 os.remove(fil)
209 except OSError:
210 continue
211 for fil in cpts:
212 try:
213 os.remove('my_%s.cpt' % fil)
214 except OSError:
215 continue
218def xy_to_latlon(source, x, y):
219 '''
220 Convert x and y relative coordinates on extended ruptures into latlon.
222 :param source:
223 Extended source class, on which the given point is located.
224 :type source:
225 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or
226 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
228 :param x:
229 Relative point coordinate along strike (range: -1:1).
230 :type x:
231 float or :py:class:`~numpy.ndarray`
233 :param y:
234 Relative downdip point coordinate (range: -1:1).
235 :type y:
236 float or :py:class:`~numpy.ndarray`
238 :returns:
239 Latitude and Longitude of the given point in [deg].
240 :rtype:
241 tuple of float
242 '''
244 s = source
245 ax, ay = map_anchor[s.anchor]
247 length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
248 strike, dip = s.strike * d2r, s.dip * d2r
250 northing = num.cos(strike)*length - num.cos(dip) * num.sin(strike)*width
251 easting = num.sin(strike)*length + num.cos(dip) * num.cos(strike)*width
253 northing += source.north_shift
254 easting += source.east_shift
256 return pod.ne_to_latlon(s.lat, s.lon, northing, easting)
259def xy_to_lw(source, x, y):
260 '''
261 Convert relative coordinates on extended ruptures into length and width.
263 :param source:
264 Extended source, on which the given points are located.
265 :type source:
266 :py:class:`~pyrocko.gf.seismosizer.RectangularSource` or
267 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
269 :param x:
270 Relative point coordinates along strike (range: -1:1).
271 :type x:
272 float or :py:class:`~numpy.ndarray`
274 :param y:
275 Relative downdip point coordinates (range: -1:1).
276 :type y:
277 float or :py:class:`~numpy.ndarray`
279 :returns:
280 Length and downdip width of the given points from the anchor in [m].
281 :rtype:
282 tuple of float
283 '''
285 length, width = source.length, source.width
287 ax, ay = map_anchor[source.anchor]
289 lengths = (x - ax) / 2. * length
290 widths = (y - ay) / 2. * width
292 return lengths, widths
295cbar_anchor = {
296 'center': 'MC',
297 'center_left': 'ML',
298 'center_right': 'MR',
299 'top': 'TC',
300 'top_left': 'TL',
301 'top_right': 'TR',
302 'bottom': 'BC',
303 'bottom_left': 'BL',
304 'bottom_right': 'BR'}
307cbar_helper = {
308 'traction': {
309 'unit': 'MPa',
310 'factor': 1e-6},
311 'tx': {
312 'unit': 'MPa',
313 'factor': 1e-6},
314 'ty': {
315 'unit': 'MPa',
316 'factor': 1e-6},
317 'tz': {
318 'unit': 'MPa',
319 'factor': 1e-6},
320 'time': {
321 'unit': 's',
322 'factor': 1.},
323 'strike': {
324 'unit': '°',
325 'factor': 1.},
326 'dip': {
327 'unit': '°',
328 'factor': 1.},
329 'vr': {
330 'unit': 'km/s',
331 'factor': 1e-3},
332 'length': {
333 'unit': 'km',
334 'factor': 1e-3},
335 'width': {
336 'unit': 'km',
337 'factor': 1e-3}
338}
341fonttype = 'Helvetica'
343c2disl = dict([('x', 0), ('y', 1), ('z', 2)])
346def _make_gmt_conf(fontcolor, size):
347 '''
348 Update different gmt parameters depending on figure size and fontcolor.
350 :param fontcolor:
351 GMT readable colorcode or colorstring for the font.
352 :type fontcolor:
353 str
355 :param size:
356 Tuple of the figure size (width, height) in [cm].
357 :type size:
358 tuple of float
360 :returns:
361 Best fitting fontsize, GMT configuration parameter
362 :rtype:
363 float, dict
364 '''
366 color = fontcolor
367 fontsize = num.max(size)
369 font = '%gp,%s' % (fontsize, fonttype)
371 pen_size = fontsize / 16.
372 tick_size = num.min(size) / 200.
374 return fontsize, dict(
375 MAP_FRAME_PEN='%s' % color,
376 MAP_TICK_PEN_PRIMARY='%gp,%s' % (pen_size, color),
377 MAP_TICK_PEN_SECONDARY='%gp,%s' % (pen_size, color),
378 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size,
379 MAP_TICK_LENGTH_SECONDARY='%gc' % (tick_size * 3),
380 FONT_ANNOT_PRIMARY='%s-Bold,%s' % (font, color),
381 FONT_LABEL='%s-Bold,%s' % (font, color),
382 FONT_TITLE='%s-Bold,%s' % (font, color),
383 PS_CHAR_ENCODING='ISOLatin1+',
384 MAP_FRAME_TYPE='fancy',
385 FORMAT_GEO_MAP='D',
386 PS_PAGE_ORIENTATION='portrait',
387 MAP_GRID_PEN_PRIMARY='thinnest,%s' % color,
388 MAP_ANNOT_OBLIQUE='6')
391class SourceError(Exception):
392 pass
395class RuptureMap(Map):
396 '''
397 Map plotting of attributes and results of the
398 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
399 '''
401 def __init__(
402 self,
403 source=None,
404 fontcolor='darkslategrey',
405 width=20.,
406 height=14.,
407 margins=None,
408 color_wet=(216, 242, 254),
409 color_dry=(238, 236, 230),
410 topo_cpt_wet='light_sea_uniform',
411 topo_cpt_dry='light_land_uniform',
412 show_cities=False,
413 *args,
414 **kwargs):
416 size = (width, height)
417 fontsize, gmt_config = _make_gmt_conf(fontcolor, size)
419 if margins is None:
420 margins = [
421 fontsize * 0.15, num.min(size) / 200.,
422 num.min(size) / 200., fontsize * 0.05]
424 Map.__init__(self, *args, margins=margins, width=width, height=height,
425 gmt_config=gmt_config,
426 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet,
427 color_wet=color_wet, color_dry=color_dry,
428 **kwargs)
430 if show_cities:
431 self.draw_cities()
433 self._source = source
434 self._fontcolor = fontcolor
435 self._fontsize = fontsize
436 self.axes_layout = 'WeSn'
438 @property
439 def size(self):
440 '''
441 Figure size in [cm].
442 '''
443 return (self.width, self.height)
445 @property
446 def font(self):
447 '''
448 Font style (size and type).
449 '''
451 return '%sp,%s' % (self._fontsize, fonttype)
453 @property
454 def source(self):
455 '''
456 PseudoDynamicRupture whose attributes are plotted.
458 Note, that source.patches attribute needs to be calculated in advance.
459 '''
461 if self._source is None:
462 raise SourceError('No source given. Please define it!')
464 if not isinstance(self._source, PseudoDynamicRupture):
465 raise SourceError('This function is only capable for a source of'
466 ' type: %s' % type(PseudoDynamicRupture()))
468 if not self._source.patches:
469 raise TypeError('No source patches are defined. Please run'
470 '"source.discretize_patches()" on your source')
472 return self._source
474 @source.setter
475 def source(self, source):
476 self._source = source
478 def _get_topotile(self):
479 if self._dems is None:
480 self._setup()
482 try:
483 t, _ = self._get_topo_tile('land')
484 except NoTopo:
485 wesn = self.wesn
487 nx = int(num.ceil(
488 self.width * cm2inch * self.topo_resolution_max))
489 ny = int(num.ceil(
490 self.height * cm2inch * self.topo_resolution_max))
492 data = num.zeros((nx, ny))
494 t = Tile(wesn[0], wesn[2],
495 (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny,
496 data)
498 return t
500 def _patches_to_lw(self):
501 '''
502 Generate regular rect. length-width grid based on the patch distance.
504 Prerequesite is a regular grid of patches (constant lengths and widths)
505 Both coordinates are given relative to the source anchor point in [m]
506 The grid is extended from the patch centres to the edges of the source.
508 :returns:
509 Lengths along strike, widths downdip in [m].
510 :rtype:
511 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray`
512 '''
514 source = self.source
515 patches = source.patches
517 patch_l, patch_w = patches[0].length, patches[0].width
519 patch_lengths = num.concatenate((
520 num.array([0.]),
521 num.array([il*patch_l+patch_l/2. for il in range(source.nx)]),
522 num.array([patch_l * source.nx])))
524 patch_widths = num.concatenate((
525 num.array([0.]),
526 num.array([iw*patch_w+patch_w/2. for iw in range(source.ny)]),
527 num.array([patch_w * source.ny])))
529 ax, ay = map_anchor[source.anchor]
531 patch_lengths -= source.length * (ax + 1.) / 2.
532 patch_widths -= source.width * (ay + 1.) / 2.
534 return patch_lengths, patch_widths
536 def _xy_to_lw(self, x, y):
537 '''
538 Generate regular rect. length-width grid based on the xy coordinates.
540 Prerequesite is a regular grid with constant dx and dy. x and y are
541 relative coordinates on the rupture plane (range -1:1) along strike and
542 downdip.
543 Length and width are obtained relative to the source anchor point
544 in [m].
546 :returns:
547 Lengths along strike, widths downdip in [m].
548 :rtype:
549 :py:class:`~numpy.ndarray`, :py:class:`~numpy.ndarray`
550 '''
552 x, y = num.unique(x), num.unique(y)
553 dx, dy = x[1] - x[0], y[1] - y[0]
555 if any(num.abs(num.diff(x) - dx) >= 1e-6):
556 raise ValueError('Regular grid with constant spacing needed.'
557 ' Please check the x coordinates.')
559 if any(num.abs(num.diff(y) - dy) >= 1e-6):
560 raise ValueError('Regular grid with constant spacing needed.'
561 ' Please check the y coordinates.')
563 return xy_to_lw(self.source, x, y)
565 def _tile_to_lw(self, ref_lat, ref_lon,
566 north_shift=0., east_shift=0., strike=0.):
568 '''
569 Coordinate transformation from the topo tile grid into length-width.
571 The topotile latlon grid is rotated into the length width grid. The
572 width is defined here as its horizontal component. The resulting grid
573 is used for interpolation of grid data.
575 :param ref_lat:
576 Reference latitude, from which length-width relative coordinates
577 grid are calculated.
578 :type ref_lat:
579 float
581 :param ref_lon:
582 Reference longitude, from which length-width relative coordinates
583 grid are calculated.
584 :type ref_lon:
585 float
587 :param north_shift:
588 North shift of the reference point in [m].
589 :type north_shift:
590 optional, float
592 :param east_shift:
593 East shift of the reference point [m].
594 :type east_shift:
595 optional, float
597 :param strike:
598 striking of the length axis compared to the North axis in [deg].
599 :type strike:
600 optional, float
602 :returns:
603 Topotile grid nodes as array of length-width coordinates.
604 :rtype:
605 :py:class:`~numpy.ndarray`, ``(n_nodes, 2)``
606 '''
608 t = self._get_topotile()
609 grid_lats = t.y()
610 grid_lons = t.x()
612 meshgrid_lons, meshgrid_lats = num.meshgrid(grid_lons, grid_lats)
614 grid_northing, grid_easting = pod.latlon_to_ne_numpy(
615 ref_lat, ref_lon, meshgrid_lats.flatten(), meshgrid_lons.flatten())
617 grid_northing -= north_shift
618 grid_easting -= east_shift
620 strike *= d2r
621 sin, cos = num.sin(strike), num.cos(strike)
622 points_out = num.zeros((grid_northing.shape[0], 2))
623 points_out[:, 0] = -sin * grid_northing + cos * grid_easting
624 points_out[:, 1] = cos * grid_northing + sin * grid_easting
626 return points_out
628 def _prep_patch_grid_data(self, data):
629 '''
630 Extend patch data from patch centres to the outer source edges.
632 Patch data is always defined in the centre of the patches. For
633 interpolation the data is extended here to the edges of the rupture
634 plane.
636 :param data:
637 Patch wise data.
638 :type data:
639 :py:class:`~numpy.ndarray`
641 :returns:
642 Extended data array.
643 :rtype:
644 :py:class:`~numpy.ndarray`
645 '''
647 shape = (self.source.nx + 2, self.source.ny + 2)
648 data = data.reshape(self.source.nx, self.source.ny).copy()
650 data_new = num.zeros(shape)
651 data_new[1:-1, 1:-1] = data
652 data_new[1:-1, 0] = data[:, 0]
653 data_new[1:-1, -1] = data[:, -1]
654 data_new[0, 1:-1] = data[0, :]
655 data_new[-1, 1:-1] = data[-1, :]
657 for i, j in zip([-1, -1, 0, 0], [-1, 0, -1, 0]):
658 data_new[i, j] = data[i, j]
660 return data_new
662 def _regular_data_to_grid(self, lengths, widths, data, filename):
663 '''
664 Interpolate regular data onto topotile grid.
666 Regular gridded data is interpolated onto the latlon grid of the
667 topotile. It is then stored as a gmt-readable .grd-file.
669 :param lengths:
670 Grid coordinates along strike relative to anchor in [m].
671 :type lengths:
672 :py:class:`~numpy.ndarray`
674 :param widths:
675 Grid coordinates downdip relative to anchor in [m].
676 :type widths:
677 :py:class:`~numpy.ndarray`
679 :param data:
680 Data grid array.
681 :type data:
682 :py:class:`~numpy.ndarray`
684 :param filename:
685 Filename, where grid is stored.
686 :type filename:
687 str
688 '''
690 source = self.source
692 interpolator = scrgi(
693 (widths * num.cos(d2r * source.dip), lengths),
694 data.T,
695 bounds_error=False,
696 method='nearest')
698 points_out = self._tile_to_lw(
699 ref_lat=source.lat,
700 ref_lon=source.lon,
701 north_shift=source.north_shift,
702 east_shift=source.east_shift,
703 strike=source.strike)
705 t = self._get_topotile()
706 t.data = num.zeros_like(t.data, dtype=float)
707 t.data[:] = num.nan
709 t.data = interpolator(points_out).reshape(t.data.shape)
711 _save_grid(t.y(), t.x(), t.data, filename=filename)
713 def patch_data_to_grid(self, data, *args, **kwargs):
714 '''
715 Generate a grid file based on regular patch wise data.
717 :param data:
718 Patchwise grid data.
719 :type data:
720 :py:class:`~numpy.ndarray`
721 '''
723 lengths, widths = self._patches_to_lw()
725 data_new = self._prep_patch_grid_data(data)
727 self._regular_data_to_grid(lengths, widths, data_new, *args, **kwargs)
729 def xy_data_to_grid(self, x, y, data, *args, **kwargs):
730 '''
731 Generate a grid file based on gridded data using xy coordinates.
733 Convert a grid based on relative fault coordinates (range -1:1) along
734 strike (x) and downdip (y) into a .grd file.
736 :param x:
737 Relative point coordinate along strike (range: -1:1).
738 :type x:
739 float or :py:class:`~numpy.ndarray`
741 :param y:
742 Relative downdip point coordinate (range: -1:1).
743 :type y:
744 float or :py:class:`~numpy.ndarray`
746 :param data:
747 Patchwise grid data.
748 :type data:
749 :py:class:`~numpy.ndarray`
750 '''
752 lengths, widths = self._xy_to_lw(x, y)
754 self._regular_data_to_grid(
755 lengths, widths, data.reshape((lengths.shape[0], widths.shape[0])),
756 *args, **kwargs)
758 def draw_image(self, gridfile, cmap, cbar=True, **kwargs):
759 '''
760 Draw grid data as image and include, if whished, a colorbar.
762 :param gridfile:
763 File of the grid which shall be plotted.
764 :type gridfile:
765 str
767 :param cmap:
768 Name of the colormap, which shall be used. A .cpt-file
769 "my_<cmap>.cpt" must exist.
770 :type cmap:
771 str
773 :param cbar:
774 If ``True``, a colorbar corresponding to the grid data is
775 added. Keyword arguments are parsed to it.
776 :type cbar:
777 optional, bool
778 '''
780 self.gmt.grdimage(
781 gridfile,
782 C='my_%s.cpt' % cmap,
783 E='200',
784 Q=True,
785 n='+t0.0',
786 *self.jxyr)
788 if cbar:
789 self.draw_colorbar(cmap=cmap, **kwargs)
791 def draw_contour(
792 self,
793 gridfile,
794 contour_int,
795 anot_int,
796 angle=None,
797 unit='',
798 color='',
799 style='',
800 **kwargs):
802 '''
803 Draw grid data as contour lines.
805 :param gridfile:
806 File of the grid which shall be plotted.
807 :type gridfile:
808 str
810 :param contour_int:
811 Interval of contour lines in units of the gridfile.
812 :type contour_int:
813 float
815 :param anot_int:
816 Interval of labelled contour lines in units of the gridfile. Must
817 be a integer multiple of contour_int.
818 :type anot_int:
819 float
821 :param angle:
822 Rotation angle of the labels in [deg].
823 :type angle:
824 optional, float
826 :param unit:
827 Name of the unit in the grid file. It will be displayed behind the
828 label on labelled contour lines.
829 :type unit:
830 optional, str
832 :param color:
833 GMT readable color code or string of the contour lines.
834 :type color:
835 optional, str
837 :param style:
838 Line style of the contour lines. If not given, solid lines are
839 plotted.
840 :type style:
841 optional, str
842 '''
844 pen_size = self._fontsize / 40.
846 if not color:
847 color = self._fontcolor
849 a_string = '%g+f%s,%s+r%gc+u%s' % (
850 anot_int, self.font, color, pen_size*4, unit)
851 if angle:
852 a_string += '+a%g' % angle
853 c_string = '%g' % contour_int
855 if kwargs:
856 kwargs['A'], kwargs['C'] = a_string, c_string
857 else:
858 kwargs = dict(A=a_string, C=c_string)
860 if style:
861 style = ',' + style
863 args = ['-Wc%gp,%s%s+s' % (pen_size, color, style)]
865 self.gmt.grdcontour(
866 gridfile,
867 S='10',
868 W='a%gp,%s%s+s' % (pen_size*4, color, style),
869 *self.jxyr + args,
870 **kwargs)
872 def draw_colorbar(self, cmap, label='', anchor='top_right', **kwargs):
873 '''
874 Draw a colorbar based on a existing colormap.
876 :param cmap:
877 Name of the colormap, which shall be used. A .cpt-file
878 "my_<cmap>.cpt" must exist.
879 :type cmap:
880 str
882 :param label:
883 Title of the colorbar.
884 :type label:
885 optional, str
887 :param anchor:
888 Placement of the colorbar. Combine 'top', 'center' and 'bottom'
889 with 'left', None for middle and 'right'.
890 :type anchor:
891 optional, str
892 '''
894 if not kwargs:
895 kwargs = {}
897 if label:
898 kwargs['B'] = 'af+l%s' % label
900 kwargs['C'] = 'my_%s.cpt' % cmap
901 a_str = cbar_anchor[anchor]
903 w = self.width / 3.
904 h = w / 10.
906 lgap = rgap = w / 10.
907 bgap, tgap = h, h / 10.
909 dx, dy = 2.5 * lgap, 2. * tgap
911 if 'bottom' in anchor:
912 dy += 4 * h
914 self.gmt.psscale(
915 D='j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy),
916 F='+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap),
917 *self.jxyr,
918 **kwargs)
920 def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs):
921 '''
922 Draw vectors based on two grid files.
924 Two grid files for vector lengths in x and y need to be given. The
925 function calls gmt.grdvector. All arguments defined for this function
926 in gmt can be passed as keyword arguments. Different standard settings
927 are applied if not defined differently.
929 :param x_gridfile:
930 File of the grid defining vector lengths in x.
931 :type x_gridfile:
932 str
934 :param y_gridfile:
935 File of the grid defining vector lengths in y.
936 :type y_gridfile:
937 str
939 :param vcolor:
940 Vector face color as defined in "G" option.
941 :type vcolor:
942 str
943 '''
945 kwargs['S'] = kwargs.get('S', 'il1.')
946 kwargs['I'] = kwargs.get('I', 'x20')
947 kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black')
948 kwargs['Q'] = kwargs.get('Q', '4c+e+n5c+h1')
950 self.gmt.grdvector(
951 x_gridfile, y_gridfile,
952 G='%s' % 'lightgrey' if not vcolor else vcolor,
953 *self.jxyr,
954 **kwargs)
956 def draw_dynamic_data(self, data, **kwargs):
957 '''
958 Draw an image of any data gridded on the patches e.g dislocation.
960 :param data:
961 Patchwise grid data.
962 :type data:
963 :py:class:`~numpy.ndarray`
964 '''
966 plot_data = data
968 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
970 clim = kwargs.pop('clim', (plot_data.min(), plot_data.max()))
972 cpt = []
973 if not op.exists('my_%s.cpt' % kwargs['cmap']):
974 make_colormap(self.gmt, clim[0], clim[1],
975 cmap=kwargs['cmap'], space=False)
976 cpt = [kwargs['cmap']]
978 tmp_grd_file = 'tmpdata.grd'
979 self.patch_data_to_grid(plot_data, tmp_grd_file)
980 self.draw_image(tmp_grd_file, **kwargs)
982 clear_temp(gridfiles=[tmp_grd_file], cpts=cpt)
984 def draw_patch_parameter(self, attribute, **kwargs):
985 '''
986 Draw an image of a chosen patch attribute e.g traction.
988 :param attribute:
989 Patch attribute, which is plotted. All patch attributes can be
990 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`)
991 and also ``traction``, ``tx``, ``ty`` or ``tz`` to display the
992 length or the single components of the traction vector.
993 :type attribute:
994 str
995 '''
997 a = attribute
998 source = self.source
1000 if a == 'traction':
1001 data = num.linalg.norm(source.get_tractions(), axis=1)
1002 elif a == 'tx':
1003 data = source.get_tractions()[:, 0]
1004 elif a == 'ty':
1005 data = source.get_tractions()[:, 1]
1006 elif a == 'tz':
1007 data = source.get_tractions()[:, 2]
1008 else:
1009 data = source.get_patch_attribute(attribute)
1011 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1012 data *= factor
1014 kwargs['label'] = kwargs.get(
1015 'label',
1016 '%s [%s]' % (a, cbar_helper[a]['unit']))
1018 self.draw_dynamic_data(data, **kwargs)
1020 def draw_time_contour(self, store, clevel=[], **kwargs):
1021 '''
1022 Draw high contour lines of the rupture front propgation time.
1024 :param store:
1025 Greens function store, which is used for time calculation.
1026 :type store:
1027 :py:class:`~pyrocko.gf.store.Store`
1028 :param clevel:
1029 List of times, for which contour lines are drawn.
1030 :type clevel:
1031 optional, list of float
1032 '''
1034 _, _, _, _, points_xy = self.source._discretize_points(store, cs='xyz')
1035 _, _, times, _, _, _ = self.source.get_vr_time_interpolators(store)
1037 scaler = AutoScaler(mode='0-max', approx_ticks=8)
1038 scale = scaler.make_scale([num.min(times), num.max(times)])
1040 if clevel:
1041 if len(clevel) > 1:
1042 kwargs['anot_int'] = num.min(num.diff(clevel))
1043 else:
1044 kwargs['anot_int'] = clevel[0]
1046 kwargs['contour_int'] = kwargs['anot_int']
1047 kwargs['L'] = '0/%g' % num.max(clevel)
1049 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1050 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1051 kwargs['unit'] = kwargs.get('unit', cbar_helper['time']['unit'])
1052 kwargs['L'] = kwargs.get('L', '0/%g' % (num.max(times) + 1.))
1053 kwargs['G'] = kwargs.get('G', 'n2/3c')
1055 tmp_grd_file = 'tmpdata.grd'
1056 self.xy_data_to_grid(points_xy[:, 0], points_xy[:, 1],
1057 times, tmp_grd_file)
1058 self.draw_contour(tmp_grd_file, **kwargs)
1060 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
1062 def draw_points(self, lats, lons, symbol='point', size=None, **kwargs):
1063 '''
1064 Draw points at given locations.
1066 :param lats:
1067 Point latitude coordinates in [deg].
1068 :type lats:
1069 iterable
1071 :param lons:
1072 Point longitude coordinates in [deg].
1073 :type lons:
1074 iterable
1076 :param symbol:
1077 Define symbol of the points (``'star', 'circle', 'point',
1078 'triangle'``) - default is ``point``.
1079 :type symbol:
1080 optional, str
1082 :param size:
1083 Size of the points in [points].
1084 :type size:
1085 optional, float
1086 '''
1088 sym_to_gmt = dict(
1089 star='a',
1090 circle='c',
1091 point='p',
1092 triangle='t')
1094 lats = num.atleast_1d(lats)
1095 lons = num.atleast_1d(lons)
1097 if lats.shape[0] != lons.shape[0]:
1098 raise IndexError('lats and lons do not have the same shape!')
1100 if size is None:
1101 size = self._fontsize / 3.
1103 kwargs['S'] = kwargs.get('S', sym_to_gmt[symbol] + '%gp' % size)
1104 kwargs['G'] = kwargs.get('G', gmtpy.color('scarletred2'))
1105 kwargs['W'] = kwargs.get('W', '2p,%s' % self._fontcolor)
1107 self.gmt.psxy(
1108 in_columns=[lons, lats],
1109 *self.jxyr,
1110 **kwargs)
1112 def draw_nucleation_point(self, **kwargs):
1113 '''
1114 Plot the nucleation point onto the map.
1115 '''
1117 nlat, nlon = xy_to_latlon(
1118 self.source, self.source.nucleation_x, self.source.nucleation_y)
1120 self.draw_points(nlat, nlon, **kwargs)
1122 def draw_dislocation(self, time=None, component='', **kwargs):
1123 '''
1124 Draw dislocation onto map at any time.
1126 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given
1127 component the patchwise dislocation is plotted onto the map.
1129 :param time:
1130 Time after origin, for which dislocation is computed. If ``None``,
1131 ``tmax`` is taken.
1132 :type time:
1133 optional, float
1135 :param component:
1136 Dislocation component, which shall be plotted: ``x`` along strike,
1137 ``y`` along updip, ``z`` normal. If ``None``, the
1138 length of the dislocation vector is plotted.
1139 :type component: optional, str
1140 '''
1142 disl = self.source.get_slip(time=time)
1144 if component:
1145 data = disl[:, c2disl[component]]
1146 else:
1147 data = num.linalg.norm(disl, axis=1)
1149 kwargs['label'] = kwargs.get(
1150 'label', 'u%s [m]' % (component))
1152 self.draw_dynamic_data(data, **kwargs)
1154 def draw_dislocation_contour(
1155 self, time=None, component=None, clevel=[], **kwargs):
1156 '''
1157 Draw dislocation contour onto map at any time.
1159 For a given time (if ``time`` is ``None``, ``tmax`` is used) and given
1160 component the patchwise dislocation is plotted as contour onto the map.
1162 :param time:
1163 Time after origin, for which dislocation is computed. If ``None``,
1164 ``tmax`` is taken.
1165 :type time:
1166 optional, float
1168 :param component:
1169 Dislocation component, which shall be plotted: ``x`` along strike,
1170 ``y`` along updip, ``z`` normal``. If ``None``, the length of the
1171 dislocation vector is plotted.
1172 :type component: optional, str
1174 :param clevel:
1175 Times, for which contour lines are drawn.
1176 :type clevel:
1177 optional, list of float
1178 '''
1180 disl = self.source.get_slip(time=time)
1182 if component:
1183 data = disl[:, c2disl[component]]
1184 else:
1185 data = num.linalg.norm(disl, axis=1)
1187 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1188 scale = scaler.make_scale([num.min(data), num.max(data)])
1190 if clevel:
1191 if len(clevel) > 1:
1192 kwargs['anot_int'] = num.min(num.diff(clevel))
1193 else:
1194 kwargs['anot_int'] = clevel[0]
1196 kwargs['contour_int'] = kwargs['anot_int']
1197 kwargs['L'] = '%g/%g' % (
1198 num.min(clevel) - kwargs['contour_int'],
1199 num.max(clevel) + kwargs['contour_int'])
1200 else:
1201 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1202 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1203 kwargs['L'] = kwargs.get('L', '%g/%g' % (
1204 num.min(data) - 1., num.max(data) + 1.))
1206 kwargs['unit'] = kwargs.get('unit', ' m')
1207 kwargs['G'] = kwargs.get('G', 'n2/3c')
1209 tmp_grd_file = 'tmpdata.grd'
1210 self.patch_data_to_grid(data, tmp_grd_file)
1211 self.draw_contour(tmp_grd_file, **kwargs)
1213 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
1215 def draw_dislocation_vector(self, time=None, **kwargs):
1216 '''
1217 Draw vector arrows onto map indicating direction of dislocation.
1219 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1220 and given component the dislocation is plotted as vectors onto the map.
1222 :param time:
1223 Time after origin [s], for which dislocation is computed. If
1224 ``None``, ``tmax`` is used.
1225 :type time:
1226 optional, float
1227 '''
1229 disl = self.source.get_slip(time=time)
1231 p_strike = self.source.get_patch_attribute('strike') * d2r
1232 p_dip = self.source.get_patch_attribute('dip') * d2r
1234 disl_dh = num.cos(p_dip) * disl[:, 1]
1235 disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh
1236 disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh
1238 tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')]
1240 self.patch_data_to_grid(disl_n, tmp_grd_files[0])
1241 self.patch_data_to_grid(disl_e, tmp_grd_files[1])
1243 self.draw_vector(
1244 tmp_grd_files[1], tmp_grd_files[0],
1245 **kwargs)
1247 clear_temp(gridfiles=tmp_grd_files, cpts=[])
1249 def draw_top_edge(self, **kwargs):
1250 '''
1251 Indicate rupture top edge on map.
1252 '''
1254 outline = self.source.outline(cs='latlondepth')
1255 top_edge = outline[:2, :]
1257 kwargs = kwargs or {}
1258 kwargs['W'] = kwargs.get(
1259 'W', '%gp,%s' % (self._fontsize / 10., gmtpy.color('orange3')))
1261 self.gmt.psxy(
1262 in_columns=[top_edge[:, 1], top_edge[:, 0]],
1263 *self.jxyr,
1264 **kwargs)
1267class RuptureView(Object):
1268 '''
1269 Plot of attributes and results of the
1270 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1271 '''
1273 _patches_to_lw = RuptureMap._patches_to_lw
1275 def __init__(self, source=None, figsize=None, fontsize=None):
1276 self._source = source
1277 self._axes = None
1279 self.figsize = figsize or mpl_papersize('halfletter', 'landscape')
1280 self.fontsize = fontsize or 10
1281 mpl_init(fontsize=self.fontsize)
1283 self._fig = None
1284 self._axes = None
1285 self._is_1d = False
1287 @property
1288 def source(self):
1289 '''
1290 PseudoDynamicRupture whose attributes are plotted.
1292 Note, that source.patches attribute needs to be calculated for
1293 :type source: :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1294 '''
1296 if self._source is None:
1297 raise SourceError('No source given. Please define it!')
1299 if not isinstance(self._source, PseudoDynamicRupture):
1300 raise SourceError('This function is only capable for a source of'
1301 ' type: %s' % type(PseudoDynamicRupture()))
1303 if not self._source.patches:
1304 raise TypeError('No source patches are defined. Please run'
1305 '\"discretize_patches\" on your source')
1307 return self._source
1309 @source.setter
1310 def source(self, source):
1311 self._source = source
1313 def _setup(
1314 self,
1315 title='',
1316 xlabel='',
1317 ylabel='',
1318 aspect=1.,
1319 spatial_plot=True,
1320 **kwargs):
1322 if self._fig is not None and self._axes is not None:
1323 return
1325 self._fig = plt.figure(figsize=self.figsize)
1327 self._axes = self._fig.add_subplot(1, 1, 1, aspect=aspect)
1328 ax = self._axes
1330 if ax is not None:
1331 ax.set_title(title)
1332 ax.grid(alpha=.3)
1333 ax.set_xlabel(xlabel)
1334 ax.set_ylabel(ylabel)
1336 if spatial_plot:
1337 ax.xaxis.set_major_formatter(FuncFormatter(km_formatter))
1338 ax.yaxis.set_major_formatter(FuncFormatter(km_formatter))
1340 def _clear_all(self):
1341 plt.cla()
1342 plt.clf()
1343 plt.close()
1345 self._fig, self._axes = None, None
1347 def _draw_scatter(self, x, y, *args, **kwargs):
1348 default_kwargs = dict(
1349 linewidth=0,
1350 marker='o',
1351 markerfacecolor=mpl_color('skyblue2'),
1352 markersize=6.,
1353 markeredgecolor=mpl_color('skyblue3'))
1354 default_kwargs.update(kwargs)
1356 if self._axes is not None:
1357 self._axes.plot(x, y, *args, **default_kwargs)
1359 def _draw_image(self, length, width, data, *args, **kwargs):
1360 if self._axes is not None:
1361 if 'extent' not in kwargs:
1362 kwargs['extent'] = [
1363 num.min(length), num.max(length),
1364 num.max(width), num.min(width)]
1366 im = self._axes.imshow(
1367 data,
1368 *args,
1369 interpolation='none',
1370 vmin=kwargs.get('clim', [None])[0],
1371 vmax=kwargs.get('clim', [None, None])[1],
1372 **kwargs)
1374 del kwargs['extent']
1375 if 'aspect' in kwargs:
1376 del kwargs['aspect']
1377 if 'clim' in kwargs:
1378 del kwargs['clim']
1379 if 'cmap' in kwargs:
1380 del kwargs['cmap']
1382 plt.colorbar(
1383 im, *args, shrink=0.9, pad=0.03, aspect=15., **kwargs)
1385 def _draw_contour(self, x, y, data, clevel=None, unit='', *args, **kwargs):
1386 setup_kwargs = dict(
1387 xlabel=kwargs.pop('xlabel', 'along strike [km]'),
1388 ylabel=kwargs.pop('ylabel', 'down dip [km]'),
1389 title=kwargs.pop('title', ''),
1390 aspect=kwargs.pop('aspect', 1))
1392 self._setup(**setup_kwargs)
1394 if self._axes is not None:
1395 if clevel is None:
1396 scaler = AutoScaler(mode='min-max')
1397 scale = scaler.make_scale([data.min(), data.max()])
1399 clevel = num.arange(scale[0], scale[1] + scale[2], scale[2])
1401 if not isinstance(clevel, num.ndarray):
1402 clevel = num.array([clevel])
1404 clevel = clevel[clevel < data.max()]
1406 cont = self._axes.contour(
1407 x, y, data, clevel, *args,
1408 linewidths=1.5, **kwargs)
1410 plt.setp(cont.collections, path_effects=[
1411 patheffects.withStroke(linewidth=2.0, foreground="beige"),
1412 patheffects.Normal()])
1414 clabels = self._axes.clabel(
1415 cont, clevel[::2], *args,
1416 inline=1, fmt='%g' + '%s' % unit,
1417 inline_spacing=15, rightside_up=True, use_clabeltext=True,
1418 **kwargs)
1420 plt.setp(clabels, path_effects=[
1421 patheffects.withStroke(linewidth=1.25, foreground="beige"),
1422 patheffects.Normal()])
1424 def draw_points(self, length, width, *args, **kwargs):
1425 '''
1426 Draw a point onto the figure.
1428 Args and kwargs can be defined according to
1429 :py:func:`matplotlib.pyplot.scatter`.
1431 :param length:
1432 Point(s) coordinate on the rupture plane along strike relative to
1433 the anchor point in [m].
1434 :type length:
1435 float, :py:class:`~numpy.ndarray`
1437 :param width:
1438 Point(s) coordinate on the rupture plane along downdip relative to
1439 the anchor point in [m].
1440 :type width:
1441 float, :py:class:`~numpy.ndarray`
1442 '''
1444 if self._axes is not None:
1445 kwargs['color'] = kwargs.get('color', mpl_color('scarletred2'))
1446 kwargs['s'] = kwargs.get('s', 100.)
1448 self._axes.scatter(length, width, *args, **kwargs)
1450 def draw_dynamic_data(self, data, **kwargs):
1451 '''
1452 Draw an image of any data gridded on the patches e.g dislocation.
1454 :param data:
1455 Patchwise grid data.
1456 :type data:
1457 :py:class:`~numpy.ndarray`
1458 '''
1460 plot_data = data
1462 anchor_x, anchor_y = map_anchor[self.source.anchor]
1464 length, width = xy_to_lw(
1465 self.source, num.array([-1., 1.]), num.array([-1., 1.]))
1467 data = plot_data.reshape(self.source.nx, self.source.ny).T
1469 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
1471 setup_kwargs = dict(
1472 xlabel='along strike [km]',
1473 ylabel='down dip [km]',
1474 title='',
1475 aspect=1)
1476 setup_kwargs.update(kwargs)
1478 kwargs = {k: v for k, v in kwargs.items() if k not in
1479 ('xlabel', 'ylabel', 'title')}
1480 self._setup(**setup_kwargs)
1481 self._draw_image(length=length, width=width, data=data, **kwargs)
1483 def draw_patch_parameter(self, attribute, **kwargs):
1484 '''
1485 Draw an image of a chosen patch attribute e.g traction.
1487 :param attribute:
1488 Patch attribute, which is plotted. All patch attributes can be
1489 taken (see doc of :py:class:`~pyrocko.modelling.okada.OkadaSource`)
1490 and also ``'traction', 'tx', 'ty', 'tz'`` to display the length
1491 or the single components of the traction vector.
1492 :type attribute:
1493 str
1494 '''
1496 a = attribute
1497 source = self.source
1499 if a == 'traction':
1500 data = num.linalg.norm(source.get_tractions(), axis=1)
1501 elif a == 'tx':
1502 data = source.get_tractions()[:, 0]
1503 elif a == 'ty':
1504 data = source.get_tractions()[:, 1]
1505 elif a == 'tz':
1506 data = source.get_tractions()[:, 2]
1507 else:
1508 data = source.get_patch_attribute(attribute)
1510 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1511 data *= factor
1513 kwargs['label'] = kwargs.get(
1514 'label',
1515 '%s [%s]' % (a, cbar_helper[a]['unit']))
1517 return self.draw_dynamic_data(data=data, **kwargs)
1519 def draw_time_contour(self, store, clevel=[], **kwargs):
1520 '''
1521 Draw high resolution contours of the rupture front propgation time
1523 :param store:
1524 Greens function store, which is used for time calculation.
1525 :type store:
1526 :py:class:`~pyrocko.gf.store.Store`
1528 :param clevel:
1529 Levels of the contour lines. If no levels are given, they are
1530 automatically computed based on ``tmin`` and ``tmax``.
1531 :type clevel:
1532 optional, list
1533 '''
1535 source = self.source
1536 default_kwargs = dict(
1537 colors='#474747ff'
1538 )
1539 default_kwargs.update(kwargs)
1541 _, _, _, _, points_xy = source._discretize_points(store, cs='xyz')
1542 _, _, _, _, time_interpolator, _ = source.get_vr_time_interpolators(
1543 store)
1545 times = time_interpolator.values
1547 scaler = AutoScaler(mode='min-max', approx_ticks=8)
1548 scale = scaler.make_scale([times.min(), times.max()])
1550 if len(clevel) == 0:
1551 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1553 points_x = time_interpolator.grid[0]
1554 points_y = time_interpolator.grid[1]
1556 self._draw_contour(points_x, points_y, data=times.T,
1557 clevel=clevel, unit='s', **default_kwargs)
1559 def draw_nucleation_point(self, **kwargs):
1560 '''
1561 Draw the nucleation point onto the map.
1562 '''
1564 nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y
1565 length, width = xy_to_lw(self.source, nuc_x, nuc_y)
1566 self.draw_points(length, width, marker='o', **kwargs)
1568 def draw_dislocation(self, time=None, component='', **kwargs):
1569 '''
1570 Draw dislocation onto map at any time.
1572 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1573 and given component the patchwise dislocation is plotted onto the map.
1575 :param time:
1576 Time after origin [s], for which dislocation is computed.
1577 If ``None``, ``tmax`` is taken.
1578 :type time:
1579 optional, float
1581 :param component:
1582 Dislocation component, which shall be plotted: ``x``
1583 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
1584 length of the dislocation vector is plotted
1585 :type component: optional, str
1586 '''
1588 disl = self.source.get_slip(time=time)
1590 if component:
1591 data = disl[:, c2disl[component]]
1592 else:
1593 data = num.linalg.norm(disl, axis=1)
1595 kwargs['label'] = kwargs.get(
1596 'label', 'u%s [m]' % (component))
1598 self.draw_dynamic_data(data, **kwargs)
1600 def draw_dislocation_contour(
1601 self, time=None, component=None, clevel=[], **kwargs):
1602 '''
1603 Draw dislocation contour onto map at any time.
1605 For a given time (if time is ``None``, ``tmax`` is used) and given
1606 component the patchwise dislocation is plotted as contour onto the map.
1608 :param time:
1609 Time after origin, for which dislocation is computed. If
1610 ``None``, ``tmax`` is taken.
1611 :type time:
1612 optional, float
1614 :param component:
1615 Dislocation component, which shall be plotted. ``x``
1616 along strike, ``y`` along updip, ``z`` - normal. If ``None``
1617 is given, the length of the dislocation vector is plotted.
1618 :type component:
1619 optional, str
1620 '''
1622 disl = self.source.get_slip(time=time)
1624 if component:
1625 data = disl[:, c2disl[component]]
1626 else:
1627 data = num.linalg.norm(disl, axis=1)
1629 data = data.reshape(self.source.ny, self.source.nx, order='F')
1631 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1632 scale = scaler.make_scale([data.min(), data.max()])
1634 if len(clevel) == 0:
1635 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1637 anchor_x, anchor_y = map_anchor[self.source.anchor]
1639 length, width = self._patches_to_lw()
1640 length, width = length[1:-1], width[1:-1]
1642 kwargs['colors'] = kwargs.get('colors', '#474747ff')
1644 self._setup(**kwargs)
1645 self._draw_contour(
1646 length, width, data=data, clevel=clevel, unit='m', **kwargs)
1648 def draw_source_dynamics(
1649 self, variable, store, deltat=None, *args, **kwargs):
1650 '''
1651 Display dynamic source parameter.
1653 Fast inspection possibility for the cumulative moment and the source
1654 time function approximation (assuming equal paths between different
1655 patches and observation point - valid for an observation point in the
1656 far field perpendicular to the source strike), so the cumulative moment
1657 rate function.
1659 :param variable:#
1660 Dynamic parameter, which shall be plotted. Choose
1661 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1662 :type variable:
1663 str
1665 :param store:
1666 Greens function store, whose store.config.deltat defines
1667 the time increment between two parameter snapshots. If ``store`` is
1668 not given, the time increment is defined is taken from ``deltat``.
1669 :type store:
1670 :py:class:`~pyrocko.gf.store.Store`
1672 :param deltat:
1673 Time increment between two parameter snapshots. If not
1674 given, store.config.deltat is used to define ``deltat``.
1675 :type deltat:
1676 optional, float
1677 '''
1679 v = variable
1681 data, times = self.source.get_moment_rate(store=store, deltat=deltat)
1682 deltat = num.concatenate([(num.diff(times)[0], ), num.diff(times)])
1684 if v in ('moment_rate', 'stf'):
1685 name, unit = 'dM/dt', 'Nm/s'
1686 elif v in ('cumulative_moment', 'moment'):
1687 data = num.cumsum(data * deltat)
1688 name, unit = 'M', 'Nm'
1689 else:
1690 raise ValueError('No dynamic data for given variable %s found' % v)
1692 self._setup(xlabel='time [s]',
1693 ylabel='%s / %.2g %s' % (name, data.max(), unit),
1694 aspect='auto',
1695 spatial_plot=False)
1696 self._draw_scatter(times, data/num.max(data), *args, **kwargs)
1697 self._is_1d = True
1699 def draw_patch_dynamics(
1700 self, variable, nx, ny, store=None, deltat=None, *args, **kwargs):
1701 '''
1702 Display dynamic boundary element / patch parameter.
1704 Fast inspection possibility for different dynamic parameter for a
1705 single patch / boundary element. The chosen parameter is plotted for
1706 the chosen patch.
1708 :param variable:
1709 Dynamic parameter, which shall be plotted. Choose
1710 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment').
1711 :type variable:
1712 str
1714 :param nx:
1715 Patch index along strike (range: 0:source.nx - 1).
1716 :type nx:
1717 int
1719 :param nx:
1720 Patch index downdip (range: 0:source.ny - 1).
1721 :type nx:
1722 int
1724 :param store:
1725 Greens function store, whose store.config.deltat defines
1726 the time increment between two parameter snapshots. If ``store`` is
1727 not given, the time increment is defined is taken from ``deltat``.
1728 :type store:
1729 optional, :py:class:`~pyrocko.gf.store.Store`
1731 :param deltat:
1732 Time increment between two parameter snapshots. If not given,
1733 store.config.deltat is used to define ``deltat``.
1734 :type deltat:
1735 optional, float
1736 '''
1738 v = variable
1739 source = self.source
1740 idx = nx * source.ny + nx
1742 m = re.match(r'dislocation_([xyz])', v)
1744 if v in ('moment_rate', 'cumulative_moment', 'moment', 'stf'):
1745 data, times = source.get_moment_rate_patches(
1746 store=store, deltat=deltat)
1747 elif 'dislocation' in v or 'slip_rate' == v:
1748 ddisloc, times = source.get_delta_slip(store=store, deltat=deltat)
1750 if v in ('moment_rate', 'stf'):
1751 data, times = source.get_moment_rate_patches(
1752 store=store, deltat=deltat)
1753 name, unit = 'dM/dt', 'Nm/s'
1754 elif v in ('cumulative_moment', 'moment'):
1755 data, times = source.get_moment_rate_patches(
1756 store=store, deltat=deltat)
1757 data = num.cumsum(data, axis=1)
1758 name, unit = 'M', 'Nm'
1759 elif v == 'slip_rate':
1760 data, times = source.get_slip_rate(store=store, deltat=deltat)
1761 name, unit = 'du/dt', 'm/s'
1762 elif v == 'dislocation':
1763 data, times = source.get_delta_slip(store=store, deltat=deltat)
1764 data = num.linalg.norm(num.cumsum(data, axis=2), axis=1)
1765 name, unit = 'du', 'm'
1766 elif m:
1767 data, times = source.get_delta_slip(store=store, deltat=deltat)
1768 data = num.cumsum(data, axis=2)[:, c2disl[m.group(1)], :]
1769 name, unit = 'du%s' % m.group(1), 'm'
1770 else:
1771 raise ValueError('No dynamic data for given variable %s found' % v)
1773 self._setup(xlabel='time [s]',
1774 ylabel='%s / %.2g %s' % (name, num.max(data), unit),
1775 aspect='auto',
1776 spatial_plot=False)
1777 self._draw_scatter(times, data[idx, :]/num.max(data),
1778 *args, **kwargs)
1779 self._is_1d = True
1781 def finalize(self):
1782 if self._is_1d:
1783 return
1785 length, width = xy_to_lw(
1786 self.source, num.array([-1., 1.]), num.array([1., -1.]))
1788 self._axes.set_xlim(length)
1789 self._axes.set_ylim(width)
1791 def gcf(self):
1792 self.finalize()
1793 return self._fig
1795 def save(self, filename, dpi=None):
1796 '''
1797 Save plot to file.
1799 :param filename:
1800 Filename and path, where the plot is stored.
1801 :type filename:
1802 str
1804 :param dpi:
1805 Resolution of the output plot in [dpi].
1806 :type dpi:
1807 int
1808 '''
1810 self.finalize()
1811 try:
1812 self._fig.savefig(fname=filename, dpi=dpi, bbox_inches='tight')
1813 except TypeError:
1814 self._fig.savefig(filename=filename, dpi=dpi, bbox_inches='tight')
1816 self._clear_all()
1818 def show_plot(self):
1819 '''
1820 Show plot.
1821 '''
1823 self.finalize()
1824 plt.show()
1825 self._clear_all()
1828def render_movie(fn_path, output_path, framerate=20):
1829 '''
1830 Generate a mp4 movie based on given png files using
1831 `ffmpeg <https://ffmpeg.org>`_.
1833 Render a movie based on a set of given .png files in fn_path. All files
1834 must have a filename specified by ``fn_path`` (e.g. giving ``fn_path`` with
1835 ``/temp/f%04.png`` a valid png filename would be ``/temp/f0001.png``). The
1836 files must have a numbering, indicating their order in the movie.
1838 :param fn_path:
1839 Path and fileformat specification of the input .png files.
1840 :type fn_path:
1841 str
1843 :param output_path:
1844 Path and filename of the output ``.mp4`` movie file.
1845 :type output_path:
1846 str
1848 :param deltat:
1849 Time between individual frames (``1 / framerate``) in [s].
1850 :type deltat:
1851 optional, float
1852 '''
1854 try:
1855 check_call(['ffmpeg', '-loglevel', 'panic'])
1856 except CalledProcessError:
1857 pass
1858 except (TypeError, IOError):
1859 logger.warning(
1860 'Package ffmpeg needed for movie rendering. Please install it '
1861 '(e.g. on linux distr. via sudo apt-get ffmpeg) and retry.')
1862 return False
1864 try:
1865 check_call([
1866 'ffmpeg', '-loglevel', 'info', '-y',
1867 '-framerate', '%g' % framerate,
1868 '-i', fn_path,
1869 '-vcodec', 'libx264',
1870 '-preset', 'medium',
1871 '-tune', 'animation',
1872 '-pix_fmt', 'yuv420p',
1873 '-movflags', '+faststart',
1874 '-filter:v', "crop='floor(in_w/2)*2:floor(in_h/2)*2'",
1875 '-crf', '15',
1876 output_path])
1878 return True
1879 except CalledProcessError as e:
1880 logger.warning(e)
1881 return False
1884def render_gif(fn, output_path, loops=-1):
1885 '''
1886 Generate a gif based on a given mp4 using ffmpeg.
1888 Render a gif based on a given .mp4 movie file in ``fn`` path.
1890 :param fn:
1891 Path and file name of the input .mp4 file.
1892 :type fn:
1893 str
1895 :param output_path:
1896 Path and filename of the output animated gif file.
1897 :type output_path:
1898 str
1899 :param loops:
1900 Number of gif repeats (loops). ``-1`` is not repetition, ``0``
1901 infinite.
1902 :type loops:
1903 optional, integer
1904 '''
1906 try:
1907 check_call(['ffmpeg', '-loglevel', 'panic'])
1908 except CalledProcessError:
1909 pass
1910 except (TypeError, IOError):
1911 logger.warning(
1912 'Package ffmpeg needed for movie rendering. Please install it '
1913 '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.')
1914 return False
1916 try:
1917 check_call([
1918 'ffmpeg', '-hide_banner', '-loglevel', 'panic', '-i',
1919 fn,
1920 '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse',
1921 '-loop', '%d' % loops,
1922 output_path])
1924 return True
1925 except CalledProcessError as e:
1926 logger.warning(e)
1927 return False
1930def rupture_movie(
1931 source,
1932 store,
1933 variable='dislocation',
1934 draw_time_contours=False,
1935 fn_path='.',
1936 prefix='',
1937 plot_type='map',
1938 deltat=None,
1939 framerate=None,
1940 store_images=False,
1941 render_as_gif=False,
1942 gif_loops=-1,
1943 **kwargs):
1944 '''
1945 Generate a movie based on a given source for dynamic parameter.
1947 Create a MPEG-4 movie or gif of one of the following dynamic parameters
1948 (``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1949 (along updip), ``dislocation_z`` (normal), ``slip_rate``, ``moment_rate``).
1950 If desired, the single snap shots can be stored as images as well.
1951 ``kwargs`` have to be given according to the chosen ``plot_type``.
1953 :param source:
1954 Pseudo dynamic rupture, for which the movie is produced.
1955 :type source:
1956 :py:class:`~pyrocko.gf.seismosizer.PseudoDynamicRupture`
1958 :param store:
1959 Greens function store, which is used for time calculation. If
1960 ``deltat`` is not given, it is taken from the store.config.deltat
1961 :type store:
1962 :py:class:`~pyrocko.gf.store.Store`
1964 :param variable:
1965 Dynamic parameter, which shall be plotted. Choose between
1966 ``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1967 (along updip), ``dislocation_z`` (normal), ``slip_rate`` and
1968 ``moment_rate``, default ``dislocation``.
1969 :type variable:
1970 optional, str
1972 :param draw_time_contours:
1973 If ``True``, corresponding isochrones are drawn on the each plots.
1974 :type draw_time_contours:
1975 optional, bool
1977 :param fn_path:
1978 Absolut or relative path, where movie (and optional images) are stored.
1979 :type fn_path:
1980 optional, str
1982 :param prefix:
1983 File prefix used for the movie (and optional image) files.
1984 :type prefix:
1985 optional, str
1987 :param plot_type:
1988 Choice of plot type: ``map``, ``view`` (map plot using
1989 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap`
1990 or plane view using
1991 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`).
1992 :type plot_type:
1993 optional, str
1995 :param deltat:
1996 Time between parameter snapshots. If not given, store.config.deltat is
1997 used to define ``deltat``.
1998 :type deltat:
1999 optional, float
2001 :param store_images:
2002 Choice to store the single .png parameter snapshots in ``fn_path`` or
2003 not.
2004 :type store_images:
2005 optional, bool
2007 :param render_as_gif:
2008 If ``True``, the movie is converted into a gif. If ``False``, the movie
2009 is returned as mp4.
2010 :type render_as_gif:
2011 optional, bool
2013 :param gif_loops:
2014 If ``render_as_gif`` is ``True``, a gif with ``gif_loops`` number of
2015 loops (repetitions) is returned. ``-1`` is no repetition, ``0``
2016 infinite.
2017 :type gif_loops:
2018 optional, integer
2019 '''
2021 v = variable
2022 assert plot_type in ('map', 'view')
2024 if not source.patches:
2025 source.discretize_patches(store, interpolation='multilinear')
2027 if source.coef_mat is None:
2028 source.calc_coef_mat()
2030 prefix = prefix or v
2031 deltat = deltat or store.config.deltat
2032 framerate = max(framerate or int(1./deltat), 1)
2034 if v == 'moment_rate':
2035 data, times = source.get_moment_rate_patches(deltat=deltat)
2036 name, unit = 'dM/dt', 'Nm/s'
2037 elif 'dislocation' in v or 'slip_rate' == v:
2038 ddisloc, times = source.get_delta_slip(deltat=deltat)
2039 else:
2040 raise ValueError('No dynamic data for given variable %s found' % v)
2042 deltat = times[1] - times[0]
2044 m = re.match(r'dislocation_([xyz])', v)
2045 if m:
2046 data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]]
2047 name, unit = 'du%s' % m.group(1), 'm'
2048 elif v == 'dislocation':
2049 data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2)
2050 name, unit = 'du', 'm'
2051 elif v == 'slip_rate':
2052 data = num.linalg.norm(ddisloc, axis=2) / deltat
2053 name, unit = 'du/dt', 'm/s'
2055 if plot_type == 'map':
2056 plt_base = RuptureMap
2057 elif plot_type == 'view':
2058 plt_base = RuptureView
2059 else:
2060 raise AttributeError('invalid type: %s' % plot_type)
2062 attrs_base = get_kwargs(plt_base)
2063 kwargs_base = dict([k for k in kwargs.items() if k[0] in attrs_base])
2064 kwargs_plt = dict([k for k in kwargs.items() if k[0] not in kwargs_base])
2066 if 'clim' in kwargs_plt:
2067 data = num.clip(data, kwargs_plt['clim'][0], kwargs_plt['clim'][1])
2068 else:
2069 kwargs_plt['clim'] = [num.min(data), num.max(data)]
2071 if 'label' not in kwargs_plt:
2072 vmax = num.max(num.abs(kwargs_plt['clim']))
2073 data /= vmax
2075 kwargs_plt['label'] = '%s / %.2g %s' % (name, vmax, unit)
2076 kwargs_plt['clim'] = [i / vmax for i in kwargs_plt['clim']]
2078 with tempfile.TemporaryDirectory(suffix='pyrocko') as temp_path:
2079 fns_temp = [op.join(temp_path, 'f%09d.png' % (it + 1))
2080 for it, _ in enumerate(times)]
2081 fn_temp_path = op.join(temp_path, 'f%09d.png')
2083 for it, (t, ft) in enumerate(zip(times, fns_temp)):
2084 plot_data = data[:, it]
2086 plt = plt_base(source=source, **kwargs_base)
2087 plt.draw_dynamic_data(plot_data, **kwargs_plt)
2088 plt.draw_nucleation_point()
2090 if draw_time_contours:
2091 plt.draw_time_contour(store, clevel=[t])
2093 plt.save(ft)
2095 fn_mp4 = op.join(temp_path, 'movie.mp4')
2096 return_code = render_movie(
2097 fn_temp_path,
2098 output_path=fn_mp4,
2099 framerate=framerate)
2101 if render_as_gif and return_code:
2102 render_gif(fn=fn_mp4, output_path=op.join(
2103 fn_path, '%s_%s_gif.gif' % (prefix, plot_type)),
2104 loops=gif_loops)
2106 elif return_code:
2107 shutil.move(fn_mp4, op.join(
2108 fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
2110 else:
2111 logger.error('ffmpeg failed. Exit')
2113 if store_images:
2114 fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t))
2115 for t in times]
2117 for ft, f in zip(fns_temp, fns):
2118 shutil.move(ft, f)
2121__all__ = [
2122 'make_colormap',
2123 'clear_temp',
2124 'xy_to_latlon',
2125 'xy_to_lw',
2126 'SourceError',
2127 'RuptureMap',
2128 'RuptureView',
2129 'rupture_movie',
2130 'render_movie']