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: Grid latitude coordinates [degree]
68 :type lats: iterable
69 :param lons: Grid longitude coordinates [degree]
70 :type lons: iterable
71 :param data: Grid data of any kind
72 :type data: :py:class:`numpy.ndarray`, ``(n_lons, n_lats)``
73 :param filename: Filename of the written grid file
74 :type filename: str
75 '''
77 gmtpy.savegrd(lons, lats, data, filename=filename, naming='lonlat')
80def _mplcmap_to_gmtcpt_code(mplcmap, steps=256):
81 '''
82 Get gmt readable R/G/A code from a given matplotlib colormap
84 :param mplcmap: Name of the demanded matplotlib colormap
85 :type mplcmap: str
87 :returns: Series of comma seperate R/G/B values for gmtpy usage
88 :rtype: str
89 '''
91 cmap = cm.get_cmap(mplcmap)
93 rgbas = [cmap(i) for i in num.linspace(0, 255, steps).astype(num.int64)]
95 return ','.join(['%d/%d/%d' % (
96 c[0] * 255, c[1] * 255, c[2] * 255) for c in rgbas])
99def make_colormap(
100 gmt,
101 vmin,
102 vmax,
103 C=None,
104 cmap=None,
105 space=False):
106 '''
107 Create gmt-readable colormap cpt file called my_<cmap>.cpt
109 :type vmin: Minimum value covered by the colormap
110 :param vmin: float
111 :type vmax: Maximum value covered by the colormap
112 :param vmax: float
113 :type C: comma seperated R/G/B values for cmap definition.
114 :param C: optional, str
115 :type cmap: Name of the colormap. Colormap is stored as "my_<cmap>.cpt".
116 If name is equivalent to a matplotlib colormap, R/G/B strings are
117 extracted from this colormap.
118 :param cmap: optional, str
119 :type space: If True, the range of the colormap is broadened below vmin and
120 above vmax.
121 :param space: optional, bool
122 '''
124 scaler = AutoScaler(mode='min-max')
125 scale = scaler.make_scale((vmin, vmax))
127 incr = scale[2]
128 margin = 0.
130 if vmin == vmax:
131 space = True
133 if space:
134 margin = incr
136 msg = ('Please give either a valid color code or a'
137 ' valid matplotlib colormap name.')
139 if C is None and cmap is None:
140 raise ValueError(msg)
142 if C is None:
143 try:
144 C = _mplcmap_to_gmtcpt_code(cmap)
145 except ValueError:
146 raise ValueError(msg)
148 if cmap is None:
149 logger.warning(
150 'No colormap name given. Uses temporary filename instead')
151 cmap = 'temp_cmap'
153 return gmt.makecpt(
154 C=C,
155 D='o',
156 T='%g/%g/%g' % (
157 vmin - margin, vmax + margin, incr),
158 Z=True,
159 out_filename='my_%s.cpt' % cmap,
160 suppress_defaults=True)
163def clear_temp(gridfiles=[], cpts=[]):
164 '''
165 Clear all temporary needed grid and colormap cpt files
167 :param gridfiles: List of all "...grd" files, which shall be deleted
168 :type gridfiles: optional, list
169 :param cpts: List of all cmaps, whose corresponding "my_<cmap>.cpt" file
170 shall be deleted
171 :type cpts: optional, list
172 '''
174 for fil in gridfiles:
175 try:
176 os.remove(fil)
177 except OSError:
178 continue
179 for fil in cpts:
180 try:
181 os.remove('my_%s.cpt' % fil)
182 except OSError:
183 continue
186def xy_to_latlon(source, x, y):
187 '''
188 Convert x and y relative coordinates on extended ruptures into latlon
190 :param source: Extended source class, on which the given point is located
191 :type source: :py:class:`pyrocko.gf.seismosizer.RectangularSource` or
192 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
193 :param x: Relative point coordinate along strike (range: -1:1)
194 :type x: float or :py:class:`numpy.ndarray`
195 :param y: Relative downdip point coordinate (range: -1:1)
196 :type y: float or :py:class:`numpy.ndarray`
198 :returns: Latitude and Longitude [degrees] of the given point
199 :rtype: tuple of float
200 '''
202 s = source
203 ax, ay = map_anchor[s.anchor]
205 length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
206 strike, dip = s.strike * d2r, s.dip * d2r
208 northing = num.cos(strike)*length - num.cos(dip) * num.sin(strike)*width
209 easting = num.sin(strike)*length + num.cos(dip) * num.cos(strike)*width
211 northing += source.north_shift
212 easting += source.east_shift
214 return pod.ne_to_latlon(s.lat, s.lon, northing, easting)
217def xy_to_lw(source, x, y):
218 '''
219 Convert x and y relative coordinates on extended ruptures into length width
221 :param source: Extended source class, on which the given points are located
222 :type source: :py:class:`pyrocko.gf.seismosizer.RectangularSource` or
223 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
224 :param x: Relative point coordinates along strike (range: -1:1)
225 :type x: float or :py:class:`numpy.ndarray`
226 :param y: Relative downdip point coordinates (range: -1:1)
227 :type y: float or :py:class:`numpy.ndarray`
229 :returns: length and downdip width [m] of the given points from the anchor
230 :rtype: tuple of float
231 '''
233 length, width = source.length, source.width
235 ax, ay = map_anchor[source.anchor]
237 lengths = (x - ax) / 2. * length
238 widths = (y - ay) / 2. * width
240 return lengths, widths
243cbar_anchor = {
244 'center': 'MC',
245 'center_left': 'ML',
246 'center_right': 'MR',
247 'top': 'TC',
248 'top_left': 'TL',
249 'top_right': 'TR',
250 'bottom': 'BC',
251 'bottom_left': 'BL',
252 'bottom_right': 'BR'}
255cbar_helper = {
256 'traction': {
257 'unit': 'MPa',
258 'factor': 1e-6},
259 'tx': {
260 'unit': 'MPa',
261 'factor': 1e-6},
262 'ty': {
263 'unit': 'MPa',
264 'factor': 1e-6},
265 'tz': {
266 'unit': 'MPa',
267 'factor': 1e-6},
268 'time': {
269 'unit': 's',
270 'factor': 1.},
271 'strike': {
272 'unit': '°',
273 'factor': 1.},
274 'dip': {
275 'unit': '°',
276 'factor': 1.},
277 'vr': {
278 'unit': 'km/s',
279 'factor': 1e-3},
280 'length': {
281 'unit': 'km',
282 'factor': 1e-3},
283 'width': {
284 'unit': 'km',
285 'factor': 1e-3}
286}
289fonttype = 'Helvetica'
291c2disl = dict([('x', 0), ('y', 1), ('z', 2)])
294def _make_gmt_conf(fontcolor, size):
295 '''
296 Update different gmt parameters depending on figure size and fontcolor
298 :param fontcolor: GMT readable colorcode / colorstring for the font
299 :type fontcolor: str
300 :param size: Tuple of the figure size (width, height) [centimetre]
301 :type size: tuple of float
303 :returns: estimate best fitting fontsize,
304 Dictionary of different gmt configuration parameter
305 :rtype: float, dict
306 '''
308 color = fontcolor
309 fontsize = num.max(size)
311 font = '%gp,%s' % (fontsize, fonttype)
313 pen_size = fontsize / 16.
314 tick_size = num.min(size) / 200.
316 return fontsize, dict(
317 MAP_FRAME_PEN='%s' % color,
318 MAP_TICK_PEN_PRIMARY='%gp,%s' % (pen_size, color),
319 MAP_TICK_PEN_SECONDARY='%gp,%s' % (pen_size, color),
320 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size,
321 MAP_TICK_LENGTH_SECONDARY='%gc' % (tick_size * 3),
322 FONT_ANNOT_PRIMARY='%s-Bold,%s' % (font, color),
323 FONT_LABEL='%s-Bold,%s' % (font, color),
324 FONT_TITLE='%s-Bold,%s' % (font, color),
325 PS_CHAR_ENCODING='ISOLatin1+',
326 MAP_FRAME_TYPE='fancy',
327 FORMAT_GEO_MAP='D',
328 PS_PAGE_ORIENTATION='portrait',
329 MAP_GRID_PEN_PRIMARY='thinnest,%s' % color,
330 MAP_ANNOT_OBLIQUE='6')
333class SourceError(Exception):
334 pass
337class RuptureMap(Map):
338 ''' Map plotting of attributes and results of the
339 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
340 '''
342 def __init__(
343 self,
344 source=None,
345 fontcolor='darkslategrey',
346 width=20.,
347 height=14.,
348 margins=None,
349 color_wet=(216, 242, 254),
350 color_dry=(238, 236, 230),
351 topo_cpt_wet='light_sea_uniform',
352 topo_cpt_dry='light_land_uniform',
353 show_cities=False,
354 *args,
355 **kwargs):
357 size = (width, height)
358 fontsize, gmt_config = _make_gmt_conf(fontcolor, size)
360 if margins is None:
361 margins = [
362 fontsize * 0.15, num.min(size) / 200.,
363 num.min(size) / 200., fontsize * 0.05]
365 Map.__init__(self, *args, margins=margins, width=width, height=height,
366 gmt_config=gmt_config,
367 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet,
368 color_wet=color_wet, color_dry=color_dry,
369 **kwargs)
371 if show_cities:
372 self.draw_cities()
374 self._source = source
375 self._fontcolor = fontcolor
376 self._fontsize = fontsize
377 self.axes_layout = 'WeSn'
379 @property
380 def size(self):
381 '''
382 Figure size [cm]
383 '''
385 return (self.width, self.height)
387 @property
388 def font(self):
389 '''
390 Font style (size and type)
391 '''
393 return '%sp,%s' % (self._fontsize, fonttype)
395 @property
396 def source(self):
397 '''
398 PseudoDynamicRupture whose attributes are plotted.
400 Note, that source.patches attribute needs to be calculated
401 :type source: :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
402 '''
404 if self._source is None:
405 raise SourceError('No source given. Please define it!')
407 if not isinstance(self._source, PseudoDynamicRupture):
408 raise SourceError('This function is only capable for a source of'
409 ' type: %s' % type(PseudoDynamicRupture()))
411 if not self._source.patches:
412 raise TypeError('No source patches are defined. Please run'
413 '"source.discretize_patches()" on your source')
415 return self._source
417 @source.setter
418 def source(self, source):
419 self._source = source
421 def _get_topotile(self):
422 if self._dems is None:
423 self._setup()
425 try:
426 t, _ = self._get_topo_tile('land')
427 except NoTopo:
428 wesn = self.wesn
430 nx = int(num.ceil(
431 self.width * cm2inch * self.topo_resolution_max))
432 ny = int(num.ceil(
433 self.height * cm2inch * self.topo_resolution_max))
435 data = num.zeros((nx, ny))
437 t = Tile(wesn[0], wesn[2],
438 (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny,
439 data)
441 return t
443 def _patches_to_lw(self):
444 '''
445 Generate regular rect. length-width grid based on the patch distance
447 Prerequesite is a regular grid of patches (constant lengths and widths)
448 Both coordinates are given relative to the source anchor point [in m]
449 The grid is extended from the patch centres to the edges of the source
451 :returns: lengths along strike, widths downdip
452 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
453 '''
455 source = self.source
456 patches = source.patches
458 patch_l, patch_w = patches[0].length, patches[0].width
460 patch_lengths = num.concatenate((
461 num.array([0.]),
462 num.array([il*patch_l+patch_l/2. for il in range(source.nx)]),
463 num.array([patch_l * source.nx])))
465 patch_widths = num.concatenate((
466 num.array([0.]),
467 num.array([iw*patch_w+patch_w/2. for iw in range(source.ny)]),
468 num.array([patch_w * source.ny])))
470 ax, ay = map_anchor[source.anchor]
472 patch_lengths -= source.length * (ax + 1.) / 2.
473 patch_widths -= source.width * (ay + 1.) / 2.
475 return patch_lengths, patch_widths
477 def _xy_to_lw(self, x, y):
478 '''
479 Generate regular rect. length-width grid based on the xy coordinates
481 Prerequesite is a regular grid with constant dx and dy. x and y are
482 relative coordinates on the rupture plane (range -1:1) along strike and
483 downdip.
484 Length and width are obtained relative to the source anchor point
485 [in m].
487 :returns: lengths along strike [m], widths downdip [m]
488 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
489 '''
491 x, y = num.unique(x), num.unique(y)
492 dx, dy = x[1] - x[0], y[1] - y[0]
494 if any(num.abs(num.diff(x) - dx) >= 1e-6):
495 raise ValueError('Regular grid with constant spacing needed.'
496 ' Please check the x coordinates.')
498 if any(num.abs(num.diff(y) - dy) >= 1e-6):
499 raise ValueError('Regular grid with constant spacing needed.'
500 ' Please check the y coordinates.')
502 return xy_to_lw(self.source, x, y)
504 def _tile_to_lw(self, ref_lat, ref_lon,
505 north_shift=0., east_shift=0., strike=0.):
507 '''
508 Coordinate transformation from the topo tile grid into length-width
510 The topotile latlon grid is rotated into the length width grid. The
511 width is defined here as its horizontal component. The resulting grid
512 is used for interpolation of grid data.
514 :param ref_lat: Reference latitude, from which length-width relative
515 coordinatesgrid are calculated
516 :type ref_lat: float
517 :param ref_lon: Reference longitude, from which length-width relative
518 coordinatesgrid are calculated
519 :type ref_lon: float
520 :param north_shift: North shift of the reference point [m]
521 :type north_shift: optional, float
522 :param east_shift: East shift of the reference point [m]
523 :type east_shift: optional, float
524 :param strike: striking of the length axis compared to the North axis
525 [degree]
526 :type strike: optional, float
528 :returns: topotile grid nodes as array of length-width coordinates
529 :rtype: :py:class:`numpy.ndarray`, ``(n_nodes, 2)``
530 '''
532 t = self._get_topotile()
533 grid_lats = t.y()
534 grid_lons = t.x()
536 meshgrid_lons, meshgrid_lats = num.meshgrid(grid_lons, grid_lats)
538 grid_northing, grid_easting = pod.latlon_to_ne_numpy(
539 ref_lat, ref_lon, meshgrid_lats.flatten(), meshgrid_lons.flatten())
541 grid_northing -= north_shift
542 grid_easting -= east_shift
544 strike *= d2r
545 sin, cos = num.sin(strike), num.cos(strike)
546 points_out = num.zeros((grid_northing.shape[0], 2))
547 points_out[:, 0] = -sin * grid_northing + cos * grid_easting
548 points_out[:, 1] = cos * grid_northing + sin * grid_easting
550 return points_out
552 def _prep_patch_grid_data(self, data):
553 '''
554 Extend patch data from patch centres to the outer source edges
556 Patch data is always defined in the centre of the patches. For
557 interpolation the data is extended here to the edges of the rupture
558 plane.
560 :param data: Patch wise data
561 :type data: :py:class:`numpy.ndarray`
563 :returns: Extended data array
564 :rtype: :py:class:`numpy.ndarray`
565 '''
567 shape = (self.source.nx + 2, self.source.ny + 2)
568 data = data.reshape(self.source.nx, self.source.ny).copy()
570 data_new = num.zeros(shape)
571 data_new[1:-1, 1:-1] = data
572 data_new[1:-1, 0] = data[:, 0]
573 data_new[1:-1, -1] = data[:, -1]
574 data_new[0, 1:-1] = data[0, :]
575 data_new[-1, 1:-1] = data[-1, :]
577 for i, j in zip([-1, -1, 0, 0], [-1, 0, -1, 0]):
578 data_new[i, j] = data[i, j]
580 return data_new
582 def _regular_data_to_grid(self, lengths, widths, data, filename):
583 '''
584 Interpolate regular data onto topotile grid
586 Regular gridded data is interpolated onto the latlon grid of the
587 topotile. It is then stored as a gmt-readable .grd-file.
589 :param lengths: Grid coordinates along strike relative to anchor [m]
590 :type lengths: :py:class:`numpy.ndarray`
591 :param widths: Grid coordinates downdip relative to anchor [m]
592 :type widths: :py:class:`numpy.ndarray`
593 :param data: Data grid array
594 :type data: :py:class:`numpy.ndarray`
595 :param filename: Filename, where grid is stored
596 :type filename: str
597 '''
599 source = self.source
601 interpolator = scrgi(
602 (widths * num.cos(d2r * source.dip), lengths),
603 data.T,
604 bounds_error=False,
605 method='nearest')
607 points_out = self._tile_to_lw(
608 ref_lat=source.lat,
609 ref_lon=source.lon,
610 north_shift=source.north_shift,
611 east_shift=source.east_shift,
612 strike=source.strike)
614 t = self._get_topotile()
615 t.data = num.zeros_like(t.data, dtype=num.float)
616 t.data[:] = num.nan
618 t.data = interpolator(points_out).reshape(t.data.shape)
620 _save_grid(t.y(), t.x(), t.data, filename=filename)
622 def patch_data_to_grid(self, data, *args, **kwargs):
623 '''
624 Generate a grid file based on regular patch wise data.
626 :param data: Patchwise data grid array
627 :type data: :py:class:`numpy.ndarray`
628 '''
630 lengths, widths = self._patches_to_lw()
632 data_new = self._prep_patch_grid_data(data)
634 self._regular_data_to_grid(lengths, widths, data_new, *args, **kwargs)
636 def xy_data_to_grid(self, x, y, data, *args, **kwargs):
637 '''
638 Generate a grid file based on regular gridded data using xy coordinates
640 Convert a grid based on relative fault coordinates (range -1:1) along
641 strike (x) and downdip (y) into a .grd file.
643 :param x: Relative point coordinate along strike (range: -1:1)
644 :type x: float or :py:class:`numpy.ndarray`
645 :param y: Relative downdip point coordinate (range: -1:1)
646 :type y: float or :py:class:`numpy.ndarray`
647 :param data: Patchwise data grid array
648 :type data: :py:class:`numpy.ndarray`
649 '''
651 lengths, widths = self._xy_to_lw(x, y)
653 self._regular_data_to_grid(
654 lengths, widths, data.reshape((lengths.shape[0], widths.shape[0])),
655 *args, **kwargs)
657 def draw_image(self, gridfile, cmap, cbar=True, **kwargs):
658 '''
659 Draw grid data as image and include, if whished, a colorbar
661 :param gridfile: File of the grid which shall be plotted
662 :type gridfile: str
663 :param cmap: Name of the colormap, which shall be used. A .cpt-file
664 "my_<cmap>.cpt" must exist
665 :type cmap: str
666 :param cbar: If True, a colorbar corresponding to the grid data is
667 added. Keywordarguments are parsed to it.
668 :type cbar: optional, bool
669 '''
671 self.gmt.grdimage(
672 gridfile,
673 C='my_%s.cpt' % cmap,
674 E='200',
675 Q=True,
676 n='+t0.0',
677 *self.jxyr)
679 if cbar:
680 self.draw_colorbar(cmap=cmap, **kwargs)
682 def draw_contour(
683 self,
684 gridfile,
685 contour_int,
686 anot_int,
687 angle=None,
688 unit='',
689 color='',
690 style='',
691 **kwargs):
693 '''
694 Draw grid data as contour lines
696 :param gridfile: File of the grid which shall be plotted
697 :type gridfile: str
698 :param contour_int: Interval of contour lines in units of the gridfile
699 :type contour_int: float
700 :param anot_int: Interval of labelled contour lines in units of the
701 gridfile. Must be a integer multiple of contour_int
702 :type anot_int: float
703 :param angle: Rotation angle of the labels [degree]
704 :type angle: optional, float
705 :param unit: Name of the unit in the grid file. It will be displayed
706 behind the label on labelled contour lines
707 :type unit: optional, str
708 :param color: GMT readable color code/str of the contour lines
709 :type color: optional, str
710 :param style: Line style of the contour lines. If not given, solid
711 lines are plotted
712 :type style: optional, str
713 '''
715 pen_size = self._fontsize / 40.
717 if not color:
718 color = self._fontcolor
720 a_string = '%g+f%s,%s+r%gc+u%s' % (
721 anot_int, self.font, color, pen_size*4, unit)
722 if angle:
723 a_string += '+a%g' % angle
724 c_string = '%g' % contour_int
726 if kwargs:
727 kwargs['A'], kwargs['C'] = a_string, c_string
728 else:
729 kwargs = dict(A=a_string, C=c_string)
731 if style:
732 style = ',' + style
734 args = ['-Wc%gp,%s%s+s' % (pen_size, color, style)]
736 self.gmt.grdcontour(
737 gridfile,
738 S='10',
739 W='a%gp,%s%s+s' % (pen_size*4, color, style),
740 *self.jxyr + args,
741 **kwargs)
743 def draw_colorbar(self, cmap, label='', anchor='top_right', **kwargs):
744 '''
745 Draw a colorbar based on a existing colormap
747 :param cmap: Name of the colormap, which shall be used. A .cpt-file
748 "my_<cmap>.cpt" must exist
749 :type cmap: str
750 :param label: Title of the colorbar
751 :type label: optional, str
752 :param anchor: Placement of the colorbar. Combine 'top', 'center' and
753 'bottom' with 'left', None for middle and 'right'
754 :type anchor: optional, str
755 '''
757 if not kwargs:
758 kwargs = {}
760 if label:
761 kwargs['B'] = 'af+l%s' % label
763 kwargs['C'] = 'my_%s.cpt' % cmap
764 a_str = cbar_anchor[anchor]
766 w = self.width / 3.
767 h = w / 10.
769 lgap = rgap = w / 10.
770 bgap, tgap = h, h / 10.
772 dx, dy = 2.5 * lgap, 2. * tgap
774 if 'bottom' in anchor:
775 dy += 4 * h
777 self.gmt.psscale(
778 D='j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy),
779 F='+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap),
780 *self.jxyr,
781 **kwargs)
783 def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs):
784 ''' Draw vectors based on two grid files
786 Two grid files for vector lengths in x and y need to be given. The
787 function calls gmt.grdvector. All arguments defined for this function
788 in gmt can be passed as keyword arguments. Different standard settings
789 are applied if not defined differently.
791 :param x_gridfile: File of the grid defining vector lengths in x
792 :type x_gridfile: str
793 :param y_gridfile: File of the grid defining vector lengths in y
794 :type y_gridfile: str
795 :param vcolor: Vector face color as defined in "G" option
796 :type vcolor: str
797 '''
799 kwargs['S'] = kwargs.get('S', 'il1.')
800 kwargs['I'] = kwargs.get('I', 'x20')
801 kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black')
802 kwargs['Q'] = kwargs.get('Q', '4c+e+n5c+h1')
804 self.gmt.grdvector(
805 x_gridfile, y_gridfile,
806 G='%s' % 'lightgrey' if not vcolor else vcolor,
807 *self.jxyr,
808 **kwargs)
810 def draw_dynamic_data(self, data, **kwargs):
811 '''
812 Draw an image of any data gridded on the patches e.g dislocation
814 :param data: Patchwise data grid array
815 :type data: :py:class:`numpy.ndarray`
816 '''
818 plot_data = data
820 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
822 clim = kwargs.pop('clim', (plot_data.min(), plot_data.max()))
824 cpt = []
825 if not op.exists('my_%s.cpt' % kwargs['cmap']):
826 make_colormap(self.gmt, clim[0], clim[1],
827 cmap=kwargs['cmap'], space=False)
828 cpt = [kwargs['cmap']]
830 tmp_grd_file = 'tmpdata.grd'
831 self.patch_data_to_grid(plot_data, tmp_grd_file)
832 self.draw_image(tmp_grd_file, **kwargs)
834 clear_temp(gridfiles=[tmp_grd_file], cpts=cpt)
836 def draw_patch_parameter(self, attribute, **kwargs):
837 '''Draw an image of a chosen patch attribute e.g traction
839 :param attribute: Patch attribute, which is plotted. All patch
840 attributes can be taken (see doc of
841 :py:class:`pyrocko.modelling.okada.OkadaSource`) and also
842 ``traction``, ``tx``, ``ty`` or ``tz`` to display the
843 length or the single components of the traction vector.
844 :type attribute: str
845 '''
847 a = attribute
848 source = self.source
850 if a == 'traction':
851 data = num.linalg.norm(source.get_tractions(), axis=1)
852 elif a == 'tx':
853 data = source.get_tractions()[:, 0]
854 elif a == 'ty':
855 data = source.get_tractions()[:, 1]
856 elif a == 'tz':
857 data = source.get_tractions()[:, 2]
858 else:
859 data = source.get_patch_attribute(attribute)
861 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
862 data *= factor
864 kwargs['label'] = kwargs.get(
865 'label',
866 '%s [%s]' % (a, cbar_helper[a]['unit']))
868 self.draw_dynamic_data(data, **kwargs)
870 def draw_time_contour(self, store, clevel=[], **kwargs):
871 '''Draw high contour lines of the rupture front propgation time
873 :param store: Greens function store, which is used for time calculation
874 :type store: :py:class:`pyrocko.gf.store.Store`
875 :param clevel: List of times, for which contour lines are drawn,
876 optional
877 :type clevel: list of float
878 '''
880 _, _, _, _, points_xy = self.source._discretize_points(store, cs='xyz')
881 _, _, times, _, _, _ = self.source.get_vr_time_interpolators(store)
883 scaler = AutoScaler(mode='0-max', approx_ticks=8)
884 scale = scaler.make_scale([num.min(times), num.max(times)])
886 if clevel:
887 if len(clevel) > 1:
888 kwargs['anot_int'] = num.min(num.diff(clevel))
889 else:
890 kwargs['anot_int'] = clevel[0]
892 kwargs['contour_int'] = kwargs['anot_int']
893 kwargs['L'] = '0/%g' % num.max(clevel)
895 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
896 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
897 kwargs['unit'] = kwargs.get('unit', cbar_helper['time']['unit'])
898 kwargs['L'] = kwargs.get('L', '0/%g' % (num.max(times) + 1.))
899 kwargs['G'] = kwargs.get('G', 'n2/3c')
901 tmp_grd_file = 'tmpdata.grd'
902 self.xy_data_to_grid(points_xy[:, 0], points_xy[:, 1],
903 times, tmp_grd_file)
904 self.draw_contour(tmp_grd_file, **kwargs)
906 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
908 def draw_points(self, lats, lons, symbol='point', size=None, **kwargs):
909 '''Draw points at given locations
911 :param lats: Point latitude coordinates [degree]
912 :type lats: iterable
913 :param lons: Point longitude coordinates [degree]
914 :type lons: iterable
915 :param symbol: Define symbol of the points
916 (``'star', 'circle', 'point', 'triangle'``) - default is ``point``
917 :type symbol: optional, str
918 :param size: Size of the points in points
919 :type size: optional, float
920 '''
922 sym_to_gmt = dict(
923 star='a',
924 circle='c',
925 point='p',
926 triangle='t')
928 lats = num.atleast_1d(lats)
929 lons = num.atleast_1d(lons)
931 if lats.shape[0] != lons.shape[0]:
932 raise IndexError('lats and lons do not have the same shape!')
934 if size is None:
935 size = self._fontsize / 3.
937 kwargs['S'] = kwargs.get('S', sym_to_gmt[symbol] + '%gp' % size)
938 kwargs['G'] = kwargs.get('G', gmtpy.color('scarletred2'))
939 kwargs['W'] = kwargs.get('W', '2p,%s' % self._fontcolor)
941 self.gmt.psxy(
942 in_columns=[lons, lats],
943 *self.jxyr,
944 **kwargs)
946 def draw_nucleation_point(self, **kwargs):
947 ''' Plot the nucleation point onto the map '''
949 nlat, nlon = xy_to_latlon(
950 self.source, self.source.nucleation_x, self.source.nucleation_y)
952 self.draw_points(nlat, nlon, **kwargs)
954 def draw_dislocation(self, time=None, component='', **kwargs):
955 ''' Draw dislocation onto map at any time
957 For a given time (if ``time`` is ``None``, ``tmax`` is used)
958 and given component the patchwise dislocation is plotted onto the map.
960 :param time: time after origin, for which dislocation is computed. If
961 ``None``, ``tmax`` is taken.
962 :type time: optional, float
963 :param component: Dislocation component, which shall be plotted: ``x``
964 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
965 length of the dislocation vector is plotted
966 '''
968 disl = self.source.get_slip(time=time)
970 if component:
971 data = disl[:, c2disl[component]]
972 else:
973 data = num.linalg.norm(disl, axis=1)
975 kwargs['label'] = kwargs.get(
976 'label', 'u%s [m]' % (component))
978 self.draw_dynamic_data(data, **kwargs)
980 def draw_dislocation_contour(
981 self, time=None, component=None, clevel=[], **kwargs):
982 ''' Draw dislocation contour onto map at any time
984 For a given time (if ``time`` is ``None``, ``tmax`` is used)
985 and given component the patchwise dislocation is plotted as contour
986 onto the map.
988 :param time: time after origin, for which dislocation is computed. If
989 ``None``, ``tmax`` is taken.
990 :type time: optional, float
991 :param component: Dislocation component, which shall be plotted: ``x``
992 along strike, ``y`` along updip, ``z`` normal``. If ``None``,
993 the length of the dislocation vector is plotted
994 :param clevel: List of times, for which contour lines are drawn
995 :type clevel: optional, list of float
996 '''
998 disl = self.source.get_slip(time=time)
1000 if component:
1001 data = disl[:, c2disl[component]]
1002 else:
1003 data = num.linalg.norm(disl, axis=1)
1005 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1006 scale = scaler.make_scale([num.min(data), num.max(data)])
1008 if clevel:
1009 if len(clevel) > 1:
1010 kwargs['anot_int'] = num.min(num.diff(clevel))
1011 else:
1012 kwargs['anot_int'] = clevel[0]
1014 kwargs['contour_int'] = kwargs['anot_int']
1015 kwargs['L'] = '%g/%g' % (
1016 num.min(clevel) - kwargs['contour_int'],
1017 num.max(clevel) + kwargs['contour_int'])
1018 else:
1019 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1020 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1021 kwargs['L'] = kwargs.get('L', '%g/%g' % (
1022 num.min(data) - 1., num.max(data) + 1.))
1024 kwargs['unit'] = kwargs.get('unit', ' m')
1025 kwargs['G'] = kwargs.get('G', 'n2/3c')
1027 tmp_grd_file = 'tmpdata.grd'
1028 self.patch_data_to_grid(data, tmp_grd_file)
1029 self.draw_contour(tmp_grd_file, **kwargs)
1031 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
1033 def draw_dislocation_vector(self, time=None, **kwargs):
1034 '''Draw vector arrows onto map indicating direction of dislocation
1036 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1037 and given component the dislocation is plotted as vectors onto the map.
1039 :param time: time after origin [s], for which dislocation is computed.
1040 If ``None``, ``tmax`` is used.
1041 :type time: optional, float
1042 '''
1044 disl = self.source.get_slip(time=time)
1046 p_strike = self.source.get_patch_attribute('strike') * d2r
1047 p_dip = self.source.get_patch_attribute('dip') * d2r
1049 disl_dh = num.cos(p_dip) * disl[:, 1]
1050 disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh
1051 disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh
1053 tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')]
1055 self.patch_data_to_grid(disl_n, tmp_grd_files[0])
1056 self.patch_data_to_grid(disl_e, tmp_grd_files[1])
1058 self.draw_vector(
1059 tmp_grd_files[1], tmp_grd_files[0],
1060 **kwargs)
1062 clear_temp(gridfiles=tmp_grd_files, cpts=[])
1064 def draw_top_edge(self, **kwargs):
1065 '''Indicate rupture top edge on map
1066 '''
1068 outline = self.source.outline(cs='latlondepth')
1069 top_edge = outline[:2, :]
1071 kwargs = kwargs or {}
1072 kwargs['W'] = kwargs.get(
1073 'W', '%gp,%s' % (self._fontsize / 10., gmtpy.color('orange3')))
1075 self.gmt.psxy(
1076 in_columns=[top_edge[:, 1], top_edge[:, 0]],
1077 *self.jxyr,
1078 **kwargs)
1081class RuptureView(Object):
1082 ''' Plot of attributes and results of the
1083 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1084 '''
1086 _patches_to_lw = RuptureMap._patches_to_lw
1088 def __init__(self, source=None, figsize=None, fontsize=None):
1089 self._source = source
1090 self._axes = None
1092 self.figsize = figsize or mpl_papersize('halfletter', 'landscape')
1093 self.fontsize = fontsize or 10
1094 mpl_init(fontsize=self.fontsize)
1096 self._fig = None
1097 self._axes = None
1098 self._is_1d = False
1100 @property
1101 def source(self):
1102 ''' PseudoDynamicRupture whose attributes are plotted.
1104 Note, that source.patches attribute needs to be calculated for
1105 :type source: :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
1106 '''
1108 if self._source is None:
1109 raise SourceError('No source given. Please define it!')
1111 if not isinstance(self._source, PseudoDynamicRupture):
1112 raise SourceError('This function is only capable for a source of'
1113 ' type: %s' % type(PseudoDynamicRupture()))
1115 if not self._source.patches:
1116 raise TypeError('No source patches are defined. Please run'
1117 '\"discretize_patches\" on your source')
1119 return self._source
1121 @source.setter
1122 def source(self, source):
1123 self._source = source
1125 def _setup(
1126 self,
1127 title='',
1128 xlabel='',
1129 ylabel='',
1130 aspect=1.,
1131 spatial_plot=True,
1132 **kwargs):
1134 if self._fig is not None and self._axes is not None:
1135 return
1137 self._fig = plt.figure(figsize=self.figsize)
1139 self._axes = self._fig.add_subplot(1, 1, 1, aspect=aspect)
1140 ax = self._axes
1142 if ax is not None:
1143 ax.set_title(title)
1144 ax.grid(alpha=.3)
1145 ax.set_xlabel(xlabel)
1146 ax.set_ylabel(ylabel)
1148 if spatial_plot:
1149 ax.xaxis.set_major_formatter(FuncFormatter(km_formatter))
1150 ax.yaxis.set_major_formatter(FuncFormatter(km_formatter))
1152 def _clear_all(self):
1153 plt.cla()
1154 plt.clf()
1155 plt.close()
1157 self._fig, self._axes = None, None
1159 def _draw_scatter(self, x, y, *args, **kwargs):
1160 default_kwargs = dict(
1161 linewidth=0,
1162 marker='o',
1163 markerfacecolor=mpl_color('skyblue2'),
1164 markersize=6.,
1165 markeredgecolor=mpl_color('skyblue3'))
1166 default_kwargs.update(kwargs)
1168 if self._axes is not None:
1169 self._axes.plot(x, y, *args, **default_kwargs)
1171 def _draw_image(self, length, width, data, *args, **kwargs):
1172 if self._axes is not None:
1173 if 'extent' not in kwargs:
1174 kwargs['extent'] = [
1175 num.min(length), num.max(length),
1176 num.max(width), num.min(width)]
1178 im = self._axes.imshow(
1179 data,
1180 *args,
1181 interpolation='none',
1182 vmin=kwargs.get('clim', [None])[0],
1183 vmax=kwargs.get('clim', [None, None])[1],
1184 **kwargs)
1186 del kwargs['extent']
1187 if 'aspect' in kwargs:
1188 del kwargs['aspect']
1189 if 'clim' in kwargs:
1190 del kwargs['clim']
1191 if 'cmap' in kwargs:
1192 del kwargs['cmap']
1194 plt.colorbar(
1195 im, *args, shrink=0.9, pad=0.03, aspect=15., **kwargs)
1197 def _draw_contour(self, x, y, data, clevel=None, unit='', *args, **kwargs):
1198 setup_kwargs = dict(
1199 xlabel=kwargs.pop('xlabel', 'along strike [km]'),
1200 ylabel=kwargs.pop('ylabel', 'down dip [km]'),
1201 title=kwargs.pop('title', ''),
1202 aspect=kwargs.pop('aspect', 1))
1204 self._setup(**setup_kwargs)
1206 if self._axes is not None:
1207 if clevel is None:
1208 scaler = AutoScaler(mode='min-max')
1209 scale = scaler.make_scale([data.min(), data.max()])
1211 clevel = num.arange(scale[0], scale[1] + scale[2], scale[2])
1213 if not isinstance(clevel, num.ndarray):
1214 clevel = num.array([clevel])
1216 clevel = clevel[clevel < data.max()]
1218 cont = self._axes.contour(
1219 x, y, data, clevel, *args,
1220 linewidths=1.5, **kwargs)
1222 plt.setp(cont.collections, path_effects=[
1223 patheffects.withStroke(linewidth=2.0, foreground="beige"),
1224 patheffects.Normal()])
1226 clabels = self._axes.clabel(
1227 cont, clevel[::2], *args,
1228 inline=1, fmt='%g' + '%s' % unit,
1229 inline_spacing=15, rightside_up=True, use_clabeltext=True,
1230 **kwargs)
1232 plt.setp(clabels, path_effects=[
1233 patheffects.withStroke(linewidth=1.25, foreground="beige"),
1234 patheffects.Normal()])
1236 def draw_points(self, length, width, *args, **kwargs):
1237 ''' Draw a point onto the figure.
1239 Args and kwargs can be defined according to
1240 :py:func:`matplotlib.pyplot.scatter`.
1242 :param length: Point(s) coordinate on the rupture plane along strike
1243 relative to the anchor point [m]
1244 :type length: float, :py:class:`numpy.ndarray`
1245 :param width: Point(s) coordinate on the rupture plane along downdip
1246 relative to the anchor point [m]
1247 :type width: float, :py:class:`numpy.ndarray`
1248 '''
1250 if self._axes is not None:
1251 kwargs['color'] = kwargs.get('color', mpl_color('scarletred2'))
1252 kwargs['s'] = kwargs.get('s', 100.)
1254 self._axes.scatter(length, width, *args, **kwargs)
1256 def draw_dynamic_data(self, data, **kwargs):
1257 ''' Draw an image of any data gridded on the patches e.g dislocation
1259 :param data: Patchwise data grid array
1260 :type data: :py:class:`numpy.ndarray`
1261 '''
1263 plot_data = data
1265 anchor_x, anchor_y = map_anchor[self.source.anchor]
1267 length, width = xy_to_lw(
1268 self.source, num.array([-1., 1.]), num.array([-1., 1.]))
1270 data = plot_data.reshape(self.source.nx, self.source.ny).T
1272 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
1274 setup_kwargs = dict(
1275 xlabel='along strike [km]',
1276 ylabel='down dip [km]',
1277 title='',
1278 aspect=1)
1279 setup_kwargs.update(kwargs)
1281 kwargs = {k: v for k, v in kwargs.items() if k not in
1282 ('xlabel', 'ylabel', 'title')}
1283 self._setup(**setup_kwargs)
1284 self._draw_image(length=length, width=width, data=data, **kwargs)
1286 def draw_patch_parameter(self, attribute, **kwargs):
1287 ''' Draw an image of a chosen patch attribute e.g traction
1289 :param attribute: Patch attribute, which is plotted. All patch
1290 attributes can be taken (see doc of
1291 :py:class:`pyrocko.modelling.okada.OkadaSource`) and also
1292 ``'traction', 'tx', 'ty', 'tz'`` to display the length
1293 or the single components of the traction vector.
1294 :type attribute: str
1295 '''
1297 a = attribute
1298 source = self.source
1300 if a == 'traction':
1301 data = num.linalg.norm(source.get_tractions(), axis=1)
1302 elif a == 'tx':
1303 data = source.get_tractions()[:, 0]
1304 elif a == 'ty':
1305 data = source.get_tractions()[:, 1]
1306 elif a == 'tz':
1307 data = source.get_tractions()[:, 2]
1308 else:
1309 data = source.get_patch_attribute(attribute)
1311 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1312 data *= factor
1314 kwargs['label'] = kwargs.get(
1315 'label',
1316 '%s [%s]' % (a, cbar_helper[a]['unit']))
1318 return self.draw_dynamic_data(data=data, **kwargs)
1320 def draw_time_contour(self, store, clevel=[], **kwargs):
1321 ''' Draw high resolution contours of the rupture front propgation time
1323 :param store: Greens function store, which is used for time calculation
1324 :type store: :py:class:`pyrocko.gf.store.Store`
1325 :param clevel: levels of the contour lines. If no levels are given,
1326 they are automatically computed based on tmin and tmax
1327 :type clevel: optional, list
1328 '''
1329 source = self.source
1330 default_kwargs = dict(
1331 colors='#474747ff'
1332 )
1333 default_kwargs.update(kwargs)
1335 _, _, _, _, points_xy = source._discretize_points(store, cs='xyz')
1336 _, _, _, _, time_interpolator, _ = source.get_vr_time_interpolators(
1337 store)
1339 times = time_interpolator.values
1341 scaler = AutoScaler(mode='min-max', approx_ticks=8)
1342 scale = scaler.make_scale([times.min(), times.max()])
1344 if len(clevel) == 0:
1345 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1347 points_x = time_interpolator.grid[0]
1348 points_y = time_interpolator.grid[1]
1350 self._draw_contour(points_x, points_y, data=times.T,
1351 clevel=clevel, unit='s', **default_kwargs)
1353 def draw_nucleation_point(self, **kwargs):
1354 ''' Draw the nucleation point onto the map '''
1356 nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y
1357 length, width = xy_to_lw(self.source, nuc_x, nuc_y)
1358 self.draw_points(length, width, marker='o', **kwargs)
1360 def draw_dislocation(self, time=None, component='', **kwargs):
1361 ''' Draw dislocation onto map at any time
1363 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1364 and given component the patchwise dislocation is plotted onto the map.
1366 :param time: time after origin [s], for which dislocation is computed.
1367 If ``None``, ``tmax`` is taken.
1368 :type time: optional, float
1369 :param component: Dislocation component, which shall be plotted: ``x``
1370 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
1371 length of the dislocation vector is plotted
1372 '''
1374 disl = self.source.get_slip(time=time)
1376 if component:
1377 data = disl[:, c2disl[component]]
1378 else:
1379 data = num.linalg.norm(disl, axis=1)
1381 kwargs['label'] = kwargs.get(
1382 'label', 'u%s [m]' % (component))
1384 self.draw_dynamic_data(data, **kwargs)
1386 def draw_dislocation_contour(
1387 self, time=None, component=None, clevel=[], **kwargs):
1388 ''' Draw dislocation contour onto map at any time
1390 For a given time (if time is ``None``, ``tmax`` is used) and given
1391 component the patchwise dislocation is plotted as contour onto the map.
1393 :param time: time after origin, for which dislocation is computed. If
1394 None, tmax is taken.
1395 :type time: optional, float
1396 :param component: Dislocation component, which shall be plotted. ``x``
1397 along strike, ``y`` along updip, ``z`` - normal. If ``None``
1398 is given, the length of the dislocation vector is plotted
1399 :type component: str
1400 '''
1402 disl = self.source.get_slip(time=time)
1404 if component:
1405 data = disl[:, c2disl[component]]
1406 else:
1407 data = num.linalg.norm(disl, axis=1)
1409 data = data.reshape(self.source.ny, self.source.nx, order='F')
1411 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1412 scale = scaler.make_scale([data.min(), data.max()])
1414 if len(clevel) == 0:
1415 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1417 anchor_x, anchor_y = map_anchor[self.source.anchor]
1419 length, width = self._patches_to_lw()
1420 length, width = length[1:-1], width[1:-1]
1422 kwargs['colors'] = kwargs.get('colors', '#474747ff')
1424 self._setup(**kwargs)
1425 self._draw_contour(
1426 length, width, data=data, clevel=clevel, unit='m', **kwargs)
1428 def draw_source_dynamics(
1429 self, variable, store, deltat=None, *args, **kwargs):
1430 ''' Display dynamic source parameter
1432 Fast inspection possibility for the cumulative moment and the source
1433 time function approximation (assuming equal paths between different
1434 patches and observation point - valid for an observation point in the
1435 far field perpendicular to the source strike), so the cumulative moment
1436 rate function.
1438 :param variable: Dynamic parameter, which shall be plotted. Choose
1439 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1440 :type variable: str
1441 :param store: Greens function store, whose store.config.deltat defines
1442 the time increment between two parameter snapshots. If store is not
1443 given, the time increment is defined is taken from deltat.
1444 :type store: :py:class:`pyrocko.gf.store.Store`
1445 :param deltat: Time increment between two parameter snapshots. If not
1446 given, store.config.deltat is used to define deltat
1447 :type deltat: optional, float
1448 '''
1450 v = variable
1452 data, times = self.source.get_moment_rate(store=store, deltat=deltat)
1453 deltat = num.concatenate([(num.diff(times)[0], ), num.diff(times)])
1455 if v in ('moment_rate', 'stf'):
1456 name, unit = 'dM/dt', 'Nm/s'
1457 elif v in ('cumulative_moment', 'moment'):
1458 data = num.cumsum(data * deltat)
1459 name, unit = 'M', 'Nm'
1460 else:
1461 raise ValueError('No dynamic data for given variable %s found' % v)
1463 self._setup(xlabel='time [s]',
1464 ylabel='%s / %.2g %s' % (name, data.max(), unit),
1465 aspect='auto',
1466 spatial_plot=False)
1467 self._draw_scatter(times, data/num.max(data), *args, **kwargs)
1468 self._is_1d = True
1470 def draw_patch_dynamics(
1471 self, variable, nx, ny, store=None, deltat=None, *args, **kwargs):
1472 ''' Display dynamic boundary element / patch parameter
1474 Fast inspection possibility for different dynamic parameter for a
1475 single patch / boundary element. The chosen parameter is plotted for
1476 the chosen patch.
1478 :param variable: Dynamic parameter, which shall be plotted. Choose
1479 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1480 :type variable: str
1481 :param nx: Patch index along strike (range: 0:self.source.nx - 1)
1482 :type nx: int
1483 :param nx: Patch index downdip (range: 0:self.source.ny - 1)
1484 :type nx: int
1485 :param store: Greens function store, whose store.config.deltat defines
1486 the time increment between two parameter snapshots. If store is not
1487 given, the time increment is defined is taken from deltat.
1488 :type store: optional, :py:class:`pyrocko.gf.store.Store`
1489 :param deltat: Time increment between two parameter snapshots. If not
1490 given, store.config.deltat is used to define deltat
1491 :type deltat: optional, float
1492 '''
1494 v = variable
1495 source = self.source
1496 idx = nx * source.ny + nx
1498 m = re.match(r'dislocation_([xyz])', v)
1500 if v in ('moment_rate', 'cumulative_moment', 'moment', 'stf'):
1501 data, times = source.get_moment_rate_patches(
1502 store=store, deltat=deltat)
1503 elif 'dislocation' in v or 'slip_rate' == v:
1504 ddisloc, times = source.get_delta_slip(store=store, deltat=deltat)
1506 if v in ('moment_rate', 'stf'):
1507 data, times = source.get_moment_rate_patches(
1508 store=store, deltat=deltat)
1509 name, unit = 'dM/dt', 'Nm/s'
1510 elif v in ('cumulative_moment', 'moment'):
1511 data, times = source.get_moment_rate_patches(
1512 store=store, deltat=deltat)
1513 data = num.cumsum(data, axis=1)
1514 name, unit = 'M', 'Nm'
1515 elif v == 'slip_rate':
1516 data, times = source.get_slip_rate(store=store, deltat=deltat)
1517 name, unit = 'du/dt', 'm/s'
1518 elif v == 'dislocation':
1519 data, times = source.get_delta_slip(store=store, deltat=deltat)
1520 data = num.linalg.norm(num.cumsum(data, axis=2), axis=1)
1521 name, unit = 'du', 'm'
1522 elif m:
1523 data, times = source.get_delta_slip(store=store, deltat=deltat)
1524 data = num.cumsum(data, axis=2)[:, c2disl[m.group(1)], :]
1525 name, unit = 'du%s' % m.group(1), 'm'
1526 else:
1527 raise ValueError('No dynamic data for given variable %s found' % v)
1529 self._setup(xlabel='time [s]',
1530 ylabel='%s / %.2g %s' % (name, num.max(data), unit),
1531 aspect='auto',
1532 spatial_plot=False)
1533 self._draw_scatter(times, data[idx, :]/num.max(data),
1534 *args, **kwargs)
1535 self._is_1d = True
1537 def finalize(self):
1538 if self._is_1d:
1539 return
1541 length, width = xy_to_lw(
1542 self.source, num.array([-1., 1.]), num.array([1., -1.]))
1544 self._axes.set_xlim(length)
1545 self._axes.set_ylim(width)
1547 def gcf(self):
1548 self.finalize()
1549 return self._fig
1551 def save(self, filename, dpi=None):
1552 ''' Save plot to file
1554 :param filename: filename and path, where the plot is stored at
1555 :type filename: str
1556 :param dpi: Resolution of the output plot [dpi]
1557 :type dpi: int
1558 '''
1559 self.finalize()
1560 try:
1561 self._fig.savefig(fname=filename, dpi=dpi, bbox_inches='tight')
1562 except TypeError:
1563 self._fig.savefig(filename=filename, dpi=dpi, bbox_inches='tight')
1565 self._clear_all()
1567 def show_plot(self):
1568 ''' Show plot '''
1569 self.finalize()
1570 plt.show()
1571 self._clear_all()
1574def render_movie(fn_path, output_path, framerate=20):
1575 ''' Generate a mp4 movie based on given png files using
1576 `ffmpeg <https://ffmpeg.org>`_.
1578 Render a movie based on a set of given .png files in fn_path. All files
1579 must have a filename specified by ``fn_path`` (e.g. giving ``fn_path`` with
1580 ``/temp/f%04.png`` a valid png filename would be ``/temp/f0001.png``). The
1581 files must have a numbering, indicating their order in the movie.
1583 :param fn_path: Path and fileformat specification of the input .png files.
1584 :type fn_path: str
1585 :param output_path: Path and filename of the output ``.mp4`` movie file
1586 :type output_path: str
1587 :param deltat: Time between individual frames (``1 / framerate``) [s]
1588 :type deltat: optional, float
1590 '''
1591 try:
1592 check_call(['ffmpeg', '-loglevel', 'panic'])
1593 except CalledProcessError:
1594 pass
1595 except (TypeError, IOError):
1596 logger.warning(
1597 'Package ffmpeg needed for movie rendering. Please install it '
1598 '(e.g. on linux distr. via sudo apt-get ffmpeg) and retry.')
1599 return False
1601 try:
1602 check_call([
1603 'ffmpeg', '-loglevel', 'info', '-y',
1604 '-framerate', '%g' % framerate,
1605 '-i', fn_path,
1606 '-vcodec', 'libx264',
1607 '-preset', 'medium',
1608 '-tune', 'animation',
1609 '-pix_fmt', 'yuv420p',
1610 '-movflags', '+faststart',
1611 '-filter:v', "crop='floor(in_w/2)*2:floor(in_h/2)*2'",
1612 '-crf', '15',
1613 output_path])
1615 return True
1616 except CalledProcessError as e:
1617 logger.warning(e)
1618 return False
1621def render_gif(fn, output_path, loops=-1):
1622 ''' Generate a gif based on a given mp4 using ffmpeg
1624 Render a gif based on a given .mp4 movie file in ``fn`` path.
1626 :param fn: Path and file name of the input .mp4 file.
1627 :type fn: str
1628 :param output_path: Path and filename of the output animated gif file
1629 :type output_path: str
1630 :param loops: Number of gif repeats (loops). ``-1`` is not repetition,
1631 ``0`` infinite
1632 :type loops: optional, integer
1633 '''
1635 try:
1636 check_call(['ffmpeg', '-loglevel', 'panic'])
1637 except CalledProcessError:
1638 pass
1639 except (TypeError, IOError):
1640 logger.warning(
1641 'Package ffmpeg needed for movie rendering. Please install it '
1642 '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.')
1643 return False
1645 try:
1646 check_call([
1647 'ffmpeg', '-hide_banner', '-loglevel', 'panic', '-i',
1648 fn,
1649 '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse',
1650 '-loop', '%d' % loops,
1651 output_path])
1653 return True
1654 except CalledProcessError as e:
1655 logger.warning(e)
1656 return False
1659def rupture_movie(
1660 source,
1661 store,
1662 variable='dislocation',
1663 draw_time_contours=False,
1664 fn_path='.',
1665 prefix='',
1666 plot_type='map',
1667 deltat=None,
1668 framerate=None,
1669 store_images=False,
1670 render_as_gif=False,
1671 gif_loops=-1,
1672 **kwargs):
1673 ''' Generate a movie based on a given source for dynamic parameter
1675 Create a MPEG-4 movie or gif of one of the following dynamic parameters
1676 (``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1677 (along updip), ``dislocation_z`` (normal), ``slip_rate``, ``moment_rate``).
1678 If desired, the single snap shots can be stored as images as well.
1679 ``kwargs`` have to be given according to the chosen ``plot_type``.
1681 :param source: Pseudo dynamic rupture, for which the movie is produced
1682 :type source: :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
1683 :param store: Greens function store, which is used for time calculation. If
1684 deltat is not given, it is taken from the store.config.deltat
1685 :type store: :py:class:`pyrocko.gf.store.Store`
1686 :param variable: Dynamic parameter, which shall be plotted. Choose between
1687 ``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1688 (along updip), ``dislocation_z`` (normal), ``slip_rate`` and
1689 ``moment_rate``, default ``dislocation``
1690 :type variable: optional, str
1691 :param draw_time_contours: If True, corresponding isochrones are drawn on
1692 the each plots
1693 :type draw_time_contours: optional, bool
1694 :param fn_path: Absolut or relative path, where movie (and optional images)
1695 are stored
1696 :type fn_path: optional, str
1697 :param prefix: File prefix used for the movie (and optional image) files
1698 :type prefix: optional, str
1699 :param plot_type: Choice of plot type: ``map``, ``view`` (map plot using
1700 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap`
1701 or plane view using
1702 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`)
1703 :type plot_type: optional, str
1704 :param deltat: Time between parameter snapshots. If not given,
1705 store.config.deltat is used to define deltat
1706 :type deltat: optional, float
1707 :param store_images: Choice to store the single .png parameter snapshots in
1708 fn_path or not.
1709 :type store_images: optional, bool
1710 :param render_as_gif: If ``True``, the movie is converted into a gif. If
1711 ``False``, the movie is returned as mp4
1712 :type render_as_gif: optional, bool
1713 :param gif_loops: If ``render_as_gif`` is ``True``, a gif with
1714 ``gif_loops`` number of loops (repetitions) is returned.
1715 ``-1`` is no repetition, ``0`` infinite.
1716 :type gif_loops: optional, integer
1717 '''
1719 v = variable
1720 assert plot_type in ('map', 'view')
1722 if not source.patches:
1723 source.discretize_patches(store, interpolation='multilinear')
1725 if source.coef_mat is None:
1726 source.calc_coef_mat()
1728 prefix = prefix or v
1729 deltat = deltat or store.config.deltat
1730 framerate = max(framerate or int(1./deltat), 1)
1732 if v == 'moment_rate':
1733 data, times = source.get_moment_rate_patches(deltat=deltat)
1734 name, unit = 'dM/dt', 'Nm/s'
1735 elif 'dislocation' in v or 'slip_rate' == v:
1736 ddisloc, times = source.get_delta_slip(deltat=deltat)
1737 else:
1738 raise ValueError('No dynamic data for given variable %s found' % v)
1740 deltat = times[1] - times[0]
1742 m = re.match(r'dislocation_([xyz])', v)
1743 if m:
1744 data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]]
1745 name, unit = 'du%s' % m.group(1), 'm'
1746 elif v == 'dislocation':
1747 data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2)
1748 name, unit = 'du', 'm'
1749 elif v == 'slip_rate':
1750 data = num.linalg.norm(ddisloc, axis=2) / deltat
1751 name, unit = 'du/dt', 'm/s'
1753 if plot_type == 'map':
1754 plt_base = RuptureMap
1755 elif plot_type == 'view':
1756 plt_base = RuptureView
1757 else:
1758 raise AttributeError('invalid type: %s' % plot_type)
1760 attrs_base = get_kwargs(plt_base)
1761 kwargs_base = dict([k for k in kwargs.items() if k[0] in attrs_base])
1762 kwargs_plt = dict([k for k in kwargs.items() if k[0] not in kwargs_base])
1764 if 'clim' in kwargs_plt:
1765 data = num.clip(data, kwargs_plt['clim'][0], kwargs_plt['clim'][1])
1766 else:
1767 kwargs_plt['clim'] = [num.min(data), num.max(data)]
1769 if 'label' not in kwargs_plt:
1770 vmax = num.max(num.abs(kwargs_plt['clim']))
1771 data /= vmax
1773 kwargs_plt['label'] = '%s / %.2g %s' % (name, vmax, unit)
1774 kwargs_plt['clim'] = [i / vmax for i in kwargs_plt['clim']]
1776 with tempfile.TemporaryDirectory(suffix='pyrocko') as temp_path:
1777 fns_temp = [op.join(temp_path, 'f%09d.png' % (it + 1))
1778 for it, _ in enumerate(times)]
1779 fn_temp_path = op.join(temp_path, 'f%09d.png')
1781 for it, (t, ft) in enumerate(zip(times, fns_temp)):
1782 plot_data = data[:, it]
1784 plt = plt_base(source=source, **kwargs_base)
1785 plt.draw_dynamic_data(plot_data, **kwargs_plt)
1786 plt.draw_nucleation_point()
1788 if draw_time_contours:
1789 plt.draw_time_contour(store, clevel=[t])
1791 plt.save(ft)
1793 fn_mp4 = op.join(temp_path, 'movie.mp4')
1794 return_code = render_movie(
1795 fn_temp_path,
1796 output_path=fn_mp4,
1797 framerate=framerate)
1799 if render_as_gif and return_code:
1800 render_gif(fn=fn_mp4, output_path=op.join(
1801 fn_path, '%s_%s_gif.gif' % (prefix, plot_type)),
1802 loops=gif_loops)
1804 elif return_code:
1805 shutil.move(fn_mp4, op.join(
1806 fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
1808 else:
1809 logger.error('ffmpeg failed. Exit')
1811 if store_images:
1812 fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t))
1813 for t in times]
1815 for ft, f in zip(fns_temp, fns):
1816 shutil.move(ft, f)
1819__all__ = [
1820 'make_colormap',
1821 'clear_temp',
1822 'xy_to_latlon',
1823 'xy_to_lw',
1824 'SourceError',
1825 'RuptureMap',
1826 'RuptureView',
1827 'rupture_movie',
1828 'render_movie']