1# -*- coding: utf-8 -*- 2# ------------------------------------------------------------------- 3# Filename: beachball.py 4# Purpose: Draws a beach ball diagram of an earthquake focal mechanism. 5# Author: Robert Barsch 6# Email: barsch@egu.eu 7# 8# Copyright (C) 2008-2012 Robert Barsch 9# --------------------------------------------------------------------- 10 11""" 12Draws a beachball diagram of an earthquake focal mechanism 13 14Most source code provided here are adopted from 15 161. MatLab script `bb.m`_ written by Andy Michael, Chen Ji and Oliver Boyd. 172. ps_meca program from the `Generic Mapping Tools (GMT)`_. 18 19:copyright: 20 The ObsPy Development Team (devs@obspy.org) 21:license: 22 GNU Lesser General Public License, Version 3 23 (https://www.gnu.org/copyleft/lesser.html) 24 25.. _`Generic Mapping Tools (GMT)`: https://gmt.soest.hawaii.edu 26.. _`bb.m`: http://www.ceri.memphis.edu/people/olboyd/Software/Software.html 27""" 28from __future__ import (absolute_import, division, print_function, 29 unicode_literals) 30from future.builtins import * # NOQA @UnusedWildImport 31 32import io 33import warnings 34 35import numpy as np 36from matplotlib import path as mplpath 37from matplotlib import collections, patches, transforms 38from decorator import decorator 39 40 41D2R = np.pi / 180 42R2D = 180 / np.pi 43EPSILON = 0.00001 44 45 46@decorator 47def mopad_fallback(func, *args, **kwargs): 48 try: 49 result = func(*args, **kwargs) 50 except IndexError: 51 msg = "Encountered an exception while plotting the beachball. " \ 52 "Falling back to the mopad wrapper which is slower but more " \ 53 "stable." 54 warnings.warn(msg) 55 56 # Could be done with the inspect module but this wrapper is only a 57 # single purpose wrapper and thus KISS. 58 arguments = ["fm", "linewidth", "facecolor", "bgcolor", "edgecolor", 59 "alpha", "xy", "width", "size", "nofill", "zorder", 60 "axes"] 61 62 final_kwargs = {} 63 for _i, arg in enumerate(args): 64 final_kwargs[arguments[_i]] = arg 65 66 final_kwargs.update(kwargs) 67 68 from .mopad_wrapper import beach as _mopad_beach 69 70 result = _mopad_beach(**final_kwargs) 71 72 return result 73 74 75@mopad_fallback 76def beach(fm, linewidth=2, facecolor='b', bgcolor='w', edgecolor='k', 77 alpha=1.0, xy=(0, 0), width=200, size=100, nofill=False, 78 zorder=100, axes=None): 79 """ 80 Return a beach ball as a collection which can be connected to an 81 current matplotlib axes instance (ax.add_collection). 82 83 S1, D1, and R1, the strike, dip and rake of one of the focal planes, can 84 be vectors of multiple focal mechanisms. 85 86 :param fm: Focal mechanism that is either number of mechanisms (NM) by 3 87 (strike, dip, and rake) or NM x 6 (M11, M22, M33, M12, M13, M23 - the 88 six independent components of the moment tensor, where the coordinate 89 system is 1,2,3 = Up,South,East which equals r,theta,phi - 90 Harvard/Global CMT convention). The relation to Aki and Richards 91 x,y,z equals North,East,Down convention is as follows: Mrr=Mzz, 92 Mtt=Mxx, Mpp=Myy, Mrt=Mxz, Mrp=-Myz, Mtp=-Mxy. 93 The strike is of the first plane, clockwise relative to north. 94 The dip is of the first plane, defined clockwise and perpendicular to 95 strike, relative to horizontal such that 0 is horizontal and 90 is 96 vertical. The rake is of the first focal plane solution. 90 moves the 97 hanging wall up-dip (thrust), 0 moves it in the strike direction 98 (left-lateral), -90 moves it down-dip (normal), and 180 moves it 99 opposite to strike (right-lateral). 100 :param facecolor: Color to use for quadrants of tension; can be a string, 101 e.g. ``'r'``, ``'b'`` or three component color vector, [R G B]. 102 Defaults to ``'b'`` (blue). 103 :param bgcolor: The background color. Defaults to ``'w'`` (white). 104 :param edgecolor: Color of the edges. Defaults to ``'k'`` (black). 105 :param alpha: The alpha level of the beach ball. Defaults to ``1.0`` 106 (opaque). 107 :param xy: Origin position of the beach ball as tuple. Defaults to 108 ``(0, 0)``. 109 :type width: int or tuple 110 :param width: Symbol size of beach ball, or tuple for elliptically 111 shaped patches. Defaults to size ``200``. 112 :param size: Controls the number of interpolation points for the 113 curves. Minimum is automatically set to ``100``. 114 :param nofill: Do not fill the beach ball, but only plot the planes. 115 :param zorder: Set zorder. Artists with lower zorder values are drawn 116 first. 117 :type axes: :class:`matplotlib.axes.Axes` 118 :param axes: Used to make beach balls circular on non-scaled axes. Also 119 maintains the aspect ratio when resizing the figure. Will not add 120 the returned collection to the axes instance. 121 """ 122 # check if one or two widths are specified (Circle or Ellipse) 123 try: 124 assert(len(width) == 2) 125 except TypeError: 126 width = (width, width) 127 mt = None 128 np1 = None 129 if isinstance(fm, MomentTensor): 130 mt = fm 131 np1 = mt2plane(mt) 132 elif isinstance(fm, NodalPlane): 133 np1 = fm 134 elif len(fm) == 6: 135 mt = MomentTensor(fm[0], fm[1], fm[2], fm[3], fm[4], fm[5], 0) 136 np1 = mt2plane(mt) 137 elif len(fm) == 3: 138 np1 = NodalPlane(fm[0], fm[1], fm[2]) 139 else: 140 raise TypeError("Wrong input value for 'fm'.") 141 142 # Only at least size 100, i.e. 100 points in the matrix are allowed 143 if size < 100: 144 size = 100 145 146 # Return as collection 147 plot_dc_used = True 148 if mt: 149 (t, n, p) = mt2axes(mt.normalized) 150 if np.fabs(n.val) < EPSILON and np.fabs(t.val + p.val) < EPSILON: 151 colors, p = plot_dc(np1, size, xy=xy, width=width) 152 else: 153 colors, p = plot_mt(t, n, p, size, 154 plot_zerotrace=True, xy=xy, width=width) 155 plot_dc_used = False 156 else: 157 colors, p = plot_dc(np1, size=size, xy=xy, width=width) 158 159 col = collections.PatchCollection(p, match_original=False) 160 if nofill: 161 col.set_facecolor('none') 162 else: 163 # Replace color dummies 'b' and 'w' by face and bgcolor 164 fc = [facecolor if c == 'b' else bgcolor for c in colors] 165 col.set_facecolors(fc) 166 167 # Use the given axes to maintain the aspect ratio of beachballs on figure 168 # resize. 169 if axes is not None: 170 # This is what holds the aspect ratio (but breaks the positioning) 171 col.set_transform(transforms.IdentityTransform()) 172 # Next is a dirty hack to fix the positioning: 173 # 1. Need to bring the all patches to the origin (0, 0). 174 for p in col._paths: 175 p.vertices -= xy 176 # 2. Then use the offset property of the collection to position the 177 # patches 178 col.set_offsets(xy) 179 col._transOffset = axes.transData 180 181 col.set_edgecolor(edgecolor) 182 col.set_alpha(alpha) 183 col.set_linewidth(linewidth) 184 col.set_zorder(zorder) 185 186 # warn about color blending bug, see #1464 187 if alpha != 1 and not nofill and not plot_dc_used: 188 msg = ("There is a known bug when plotting semi-transparent patches " 189 "for non-DC sources, which leads to blending of pressure and " 190 "tension color, see issue #1464.") 191 warnings.warn(msg) 192 193 return col 194 195 196def beachball(fm, linewidth=2, facecolor='b', bgcolor='w', edgecolor='k', 197 alpha=1.0, xy=(0, 0), width=200, size=100, nofill=False, 198 zorder=100, outfile=None, format=None, fig=None): 199 """ 200 Draws a beach ball diagram of an earthquake focal mechanism. 201 202 S1, D1, and R1, the strike, dip and rake of one of the focal planes, can 203 be vectors of multiple focal mechanisms. 204 205 :param fm: Focal mechanism that is either number of mechanisms (NM) by 3 206 (strike, dip, and rake) or NM x 6 (M11, M22, M33, M12, M13, M23 - the 207 six independent components of the moment tensor, where the coordinate 208 system is 1,2,3 = Up,South,East which equals r,theta,phi). The strike 209 is of the first plane, clockwise relative to north. 210 The dip is of the first plane, defined clockwise and perpendicular to 211 strike, relative to horizontal such that 0 is horizontal and 90 is 212 vertical. The rake is of the first focal plane solution. 90 moves the 213 hanging wall up-dip (thrust), 0 moves it in the strike direction 214 (left-lateral), -90 moves it down-dip (normal), and 180 moves it 215 opposite to strike (right-lateral). 216 :param facecolor: Color to use for quadrants of tension; can be a string, 217 e.g. ``'r'``, ``'b'`` or three component color vector, [R G B]. 218 Defaults to ``'b'`` (blue). 219 :param bgcolor: The background color. Defaults to ``'w'`` (white). 220 :param edgecolor: Color of the edges. Defaults to ``'k'`` (black). 221 :param alpha: The alpha level of the beach ball. Defaults to ``1.0`` 222 (opaque). 223 :param xy: Origin position of the beach ball as tuple. Defaults to 224 ``(0, 0)``. 225 :type width: int 226 :param width: Symbol size of beach ball. Defaults to ``200``. 227 :param size: Controls the number of interpolation points for the 228 curves. Minimum is automatically set to ``100``. 229 :param nofill: Do not fill the beach ball, but only plot the planes. 230 :param zorder: Set zorder. Artists with lower zorder values are drawn 231 first. 232 :param outfile: Output file string. Also used to automatically 233 determine the output format. Supported file formats depend on your 234 matplotlib backend. Most backends support png, pdf, ps, eps and 235 svg. Defaults to ``None``. 236 :param format: Format of the graph picture. If no format is given the 237 outfile parameter will be used to try to automatically determine 238 the output format. If no format is found it defaults to png output. 239 If no outfile is specified but a format is, than a binary 240 imagestring will be returned. 241 Defaults to ``None``. 242 :param fig: Give an existing figure instance to plot into. New Figure if 243 set to ``None``. 244 """ 245 import matplotlib.pyplot as plt 246 plot_width = width * 0.95 247 248 # plot the figure 249 if not fig: 250 fig = plt.figure(figsize=(3, 3), dpi=100) 251 fig.subplots_adjust(left=0, bottom=0, right=1, top=1) 252 fig.set_figheight(width // 100) 253 fig.set_figwidth(width // 100) 254 ax = fig.add_subplot(111, aspect='equal') 255 256 # hide axes + ticks 257 ax.axison = False 258 259 # plot the collection 260 collection = beach(fm, linewidth=linewidth, facecolor=facecolor, 261 edgecolor=edgecolor, bgcolor=bgcolor, 262 alpha=alpha, nofill=nofill, xy=xy, 263 width=plot_width, size=size, zorder=zorder) 264 ax.add_collection(collection) 265 266 ax.autoscale_view(tight=False, scalex=True, scaley=True) 267 # export 268 if outfile: 269 if format: 270 fig.savefig(outfile, dpi=100, transparent=True, format=format) 271 else: 272 fig.savefig(outfile, dpi=100, transparent=True) 273 elif format and not outfile: 274 imgdata = io.BytesIO() 275 fig.savefig(imgdata, format=format, dpi=100, transparent=True) 276 imgdata.seek(0) 277 return imgdata.read() 278 else: 279 plt.show() 280 return fig 281 282 283def plot_mt(T, N, P, size=200, plot_zerotrace=True, # noqa 284 x0=0, y0=0, xy=(0, 0), width=200): 285 """ 286 Uses a principal axis T, N and P to draw a beach ball plot. 287 288 :param ax: axis object of a matplotlib figure 289 :param T: :class:`~PrincipalAxis` 290 :param N: :class:`~PrincipalAxis` 291 :param P: :class:`~PrincipalAxis` 292 293 Adapted from ps_tensor / utilmeca.c / `Generic Mapping Tools (GMT)`_. 294 295 .. _`Generic Mapping Tools (GMT)`: https://gmt.soest.hawaii.edu 296 """ 297 # check if one or two widths are specified (Circle or Ellipse) 298 try: 299 assert(len(width) == 2) 300 except TypeError: 301 width = (width, width) 302 collect = [] 303 colors = [] 304 res = [value / float(size) for value in width] 305 b = 1 306 big_iso = 0 307 j = 1 308 j2 = 0 309 j3 = 0 310 n = 0 311 azi = np.zeros((3, 2)) 312 x = np.zeros(400) 313 y = np.zeros(400) 314 x2 = np.zeros(400) 315 y2 = np.zeros(400) 316 x3 = np.zeros(400) 317 y3 = np.zeros(400) 318 xp1 = np.zeros(800) 319 yp1 = np.zeros(800) 320 xp2 = np.zeros(400) 321 yp2 = np.zeros(400) 322 323 a = np.zeros(3) 324 p = np.zeros(3) 325 v = np.zeros(3) 326 a[0] = T.strike 327 a[1] = N.strike 328 a[2] = P.strike 329 p[0] = T.dip 330 p[1] = N.dip 331 p[2] = P.dip 332 v[0] = T.val 333 v[1] = N.val 334 v[2] = P.val 335 336 vi = (v[0] + v[1] + v[2]) / 3. 337 for i in range(0, 3): 338 v[i] = v[i] - vi 339 340 radius_size = size * 0.5 341 342 if np.fabs(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) < EPSILON: 343 # pure implosion-explosion 344 if vi > 0.: 345 cir = patches.Ellipse(xy, width=width[0], height=width[1]) 346 collect.append(cir) 347 colors.append('b') 348 if vi < 0.: 349 cir = patches.Ellipse(xy, width=width[0], height=width[1]) 350 collect.append(cir) 351 colors.append('w') 352 return colors, collect 353 354 if np.fabs(v[0]) >= np.fabs(v[2]): 355 d = 0 356 m = 2 357 else: 358 d = 2 359 m = 0 360 361 if (plot_zerotrace): 362 vi = 0. 363 364 f = -v[1] / float(v[d]) 365 iso = vi / float(v[d]) 366 367 # Cliff Frohlich, Seismological Research letters, 368 # Vol 7, Number 1, January-February, 1996 369 # Unless the isotropic parameter lies in the range 370 # between -1 and 1 - f there will be no nodes whatsoever 371 372 if iso < -1: 373 cir = patches.Ellipse(xy, width=width[0], height=width[1]) 374 collect.append(cir) 375 colors.append('w') 376 return colors, collect 377 elif iso > 1 - f: 378 cir = patches.Ellipse(xy, width=width[0], height=width[1]) 379 collect.append(cir) 380 colors.append('b') 381 return colors, collect 382 383 spd = np.sin(p[d] * D2R) 384 cpd = np.cos(p[d] * D2R) 385 spb = np.sin(p[b] * D2R) 386 cpb = np.cos(p[b] * D2R) 387 spm = np.sin(p[m] * D2R) 388 cpm = np.cos(p[m] * D2R) 389 sad = np.sin(a[d] * D2R) 390 cad = np.cos(a[d] * D2R) 391 sab = np.sin(a[b] * D2R) 392 cab = np.cos(a[b] * D2R) 393 sam = np.sin(a[m] * D2R) 394 cam = np.cos(a[m] * D2R) 395 396 for i in range(0, 360): 397 fir = i * D2R 398 s2alphan = (2. + 2. * iso) / \ 399 float(3. + (1. - 2. * f) * np.cos(2. * fir)) 400 if s2alphan > 1.: 401 big_iso += 1 402 else: 403 alphan = np.arcsin(np.sqrt(s2alphan)) 404 sfi = np.sin(fir) 405 cfi = np.cos(fir) 406 san = np.sin(alphan) 407 can = np.cos(alphan) 408 409 xz = can * spd + san * sfi * spb + san * cfi * spm 410 xn = can * cpd * cad + san * sfi * cpb * cab + \ 411 san * cfi * cpm * cam 412 xe = can * cpd * sad + san * sfi * cpb * sab + \ 413 san * cfi * cpm * sam 414 415 if np.fabs(xn) < EPSILON and np.fabs(xe) < EPSILON: 416 takeoff = 0. 417 az = 0. 418 else: 419 az = np.arctan2(xe, xn) 420 if az < 0.: 421 az += np.pi * 2. 422 takeoff = np.arccos(xz / float(np.sqrt(xz * xz + xn * xn + 423 xe * xe))) 424 if takeoff > np.pi / 2.: 425 takeoff = np.pi - takeoff 426 az += np.pi 427 if az > np.pi * 2.: 428 az -= np.pi * 2. 429 r = np.sqrt(2) * np.sin(takeoff / 2.) 430 si = np.sin(az) 431 co = np.cos(az) 432 if i == 0: 433 azi[i][0] = az 434 x[i] = x0 + radius_size * r * si 435 y[i] = y0 + radius_size * r * co 436 azp = az 437 else: 438 if np.fabs(np.fabs(az - azp) - np.pi) < D2R * 10.: 439 azi[n][1] = azp 440 n += 1 441 azi[n][0] = az 442 if np.fabs(np.fabs(az - azp) - np.pi * 2.) < D2R * 2.: 443 if azp < az: 444 azi[n][0] += np.pi * 2. 445 else: 446 azi[n][0] -= np.pi * 2. 447 if n == 0: 448 x[j] = x0 + radius_size * r * si 449 y[j] = y0 + radius_size * r * co 450 j += 1 451 elif n == 1: 452 x2[j2] = x0 + radius_size * r * si 453 y2[j2] = y0 + radius_size * r * co 454 j2 += 1 455 elif n == 2: 456 x3[j3] = x0 + radius_size * r * si 457 y3[j3] = y0 + radius_size * r * co 458 j3 += 1 459 azp = az 460 azi[n][1] = az 461 462 if v[1] < 0.: 463 rgb1 = 'b' 464 rgb2 = 'w' 465 else: 466 rgb1 = 'w' 467 rgb2 = 'b' 468 469 cir = patches.Ellipse(xy, width=width[0], height=width[1]) 470 collect.append(cir) 471 colors.append(rgb2) 472 if n == 0: 473 collect.append(xy2patch(x[0:360], y[0:360], res, xy)) 474 colors.append(rgb1) 475 return colors, collect 476 elif n == 1: 477 for i in range(0, j): 478 xp1[i] = x[i] 479 yp1[i] = y[i] 480 if azi[0][0] - azi[0][1] > np.pi: 481 azi[0][0] -= np.pi * 2. 482 elif azi[0][1] - azi[0][0] > np.pi: 483 azi[0][0] += np.pi * 2. 484 if azi[0][0] < azi[0][1]: 485 az = azi[0][1] - D2R 486 while az > azi[0][0]: 487 si = np.sin(az) 488 co = np.cos(az) 489 xp1[i] = x0 + radius_size * si 490 yp1[i] = y0 + radius_size * co 491 i += 1 492 az -= D2R 493 else: 494 az = azi[0][1] + D2R 495 while az < azi[0][0]: 496 si = np.sin(az) 497 co = np.cos(az) 498 xp1[i] = x0 + radius_size * si 499 yp1[i] = y0 + radius_size * co 500 i += 1 501 az += D2R 502 collect.append(xy2patch(xp1[0:i], yp1[0:i], res, xy)) 503 colors.append(rgb1) 504 for i in range(0, j2): 505 xp2[i] = x2[i] 506 yp2[i] = y2[i] 507 if azi[1][0] - azi[1][1] > np.pi: 508 azi[1][0] -= np.pi * 2. 509 elif azi[1][1] - azi[1][0] > np.pi: 510 azi[1][0] += np.pi * 2. 511 if azi[1][0] < azi[1][1]: 512 az = azi[1][1] - D2R 513 while az > azi[1][0]: 514 si = np.sin(az) 515 co = np.cos(az) 516 xp2[i] = x0 + radius_size * si 517 i += 1 518 yp2[i] = y0 + radius_size * co 519 az -= D2R 520 else: 521 az = azi[1][1] + D2R 522 while az < azi[1][0]: 523 si = np.sin(az) 524 co = np.cos(az) 525 xp2[i] = x0 + radius_size * si 526 i += 1 527 yp2[i] = y0 + radius_size * co 528 az += D2R 529 collect.append(xy2patch(xp2[0:i], yp2[0:i], res, xy)) 530 colors.append(rgb1) 531 return colors, collect 532 elif n == 2: 533 for i in range(0, j3): 534 xp1[i] = x3[i] 535 yp1[i] = y3[i] 536 for ii in range(0, j): 537 xp1[i] = x[ii] 538 i += 1 539 yp1[i] = y[ii] 540 if big_iso: 541 ii = j2 - 1 542 while ii >= 0: 543 xp1[i] = x2[ii] 544 i += 1 545 yp1[i] = y2[ii] 546 ii -= 1 547 collect.append(xy2patch(xp1[0:i], yp1[0:i], res, xy)) 548 colors.append(rgb1) 549 return colors, collect 550 551 if azi[2][0] - azi[0][1] > np.pi: 552 azi[2][0] -= np.pi * 2. 553 elif azi[0][1] - azi[2][0] > np.pi: 554 azi[2][0] += np.pi * 2. 555 if azi[2][0] < azi[0][1]: 556 az = azi[0][1] - D2R 557 while az > azi[2][0]: 558 si = np.sin(az) 559 co = np.cos(az) 560 xp1[i] = x0 + radius_size * si 561 i += 1 562 yp1[i] = y0 + radius_size * co 563 az -= D2R 564 else: 565 az = azi[0][1] + D2R 566 while az < azi[2][0]: 567 si = np.sin(az) 568 co = np.cos(az) 569 xp1[i] = x0 + radius_size * si 570 i += 1 571 yp1[i] = y0 + radius_size * co 572 az += D2R 573 collect.append(xy2patch(xp1[0:i], yp1[0:i], res, xy)) 574 colors.append(rgb1) 575 576 for i in range(0, j2): 577 xp2[i] = x2[i] 578 yp2[i] = y2[i] 579 if azi[1][0] - azi[1][1] > np.pi: 580 azi[1][0] -= np.pi * 2. 581 elif azi[1][1] - azi[1][0] > np.pi: 582 azi[1][0] += np.pi * 2. 583 if azi[1][0] < azi[1][1]: 584 az = azi[1][1] - D2R 585 while az > azi[1][0]: 586 si = np.sin(az) 587 co = np.cos(az) 588 xp2[i] = x0 + radius_size * si 589 i += 1 590 yp2[i] = y0 + radius_size * co 591 az -= D2R 592 else: 593 az = azi[1][1] + D2R 594 while az < azi[1][0]: 595 si = np.sin(az) 596 co = np.cos(az) 597 xp2[i] = x0 + radius_size * si 598 i += 1 599 yp2[i] = y0 + radius_size * co 600 az += D2R 601 collect.append(xy2patch(xp2[0:i], yp2[0:i], res, xy)) 602 colors.append(rgb1) 603 return colors, collect 604 605 606def plot_dc(np1, size=200, xy=(0, 0), width=200): 607 """ 608 Uses one nodal plane of a double couple to draw a beach ball plot. 609 610 :param ax: axis object of a matplotlib figure 611 :param np1: :class:`~NodalPlane` 612 613 Adapted from MATLAB script 614 `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ 615 written by Andy Michael, Chen Ji and Oliver Boyd. 616 """ 617 # check if one or two widths are specified (Circle or Ellipse) 618 try: 619 assert(len(width) == 2) 620 except TypeError: 621 width = (width, width) 622 s_1 = np1.strike 623 d_1 = np1.dip 624 r_1 = np1.rake 625 626 m = 0 627 if r_1 > 180: 628 r_1 -= 180 629 m = 1 630 if r_1 < 0: 631 r_1 += 180 632 m = 1 633 634 # Get azimuth and dip of second plane 635 (s_2, d_2, _r_2) = aux_plane(s_1, d_1, r_1) 636 637 d = size / 2 638 639 if d_1 >= 90: 640 d_1 = 89.9999 641 if d_2 >= 90: 642 d_2 = 89.9999 643 644 # arange checked for numerical stability, np.pi is not multiple of 0.1 645 phi = np.arange(0, np.pi, .01) 646 l1 = np.sqrt( 647 np.power(90 - d_1, 2) / ( 648 np.power(np.sin(phi), 2) + 649 np.power(np.cos(phi), 2) * 650 np.power(90 - d_1, 2) / np.power(90, 2))) 651 l2 = np.sqrt( 652 np.power(90 - d_2, 2) / ( 653 np.power(np.sin(phi), 2) + np.power(np.cos(phi), 2) * 654 np.power(90 - d_2, 2) / np.power(90, 2))) 655 656 collect = [] 657 # plot paths, once for tension areas and once for pressure areas 658 for m_ in ((m + 1) % 2, m): 659 inc = 1 660 (x_1, y_1) = pol2cart(phi + s_1 * D2R, l1) 661 662 if m_ == 1: 663 lo = s_1 - 180 664 hi = s_2 665 if lo > hi: 666 inc = -1 667 th1 = np.arange(s_1 - 180, s_2, inc) 668 (xs_1, ys_1) = pol2cart(th1 * D2R, 90 * np.ones((1, len(th1)))) 669 (x_2, y_2) = pol2cart(phi + s_2 * D2R, l2) 670 th2 = np.arange(s_2 + 180, s_1, -inc) 671 else: 672 hi = s_1 - 180 673 lo = s_2 - 180 674 if lo > hi: 675 inc = -1 676 th1 = np.arange(hi, lo, -inc) 677 (xs_1, ys_1) = pol2cart(th1 * D2R, 90 * np.ones((1, len(th1)))) 678 (x_2, y_2) = pol2cart(phi + s_2 * D2R, l2) 679 x_2 = x_2[::-1] 680 y_2 = y_2[::-1] 681 th2 = np.arange(s_2, s_1, inc) 682 (xs_2, ys_2) = pol2cart(th2 * D2R, 90 * np.ones((1, len(th2)))) 683 x = np.concatenate((x_1, xs_1[0], x_2, xs_2[0])) 684 y = np.concatenate((y_1, ys_1[0], y_2, ys_2[0])) 685 686 x = x * d / 90 687 y = y * d / 90 688 689 # calculate resolution 690 res = [value / float(size) for value in width] 691 692 # construct the patch 693 collect.append(xy2patch(y, x, res, xy)) 694 return ['b', 'w'], collect 695 696 697def xy2patch(x, y, res, xy): 698 # check if one or two resolutions are specified (Circle or Ellipse) 699 try: 700 assert(len(res) == 2) 701 except TypeError: 702 res = (res, res) 703 # transform into the Path coordinate system 704 x = x * res[0] + xy[0] 705 y = y * res[1] + xy[1] 706 verts = list(zip(x.tolist(), y.tolist())) 707 codes = [mplpath.Path.MOVETO] 708 codes.extend([mplpath.Path.LINETO] * (len(x) - 2)) 709 codes.append(mplpath.Path.CLOSEPOLY) 710 path = mplpath.Path(verts, codes) 711 return patches.PathPatch(path) 712 713 714def pol2cart(th, r): 715 """ 716 """ 717 x = r * np.cos(th) 718 y = r * np.sin(th) 719 return (x, y) 720 721 722def strike_dip(n, e, u): 723 """ 724 Finds strike and dip of plane given normal vector having components n, e, 725 and u. 726 727 Adapted from MATLAB script 728 `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ 729 written by Andy Michael, Chen Ji and Oliver Boyd. 730 """ 731 r2d = 180 / np.pi 732 if u < 0: 733 n = -n 734 e = -e 735 u = -u 736 737 strike = np.arctan2(e, n) * r2d 738 strike = strike - 90 739 while strike >= 360: 740 strike = strike - 360 741 while strike < 0: 742 strike = strike + 360 743 x = np.sqrt(np.power(n, 2) + np.power(e, 2)) 744 dip = np.arctan2(x, u) * r2d 745 return (strike, dip) 746 747 748def aux_plane(s1, d1, r1): 749 """ 750 Get Strike and dip of second plane. 751 752 Adapted from MATLAB script 753 `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ 754 written by Andy Michael, Chen Ji and Oliver Boyd. 755 """ 756 r2d = 180 / np.pi 757 758 z = (s1 + 90) / r2d 759 z2 = d1 / r2d 760 z3 = r1 / r2d 761 # slick vector in plane 1 762 sl1 = -np.cos(z3) * np.cos(z) - np.sin(z3) * np.sin(z) * np.cos(z2) 763 sl2 = np.cos(z3) * np.sin(z) - np.sin(z3) * np.cos(z) * np.cos(z2) 764 sl3 = np.sin(z3) * np.sin(z2) 765 (strike, dip) = strike_dip(sl2, sl1, sl3) 766 767 n1 = np.sin(z) * np.sin(z2) # normal vector to plane 1 768 n2 = np.cos(z) * np.sin(z2) 769 h1 = -sl2 # strike vector of plane 2 770 h2 = sl1 771 # note h3=0 always so we leave it out 772 # n3 = np.cos(z2) 773 774 z = h1 * n1 + h2 * n2 775 z = z / np.sqrt(h1 * h1 + h2 * h2) 776 # we might get above 1.0 only due to floating point 777 # precision. Clip for those cases. 778 float64epsilon = 2.2204460492503131e-16 779 if 1.0 < abs(z) < 1.0 + 100 * float64epsilon: 780 z = np.copysign(1.0, z) 781 z = np.arccos(z) 782 rake = 0 783 if sl3 > 0: 784 rake = z * r2d 785 if sl3 <= 0: 786 rake = -z * r2d 787 return (strike, dip, rake) 788 789 790def mt2plane(mt): 791 """ 792 Calculates a nodal plane of a given moment tensor. 793 794 :param mt: :class:`~MomentTensor` 795 :return: :class:`~NodalPlane` 796 797 Adapted from MATLAB script 798 `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ 799 written by Andy Michael, Chen Ji and Oliver Boyd. 800 """ 801 (d, v) = np.linalg.eig(mt.mt) 802 d = np.array([d[1], d[0], d[2]]) 803 v = np.array([[v[1, 1], -v[1, 0], -v[1, 2]], 804 [v[2, 1], -v[2, 0], -v[2, 2]], 805 [-v[0, 1], v[0, 0], v[0, 2]]]) 806 imax = d.argmax() 807 imin = d.argmin() 808 ae = (v[:, imax] + v[:, imin]) / np.sqrt(2.0) 809 an = (v[:, imax] - v[:, imin]) / np.sqrt(2.0) 810 aer = np.sqrt(np.power(ae[0], 2) + np.power(ae[1], 2) + np.power(ae[2], 2)) 811 anr = np.sqrt(np.power(an[0], 2) + np.power(an[1], 2) + np.power(an[2], 2)) 812 ae = ae / aer 813 if not anr: 814 an = np.array([np.nan, np.nan, np.nan]) 815 else: 816 an = an / anr 817 if an[2] <= 0.: 818 an1 = an 819 ae1 = ae 820 else: 821 an1 = -an 822 ae1 = -ae 823 (ft, fd, fl) = tdl(an1, ae1) 824 return NodalPlane(360 - ft, fd, 180 - fl) 825 826 827def tdl(an, bn): 828 """ 829 Helper function for mt2plane. 830 831 Adapted from MATLAB script 832 `bb.m <http://www.ceri.memphis.edu/people/olboyd/Software/Software.html>`_ 833 written by Andy Michael, Chen Ji and Oliver Boyd. 834 """ 835 xn = an[0] 836 yn = an[1] 837 zn = an[2] 838 xe = bn[0] 839 ye = bn[1] 840 ze = bn[2] 841 aaa = 1.0 / (1000000) 842 con = 57.2957795 843 if np.fabs(zn) < aaa: 844 fd = 90. 845 axn = np.fabs(xn) 846 if axn > 1.0: 847 axn = 1.0 848 ft = np.arcsin(axn) * con 849 st = -xn 850 ct = yn 851 if st >= 0. and ct < 0: 852 ft = 180. - ft 853 if st < 0. and ct <= 0: 854 ft = 180. + ft 855 if st < 0. and ct > 0: 856 ft = 360. - ft 857 fl = np.arcsin(abs(ze)) * con 858 sl = -ze 859 if np.fabs(xn) < aaa: 860 cl = xe / yn 861 else: 862 cl = -ye / xn 863 if sl >= 0. and cl < 0: 864 fl = 180. - fl 865 if sl < 0. and cl <= 0: 866 fl = fl - 180. 867 if sl < 0. and cl > 0: 868 fl = -fl 869 else: 870 if -zn > 1.0: 871 zn = -1.0 872 fdh = np.arccos(-zn) 873 fd = fdh * con 874 sd = np.sin(fdh) 875 if sd == 0: 876 return 877 st = -xn / sd 878 ct = yn / sd 879 sx = np.fabs(st) 880 if sx > 1.0: 881 sx = 1.0 882 ft = np.arcsin(sx) * con 883 if st >= 0. and ct < 0: 884 ft = 180. - ft 885 if st < 0. and ct <= 0: 886 ft = 180. + ft 887 if st < 0. and ct > 0: 888 ft = 360. - ft 889 sl = -ze / sd 890 sx = np.fabs(sl) 891 if sx > 1.0: 892 sx = 1.0 893 fl = np.arcsin(sx) * con 894 if st == 0: 895 cl = xe / ct 896 else: 897 xxx = yn * zn * ze / sd / sd + ye 898 cl = -sd * xxx / xn 899 if ct == 0: 900 cl = ye / st 901 if sl >= 0. and cl < 0: 902 fl = 180. - fl 903 if sl < 0. and cl <= 0: 904 fl = fl - 180. 905 if sl < 0. and cl > 0: 906 fl = -fl 907 return (ft, fd, fl) 908 909 910def mt2axes(mt): 911 """ 912 Calculates the principal axes of a given moment tensor. 913 914 :param mt: :class:`~MomentTensor` 915 :return: tuple of :class:`~PrincipalAxis` T, N and P 916 917 Adapted from ps_tensor / utilmeca.c / 918 `Generic Mapping Tools (GMT) <https://gmt.soest.hawaii.edu>`_. 919 """ 920 (d, v) = np.linalg.eigh(mt.mt) 921 pl = np.arcsin(-v[0]) 922 az = np.arctan2(v[2], -v[1]) 923 for i in range(0, 3): 924 if pl[i] <= 0: 925 pl[i] = -pl[i] 926 az[i] += np.pi 927 if az[i] < 0: 928 az[i] += 2 * np.pi 929 if az[i] > 2 * np.pi: 930 az[i] -= 2 * np.pi 931 pl *= R2D 932 az *= R2D 933 934 t = PrincipalAxis(d[2], az[2], pl[2]) 935 n = PrincipalAxis(d[1], az[1], pl[1]) 936 p = PrincipalAxis(d[0], az[0], pl[0]) 937 return (t, n, p) 938 939 940class PrincipalAxis(object): 941 """ 942 A principal axis. 943 944 Strike and dip values are in degrees. 945 946 >>> a = PrincipalAxis(1.3, 20, 50) 947 >>> a.dip 948 50 949 >>> a.strike 950 20 951 >>> a.val 952 1.3 953 """ 954 def __init__(self, val=0, strike=0, dip=0): 955 self.val = val 956 self.strike = strike 957 self.dip = dip 958 959 960class NodalPlane(object): 961 """ 962 A nodal plane. 963 964 All values are in degrees. 965 966 >>> a = NodalPlane(13, 20, 50) 967 >>> a.strike 968 13 969 >>> a.dip 970 20 971 >>> a.rake 972 50 973 """ 974 def __init__(self, strike=0, dip=0, rake=0): 975 self.strike = strike 976 self.dip = dip 977 self.rake = rake 978 979 980class MomentTensor(object): 981 """ 982 A moment tensor. 983 984 >>> a = MomentTensor(1, 1, 0, 0, 0, -1, 26) 985 >>> b = MomentTensor(np.array([1, 1, 0, 0, 0, -1]), 26) 986 >>> c = MomentTensor(np.array([[1, 0, 0], [0, 1, -1], [0, -1, 0]]), 26) 987 >>> a.mt 988 array([[ 1, 0, 0], 989 [ 0, 1, -1], 990 [ 0, -1, 0]]) 991 >>> b.yz 992 -1 993 >>> a.expo 994 26 995 """ 996 def __init__(self, *args): 997 if len(args) == 2: 998 a = args[0] 999 self.expo = args[1] 1000 if len(a) == 6: 1001 # six independent components 1002 self.mt = np.array([[a[0], a[3], a[4]], 1003 [a[3], a[1], a[5]], 1004 [a[4], a[5], a[2]]]) 1005 elif isinstance(a, np.ndarray) and a.shape == (3, 3): 1006 # full matrix 1007 self.mt = a 1008 else: 1009 raise TypeError("Wrong size of input parameter.") 1010 elif len(args) == 7: 1011 # six independent components 1012 self.mt = np.array([[args[0], args[3], args[4]], 1013 [args[3], args[1], args[5]], 1014 [args[4], args[5], args[2]]]) 1015 self.expo = args[6] 1016 else: 1017 raise TypeError("Wrong size of input parameter.") 1018 1019 @property 1020 def normalized(self): 1021 return MomentTensor(self.mt_normalized, self.expo) 1022 1023 @property 1024 def mt_normalized(self): 1025 return self.mt / np.linalg.norm(self.mt) 1026 1027 @property 1028 def xx(self): 1029 return self.mt[0][0] 1030 1031 @property 1032 def xy(self): 1033 return self.mt[0][1] 1034 1035 @property 1036 def xz(self): 1037 return self.mt[0][2] 1038 1039 @property 1040 def yz(self): 1041 return self.mt[1][2] 1042 1043 @property 1044 def yy(self): 1045 return self.mt[1][1] 1046 1047 @property 1048 def zz(self): 1049 return self.mt[2][2] 1050 1051 1052if __name__ == '__main__': 1053 import doctest 1054 doctest.testmod() 1055