1# -*- coding: utf-8 -*-
2# -----------------------------------------------
3# Filename: mopad_wrapper.py
4#  Purpose: Wrapper for mopad
5#   Author: Tobias Megies, Moritz Beyreuther
6#    Email: megies@geophysik.uni-muenchen.de
7#
8# Copyright (C) 2008-2012 ObsPy Development Team
9# -----------------------------------------------
10"""
11ObsPy wrapper to the *Moment tensor Plotting and Decomposition tool* (MoPaD)
12written by Lars Krieger and Sebastian Heimann.
13
14.. seealso:: [Krieger2012]_
15
16.. warning:: The MoPaD wrapper does not yet provide the full functionality of
17    MoPaD. Please consider using the command line script ``obspy-mopad`` for
18    now if you need the full power of MoPaD.
19
20:copyright:
21    The ObsPy Development Team (devs@obspy.org)
22:license:
23    GNU Lesser General Public License, Version 3
24    (https://www.gnu.org/copyleft/lesser.html)
25"""
26from __future__ import (absolute_import, division, print_function,
27                        unicode_literals)
28from future.builtins import *  # NOQA @UnusedWildImport
29
30import numpy as np
31import matplotlib.collections as mpl_collections
32from matplotlib import patches, transforms
33
34from obspy.imaging.beachball import xy2patch
35from obspy.imaging.scripts.mopad import BeachBall as mopad_BeachBall
36from obspy.imaging.scripts.mopad import MomentTensor as mopad_MomentTensor
37from obspy.imaging.scripts.mopad import epsilon
38
39
40# seems the base system we (gmt) are using is called "USE" in mopad
41KWARG_MAP = {
42    'size': ['plot_size', 'plot_aux_plot_size'],
43    'linewidth': ['plot_nodalline_width', 'plot_outerline_width'],
44    'facecolor': ['plot_tension_colour'],
45    'edgecolor': ['plot_outerline_colour'],
46    'bgcolor': [],
47    'alpha': ['plot_total_alpha'],
48    'width': [],
49    'outfile': ['plot_outfile'],
50    'format': ['plot_outfile_format'],
51    'nofill': ['plot_only_lines']
52}
53
54
55def _normalize_focmec(fm):
56    """
57    Improve stability of plots by normalizing the moment tensors. The scale
58    does not matter for the beachballs.
59    """
60    # Only normalize 6 component tensors.
61    if len(fm) != 6:
62        return fm
63    fm = np.array(fm, dtype=np.float64)
64    fm /= np.linalg.norm(fm)
65    return fm
66
67
68def beach(fm, linewidth=2, facecolor='b', bgcolor='w', edgecolor='k',
69          alpha=1.0, xy=(0, 0), width=200, size=100, nofill=False,
70          zorder=100, mopad_basis='USE', axes=None):
71    """
72    Return a beach ball as a collection which can be connected to an
73    current matplotlib axes instance (ax.add_collection). Based on MoPaD.
74
75    S1, D1, and R1, the strike, dip and rake of one of the focal planes, can
76    be vectors of multiple focal mechanisms.
77
78    :param fm: Focal mechanism that is either number of mechanisms (NM) by 3
79        (strike, dip, and rake) or NM x 6 (M11, M22, M33, M12, M13, M23 - the
80        six independent components of the moment tensor, where the coordinate
81        system is 1,2,3 = Up,South,East which equals r,theta,phi -
82        Harvard/Global CMT convention). The relation to Aki and Richards
83        x,y,z equals North,East,Down convention is as follows: Mrr=Mzz,
84        Mtt=Mxx, Mpp=Myy, Mrt=Mxz, Mrp=-Myz, Mtp=-Mxy.
85        The strike is of the first plane, clockwise relative to north.
86        The dip is of the first plane, defined clockwise and perpendicular to
87        strike, relative to horizontal such that 0 is horizontal and 90 is
88        vertical. The rake is of the first focal plane solution. 90 moves the
89        hanging wall up-dip (thrust), 0 moves it in the strike direction
90        (left-lateral), -90 moves it down-dip (normal), and 180 moves it
91        opposite to strike (right-lateral).
92    :param facecolor: Color to use for quadrants of tension; can be a string,
93        e.g. ``'r'``, ``'b'`` or three component color vector, [R G B].
94        Defaults to ``'b'`` (blue).
95    :param bgcolor: The background color. Defaults to ``'w'`` (white).
96    :param edgecolor: Color of the edges. Defaults to ``'k'`` (black).
97    :param alpha: The alpha level of the beach ball. Defaults to ``1.0``
98        (opaque).
99    :param xy: Origin position of the beach ball as tuple. Defaults to
100        ``(0, 0)``.
101    :type width: int
102    :param width: Symbol size of beach ball. Defaults to ``200``.
103    :param size: Controls the number of interpolation points for the
104        curves. Minimum is automatically set to ``100``.
105    :param nofill: Do not fill the beach ball, but only plot the planes.
106    :param zorder: Set zorder. Artists with lower zorder values are drawn
107        first.
108    :param mopad_basis: The basis system. Defaults to ``'USE'``. See the
109        `Supported Basis Systems`_ section below for a full list of supported
110        systems.
111    :type axes: :class:`matplotlib.axes.Axes`
112    :param axes: Used to make beach balls circular on non-scaled axes. Also
113        maintains the aspect ratio when resizing the figure. Will not add
114        the returned collection to the axes instance.
115
116    .. rubric:: _`Supported Basis Systems`
117
118    ========= =================== =============================================
119    Short     Basis vectors       Usage
120    ========= =================== =============================================
121    ``'NED'`` North, East, Down   Jost and Herrmann 1989
122    ``'USE'`` Up, South, East     Global CMT Catalog, Larson et al. 2010
123    ``'XYZ'`` East, North, Up     General formulation, Jost and Herrmann 1989
124    ``'RT'``  Radial, Transverse, psmeca (GMT), Wessel and Smith 1999
125              Tangential
126    ``'NWU'`` North, West, Up     Stein and Wysession 2003
127    ========= =================== =============================================
128    """
129    fm = _normalize_focmec(fm)
130    # initialize beachball
131    mt = mopad_MomentTensor(fm, system=mopad_basis)
132    bb = mopad_BeachBall(mt, npoints=size)
133    bb._setup_BB(unit_circle=False)
134
135    # extract the coordinates and colors of the lines
136    radius = width / 2.0
137    neg_nodalline = bb._nodalline_negative_final_US
138    pos_nodalline = bb._nodalline_positive_final_US
139    tension_colour = facecolor
140    pressure_colour = bgcolor
141
142    if nofill:
143        tension_colour = 'none'
144        pressure_colour = 'none'
145
146    # based on mopads _setup_plot_US() function
147    # collect patches for the selection
148    coll = [None, None, None]
149    coll[0] = patches.Circle(xy, radius=radius)
150    coll[1] = xy2patch(neg_nodalline[0, :], neg_nodalline[1, :], radius, xy)
151    coll[2] = xy2patch(pos_nodalline[0, :], pos_nodalline[1, :], radius, xy)
152
153    # set the color of the three parts
154    fc = [None, None, None]
155    if bb._plot_clr_order > 0:
156        fc[0] = pressure_colour
157        fc[1] = tension_colour
158        fc[2] = tension_colour
159        if bb._plot_curve_in_curve != 0:
160            fc[0] = tension_colour
161            if bb._plot_curve_in_curve < 1:
162                fc[1] = pressure_colour
163                fc[2] = tension_colour
164            else:
165                coll = [coll[i] for i in (0, 2, 1)]
166                fc[1] = pressure_colour
167                fc[2] = tension_colour
168    else:
169        fc[0] = tension_colour
170        fc[1] = pressure_colour
171        fc[2] = pressure_colour
172        if bb._plot_curve_in_curve != 0:
173            fc[0] = pressure_colour
174            if bb._plot_curve_in_curve < 1:
175                fc[1] = tension_colour
176                fc[2] = pressure_colour
177            else:
178                coll = [coll[i] for i in (0, 2, 1)]
179                fc[1] = tension_colour
180                fc[2] = pressure_colour
181
182    if bb._pure_isotropic:
183        if abs(np.trace(bb._M)) > epsilon:
184            # use the circle as the most upper layer
185            coll = [coll[0]]
186            if bb._plot_clr_order < 0:
187                fc = [tension_colour]
188            else:
189                fc = [pressure_colour]
190
191    # transform the patches to a path collection and set
192    # the appropriate attributes
193    collection = mpl_collections.PatchCollection(coll, match_original=False)
194    collection.set_facecolors(fc)
195    # Use the given axes to maintain the aspect ratio of beachballs on figure
196    # resize.
197    if axes is not None:
198        # This is what holds the aspect ratio (but breaks the positioning)
199        collection.set_transform(transforms.IdentityTransform())
200        # Next is a dirty hack to fix the positioning:
201        # 1. Need to bring the all patches to the origin (0, 0).
202        for p in collection._paths:
203            p.vertices -= xy
204        # 2. Then use the offset property of the collection to position the
205        # patches
206        collection.set_offsets(xy)
207        collection._transOffset = axes.transData
208    collection.set_edgecolors(edgecolor)
209    collection.set_alpha(alpha)
210    collection.set_linewidth(linewidth)
211    collection.set_zorder(zorder)
212    return collection
213
214
215def beachball(fm, linewidth=2, facecolor='b', bgcolor='w', edgecolor='k',
216              alpha=1.0, xy=(0, 0), width=200, size=100, nofill=False,
217              zorder=100, mopad_basis='USE', outfile=None, format=None,
218              fig=None):
219    """
220    Draws a beach ball diagram of an earthquake focal mechanism. Based on
221    MoPaD.
222
223    S1, D1, and R1, the strike, dip and rake of one of the focal planes, can
224    be vectors of multiple focal mechanisms.
225
226    :param fm: Focal mechanism that is either number of mechanisms (NM) by 3
227        (strike, dip, and rake) or NM x 6 (M11, M22, M33, M12, M13, M23 - the
228        six independent components of the moment tensor, where the coordinate
229        system is 1,2,3 = Up,South,East which equals r,theta,phi). The strike
230        is of the first plane, clockwise relative to north.
231        The dip is of the first plane, defined clockwise and perpendicular to
232        strike, relative to horizontal such that 0 is horizontal and 90 is
233        vertical. The rake is of the first focal plane solution. 90 moves the
234        hanging wall up-dip (thrust), 0 moves it in the strike direction
235        (left-lateral), -90 moves it down-dip (normal), and 180 moves it
236        opposite to strike (right-lateral).
237    :param facecolor: Color to use for quadrants of tension; can be a string,
238        e.g. ``'r'``, ``'b'`` or three component color vector, [R G B].
239        Defaults to ``'b'`` (blue).
240    :param bgcolor: The background color. Defaults to ``'w'`` (white).
241    :param edgecolor: Color of the edges. Defaults to ``'k'`` (black).
242    :param alpha: The alpha level of the beach ball. Defaults to ``1.0``
243        (opaque).
244    :param xy: Origin position of the beach ball as tuple. Defaults to
245        ``(0, 0)``.
246    :type width: int
247    :param width: Symbol size of beach ball. Defaults to ``200``.
248    :param size: Controls the number of interpolation points for the
249        curves. Minimum is automatically set to ``100``.
250    :param nofill: Do not fill the beach ball, but only plot the planes.
251    :param zorder: Set zorder. Artists with lower zorder values are drawn
252        first.
253    :param mopad_basis: The basis system. Defaults to ``'USE'``. See the
254        `Supported Basis Systems`_ section below for a full list of supported
255        systems.
256    :param outfile: Output file string. Also used to automatically
257        determine the output format. Supported file formats depend on your
258        matplotlib backend. Most backends support png, pdf, ps, eps and
259        svg. Defaults to ``None``.
260    :param format: Format of the graph picture. If no format is given the
261        outfile parameter will be used to try to automatically determine
262        the output format. If no format is found it defaults to png output.
263        If no outfile is specified but a format is, than a binary
264        imagestring will be returned.
265        Defaults to ``None``.
266    :param fig: Give an existing figure instance to plot into. New Figure if
267        set to ``None``.
268
269    .. rubric:: _`Supported Basis Systems`
270
271    ========= =================== =============================================
272    Short     Basis vectors       Usage
273    ========= =================== =============================================
274    ``'NED'`` North, East, Down   Jost and Herrmann 1989
275    ``'USE'`` Up, South, East     Global CMT Catalog, Larson et al. 2010
276    ``'XYZ'`` East, North, Up     General formulation, Jost and Herrmann 1989
277    ``'RT'``  Radial, Transverse, psmeca (GMT), Wessel and Smith 1999
278              Tangential
279    ``'NWU'`` North, West, Up     Stein and Wysession 2003
280    ========= =================== =============================================
281
282    .. rubric:: Examples
283
284    (1) Using basis system ``'NED'``.
285
286        >>> from obspy.imaging.mopad_wrapper import beachball
287        >>> mt = [1, 2, 3, -4, -5, -10]
288        >>> beachball(mt, mopad_basis='NED') #doctest: +SKIP
289
290        .. plot::
291
292            from obspy.imaging.mopad_wrapper import beachball
293            mt = [1, 2, 3, -4, -5, -10]
294            beachball(mt, mopad_basis='NED')
295    """
296    fm = _normalize_focmec(fm)
297    mopad_kwargs = {}
298    loc = locals()
299    # map to kwargs used in mopad
300    for key in KWARG_MAP:
301        value = loc[key]
302        for mopad_key in KWARG_MAP[key]:
303            mopad_kwargs[mopad_key] = value
304    # convert from points to size in cm
305    for key in ['plot_aux_plot_size', 'plot_size']:
306        # 100.0 is matplotlib's default DPI for savefig
307        mopad_kwargs[key] = mopad_kwargs[key] / 100.0 * 2.54
308    # use nofill kwarg
309
310    mt = mopad_MomentTensor(fm, system=mopad_basis)
311    bb = mopad_BeachBall(mt, npoints=size)
312
313    # show plot in a window
314    if outfile is None:
315        bb.ploBB(mopad_kwargs)
316    # save plot to file
317    else:
318        # no format specified, parse it from outfile name
319        if mopad_kwargs['plot_outfile_format'] is None:
320            mopad_kwargs['plot_outfile_format'] = \
321                mopad_kwargs['plot_outfile'].split(".")[-1]
322        else:
323            # append file format if not already at end of outfile
324            if not mopad_kwargs['plot_outfile'].endswith(
325               mopad_kwargs['plot_outfile_format']):
326                mopad_kwargs['plot_outfile'] += "." + \
327                    mopad_kwargs['plot_outfile_format']
328        bb.save_BB(mopad_kwargs)
329