1#!/usr/bin/python3
2from subprocess import check_call, CalledProcessError
4import os
5import os.path as op
6import re
7import logging
8import tempfile
9import shutil
10import inspect
12import numpy as num
13from scipy.interpolate import RegularGridInterpolator as scrgi
15from matplotlib import cm, pyplot as plt, patheffects
16from matplotlib.ticker import FuncFormatter
18from pyrocko import orthodrome as pod
19from pyrocko.guts import Object
20from pyrocko.plot import mpl_init, mpl_papersize, mpl_color, AutoScaler, gmtpy
21from pyrocko.plot.automap import Map, NoTopo
22from pyrocko.gf import PseudoDynamicRupture
23from pyrocko.gf.seismosizer import map_anchor
24from pyrocko.dataset.topo.tile import Tile
26logger = logging.getLogger(__name__)
28gmtpy.check_have_gmt()
29gmt = gmtpy.GMT()
31km = 1e3
33d2m = 111180.
34m2d = 1. / d2m
36cm2inch = gmtpy.cm/gmtpy.inch
38d2r = num.pi / 180.
39r2d = 1. / d2r
42def km_formatter(v, p):
43 return '%g' % (v / km)
46def get_kwargs(cls):
47 sig = inspect.signature(cls.__init__)
48 kwargs = {}
50 if cls.T.cls == RuptureMap:
51 for p in Map.T.properties:
52 kwargs[p.name] = p._default
54 for par in sig.parameters.values():
55 if par.default is inspect._empty:
56 continue
57 kwargs[par.name] = par.default
59 return kwargs
62def _save_grid(lats, lons, data, filename):
63 '''
64 Save lat-lon gridded data as gmt .grd file
66 :param lats: Grid latitude coordinates [degree]
67 :type lats: iterable
68 :param lons: Grid longitude coordinates [degree]
69 :type lons: iterable
70 :param data: Grid data of any kind
71 :type data: :py:class:`numpy.ndarray`, ``(n_lons, n_lats)``
72 :param filename: Filename of the written grid file
73 :type filename: str
74 '''
76 gmtpy.savegrd(lons, lats, data, filename=filename, naming='lonlat')
79def _mplcmap_to_gmtcpt_code(mplcmap, steps=256):
80 '''
81 Get gmt readable R/G/A code from a given matplotlib colormap
83 :param mplcmap: Name of the demanded matplotlib colormap
84 :type mplcmap: str
86 :returns: Series of comma seperate R/G/B values for gmtpy usage
87 :rtype: str
88 '''
90 cmap = cm.get_cmap(mplcmap)
92 rgbas = [cmap(i) for i in num.linspace(0, 255, steps).astype(num.int64)]
94 return ','.join(['%d/%d/%d' % (
95 c[0] * 255, c[1] * 255, c[2] * 255) for c in rgbas])
98def make_colormap(
99 vmin,
100 vmax,
101 C=None,
102 cmap=None,
103 space=False):
104 '''
105 Create gmt-readable colormap cpt file called my_<cmap>.cpt
107 :type vmin: Minimum value covered by the colormap
108 :param vmin: float
109 :type vmax: Maximum value covered by the colormap
110 :param vmax: float
111 :type C: comma seperated R/G/B values for cmap definition.
112 :param C: optional, str
113 :type cmap: Name of the colormap. Colormap is stored as "my_<cmap>.cpt".
114 If name is equivalent to a matplotlib colormap, R/G/B strings are
115 extracted from this colormap.
116 :param cmap: optional, str
117 :type space: If True, the range of the colormap is broadened below vmin and
118 above vmax.
119 :param space: optional, bool
120 '''
122 scaler = AutoScaler(mode='min-max')
123 scale = scaler.make_scale((vmin, vmax))
125 incr = scale[2]
126 margin = 0.
128 if vmin == vmax:
129 space = True
131 if space:
132 margin = incr
134 msg = ('Please give either a valid color code or a'
135 ' valid matplotlib colormap name.')
137 if C is None and cmap is None:
138 raise ValueError(msg)
140 if C is None:
141 try:
142 C = _mplcmap_to_gmtcpt_code(cmap)
143 except ValueError:
144 raise ValueError(msg)
146 if cmap is None:
147 logger.warning(
148 'No colormap name given. Uses temporary filename instead')
149 cmap = 'temp_cmap'
151 return gmt.makecpt(
152 C=C,
153 D='o',
154 T='%g/%g/%g' % (
155 vmin - margin, vmax + margin, incr),
156 Z=True,
157 out_filename='my_%s.cpt' % cmap,
158 suppress_defaults=True)
161def clear_temp(gridfiles=[], cpts=[]):
162 '''
163 Clear all temporary needed grid and colormap cpt files
165 :param gridfiles: List of all "...grd" files, which shall be deleted
166 :type gridfiles: optional, list
167 :param cpts: List of all cmaps, whose corresponding "my_<cmap>.cpt" file
168 shall be deleted
169 :type cpts: optional, list
170 '''
172 for fil in gridfiles:
173 try:
174 os.remove(fil)
175 except OSError:
176 continue
177 for fil in cpts:
178 try:
179 os.remove('my_%s.cpt' % fil)
180 except OSError:
181 continue
184def xy_to_latlon(source, x, y):
185 '''
186 Convert x and y relative coordinates on extended ruptures into latlon
188 :param source: Extended source class, on which the given point is located
189 :type source: :py:class:`pyrocko.gf.seismosizer.RectangularSource` or
190 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
191 :param x: Relative point coordinate along strike (range: -1:1)
192 :type x: float or :py:class:`numpy.ndarray`
193 :param y: Relative downdip point coordinate (range: -1:1)
194 :type y: float or :py:class:`numpy.ndarray`
196 :returns: Latitude and Longitude [degrees] of the given point
197 :rtype: tuple of float
198 '''
200 s = source
201 ax, ay = map_anchor[s.anchor]
203 length, width = (x - ax) / 2. * s.length, (y - ay) / 2. * s.width
204 strike, dip = s.strike * d2r, s.dip * d2r
206 northing = num.cos(strike)*length - num.cos(dip) * num.sin(strike)*width
207 easting = num.sin(strike)*length + num.cos(dip) * num.cos(strike)*width
209 northing += source.north_shift
210 easting += source.east_shift
212 return pod.ne_to_latlon(s.lat, s.lon, northing, easting)
215def xy_to_lw(source, x, y):
216 '''
217 Convert x and y relative coordinates on extended ruptures into length width
219 :param source: Extended source class, on which the given points are located
220 :type source: :py:class:`pyrocko.gf.seismosizer.RectangularSource` or
221 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
222 :param x: Relative point coordinates along strike (range: -1:1)
223 :type x: float or :py:class:`numpy.ndarray`
224 :param y: Relative downdip point coordinates (range: -1:1)
225 :type y: float or :py:class:`numpy.ndarray`
227 :returns: length and downdip width [m] of the given points from the anchor
228 :rtype: tuple of float
229 '''
231 length, width = source.length, source.width
233 ax, ay = map_anchor[source.anchor]
235 lengths = (x - ax) / 2. * length
236 widths = (y - ay) / 2. * width
238 return lengths, widths
241cbar_anchor = {
242 'center': 'MC',
243 'center_left': 'ML',
244 'center_right': 'MR',
245 'top': 'TC',
246 'top_left': 'TL',
247 'top_right': 'TR',
248 'bottom': 'BC',
249 'bottom_left': 'BL',
250 'bottom_right': 'BR'}
253cbar_helper = {
254 'traction': {
255 'unit': 'MPa',
256 'factor': 1e-6},
257 'tx': {
258 'unit': 'MPa',
259 'factor': 1e-6},
260 'ty': {
261 'unit': 'MPa',
262 'factor': 1e-6},
263 'tz': {
264 'unit': 'MPa',
265 'factor': 1e-6},
266 'time': {
267 'unit': 's',
268 'factor': 1.},
269 'strike': {
270 'unit': '°',
271 'factor': 1.},
272 'dip': {
273 'unit': '°',
274 'factor': 1.},
275 'vr': {
276 'unit': 'km/s',
277 'factor': 1e-3},
278 'length': {
279 'unit': 'km',
280 'factor': 1e-3},
281 'width': {
282 'unit': 'km',
283 'factor': 1e-3}
284}
287fonttype = 'Helvetica'
289c2disl = dict([('x', 0), ('y', 1), ('z', 2)])
292def _make_gmt_conf(fontcolor, size):
293 '''
294 Update different gmt parameters depending on figure size and fontcolor
296 :param fontcolor: GMT readable colorcode / colorstring for the font
297 :type fontcolor: str
298 :param size: Tuple of the figure size (width, height) [centimetre]
299 :type size: tuple of float
301 :returns: estimate best fitting fontsize,
302 Dictionary of different gmt configuration parameter
303 :rtype: float, dict
304 '''
306 color = fontcolor
307 fontsize = num.max(size)
309 font = '%gp,%s' % (fontsize, fonttype)
311 pen_size = fontsize / 16.
312 tick_size = num.min(size) / 200.
314 return fontsize, dict(
315 MAP_FRAME_PEN='%s' % color,
316 MAP_TICK_PEN_PRIMARY='%gp,%s' % (pen_size, color),
317 MAP_TICK_PEN_SECONDARY='%gp,%s' % (pen_size, color),
318 MAP_TICK_LENGTH_PRIMARY='%gc' % tick_size,
319 MAP_TICK_LENGTH_SECONDARY='%gc' % (tick_size * 3),
320 FONT_ANNOT_PRIMARY='%s-Bold,%s' % (font, color),
321 FONT_LABEL='%s-Bold,%s' % (font, color),
322 FONT_TITLE='%s-Bold,%s' % (font, color),
323 PS_CHAR_ENCODING='ISOLatin1+',
324 MAP_FRAME_TYPE='fancy',
325 FORMAT_GEO_MAP='D',
326 PS_PAGE_ORIENTATION='portrait',
327 MAP_GRID_PEN_PRIMARY='thinnest,%s' % color,
328 MAP_ANNOT_OBLIQUE='6')
331class SourceError(Exception):
332 pass
335class RuptureMap(Map):
336 ''' Map plotting of attributes and results of the
337 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
338 '''
340 def __init__(
341 self,
342 source=None,
343 fontcolor='darkslategrey',
344 width=20.,
345 height=14.,
346 margins=None,
347 color_wet=(216, 242, 254),
348 color_dry=(238, 236, 230),
349 topo_cpt_wet='light_sea_uniform',
350 topo_cpt_dry='light_land_uniform',
351 show_cities=False,
352 *args,
353 **kwargs):
355 size = (width, height)
356 fontsize, gmt_config = _make_gmt_conf(fontcolor, size)
358 if margins is None:
359 margins = [
360 fontsize * 0.15, num.min(size) / 200.,
361 num.min(size) / 200., fontsize * 0.05]
363 Map.__init__(self, *args, margins=margins, width=width, height=height,
364 gmt_config=gmt_config,
365 topo_cpt_dry=topo_cpt_dry, topo_cpt_wet=topo_cpt_wet,
366 color_wet=color_wet, color_dry=color_dry,
367 **kwargs)
369 if show_cities:
370 self.draw_cities()
372 self._source = source
373 self._fontcolor = fontcolor
374 self._fontsize = fontsize
375 self.axes_layout = 'WeSn'
377 @property
378 def size(self):
379 '''
380 Figure size [cm]
381 '''
383 return (self.width, self.height)
385 @property
386 def font(self):
387 '''
388 Font style (size and type)
389 '''
391 return '%sp,%s' % (self._fontsize, fonttype)
393 @property
394 def source(self):
395 '''
396 PseudoDynamicRupture whose attributes are plotted.
398 Note, that source.patches attribute needs to be calculated
399 :type source: :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
400 '''
402 if self._source is None:
403 raise SourceError('No source given. Please define it!')
405 if not isinstance(self._source, PseudoDynamicRupture):
406 raise SourceError('This function is only capable for a source of'
407 ' type: %s' % type(PseudoDynamicRupture()))
409 if not self._source.patches:
410 raise TypeError('No source patches are defined. Please run'
411 '"source.discretize_patches()" on your source')
413 return self._source
415 @source.setter
416 def source(self, source):
417 self._source = source
419 def _get_topotile(self):
420 if self._dems is None:
421 self._setup()
423 try:
424 t, _ = self._get_topo_tile('land')
425 except NoTopo:
426 wesn = self.wesn
428 nx = int(num.ceil(
429 self.width * cm2inch * self.topo_resolution_max))
430 ny = int(num.ceil(
431 self.height * cm2inch * self.topo_resolution_max))
433 data = num.zeros((nx, ny))
435 t = Tile(wesn[0], wesn[2],
436 (wesn[1] - wesn[0]) / nx, (wesn[3] - wesn[2]) / ny,
437 data)
439 return t
441 def _patches_to_lw(self):
442 '''
443 Generate regular rect. length-width grid based on the patch distance
445 Prerequesite is a regular grid of patches (constant lengths and widths)
446 Both coordinates are given relative to the source anchor point [in m]
447 The grid is extended from the patch centres to the edges of the source
449 :returns: lengths along strike, widths downdip
450 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
451 '''
453 source = self.source
454 patches = source.patches
456 patch_l, patch_w = patches[0].length, patches[0].width
458 patch_lengths = num.concatenate((
459 num.array([0.]),
460 num.array([il*patch_l+patch_l/2. for il in range(source.nx)]),
461 num.array([patch_l * source.nx])))
463 patch_widths = num.concatenate((
464 num.array([0.]),
465 num.array([iw*patch_w+patch_w/2. for iw in range(source.ny)]),
466 num.array([patch_w * source.ny])))
468 ax, ay = map_anchor[source.anchor]
470 patch_lengths -= source.length * (ax + 1.) / 2.
471 patch_widths -= source.width * (ay + 1.) / 2.
473 return patch_lengths, patch_widths
475 def _xy_to_lw(self, x, y):
476 '''
477 Generate regular rect. length-width grid based on the xy coordinates
479 Prerequesite is a regular grid with constant dx and dy. x and y are
480 relative coordinates on the rupture plane (range -1:1) along strike and
481 downdip.
482 Length and width are obtained relative to the source anchor point
483 [in m].
485 :returns: lengths along strike [m], widths downdip [m]
486 :rtype: :py:class:`numpy.ndarray`, :py:class:`numpy.ndarray`
487 '''
489 x, y = num.unique(x), num.unique(y)
490 dx, dy = x[1] - x[0], y[1] - y[0]
492 if any(num.abs(num.diff(x) - dx) >= 1e-6):
493 raise ValueError('Regular grid with constant spacing needed.'
494 ' Please check the x coordinates.')
496 if any(num.abs(num.diff(y) - dy) >= 1e-6):
497 raise ValueError('Regular grid with constant spacing needed.'
498 ' Please check the y coordinates.')
500 return xy_to_lw(self.source, x, y)
502 def _tile_to_lw(self, ref_lat, ref_lon,
503 north_shift=0., east_shift=0., strike=0.):
505 '''
506 Coordinate transformation from the topo tile grid into length-width
508 The topotile latlon grid is rotated into the length width grid. The
509 width is defined here as its horizontal component. The resulting grid
510 is used for interpolation of grid data.
512 :param ref_lat: Reference latitude, from which length-width relative
513 coordinatesgrid are calculated
514 :type ref_lat: float
515 :param ref_lon: Reference longitude, from which length-width relative
516 coordinatesgrid are calculated
517 :type ref_lon: float
518 :param north_shift: North shift of the reference point [m]
519 :type north_shift: optional, float
520 :param east_shift: East shift of the reference point [m]
521 :type east_shift: optional, float
522 :param strike: striking of the length axis compared to the North axis
523 [degree]
524 :type strike: optional, float
526 :returns: topotile grid nodes as array of length-width coordinates
527 :rtype: :py:class:`numpy.ndarray`, ``(n_nodes, 2)``
528 '''
530 t = self._get_topotile()
531 grid_lats = t.y()
532 grid_lons = t.x()
534 meshgrid_lons, meshgrid_lats = num.meshgrid(grid_lons, grid_lats)
536 grid_northing, grid_easting = pod.latlon_to_ne_numpy(
537 ref_lat, ref_lon, meshgrid_lats.flatten(), meshgrid_lons.flatten())
539 grid_northing -= north_shift
540 grid_easting -= east_shift
542 strike *= d2r
543 sin, cos = num.sin(strike), num.cos(strike)
544 points_out = num.zeros((grid_northing.shape[0], 2))
545 points_out[:, 0] = -sin * grid_northing + cos * grid_easting
546 points_out[:, 1] = cos * grid_northing + sin * grid_easting
548 return points_out
550 def _prep_patch_grid_data(self, data):
551 '''
552 Extend patch data from patch centres to the outer source edges
554 Patch data is always defined in the centre of the patches. For
555 interpolation the data is extended here to the edges of the rupture
556 plane.
558 :param data: Patch wise data
559 :type data: :py:class:`numpy.ndarray`
561 :returns: Extended data array
562 :rtype: :py:class:`numpy.ndarray`
563 '''
565 shape = (self.source.nx + 2, self.source.ny + 2)
566 data = data.reshape(self.source.nx, self.source.ny).copy()
568 data_new = num.zeros(shape)
569 data_new[1:-1, 1:-1] = data
570 data_new[1:-1, 0] = data[:, 0]
571 data_new[1:-1, -1] = data[:, -1]
572 data_new[0, 1:-1] = data[0, :]
573 data_new[-1, 1:-1] = data[-1, :]
575 for i, j in zip([-1, -1, 0, 0], [-1, 0, -1, 0]):
576 data_new[i, j] = data[i, j]
578 return data_new
580 def _regular_data_to_grid(self, lengths, widths, data, filename):
581 '''
582 Interpolate regular data onto topotile grid
584 Regular gridded data is interpolated onto the latlon grid of the
585 topotile. It is then stored as a gmt-readable .grd-file.
587 :param lengths: Grid coordinates along strike relative to anchor [m]
588 :type lengths: :py:class:`numpy.ndarray`
589 :param widths: Grid coordinates downdip relative to anchor [m]
590 :type widths: :py:class:`numpy.ndarray`
591 :param data: Data grid array
592 :type data: :py:class:`numpy.ndarray`
593 :param filename: Filename, where grid is stored
594 :type filename: str
595 '''
597 source = self.source
599 interpolator = scrgi(
600 (widths * num.cos(d2r * source.dip), lengths),
601 data.T,
602 bounds_error=False,
603 method='nearest')
605 points_out = self._tile_to_lw(
606 ref_lat=source.lat,
607 ref_lon=source.lon,
608 north_shift=source.north_shift,
609 east_shift=source.east_shift,
610 strike=source.strike)
612 t = self._get_topotile()
613 t.data = num.zeros_like(t.data, dtype=num.float)
614 t.data[:] = num.nan
616 t.data = interpolator(points_out).reshape(t.data.shape)
618 _save_grid(t.y(), t.x(), t.data, filename=filename)
620 def patch_data_to_grid(self, data, *args, **kwargs):
621 '''
622 Generate a grid file based on regular patch wise data.
624 :param data: Patchwise data grid array
625 :type data: :py:class:`numpy.ndarray`
626 '''
628 lengths, widths = self._patches_to_lw()
630 data_new = self._prep_patch_grid_data(data)
632 self._regular_data_to_grid(lengths, widths, data_new, *args, **kwargs)
634 def xy_data_to_grid(self, x, y, data, *args, **kwargs):
635 '''
636 Generate a grid file based on regular gridded data using xy coordinates
638 Convert a grid based on relative fault coordinates (range -1:1) along
639 strike (x) and downdip (y) into a .grd file.
641 :param x: Relative point coordinate along strike (range: -1:1)
642 :type x: float or :py:class:`numpy.ndarray`
643 :param y: Relative downdip point coordinate (range: -1:1)
644 :type y: float or :py:class:`numpy.ndarray`
645 :param data: Patchwise data grid array
646 :type data: :py:class:`numpy.ndarray`
647 '''
649 lengths, widths = self._xy_to_lw(x, y)
651 self._regular_data_to_grid(
652 lengths, widths, data.reshape((lengths.shape[0], widths.shape[0])),
653 *args, **kwargs)
655 def draw_image(self, gridfile, cmap, cbar=True, **kwargs):
656 '''
657 Draw grid data as image and include, if whished, a colorbar
659 :param gridfile: File of the grid which shall be plotted
660 :type gridfile: str
661 :param cmap: Name of the colormap, which shall be used. A .cpt-file
662 "my_<cmap>.cpt" must exist
663 :type cmap: str
664 :param cbar: If True, a colorbar corresponding to the grid data is
665 added. Keywordarguments are parsed to it.
666 :type cbar: optional, bool
667 '''
669 self.gmt.grdimage(
670 gridfile,
671 C='my_%s.cpt' % cmap,
672 E='200',
673 Q=True,
674 n='+t0.0',
675 *self.jxyr)
677 if cbar:
678 self.draw_colorbar(cmap=cmap, **kwargs)
680 def draw_contour(
681 self,
682 gridfile,
683 contour_int,
684 anot_int,
685 angle=None,
686 unit='',
687 color='',
688 style='',
689 **kwargs):
691 '''
692 Draw grid data as contour lines
694 :param gridfile: File of the grid which shall be plotted
695 :type gridfile: str
696 :param contour_int: Interval of contour lines in units of the gridfile
697 :type contour_int: float
698 :param anot_int: Interval of labelled contour lines in units of the
699 gridfile. Must be a integer multiple of contour_int
700 :type anot_int: float
701 :param angle: Rotation angle of the labels [degree]
702 :type angle: optional, float
703 :param unit: Name of the unit in the grid file. It will be displayed
704 behind the label on labelled contour lines
705 :type unit: optional, str
706 :param color: GMT readable color code/str of the contour lines
707 :type color: optional, str
708 :param style: Line style of the contour lines. If not given, solid
709 lines are plotted
710 :type style: optional, str
711 '''
713 pen_size = self._fontsize / 40.
715 if not color:
716 color = self._fontcolor
718 a_string = '%g+f%s,%s+r%gc+u%s' % (
719 anot_int, self.font, color, pen_size*4, unit)
720 if angle:
721 a_string += '+a%g' % angle
722 c_string = '%g' % contour_int
724 if kwargs:
725 kwargs['A'], kwargs['C'] = a_string, c_string
726 else:
727 kwargs = dict(A=a_string, C=c_string)
729 if style:
730 style = ',' + style
732 args = ['-Wc%gp,%s%s+s' % (pen_size, color, style)]
734 self.gmt.grdcontour(
735 gridfile,
736 S='10',
737 W='a%gp,%s%s+s' % (pen_size*4, color, style),
738 *self.jxyr + args,
739 **kwargs)
741 def draw_colorbar(self, cmap, label='', anchor='top_right', **kwargs):
742 '''
743 Draw a colorbar based on a existing colormap
745 :param cmap: Name of the colormap, which shall be used. A .cpt-file
746 "my_<cmap>.cpt" must exist
747 :type cmap: str
748 :param label: Title of the colorbar
749 :type label: optional, str
750 :param anchor: Placement of the colorbar. Combine 'top', 'center' and
751 'bottom' with 'left', None for middle and 'right'
752 :type anchor: optional, str
753 '''
755 if not kwargs:
756 kwargs = {}
758 if label:
759 kwargs['B'] = 'af+l%s' % label
761 kwargs['C'] = 'my_%s.cpt' % cmap
762 a_str = cbar_anchor[anchor]
764 w = self.width / 3.
765 h = w / 10.
767 lgap = rgap = w / 10.
768 bgap, tgap = h, h / 10.
770 dx, dy = 2.5 * lgap, 2. * tgap
772 if 'bottom' in anchor:
773 dy += 4 * h
775 self.gmt.psscale(
776 D='j%s+w%gc/%gc+h+o%gc/%gc' % (a_str, w, h, dx, dy),
777 F='+g238/236/230+c%g/%g/%g/%g' % (lgap, rgap, bgap, tgap),
778 *self.jxyr,
779 **kwargs)
781 def draw_vector(self, x_gridfile, y_gridfile, vcolor='', **kwargs):
782 ''' Draw vectors based on two grid files
784 Two grid files for vector lengths in x and y need to be given. The
785 function calls gmt.grdvector. All arguments defined for this function
786 in gmt can be passed as keyword arguments. Different standard settings
787 are applied if not defined differently.
789 :param x_gridfile: File of the grid defining vector lengths in x
790 :type x_gridfile: str
791 :param y_gridfile: File of the grid defining vector lengths in y
792 :type y_gridfile: str
793 :param vcolor: Vector face color as defined in "G" option
794 :type vcolor: str
795 '''
797 kwargs['S'] = kwargs.get('S', 'il1.')
798 kwargs['I'] = kwargs.get('I', 'x20')
799 kwargs['W'] = kwargs.get('W', '0.3p,%s' % 'black')
800 kwargs['Q'] = kwargs.get('Q', '4c+e+n5c+h1')
802 self.gmt.grdvector(
803 x_gridfile, y_gridfile,
804 G='%s' % 'lightgrey' if not vcolor else vcolor,
805 *self.jxyr,
806 **kwargs)
808 def draw_dynamic_data(self, data, **kwargs):
809 '''
810 Draw an image of any data gridded on the patches e.g dislocation
812 :param data: Patchwise data grid array
813 :type data: :py:class:`numpy.ndarray`
814 '''
816 plot_data = data
818 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
820 clim = kwargs.pop('clim', (plot_data.min(), plot_data.max()))
822 cpt = []
823 if not op.exists('my_%s.cpt' % kwargs['cmap']):
824 make_colormap(clim[0], clim[1],
825 cmap=kwargs['cmap'], space=False)
826 cpt = [kwargs['cmap']]
828 tmp_grd_file = 'tmpdata.grd'
829 self.patch_data_to_grid(plot_data, tmp_grd_file)
830 self.draw_image(tmp_grd_file, **kwargs)
832 clear_temp(gridfiles=[tmp_grd_file], cpts=cpt)
834 def draw_patch_parameter(self, attribute, **kwargs):
835 '''Draw an image of a chosen patch attribute e.g traction
837 :param attribute: Patch attribute, which is plotted. All patch
838 attributes can be taken (see doc of
839 :py:class:`pyrocko.modelling.okada.OkadaSource`) and also
840 ``traction``, ``tx``, ``ty`` or ``tz`` to display the
841 length or the single components of the traction vector.
842 :type attribute: str
843 '''
845 a = attribute
846 source = self.source
848 if a == 'traction':
849 data = num.linalg.norm(source.get_tractions(), axis=1)
850 elif a == 'tx':
851 data = source.get_tractions()[:, 0]
852 elif a == 'ty':
853 data = source.get_tractions()[:, 1]
854 elif a == 'tz':
855 data = source.get_tractions()[:, 2]
856 else:
857 data = source.get_patch_attribute(attribute)
859 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
860 data *= factor
862 kwargs['label'] = kwargs.get(
863 'label',
864 '%s [%s]' % (a, cbar_helper[a]['unit']))
866 self.draw_dynamic_data(data, **kwargs)
868 def draw_time_contour(self, store, clevel=[], **kwargs):
869 '''Draw high contour lines of the rupture front propgation time
871 :param store: Greens function store, which is used for time calculation
872 :type store: :py:class:`pyrocko.gf.store.Store`
873 :param clevel: List of times, for which contour lines are drawn,
874 optional
875 :type clevel: list of float
876 '''
878 _, _, _, _, points_xy = self.source._discretize_points(store, cs='xyz')
879 _, _, times, _, _, _ = self.source.get_vr_time_interpolators(store)
881 scaler = AutoScaler(mode='0-max', approx_ticks=8)
882 scale = scaler.make_scale([num.min(times), num.max(times)])
884 if clevel:
885 if len(clevel) > 1:
886 kwargs['anot_int'] = num.min(num.diff(clevel))
887 else:
888 kwargs['anot_int'] = clevel[0]
890 kwargs['contour_int'] = kwargs['anot_int']
891 kwargs['L'] = '0/%g' % num.max(clevel)
893 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
894 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
895 kwargs['unit'] = kwargs.get('unit', cbar_helper['time']['unit'])
896 kwargs['L'] = kwargs.get('L', '0/%g' % (num.max(times) + 1.))
897 kwargs['G'] = kwargs.get('G', 'n2/3c')
899 tmp_grd_file = 'tmpdata.grd'
900 self.xy_data_to_grid(points_xy[:, 0], points_xy[:, 1],
901 times, tmp_grd_file)
902 self.draw_contour(tmp_grd_file, **kwargs)
904 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
906 def draw_points(self, lats, lons, symbol='point', size=None, **kwargs):
907 '''Draw points at given locations
909 :param lats: Point latitude coordinates [degree]
910 :type lats: iterable
911 :param lons: Point longitude coordinates [degree]
912 :type lons: iterable
913 :param symbol: Define symbol of the points
914 (``'star', 'circle', 'point', 'triangle'``) - default is ``point``
915 :type symbol: optional, str
916 :param size: Size of the points in points
917 :type size: optional, float
918 '''
920 sym_to_gmt = dict(
921 star='a',
922 circle='c',
923 point='p',
924 triangle='t')
926 lats = num.atleast_1d(lats)
927 lons = num.atleast_1d(lons)
929 if lats.shape[0] != lons.shape[0]:
930 raise IndexError('lats and lons do not have the same shape!')
932 if size is None:
933 size = self._fontsize / 3.
935 kwargs['S'] = kwargs.get('S', sym_to_gmt[symbol] + '%gp' % size)
936 kwargs['G'] = kwargs.get('G', gmtpy.color('scarletred2'))
937 kwargs['W'] = kwargs.get('W', '2p,%s' % self._fontcolor)
939 self.gmt.psxy(
940 in_columns=[lons, lats],
941 *self.jxyr,
942 **kwargs)
944 def draw_nucleation_point(self, **kwargs):
945 ''' Plot the nucleation point onto the map '''
947 nlat, nlon = xy_to_latlon(
948 self.source, self.source.nucleation_x, self.source.nucleation_y)
950 self.draw_points(nlat, nlon, **kwargs)
952 def draw_dislocation(self, time=None, component='', **kwargs):
953 ''' Draw dislocation onto map at any time
955 For a given time (if ``time`` is ``None``, ``tmax`` is used)
956 and given component the patchwise dislocation is plotted onto the map.
958 :param time: time after origin, for which dislocation is computed. If
959 ``None``, ``tmax`` is taken.
960 :type time: optional, float
961 :param component: Dislocation component, which shall be plotted: ``x``
962 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
963 length of the dislocation vector is plotted
964 '''
966 disl = self.source.get_slip(time=time)
968 if component:
969 data = disl[:, c2disl[component]]
970 else:
971 data = num.linalg.norm(disl, axis=1)
973 kwargs['label'] = kwargs.get(
974 'label', 'u%s [m]' % (component))
976 self.draw_dynamic_data(data, **kwargs)
978 def draw_dislocation_contour(
979 self, time=None, component=None, clevel=[], **kwargs):
980 ''' Draw dislocation contour onto map at any time
982 For a given time (if ``time`` is ``None``, ``tmax`` is used)
983 and given component the patchwise dislocation is plotted as contour
984 onto the map.
986 :param time: time after origin, for which dislocation is computed. If
987 ``None``, ``tmax`` is taken.
988 :type time: optional, float
989 :param component: Dislocation component, which shall be plotted: ``x``
990 along strike, ``y`` along updip, ``z`` normal``. If ``None``,
991 the length of the dislocation vector is plotted
992 :param clevel: List of times, for which contour lines are drawn
993 :type clevel: optional, list of float
994 '''
996 disl = self.source.get_slip(time=time)
998 if component:
999 data = disl[:, c2disl[component]]
1000 else:
1001 data = num.linalg.norm(disl, axis=1)
1003 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1004 scale = scaler.make_scale([num.min(data), num.max(data)])
1006 if clevel:
1007 if len(clevel) > 1:
1008 kwargs['anot_int'] = num.min(num.diff(clevel))
1009 else:
1010 kwargs['anot_int'] = clevel[0]
1012 kwargs['contour_int'] = kwargs['anot_int']
1013 kwargs['L'] = '%g/%g' % (
1014 num.min(clevel) - kwargs['contour_int'],
1015 num.max(clevel) + kwargs['contour_int'])
1016 else:
1017 kwargs['anot_int'] = kwargs.get('anot_int', scale[2] * 2.)
1018 kwargs['contour_int'] = kwargs.get('contour_int', scale[2])
1019 kwargs['L'] = kwargs.get('L', '%g/%g' % (
1020 num.min(data) - 1., num.max(data) + 1.))
1022 kwargs['unit'] = kwargs.get('unit', ' m')
1023 kwargs['G'] = kwargs.get('G', 'n2/3c')
1025 tmp_grd_file = 'tmpdata.grd'
1026 self.patch_data_to_grid(data, tmp_grd_file)
1027 self.draw_contour(tmp_grd_file, **kwargs)
1029 clear_temp(gridfiles=[tmp_grd_file], cpts=[])
1031 def draw_dislocation_vector(self, time=None, **kwargs):
1032 '''Draw vector arrows onto map indicating direction of dislocation
1034 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1035 and given component the dislocation is plotted as vectors onto the map.
1037 :param time: time after origin [s], for which dislocation is computed.
1038 If ``None``, ``tmax`` is used.
1039 :type time: optional, float
1040 '''
1042 disl = self.source.get_slip(time=time)
1044 p_strike = self.source.get_patch_attribute('strike') * d2r
1045 p_dip = self.source.get_patch_attribute('dip') * d2r
1047 disl_dh = num.cos(p_dip) * disl[:, 1]
1048 disl_n = num.cos(p_strike) * disl[:, 0] + num.sin(p_strike) * disl_dh
1049 disl_e = num.sin(p_strike) * disl[:, 0] - num.cos(p_strike) * disl_dh
1051 tmp_grd_files = ['tmpdata_%s.grd' % c for c in ('n', 'e')]
1053 self.patch_data_to_grid(disl_n, tmp_grd_files[0])
1054 self.patch_data_to_grid(disl_e, tmp_grd_files[1])
1056 self.draw_vector(
1057 tmp_grd_files[1], tmp_grd_files[0],
1058 **kwargs)
1060 clear_temp(gridfiles=tmp_grd_files, cpts=[])
1062 def draw_top_edge(self, **kwargs):
1063 '''Indicate rupture top edge on map
1064 '''
1066 outline = self.source.outline(cs='latlondepth')
1067 top_edge = outline[:2, :]
1069 kwargs = kwargs or {}
1070 kwargs['W'] = kwargs.get(
1071 'W', '%gp,%s' % (self._fontsize / 10., gmtpy.color('orange3')))
1073 self.gmt.psxy(
1074 in_columns=[top_edge[:, 1], top_edge[:, 0]],
1075 *self.jxyr,
1076 **kwargs)
1079class RuptureView(Object):
1080 ''' Plot of attributes and results of the
1081 :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`.
1082 '''
1084 _patches_to_lw = RuptureMap._patches_to_lw
1086 def __init__(self, source=None, figsize=None, fontsize=None):
1087 self._source = source
1088 self._axes = None
1090 self.figsize = figsize or mpl_papersize('halfletter', 'landscape')
1091 self.fontsize = fontsize or 10
1092 mpl_init(fontsize=self.fontsize)
1094 self._fig = None
1095 self._axes = None
1096 self._is_1d = False
1098 @property
1099 def source(self):
1100 ''' PseudoDynamicRupture whose attributes are plotted.
1102 Note, that source.patches attribute needs to be calculated for
1103 :type source: :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
1104 '''
1106 if self._source is None:
1107 raise SourceError('No source given. Please define it!')
1109 if not isinstance(self._source, PseudoDynamicRupture):
1110 raise SourceError('This function is only capable for a source of'
1111 ' type: %s' % type(PseudoDynamicRupture()))
1113 if not self._source.patches:
1114 raise TypeError('No source patches are defined. Please run'
1115 '\"discretize_patches\" on your source')
1117 return self._source
1119 @source.setter
1120 def source(self, source):
1121 self._source = source
1123 def _setup(
1124 self,
1125 title='',
1126 xlabel='',
1127 ylabel='',
1128 aspect=1.,
1129 spatial_plot=True,
1130 **kwargs):
1132 if self._fig is not None and self._axes is not None:
1133 return
1135 self._fig = plt.figure(figsize=self.figsize)
1137 self._axes = self._fig.add_subplot(1, 1, 1, aspect=aspect)
1138 ax = self._axes
1140 if ax is not None:
1141 ax.set_title(title)
1142 ax.grid(alpha=.3)
1143 ax.set_xlabel(xlabel)
1144 ax.set_ylabel(ylabel)
1146 if spatial_plot:
1147 ax.xaxis.set_major_formatter(FuncFormatter(km_formatter))
1148 ax.yaxis.set_major_formatter(FuncFormatter(km_formatter))
1150 def _clear_all(self):
1151 plt.cla()
1152 plt.clf()
1153 plt.close()
1155 self._fig, self._axes = None, None
1157 def _draw_scatter(self, x, y, *args, **kwargs):
1158 default_kwargs = dict(
1159 linewidth=0,
1160 marker='o',
1161 markerfacecolor=mpl_color('skyblue2'),
1162 markersize=6.,
1163 markeredgecolor=mpl_color('skyblue3'))
1164 default_kwargs.update(kwargs)
1166 if self._axes is not None:
1167 self._axes.plot(x, y, *args, **default_kwargs)
1169 def _draw_image(self, length, width, data, *args, **kwargs):
1170 if self._axes is not None:
1171 if 'extent' not in kwargs:
1172 kwargs['extent'] = [
1173 num.min(length), num.max(length),
1174 num.max(width), num.min(width)]
1176 im = self._axes.imshow(
1177 data,
1178 *args,
1179 interpolation='none',
1180 vmin=kwargs.get('clim', [None])[0],
1181 vmax=kwargs.get('clim', [None, None])[1],
1182 **kwargs)
1184 del kwargs['extent']
1185 if 'aspect' in kwargs:
1186 del kwargs['aspect']
1187 if 'clim' in kwargs:
1188 del kwargs['clim']
1189 if 'cmap' in kwargs:
1190 del kwargs['cmap']
1192 plt.colorbar(
1193 im, *args, shrink=0.9, pad=0.03, aspect=15., **kwargs)
1195 def _draw_contour(self, x, y, data, clevel=None, unit='', *args, **kwargs):
1196 setup_kwargs = dict(
1197 xlabel=kwargs.pop('xlabel', 'along strike [km]'),
1198 ylabel=kwargs.pop('ylabel', 'down dip [km]'),
1199 title=kwargs.pop('title', ''),
1200 aspect=kwargs.pop('aspect', 1))
1202 self._setup(**setup_kwargs)
1204 if self._axes is not None:
1205 if clevel is None:
1206 scaler = AutoScaler(mode='min-max')
1207 scale = scaler.make_scale([data.min(), data.max()])
1209 clevel = num.arange(scale[0], scale[1] + scale[2], scale[2])
1211 if not isinstance(clevel, num.ndarray):
1212 clevel = num.array([clevel])
1214 clevel = clevel[clevel < data.max()]
1216 cont = self._axes.contour(
1217 x, y, data, clevel, *args,
1218 linewidths=1.5, **kwargs)
1220 plt.setp(cont.collections, path_effects=[
1221 patheffects.withStroke(linewidth=2.0, foreground="beige"),
1222 patheffects.Normal()])
1224 clabels = self._axes.clabel(
1225 cont, clevel[::2], *args,
1226 inline=1, fmt='%g' + '%s' % unit,
1227 inline_spacing=15, rightside_up=True, use_clabeltext=True,
1228 **kwargs)
1230 plt.setp(clabels, path_effects=[
1231 patheffects.withStroke(linewidth=1.25, foreground="beige"),
1232 patheffects.Normal()])
1234 def draw_points(self, length, width, *args, **kwargs):
1235 ''' Draw a point onto the figure.
1237 Args and kwargs can be defined according to
1238 :py:func:`matplotlib.pyplot.scatter`.
1240 :param length: Point(s) coordinate on the rupture plane along strike
1241 relative to the anchor point [m]
1242 :type length: float, :py:class:`numpy.ndarray`
1243 :param width: Point(s) coordinate on the rupture plane along downdip
1244 relative to the anchor point [m]
1245 :type width: float, :py:class:`numpy.ndarray`
1246 '''
1248 if self._axes is not None:
1249 kwargs['color'] = kwargs.get('color', mpl_color('scarletred2'))
1250 kwargs['s'] = kwargs.get('s', 100.)
1252 self._axes.scatter(length, width, *args, **kwargs)
1254 def draw_dynamic_data(self, data, **kwargs):
1255 ''' Draw an image of any data gridded on the patches e.g dislocation
1257 :param data: Patchwise data grid array
1258 :type data: :py:class:`numpy.ndarray`
1259 '''
1261 plot_data = data
1263 anchor_x, anchor_y = map_anchor[self.source.anchor]
1265 length, width = xy_to_lw(
1266 self.source, num.array([-1., 1.]), num.array([-1., 1.]))
1268 data = plot_data.reshape(self.source.nx, self.source.ny).T
1270 kwargs['cmap'] = kwargs.get('cmap', 'afmhot_r')
1272 setup_kwargs = dict(
1273 xlabel='along strike [km]',
1274 ylabel='down dip [km]',
1275 title='',
1276 aspect=1)
1277 setup_kwargs.update(kwargs)
1279 kwargs = {k: v for k, v in kwargs.items() if k not in
1280 ('xlabel', 'ylabel', 'title')}
1281 self._setup(**setup_kwargs)
1282 self._draw_image(length=length, width=width, data=data, **kwargs)
1284 def draw_patch_parameter(self, attribute, **kwargs):
1285 ''' Draw an image of a chosen patch attribute e.g traction
1287 :param attribute: Patch attribute, which is plotted. All patch
1288 attributes can be taken (see doc of
1289 :py:class:`pyrocko.modelling.okada.OkadaSource`) and also
1290 ``'traction', 'tx', 'ty', 'tz'`` to display the length
1291 or the single components of the traction vector.
1292 :type attribute: str
1293 '''
1295 a = attribute
1296 source = self.source
1298 if a == 'traction':
1299 data = num.linalg.norm(source.get_tractions(), axis=1)
1300 elif a == 'tx':
1301 data = source.get_tractions()[:, 0]
1302 elif a == 'ty':
1303 data = source.get_tractions()[:, 1]
1304 elif a == 'tz':
1305 data = source.get_tractions()[:, 2]
1306 else:
1307 data = source.get_patch_attribute(attribute)
1309 factor = 1. if 'label' in kwargs else cbar_helper[a]['factor']
1310 data *= factor
1312 kwargs['label'] = kwargs.get(
1313 'label',
1314 '%s [%s]' % (a, cbar_helper[a]['unit']))
1316 return self.draw_dynamic_data(data=data, **kwargs)
1318 def draw_time_contour(self, store, clevel=[], **kwargs):
1319 ''' Draw high resolution contours of the rupture front propgation time
1321 :param store: Greens function store, which is used for time calculation
1322 :type store: :py:class:`pyrocko.gf.store.Store`
1323 :param clevel: levels of the contour lines. If no levels are given,
1324 they are automatically computed based on tmin and tmax
1325 :type clevel: optional, list
1326 '''
1327 source = self.source
1328 default_kwargs = dict(
1329 colors='#474747ff'
1330 )
1331 default_kwargs.update(kwargs)
1333 _, _, _, _, points_xy = source._discretize_points(store, cs='xyz')
1334 _, _, _, _, time_interpolator, _ = source.get_vr_time_interpolators(
1335 store)
1337 times = time_interpolator.values
1339 scaler = AutoScaler(mode='min-max', approx_ticks=8)
1340 scale = scaler.make_scale([times.min(), times.max()])
1342 if len(clevel) == 0:
1343 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1345 points_x = time_interpolator.grid[0]
1346 points_y = time_interpolator.grid[1]
1348 self._draw_contour(points_x, points_y, data=times.T,
1349 clevel=clevel, unit='s', **default_kwargs)
1351 def draw_nucleation_point(self, **kwargs):
1352 ''' Draw the nucleation point onto the map '''
1354 nuc_x, nuc_y = self.source.nucleation_x, self.source.nucleation_y
1355 length, width = xy_to_lw(self.source, nuc_x, nuc_y)
1356 self.draw_points(length, width, marker='o', **kwargs)
1358 def draw_dislocation(self, time=None, component='', **kwargs):
1359 ''' Draw dislocation onto map at any time
1361 For a given time (if ``time`` is ``None``, ``tmax`` is used)
1362 and given component the patchwise dislocation is plotted onto the map.
1364 :param time: time after origin [s], for which dislocation is computed.
1365 If ``None``, ``tmax`` is taken.
1366 :type time: optional, float
1367 :param component: Dislocation component, which shall be plotted: ``x``
1368 along strike, ``y`` along updip, ``z`` normal. If ``None``, the
1369 length of the dislocation vector is plotted
1370 '''
1372 disl = self.source.get_slip(time=time)
1374 if component:
1375 data = disl[:, c2disl[component]]
1376 else:
1377 data = num.linalg.norm(disl, axis=1)
1379 kwargs['label'] = kwargs.get(
1380 'label', 'u%s [m]' % (component))
1382 self.draw_dynamic_data(data, **kwargs)
1384 def draw_dislocation_contour(
1385 self, time=None, component=None, clevel=[], **kwargs):
1386 ''' Draw dislocation contour onto map at any time
1388 For a given time (if time is ``None``, ``tmax`` is used) and given
1389 component the patchwise dislocation is plotted as contour onto the map.
1391 :param time: time after origin, for which dislocation is computed. If
1392 None, tmax is taken.
1393 :type time: optional, float
1394 :param component: Dislocation component, which shall be plotted. ``x``
1395 along strike, ``y`` along updip, ``z`` - normal. If ``None``
1396 is given, the length of the dislocation vector is plotted
1397 :type component: str
1398 '''
1400 disl = self.source.get_slip(time=time)
1402 if component:
1403 data = disl[:, c2disl[component]]
1404 else:
1405 data = num.linalg.norm(disl, axis=1)
1407 data = data.reshape(self.source.ny, self.source.nx, order='F')
1409 scaler = AutoScaler(mode='min-max', approx_ticks=7)
1410 scale = scaler.make_scale([data.min(), data.max()])
1412 if len(clevel) == 0:
1413 clevel = num.arange(scale[0] + scale[2], scale[1], scale[2])
1415 anchor_x, anchor_y = map_anchor[self.source.anchor]
1417 length, width = self._patches_to_lw()
1418 length, width = length[1:-1], width[1:-1]
1420 kwargs['colors'] = kwargs.get('colors', '#474747ff')
1422 self._setup(**kwargs)
1423 self._draw_contour(
1424 length, width, data=data, clevel=clevel, unit='m', **kwargs)
1426 def draw_source_dynamics(
1427 self, variable, store, deltat=None, *args, **kwargs):
1428 ''' Display dynamic source parameter
1430 Fast inspection possibility for the cumulative moment and the source
1431 time function approximation (assuming equal paths between different
1432 patches and observation point - valid for an observation point in the
1433 far field perpendicular to the source strike), so the cumulative moment
1434 rate function.
1436 :param variable: Dynamic parameter, which shall be plotted. Choose
1437 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1438 :type variable: str
1439 :param store: Greens function store, whose store.config.deltat defines
1440 the time increment between two parameter snapshots. If store is not
1441 given, the time increment is defined is taken from deltat.
1442 :type store: :py:class:`pyrocko.gf.store.Store`
1443 :param deltat: Time increment between two parameter snapshots. If not
1444 given, store.config.deltat is used to define deltat
1445 :type deltat: optional, float
1446 '''
1448 v = variable
1450 data, times = self.source.get_moment_rate(store=store, deltat=deltat)
1451 deltat = num.concatenate([(num.diff(times)[0], ), num.diff(times)])
1453 if v in ('moment_rate', 'stf'):
1454 name, unit = 'dM/dt', 'Nm/s'
1455 elif v in ('cumulative_moment', 'moment'):
1456 data = num.cumsum(data * deltat)
1457 name, unit = 'M', 'Nm'
1458 else:
1459 raise ValueError('No dynamic data for given variable %s found' % v)
1461 self._setup(xlabel='time [s]',
1462 ylabel='%s / %.2g %s' % (name, data.max(), unit),
1463 aspect='auto',
1464 spatial_plot=False)
1465 self._draw_scatter(times, data/num.max(data), *args, **kwargs)
1466 self._is_1d = True
1468 def draw_patch_dynamics(
1469 self, variable, nx, ny, store=None, deltat=None, *args, **kwargs):
1470 ''' Display dynamic boundary element / patch parameter
1472 Fast inspection possibility for different dynamic parameter for a
1473 single patch / boundary element. The chosen parameter is plotted for
1474 the chosen patch.
1476 :param variable: Dynamic parameter, which shall be plotted. Choose
1477 between 'moment_rate' ('stf') or 'cumulative_moment' ('moment')
1478 :type variable: str
1479 :param nx: Patch index along strike (range: 0:self.source.nx - 1)
1480 :type nx: int
1481 :param nx: Patch index downdip (range: 0:self.source.ny - 1)
1482 :type nx: int
1483 :param store: Greens function store, whose store.config.deltat defines
1484 the time increment between two parameter snapshots. If store is not
1485 given, the time increment is defined is taken from deltat.
1486 :type store: optional, :py:class:`pyrocko.gf.store.Store`
1487 :param deltat: Time increment between two parameter snapshots. If not
1488 given, store.config.deltat is used to define deltat
1489 :type deltat: optional, float
1490 '''
1492 v = variable
1493 source = self.source
1494 idx = nx * source.ny + nx
1496 m = re.match(r'dislocation_([xyz])', v)
1498 if v in ('moment_rate', 'cumulative_moment', 'moment', 'stf'):
1499 data, times = source.get_moment_rate_patches(
1500 store=store, deltat=deltat)
1501 elif 'dislocation' in v or 'slip_rate' == v:
1502 ddisloc, times = source.get_delta_slip(store=store, deltat=deltat)
1504 if v in ('moment_rate', 'stf'):
1505 data, times = source.get_moment_rate_patches(
1506 store=store, deltat=deltat)
1507 name, unit = 'dM/dt', 'Nm/s'
1508 elif v in ('cumulative_moment', 'moment'):
1509 data, times = source.get_moment_rate_patches(
1510 store=store, deltat=deltat)
1511 data = num.cumsum(data, axis=1)
1512 name, unit = 'M', 'Nm'
1513 elif v == 'slip_rate':
1514 data, times = source.get_slip_rate(store=store, deltat=deltat)
1515 name, unit = 'du/dt', 'm/s'
1516 elif v == 'dislocation':
1517 data, times = source.get_delta_slip(store=store, deltat=deltat)
1518 data = num.linalg.norm(num.cumsum(data, axis=2), axis=1)
1519 name, unit = 'du', 'm'
1520 elif m:
1521 data, times = source.get_delta_slip(store=store, deltat=deltat)
1522 data = num.cumsum(data, axis=2)[:, c2disl[m.group(1)], :]
1523 name, unit = 'du%s' % m.group(1), 'm'
1524 else:
1525 raise ValueError('No dynamic data for given variable %s found' % v)
1527 self._setup(xlabel='time [s]',
1528 ylabel='%s / %.2g %s' % (name, num.max(data), unit),
1529 aspect='auto',
1530 spatial_plot=False)
1531 self._draw_scatter(times, data[idx, :]/num.max(data),
1532 *args, **kwargs)
1533 self._is_1d = True
1535 def finalize(self):
1536 if self._is_1d:
1537 return
1539 length, width = xy_to_lw(
1540 self.source, num.array([-1., 1.]), num.array([1., -1.]))
1542 self._axes.set_xlim(length)
1543 self._axes.set_ylim(width)
1545 def gcf(self):
1546 self.finalize()
1547 return self._fig
1549 def save(self, filename, dpi=None):
1550 ''' Save plot to file
1552 :param filename: filename and path, where the plot is stored at
1553 :type filename: str
1554 :param dpi: Resolution of the output plot [dpi]
1555 :type dpi: int
1556 '''
1557 self.finalize()
1558 try:
1559 self._fig.savefig(fname=filename, dpi=dpi, bbox_inches='tight')
1560 except TypeError:
1561 self._fig.savefig(filename=filename, dpi=dpi, bbox_inches='tight')
1563 self._clear_all()
1565 def show_plot(self):
1566 ''' Show plot '''
1567 self.finalize()
1568 plt.show()
1569 self._clear_all()
1572def render_movie(fn_path, output_path, framerate=20):
1573 ''' Generate a mp4 movie based on given png files using
1574 `ffmpeg <https://ffmpeg.org>`_.
1576 Render a movie based on a set of given .png files in fn_path. All files
1577 must have a filename specified by ``fn_path`` (e.g. giving ``fn_path`` with
1578 ``/temp/f%04.png`` a valid png filename would be ``/temp/f0001.png``). The
1579 files must have a numbering, indicating their order in the movie.
1581 :param fn_path: Path and fileformat specification of the input .png files.
1582 :type fn_path: str
1583 :param output_path: Path and filename of the output ``.mp4`` movie file
1584 :type output_path: str
1585 :param deltat: Time between individual frames (``1 / framerate``) [s]
1586 :type deltat: optional, float
1588 '''
1589 try:
1590 check_call(['ffmpeg', '-loglevel', 'panic'])
1591 except CalledProcessError:
1592 pass
1593 except (TypeError, IOError):
1594 logger.warning(
1595 'Package ffmpeg needed for movie rendering. Please install it '
1596 '(e.g. on linux distr. via sudo apt-get ffmpeg) and retry.')
1597 return False
1599 try:
1600 check_call([
1601 'ffmpeg', '-loglevel', 'info', '-y',
1602 '-framerate', '%g' % framerate,
1603 '-i', fn_path,
1604 '-vcodec', 'libx264',
1605 '-preset', 'medium',
1606 '-tune', 'animation',
1607 '-pix_fmt', 'yuv420p',
1608 '-movflags', '+faststart',
1609 '-filter:v', "crop='floor(in_w/2)*2:floor(in_h/2)*2'",
1610 '-crf', '15',
1611 output_path])
1613 return True
1614 except CalledProcessError as e:
1615 logger.warning(e)
1616 return False
1619def render_gif(fn, output_path, loops=-1):
1620 ''' Generate a gif based on a given mp4 using ffmpeg
1622 Render a gif based on a given .mp4 movie file in ``fn`` path.
1624 :param fn: Path and file name of the input .mp4 file.
1625 :type fn: str
1626 :param output_path: Path and filename of the output animated gif file
1627 :type output_path: str
1628 :param loops: Number of gif repeats (loops). ``-1`` is not repetition,
1629 ``0`` infinite
1630 :type loops: optional, integer
1631 '''
1633 try:
1634 check_call(['ffmpeg', '-loglevel', 'panic'])
1635 except CalledProcessError:
1636 pass
1637 except (TypeError, IOError):
1638 logger.warning(
1639 'Package ffmpeg needed for movie rendering. Please install it '
1640 '(e.g. on linux distr. via sudo apt-get ffmpeg.) and retry.')
1641 return False
1643 try:
1644 check_call([
1645 'ffmpeg', '-hide_banner', '-loglevel', 'panic', '-i',
1646 fn,
1647 '-filter_complex', 'palettegen[v1];[0:v][v1]paletteuse',
1648 '-loop', '%d' % loops,
1649 output_path])
1651 return True
1652 except CalledProcessError as e:
1653 logger.warning(e)
1654 return False
1657def rupture_movie(
1658 source,
1659 store,
1660 variable='dislocation',
1661 draw_time_contours=False,
1662 fn_path='.',
1663 prefix='',
1664 plot_type='map',
1665 deltat=None,
1666 framerate=None,
1667 store_images=False,
1668 render_as_gif=False,
1669 gif_loops=-1,
1670 **kwargs):
1671 ''' Generate a movie based on a given source for dynamic parameter
1673 Create a MPEG-4 movie or gif of one of the following dynamic parameters
1674 (``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1675 (along updip), ``dislocation_z`` (normal), ``slip_rate``, ``moment_rate``).
1676 If desired, the single snap shots can be stored as images as well.
1677 ``kwargs`` have to be given according to the chosen ``plot_type``.
1679 :param source: Pseudo dynamic rupture, for which the movie is produced
1680 :type source: :py:class:`pyrocko.gf.seismosizer.PseudoDynamicRupture`
1681 :param store: Greens function store, which is used for time calculation. If
1682 deltat is not given, it is taken from the store.config.deltat
1683 :type store: :py:class:`pyrocko.gf.store.Store`
1684 :param variable: Dynamic parameter, which shall be plotted. Choose between
1685 ``dislocation``, ``dislocation_x`` (along strike), ``dislocation_y``
1686 (along updip), ``dislocation_z`` (normal), ``slip_rate`` and
1687 ``moment_rate``, default ``dislocation``
1688 :type variable: optional, str
1689 :param draw_time_contours: If True, corresponding isochrones are drawn on
1690 the each plots
1691 :type draw_time_contours: optional, bool
1692 :param fn_path: Absolut or relative path, where movie (and optional images)
1693 are stored
1694 :type fn_path: optional, str
1695 :param prefix: File prefix used for the movie (and optional image) files
1696 :type prefix: optional, str
1697 :param plot_type: Choice of plot type: ``map``, ``view`` (map plot using
1698 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureMap`
1699 or plane view using
1700 :py:class:`~pyrocko.plot.dynamic_rupture.RuptureView`)
1701 :type plot_type: optional, str
1702 :param deltat: Time between parameter snapshots. If not given,
1703 store.config.deltat is used to define deltat
1704 :type deltat: optional, float
1705 :param store_images: Choice to store the single .png parameter snapshots in
1706 fn_path or not.
1707 :type store_images: optional, bool
1708 :param render_as_gif: If ``True``, the movie is converted into a gif. If
1709 ``False``, the movie is returned as mp4
1710 :type render_as_gif: optional, bool
1711 :param gif_loops: If ``render_as_gif`` is ``True``, a gif with
1712 ``gif_loops`` number of loops (repetitions) is returned.
1713 ``-1`` is no repetition, ``0`` infinite.
1714 :type gif_loops: optional, integer
1715 '''
1717 v = variable
1718 assert plot_type in ('map', 'view')
1720 if not source.patches:
1721 source.discretize_patches(store, interpolation='multilinear')
1723 if source.coef_mat is None:
1724 source.calc_coef_mat()
1726 prefix = prefix or v
1727 deltat = deltat or store.config.deltat
1728 framerate = max(framerate or int(1./deltat), 1)
1730 if v == 'moment_rate':
1731 data, times = source.get_moment_rate_patches(deltat=deltat)
1732 name, unit = 'dM/dt', 'Nm/s'
1733 elif 'dislocation' in v or 'slip_rate' == v:
1734 ddisloc, times = source.get_delta_slip(deltat=deltat)
1735 else:
1736 raise ValueError('No dynamic data for given variable %s found' % v)
1738 deltat = times[1] - times[0]
1740 m = re.match(r'dislocation_([xyz])', v)
1741 if m:
1742 data = num.cumsum(ddisloc, axis=1)[:, :, c2disl[m.group(1)]]
1743 name, unit = 'du%s' % m.group(1), 'm'
1744 elif v == 'dislocation':
1745 data = num.linalg.norm(num.cumsum(ddisloc, axis=1), axis=2)
1746 name, unit = 'du', 'm'
1747 elif v == 'slip_rate':
1748 data = num.linalg.norm(ddisloc, axis=2) / deltat
1749 name, unit = 'du/dt', 'm/s'
1751 if plot_type == 'map':
1752 plt_base = RuptureMap
1753 elif plot_type == 'view':
1754 plt_base = RuptureView
1755 else:
1756 raise AttributeError('invalid type: %s' % plot_type)
1758 attrs_base = get_kwargs(plt_base)
1759 kwargs_base = dict([k for k in kwargs.items() if k[0] in attrs_base])
1760 kwargs_plt = dict([k for k in kwargs.items() if k[0] not in kwargs_base])
1762 if 'clim' in kwargs_plt:
1763 data = num.clip(data, kwargs_plt['clim'][0], kwargs_plt['clim'][1])
1764 else:
1765 kwargs_plt['clim'] = [num.min(data), num.max(data)]
1767 if 'label' not in kwargs_plt:
1768 vmax = num.max(num.abs(kwargs_plt['clim']))
1769 data /= vmax
1771 kwargs_plt['label'] = '%s / %.2g %s' % (name, vmax, unit)
1772 kwargs_plt['clim'] = [i / vmax for i in kwargs_plt['clim']]
1774 with tempfile.TemporaryDirectory(suffix='pyrocko') as temp_path:
1775 fns_temp = [op.join(temp_path, 'f%09d.png' % (it + 1))
1776 for it, _ in enumerate(times)]
1777 fn_temp_path = op.join(temp_path, 'f%09d.png')
1779 for it, (t, ft) in enumerate(zip(times, fns_temp)):
1780 plot_data = data[:, it]
1782 plt = plt_base(source=source, **kwargs_base)
1783 plt.draw_dynamic_data(plot_data, **kwargs_plt)
1784 plt.draw_nucleation_point()
1786 if draw_time_contours:
1787 plt.draw_time_contour(store, clevel=[t])
1789 plt.save(ft)
1791 fn_mp4 = op.join(temp_path, 'movie.mp4')
1792 return_code = render_movie(
1793 fn_temp_path,
1794 output_path=fn_mp4,
1795 framerate=framerate)
1797 if render_as_gif and return_code:
1798 render_gif(fn=fn_mp4, output_path=op.join(
1799 fn_path, '%s_%s_gif.gif' % (prefix, plot_type)),
1800 loops=gif_loops)
1802 elif return_code:
1803 shutil.move(fn_mp4, op.join(
1804 fn_path, '%s_%s_movie.mp4' % (prefix, plot_type)))
1806 else:
1807 logger.error('ffmpeg failed. Exit')
1809 if store_images:
1810 fns = [op.join(fn_path, '%s_%s_%g.png' % (prefix, plot_type, t))
1811 for t in times]
1813 for ft, f in zip(fns_temp, fns):
1814 shutil.move(ft, f)
1817__all__ = [
1818 'make_colormap',
1819 'clear_temp',
1820 'xy_to_latlon',
1821 'xy_to_lw',
1822 'SourceError',
1823 'RuptureMap',
1824 'RuptureView',
1825 'rupture_movie',
1826 'render_movie']