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