1#!/usr/bin/env python 2# -*- coding: utf-8 -*- 3# ------------------------------------------------------------------ 4# Filename: mopad.py 5# Purpose: Moment tensor Plotting and Decomposition tool 6# Author: Lars Krieger, Sebastian Heimann 7# Email: lars.krieger@zmaw.de, sebastian.heimann@zmaw.de 8# 9# Copyright (C) 2010 Lars Krieger, Sebastian Heimann 10# -------------------------------------------------------------------- 11# 12# Completely skip this file for flake8 testing - it is not our code. 13# 14# flake8: noqa 15""" 16MoPaD command line utility. 17 18USAGE: obspy-mopad [plot,decompose,gmt,convert] SOURCE_MECHANISM [OPTIONS] 19 20:: 21 22 ####################################################################### 23 ######################### MoPaD #################################### 24 ######### Moment tensor Plotting and Decomposition tool ############### 25 ####################################################################### 26 27 Multi method tool for: 28 29 - Plotting and saving of focal sphere diagrams ('Beachballs'). 30 31 - Decomposition and Conversion of seismic moment tensors. 32 33 - Generating coordinates, describing a focal sphere diagram, to be 34 piped into GMT's psxy (Useful where psmeca or pscoupe fail.) 35 36 For more help, please run ``python mopad.py --help``. 37 38 39 ####################################################################### 40 41 Version 0.7 42 43 ####################################################################### 44 45 Copyright (C) 2010 46 Lars Krieger & Sebastian Heimann 47 48 Contact 49 lars.krieger@zmaw.de & sebastian.heimann@zmaw.de 50 51 ####################################################################### 52 53 License: 54 55 GNU Lesser General Public License, Version 3 56 57 This program is free software; you can redistribute it and/or 58 modify it under the terms of the GNU Lesser General Public License 59 as published by the Free Software Foundation; either version 3 60 of the License, or (at your option) any later version. 61 62 This program is distributed in the hope that it will be useful, 63 but WITHOUT ANY WARRANTY; without even the implied warranty of 64 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 65 GNU Lesser General Public License for more details. 66 67 You should have received a copy of the GNU Lesser General Public License 68 along with this program; if not, write to the Free Software 69 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 70 02110-1301, USA. 71""" 72from __future__ import (absolute_import, division, print_function, 73 unicode_literals) 74from future.builtins import * # NOQA 75 76import io 77import math 78import os 79import os.path 80import sys 81 82import numpy as np 83 84from obspy import __version__ 85 86 87MOPAD_VERSION = 0.7 88 89 90# constants: 91dynecm = 1e-7 92pi = np.pi 93 94epsilon = 1e-13 95 96rad2deg = 180. / pi 97 98 99class MTError(Exception): 100 pass 101 102 103class MomentTensor: 104 """ 105 """ 106 def __init__(self, M=None, system='NED', debug=0): 107 """ 108 Creates a moment tensor object on the basis of a provided mechanism M. 109 110 If M is a non symmetric 3x3-matrix, the upper right triangle 111 of the matrix is taken as reference. M is symmetrisised 112 w.r.t. these entries. If M is provided as a 3-,4-,6-,7-tuple 113 or array, it is converted into a matrix internally according 114 to standard conventions (Aki & Richards). 115 116 'system' may be chosen as 'NED','USE','NWU', or 'XYZ'. 117 118 'debug' enables output on the shell at the intermediate steps. 119 """ 120 121 source_mechanism = M 122 self._original_M = M[:] 123 124 # basis system: 125 self._input_basis = 'NED' 126 self._list_of_possible_input_bases = ['NED', 'USE', 'XYZ', 'NWU'] 127 self._list_of_possible_output_bases = ['NED', 'USE', 'XYZ', 'NWU'] 128 129 self._input_basis = system.upper() 130 131 # bring M to symmetric matrix form 132 self._M = self._setup_M(source_mechanism) 133 134 # transform M into NED system for internal calculations 135 self._rotate_2_NED() 136 137 # 138 # set attributes list 139 # 140 141 # decomposition: 142 self._decomposition_key = 20 143 144 # eigenvector / principal-axes system: 145 self._eigenvalues = None 146 self._eigenvectors = None 147 self._null_axis = None 148 self._t_axis = None 149 self._p_axis = None 150 self._rotation_matrix = None 151 152 # optional - maybe set afterwards by external application - for later 153 # plotting: 154 self._best_faultplane = None 155 self._auxiliary_plane = None 156 157 # 158 # RUN: 159 # 160 161 # carry out the MT decomposition - results are in basis NED 162 self._decompose_M() 163 164 # set the appropriate principal axis system: 165 self._M_to_principal_axis_system() 166 167 def _setup_M(self, mech): 168 """ 169 Brings the provided mechanism into symmetric 3x3 matrix form. 170 171 The source mechanism may be provided in different forms: 172 173 * as 3x3 matrix - symmetry is checked - one basis system has to be 174 chosen, or NED as default is taken 175 * as 3-element tuple or array - interpreted as strike, dip, slip-rake 176 angles in degree 177 * as 4-element tuple or array - interpreted as strike, dip, slip-rake 178 angles in degree + seismic scalar moment in Nm 179 * as 6-element tuple or array - interpreted as the 6 independent 180 entries of the moment tensor 181 * as 7-element tuple or array - interpreted as the 6 independent 182 entries of the moment tensor + seismic scalar moment in Nm 183 * as 9-element tuple or array - interpreted as the 9 entries of the 184 moment tensor - checked for symmetry 185 * as a nesting of one of the upper types (e.g. a list of n-tuples); 186 first element of outer nesting is taken 187 """ 188 # set source mechanism to matrix form 189 190 if mech is None: 191 print('\n ERROR !! - Please provide a mechanism !!\n') 192 raise MTError(' !! ') 193 194 # if some stupid nesting occurs 195 if len(mech) == 1: 196 mech = mech[0] 197 198 # all 9 elements are given 199 if np.prod(np.shape(mech)) == 9: 200 if np.shape(mech)[0] == 3: 201 # assure symmetry: 202 mech[1, 0] = mech[0, 1] 203 mech[2, 0] = mech[0, 2] 204 mech[2, 1] = mech[1, 2] 205 new_M = mech 206 else: 207 new_M = np.array(mech).reshape(3, 3).copy() 208 new_M[1, 0] = new_M[0, 1] 209 new_M[2, 0] = new_M[0, 2] 210 new_M[2, 1] = new_M[1, 2] 211 212 # mechanism given as 6- or 7-tuple, list or array 213 if len(mech) == 6 or len(mech) == 7: 214 M = mech 215 new_M = np.matrix([M[0], M[3], M[4], 216 M[3], M[1], M[5], 217 M[4], M[5], M[2]]).reshape(3, 3) 218 219 if len(mech) == 7: 220 new_M *= M[6] 221 222 # if given as strike, dip, rake, conventions from Jost & Herrmann hold 223 # resulting matrix is in NED-basis: 224 if len(mech) == 3 or len(mech) == 4: 225 try: 226 [float(val) for val in mech] 227 except Exception: 228 msg = "angles must be given as floats, separated by commas" 229 sys.exit('\n ERROR - %s\n ' % msg) 230 231 strike = mech[0] 232 if not 0 <= strike <= 360: 233 msg = "strike angle must be between 0° and 360°" 234 sys.exit('\n ERROR - %s\n ' % msg) 235 dip = mech[1] 236 if not -90 <= dip <= 90: 237 msg = "dip angle must be between -90° and 90°" 238 sys.exit('\n ERROR - %s\n ' % msg) 239 rake = mech[2] 240 if not -180 <= rake <= 180: 241 msg = "slip-rake angle must be between -180° and 180°" 242 sys.exit('\n ERROR - %s\n ' % msg) 243 244 moms = strikediprake_2_moments(strike, dip, rake) 245 246 new_M = np.matrix([moms[0], moms[3], moms[4], 247 moms[3], moms[1], moms[5], 248 moms[4], moms[5], moms[2]]).reshape(3, 3) 249 250 if len(mech) == 4: 251 new_M *= mech[3] 252 253 # to assure right basis system - others are meaningless, provided 254 # these angles 255 self._input_basis = 'NED' 256 257 return np.asmatrix(new_M) 258 259 def _rotate_2_NED(self): 260 """ 261 Rotates the mechanism to the basis NED. 262 263 All internal calculations are carried out within the NED space. 264 """ 265 if self._input_basis not in self._list_of_possible_input_bases: 266 print('provided input basis not implemented - please specify one', 267 end=' ') 268 print('of the following bases:', 269 self._list_of_possible_input_bases) 270 raise MTError(' !! ') 271 272 NED_2_NED = np.asmatrix(np.diag([1, 1, 1])) 273 274 rotmat_USE_2_NED = NED_2_NED.copy() 275 rotmat_USE_2_NED[:] = 0 276 rotmat_USE_2_NED[0, 1] = -1 277 rotmat_USE_2_NED[1, 2] = 1 278 rotmat_USE_2_NED[2, 0] = -1 279 280 rotmat_XYZ_2_NED = NED_2_NED.copy() 281 rotmat_XYZ_2_NED[:] = 0 282 rotmat_XYZ_2_NED[0, 1] = 1 283 rotmat_XYZ_2_NED[1, 0] = 1 284 rotmat_XYZ_2_NED[2, 2] = -1 285 286 rotmat_NWU_2_NED = NED_2_NED.copy() 287 rotmat_NWU_2_NED[1, 1] = -1 288 rotmat_NWU_2_NED[2, 2] = -1 289 290 if self._input_basis == 'NED': 291 pass 292 elif self._input_basis == 'USE': 293 self._M = np.dot(rotmat_USE_2_NED, 294 np.dot(self._M, rotmat_USE_2_NED.T)) 295 elif self._input_basis == 'XYZ': 296 self._M = np.dot(rotmat_XYZ_2_NED, 297 np.dot(self._M, rotmat_XYZ_2_NED.T)) 298 elif self._input_basis == 'NWU': 299 self._M = np.dot(rotmat_NWU_2_NED, 300 np.dot(self._M, rotmat_NWU_2_NED.T)) 301 302 def _decompose_M(self): 303 """ 304 Running the decomposition of the moment tensor object. 305 306 the standard decompositions M = Isotropic + DC + (CLVD or 2nd DC) are 307 supported (C.f. Jost & Herrmann, Aki & Richards) 308 """ 309 if self._decomposition_key == 20: 310 self._standard_decomposition() 311 elif self._decomposition_key == 21: 312 self._decomposition_w_2DC() 313 elif self._decomposition_key == 31: 314 self._decomposition_w_3DC() 315 else: 316 raise MTError(' only standard decompositions supported ') 317 318 def _standard_decomposition(self): 319 """ 320 Decomposition according Aki & Richards and Jost & Herrmann into 321 322 isotropic + deviatoric 323 = isotropic + DC + CLVD 324 325 parts of the input moment tensor. 326 327 results are given as attributes, callable via the get_* function: 328 329 DC, CLVD, DC_percentage, seismic_moment, moment_magnitude 330 """ 331 M = self._M 332 333 # isotropic part 334 M_iso = np.diag(np.array([1. / 3 * np.trace(M), 335 1. / 3 * np.trace(M), 336 1. / 3 * np.trace(M)])) 337 M0_iso = abs(1. / 3 * np.trace(M)) 338 339 # deviatoric part 340 M_devi = M - M_iso 341 342 self._isotropic = M_iso 343 self._deviatoric = M_devi 344 345 # eigenvalues and -vectors 346 eigenwtot, eigenvtot = np.linalg.eig(M) 347 348 # eigenvalues and -vectors of the deviatoric part 349 eigenw1, eigenv1 = np.linalg.eig(M_devi) 350 351 # eigenvalues in ascending order: 352 eigenw = np.real(np.take(eigenw1, np.argsort(abs(eigenwtot)))) 353 eigenv = np.real(np.take(eigenv1, np.argsort(abs(eigenwtot)), 1)) 354 355 # eigenvalues in ascending order in absolute value!!: 356 eigenw_devi = np.real(np.take(eigenw1, np.argsort(abs(eigenw1)))) 357 358 M0_devi = max(abs(eigenw_devi)) 359 360 # named according to Jost & Herrmann: 361 a1 = eigenv[:, 0] 362 a2 = eigenv[:, 1] 363 a3 = eigenv[:, 2] 364 365 # eigen values can be zero in some cases. this is handled in the 366 # following try/except. 367 with np.errstate(all='ignore'): 368 F = -eigenw_devi[0] / eigenw_devi[2] 369 370 M_DC = \ 371 eigenw[2] * (1 - 2 * F) * (np.outer(a3, a3) - np.outer(a2, a2)) 372 M_CLVD = eigenw[2] * F * (2 * np.outer(a3, a3) - np.outer(a2, a2) - 373 np.outer(a1, a1)) 374 375 try: 376 M_DC_percentage = int(round((1 - 2 * abs(F)) * 100, 6)) 377 except ValueError: 378 # this should only occur in the pure isotropic case 379 M_DC_percentage = 0. 380 381 # according to Bowers & Hudson: 382 M0 = M0_iso + M0_devi 383 384 M_iso_percentage = int(round(M0_iso / M0 * 100, 6)) 385 self._iso_percentage = M_iso_percentage 386 387 self._DC = M_DC 388 self._CLVD = M_CLVD 389 self._DC_percentage = int(round((100 - M_iso_percentage) * 390 M_DC_percentage / 100.)) 391 392 self._seismic_moment = M0 393 self._moment_magnitude = \ 394 np.log10(self._seismic_moment * 1.0e7) / 1.5 - 10.7 395 396 def _decomposition_w_2DC(self): 397 """ 398 Decomposition according Aki & Richards and Jost & Herrmann into 399 400 isotropic + deviatoric 401 = isotropic + DC + DC2 402 403 parts of the input moment tensor. 404 405 results are given as attributes, callable via the get_* function: 406 407 DC1, DC2, DC_percentage, seismic_moment, moment_magnitude 408 """ 409 M = self._M 410 411 # isotropic part 412 M_iso = np.diag(np.array([1. / 3 * np.trace(M), 413 1. / 3 * np.trace(M), 414 1. / 3 * np.trace(M)])) 415 M0_iso = abs(1. / 3 * np.trace(M)) 416 417 # deviatoric part 418 M_devi = M - M_iso 419 420 self._isotropic = M_iso 421 self._deviatoric = M_devi 422 423 # eigenvalues and -vectors of the deviatoric part 424 eigenw1, eigenv1 = np.linalg.eig(M_devi) 425 426 # eigenvalues in ascending order of their absolute values: 427 eigenw = np.real(np.take(eigenw1, np.argsort(abs(eigenw1)))) 428 eigenv = np.real(np.take(eigenv1, np.argsort(abs(eigenw1)), 1)) 429 430 M0_devi = max(abs(eigenw)) 431 432 # named according to Jost & Herrmann: 433 a1 = eigenv[:, 0] 434 a2 = eigenv[:, 1] 435 a3 = eigenv[:, 2] 436 437 M_DC = eigenw[2] * (np.outer(a3, a3) - np.outer(a2, a2)) 438 M_DC2 = eigenw[0] * (np.outer(a1, a1) - np.outer(a2, a2)) 439 440 M_DC_percentage = int(round(abs(eigenw[2] / (abs(eigenw[2]) + 441 abs(eigenw[0]))) * 100.)) 442 443 # according to Bowers & Hudson: 444 M0 = M0_iso + M0_devi 445 446 M_iso_percentage = int(round(M0_iso / M0 * 100)) 447 self._iso_percentage = M_iso_percentage 448 449 self._DC = M_DC 450 self._DC2 = M_DC2 451 # self._DC_percentage = M_DC_percentage 452 self._DC_percentage = int(round((100 - M_iso_percentage) * 453 M_DC_percentage / 100.)) 454 # and M_DC2_percentage? 455 456 self._seismic_moment = M0 457 self._moment_magnitude = \ 458 np.log10(self._seismic_moment * 1.0e7) / 1.5 - 10.7 459 460 def _decomposition_w_3DC(self): 461 """ 462 Decomposition according Aki & Richards and Jost & Herrmann into 463 464 - isotropic 465 - deviatoric 466 - 3 DC 467 468 parts of the input moment tensor. 469 470 results are given as attributes, callable via the get_* function: 471 472 DC1, DC2, DC3, DC_percentage, seismic_moment, moment_magnitude 473 """ 474 M = self._M 475 476 # isotropic part 477 M_iso = np.diag(np.array([1. / 3 * np.trace(M), 478 1. / 3 * np.trace(M), 479 1. / 3 * np.trace(M)])) 480 M0_iso = abs(1. / 3 * np.trace(M)) 481 482 # deviatoric part 483 M_devi = M - M_iso 484 485 self._isotropic = M_iso 486 self._deviatoric = M_devi 487 488 # eigenvalues and -vectors of the deviatoric part 489 eigenw1, eigenv1 = np.linalg.eig(M_devi) 490 M0_devi = max(abs(eigenw1)) 491 492 # eigenvalues and -vectors of the full M !!!!!!!! 493 eigenw1, eigenv1 = np.linalg.eig(M) 494 495 # eigenvalues in ascending order of their absolute values: 496 eigenw = np.real(np.take(eigenw1, np.argsort(abs(eigenw1)))) 497 eigenv = np.real(np.take(eigenv1, np.argsort(abs(eigenw1)), 1)) 498 499 # named according to Jost & Herrmann: 500 a1 = eigenv[:, 0] 501 a2 = eigenv[:, 1] 502 a3 = eigenv[:, 2] 503 504 M_DC1 = 1. / 3. * (eigenw[0] - eigenw[1]) * (np.outer(a1, a1) - 505 np.outer(a2, a2)) 506 M_DC2 = 1. / 3. * (eigenw[1] - eigenw[2]) * (np.outer(a2, a2) - 507 np.outer(a3, a3)) 508 M_DC3 = 1. / 3. * (eigenw[2] - eigenw[0]) * (np.outer(a3, a3) - 509 np.outer(a1, a1)) 510 511 M_DC1_perc = int(100 * abs((eigenw[0] - eigenw[1])) / 512 (abs((eigenw[1] - eigenw[2])) + 513 abs((eigenw[1] - eigenw[2])) + 514 abs((eigenw[2] - eigenw[0])))) 515 M_DC2_perc = int(100 * abs((eigenw[1] - eigenw[2])) / 516 (abs((eigenw[1] - eigenw[2])) + 517 abs((eigenw[1] - eigenw[2])) + 518 abs((eigenw[2] - eigenw[0])))) 519 520 self._DC = M_DC1 521 self._DC2 = M_DC2 522 self._DC3 = M_DC3 523 524 self._DC_percentage = M_DC1_perc 525 self._DC2_percentage = M_DC2_perc 526 527 # according to Bowers & Hudson: 528 M0 = M0_iso + M0_devi 529 530 M_iso_percentage = int(M0_iso / M0 * 100) 531 self._iso_percentage = M_iso_percentage 532 533 # self._seismic_moment = np.sqrt(1./2*nnp.sum(eigenw**2) ) 534 self._seismic_moment = M0 535 self._moment_magnitude = \ 536 np.log10(self._seismic_moment * 1.0e7) / 1.5 - 10.7 537 538 def _M_to_principal_axis_system(self): 539 """ 540 Read in Matrix M and set up eigenvalues (EW) and eigenvectors 541 (EV) for setting up the principal axis system. 542 543 The internal convention is the 'HNS'-system: H is the 544 eigenvector for the smallest absolute eigenvalue, S is the 545 eigenvector for the largest absolute eigenvalue, N is the null 546 axis. 547 548 Naming due to the geometry: a CLVD is 549 Symmetric to the S-axis, 550 Null-axis is common sense, and the third (auxiliary) axis 551 Helps to construct the R³. 552 553 Additionally builds matrix for basis transformation back to NED system. 554 555 The eigensystem setup defines the colouring order for a later 556 plotting in the BeachBall class. This order is set by the 557 '_plot_clr_order' attribute. 558 """ 559 M = self._M 560 M_devi = self._deviatoric 561 562 # working in framework of 3 principal axes: 563 # eigenvalues (EW) are in order from high to low 564 # - neutral axis N, belongs to middle EW 565 # - symmetry axis S ('sigma') belongs to EW with largest absolute value 566 # (P- or T-axis) 567 # - auxiliary axis H ('help') belongs to remaining EW (T- or P-axis) 568 # EW sorting from lowest to highest value 569 EW_devi, EV_devi = np.linalg.eigh(M_devi) 570 EW_order = np.argsort(EW_devi) 571 572 # print('order', EW_order) 573 574 if 1: # self._plot_isotropic_part: 575 trace_M = np.trace(M) 576 if abs(trace_M) < epsilon: 577 trace_M = 0 578 EW, EV = np.linalg.eigh(M) 579 for i, ew in enumerate(EW): 580 if abs(EW[i]) < epsilon: 581 EW[i] = 0 582 else: 583 trace_M = np.trace(M_devi) 584 if abs(trace_M) < epsilon: 585 trace_M = 0 586 587 EW, EV = np.linalg.eigh(M_devi) 588 for i, ew in enumerate(EW): 589 if abs(EW[i]) < epsilon: 590 EW[i] = 0 591 592 EW1_devi = EW_devi[EW_order[0]] 593 EW2_devi = EW_devi[EW_order[1]] 594 EW3_devi = EW_devi[EW_order[2]] 595 EV1_devi = EV_devi[:, EW_order[0]] 596 EV2_devi = EV_devi[:, EW_order[1]] 597 EV3_devi = EV_devi[:, EW_order[2]] 598 599 EW1 = EW[EW_order[0]] 600 EW2 = EW[EW_order[1]] 601 EW3 = EW[EW_order[2]] 602 EV1 = EV[:, EW_order[0]] 603 EV2 = EV[:, EW_order[1]] 604 EV3 = EV[:, EW_order[2]] 605 606 chng_basis_tmp = np.asmatrix(np.zeros((3, 3))) 607 chng_basis_tmp[:, 0] = EV1_devi 608 chng_basis_tmp[:, 1] = EV2_devi 609 chng_basis_tmp[:, 2] = EV3_devi 610 611 symmetry_around_tension = 1 612 clr = 1 613 614 if abs(EW2_devi) < epsilon: 615 EW2_devi = 0 616 617 # implosion 618 if EW1 < 0 and EW2 < 0 and EW3 < 0: 619 symmetry_around_tension = 0 620 # logger.debug( 'IMPLOSION - symmetry around pressure axis \n\n') 621 clr = 1 622 # explosion 623 elif EW1 > 0 and EW2 > 0 and EW3 > 0: 624 symmetry_around_tension = 1 625 if abs(EW1_devi) > abs(EW3_devi): 626 symmetry_around_tension = 0 627 # logger.debug( 'EXPLOSION - symmetry around tension axis \n\n') 628 clr = -1 629 # net-implosion 630 elif EW2 < 0 and sum([EW1, EW2, EW3]) < 0: 631 if abs(EW1_devi) < abs(EW3_devi): 632 symmetry_around_tension = 1 633 clr = 1 634 else: 635 symmetry_around_tension = 1 636 clr = 1 637 # net-implosion 638 elif EW2_devi >= 0 and sum([EW1, EW2, EW3]) < 0: 639 symmetry_around_tension = 0 640 clr = -1 641 if abs(EW1_devi) < abs(EW3_devi): 642 symmetry_around_tension = 1 643 clr = 1 644 # net-explosion 645 elif EW2_devi < 0 and sum([EW1, EW2, EW3]) > 0: 646 symmetry_around_tension = 1 647 clr = 1 648 if abs(EW1_devi) > abs(EW3_devi): 649 symmetry_around_tension = 0 650 clr = -1 651 # net-explosion 652 elif EW2_devi >= 0 and sum([EW1, EW2, EW3]) > 0: 653 symmetry_around_tension = 0 654 clr = -1 655 else: 656 pass 657 if abs(EW1_devi) < abs(EW3_devi): 658 symmetry_around_tension = 1 659 clr = 1 660 if 0: # EW2 > 0 :#or (EW2 > 0 and EW2_devi > 0) : 661 symmetry_around_tension = 0 662 clr = -1 663 664 if abs(EW1_devi) >= abs(EW3_devi): 665 symmetry_around_tension = 0 666 clr = -1 667 if 0: # EW2 < 0 : 668 symmetry_around_tension = 1 669 clr = 1 670 if (EW3 < 0 and np.trace(self._M) >= 0): 671 print('Houston, we have had a problem - check M !!!!!!') 672 raise MTError(' !! ') 673 674 if trace_M == 0: 675 if EW2 == 0: 676 symmetry_around_tension = 1 677 clr = 1 678 elif 2 * abs(EW2) == abs(EW1) or 2 * abs(EW2) == abs(EW3): 679 if abs(EW1) < EW3: 680 symmetry_around_tension = 1 681 clr = 1 682 else: 683 symmetry_around_tension = 0 684 clr = -1 685 else: 686 if abs(EW1) < EW3: 687 symmetry_around_tension = 1 688 clr = 1 689 else: 690 symmetry_around_tension = 0 691 clr = -1 692 693 if symmetry_around_tension == 1: 694 EWs = EW3.copy() 695 EVs = EV3.copy() 696 EWh = EW1.copy() 697 EVh = EV1.copy() 698 699 else: 700 EWs = EW1.copy() 701 EVs = EV1.copy() 702 EWh = EW3.copy() 703 EVh = EV3.copy() 704 705 EWn = EW2 706 EVn = EV2 707 708 # build the basis system change matrix: 709 chng_basis = np.asmatrix(np.zeros((3, 3))) 710 711 # order of eigenvector's basis: (H,N,S) 712 chng_basis[:, 0] = EVh 713 chng_basis[:, 1] = EVn 714 chng_basis[:, 2] = EVs 715 716 # matrix for basis transformation 717 self._rotation_matrix = chng_basis 718 719 # collections of eigenvectors and eigenvalues 720 self._eigenvectors = [EVh, EVn, EVs] 721 self._eigenvalues = [EWh, EWn, EWs] 722 723 # principal axes 724 self._null_axis = EVn 725 self._t_axis = EV1 726 self._p_axis = EV3 727 728 # plotting order flag - important for plot in BeachBall class 729 self._plot_clr_order = clr 730 731 # collection of the faultplanes, given in strike, dip, slip-rake 732 self._faultplanes = self._find_faultplanes() 733 734 def _find_faultplanes(self): 735 """ 736 Sets the two angle-triples, describing the faultplanes of the 737 Double Couple, defined by the eigenvectors P and T of the 738 moment tensor object. 739 740 Defining a reference Double Couple with strike = dip = 741 slip-rake = 0, the moment tensor object's DC is transformed 742 (rotated) w.r.t. this orientation. The respective rotation 743 matrix yields the first fault plane angles as the Euler 744 angles. After flipping the first reference plane by 745 multiplying the appropriate flip-matrix, one gets the second fault 746 plane's geometry. 747 748 All output angles are in degree 749 750 ( 751 to check: 752 mit Sebastians Konventionen: 753 754 rotationsmatrix1 = EV Matrix von M, allerdings in der Reihenfolge TNP 755 (nicht, wie hier PNT!!!) 756 757 referenz-DC mit strike, dip, rake = 0,0,0 in NED - Darstellung: 758 M = 0,0,0,0,-1,0 759 760 davon die EV ebenfalls in eine Matrix: 761 762 trafo-matrix2 = EV Matrix von Referenz-DC in der REihenfolge TNP 763 764 effektive Rotationsmatrix = (rotationsmatrix1 * trafo-matrix2.T).T 765 766 durch check, ob det <0, schauen, ob die Matrix mit -1 multipliziert 767 werden muss 768 769 flip_matrix = 0,0,-1,0,-1,0,-1,0,0 770 771 andere DC Orientierung wird durch flip * effektive Rotationsmatrix 772 erhalten 773 774 beide Rotataionmatrizen in matrix_2_euler 775 ) 776 """ 777 # reference Double Couple (in NED basis) 778 # it has strike, dip, slip-rake = 0,0,0 779 refDC = np.matrix([[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], 780 dtype=np.float) 781 refDC_evals, refDC_evecs = np.linalg.eigh(refDC) 782 783 # matrix which is turning from one fault plane to the other 784 flip_dc = np.matrix([[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], 785 dtype=np.float) 786 787 # euler-tools need matrices of EV sorted in PNT: 788 pnt_sorted_EV_matrix = self._rotation_matrix.copy() 789 790 # resort only necessary, if abs(p) <= abs(t) 791 # print(self._plot_clr_order) 792 if self._plot_clr_order < 0: 793 pnt_sorted_EV_matrix[:, 0] = self._rotation_matrix[:, 2] 794 pnt_sorted_EV_matrix[:, 2] = self._rotation_matrix[:, 0] 795 796 # rotation matrix, describing the rotation of the eigenvector 797 # system of the input moment tensor into the eigenvector 798 # system of the reference Double Couple 799 rot_matrix_fp1 = (np.dot(pnt_sorted_EV_matrix, refDC_evecs.T)).T 800 801 # check, if rotation has right orientation 802 if np.linalg.det(rot_matrix_fp1) < 0.: 803 rot_matrix_fp1 *= -1. 804 805 # adding a rotation into the ambiguous system of the second fault plane 806 rot_matrix_fp2 = np.dot(flip_dc, rot_matrix_fp1) 807 808 fp1 = self._find_strike_dip_rake(rot_matrix_fp1) 809 fp2 = self._find_strike_dip_rake(rot_matrix_fp2) 810 811 return [fp1, fp2] 812 813 def _find_strike_dip_rake(self, rotation_matrix): 814 """ 815 Returns angles strike, dip, slip-rake in degrees, describing the fault 816 plane. 817 """ 818 (alpha, beta, gamma) = self._matrix_to_euler(rotation_matrix) 819 return (beta * rad2deg, alpha * rad2deg, -gamma * rad2deg) 820 821 def _cvec(self, x, y, z): 822 """ 823 Builds a column vector (matrix type) from a 3 tuple. 824 """ 825 return np.matrix([[x, y, z]], dtype=np.float).T 826 827 def _matrix_to_euler(self, rotmat): 828 """ 829 Returns three Euler angles alpha, beta, gamma (in radians) from a 830 rotation matrix. 831 """ 832 ex = self._cvec(1., 0., 0.) 833 ez = self._cvec(0., 0., 1.) 834 exs = rotmat.T * ex 835 ezs = rotmat.T * ez 836 enodes = np.cross(ez.T, ezs.T).T 837 if np.linalg.norm(enodes) < 1e-10: 838 enodes = exs 839 enodess = rotmat * enodes 840 cos_alpha = float((ez.T * ezs)) 841 if cos_alpha > 1.: 842 cos_alpha = 1. 843 if cos_alpha < -1.: 844 cos_alpha = -1. 845 alpha = np.arccos(cos_alpha) 846 beta = np.mod(np.arctan2(enodes[1, 0], enodes[0, 0]), np.pi * 2.) 847 gamma = np.mod(-np.arctan2(enodess[1, 0], enodess[0, 0]), np.pi * 2.) 848 return self._unique_euler(alpha, beta, gamma) 849 850 def _unique_euler(self, alpha, beta, gamma): 851 """ 852 Uniquify euler angle triplet. 853 854 Puts euler angles into ranges compatible with (dip,strike,-rake) in 855 seismology: 856 857 alpha (dip) : [0, pi/2] 858 beta (strike) : [0, 2*pi) 859 gamma (-rake) : [-pi, pi) 860 861 If alpha is near to zero, beta is replaced by beta+gamma and gamma is 862 set to zero, to prevent that additional ambiguity. 863 864 If alpha is near to pi/2, beta is put into the range [0,pi). 865 """ 866 alpha = np.mod(alpha, 2.0 * pi) 867 868 if 0.5 * pi < alpha and alpha <= pi: 869 alpha = pi - alpha 870 beta = beta + pi 871 gamma = 2.0 * pi - gamma 872 elif pi < alpha and alpha <= 1.5 * pi: 873 alpha = alpha - pi 874 gamma = pi - gamma 875 elif 1.5 * pi < alpha and alpha <= 2.0 * pi: 876 alpha = 2.0 * pi - alpha 877 beta = beta + pi 878 gamma = pi + gamma 879 880 alpha = np.mod(alpha, 2.0 * pi) 881 beta = np.mod(beta, 2.0 * pi) 882 gamma = np.mod(gamma + pi, 2.0 * pi) - pi 883 884 # If dip is exactly 90 degrees, one is still 885 # free to choose between looking at the plane from either side. 886 # Choose to look at such that beta is in the range [0,180) 887 888 # This should prevent some problems, when dip is close to 90 degrees: 889 if abs(alpha - 0.5 * pi) < 1e-10: 890 alpha = 0.5 * pi 891 if abs(beta - pi) < 1e-10: 892 beta = pi 893 if abs(beta - 2. * pi) < 1e-10: 894 beta = 0. 895 if abs(beta) < 1e-10: 896 beta = 0. 897 898 if alpha == 0.5 * pi and beta >= pi: 899 gamma = -gamma 900 beta = np.mod(beta - pi, 2.0 * pi) 901 gamma = np.mod(gamma + pi, 2.0 * pi) - pi 902 assert 0. <= beta < pi 903 assert -pi <= gamma < pi 904 905 if alpha < 1e-7: 906 beta = np.mod(beta + gamma, 2.0 * pi) 907 gamma = 0. 908 909 return (alpha, beta, gamma) 910 911 def _matrix_w_style_and_system(self, M2return, system, style): 912 """ 913 Gives the provided matrix in the desired basis system. 914 915 If the argument 'style' is set to 'fancy', a 'print' of the return 916 value yields a nice shell output of the matrix for better 917 visual control. 918 """ 919 if not system.upper() in self._list_of_possible_output_bases: 920 print('\nprovided output basis not supported - please specify', 921 end=' ') 922 print('one of the following bases: (default=NED)\n', end=' ') 923 print(self._list_of_possible_input_bases, '\n') 924 raise MTError(' !! ') 925 926 fancy = 0 927 if style.lower() in ['f', 'fan', 'fancy']: 928 fancy = 1 929 930 if system.upper() == 'NED': 931 if fancy: 932 return fancy_matrix(M2return) 933 else: 934 return M2return 935 936 elif system.upper() == 'USE': 937 if fancy: 938 return fancy_matrix(NED2USE(M2return)) 939 else: 940 return NED2USE(M2return) 941 942 elif system.upper() == 'XYZ': 943 if fancy: 944 return fancy_matrix(NED2XYZ(M2return)) 945 else: 946 return NED2XYZ(M2return) 947 948 elif system.upper() == 'NWU': 949 if fancy: 950 return fancy_matrix(NED2NWU(M2return)) 951 else: 952 return NED2NWU(M2return) 953 954 def _vector_w_style_and_system(self, vectors, system, style): 955 """ 956 Gives the provided vector(s) in the desired basis system. 957 958 If the argument 'style' is set to 'fancy', a 'print' of the return 959 value yields a nice shell output of the vector(s) for better 960 visual control. 961 962 'vectors' can be either a single array, tuple, matrix or a collection 963 in form of a list, array or matrix. 964 If it's a list, each entry will be checked, if it's 3D - if not, an 965 exception is raised. 966 If it's a matrix or array with column-length 3, the columns are 967 interpreted as vectors, otherwise, its transposed is used. 968 """ 969 if not system.upper() in self._list_of_possible_output_bases: 970 print('\n provided output basis not supported - please specify', 971 end=' ') 972 print('one of the following bases: (default=NED)\n', end=' ') 973 print(self._list_of_possible_input_bases, '\n') 974 raise MTError(' !! ') 975 976 fancy = 0 977 if style.lower() in ['f', 'fan', 'fancy']: 978 fancy = 1 979 980 lo_vectors = [] 981 982 # if list of vectors 983 if isinstance(vectors, list): 984 for vec in vectors: 985 if np.prod(np.shape(vec)) != 3: 986 print('\n please provide vector(s) from R³ \n ') 987 raise MTError(' !! ') 988 lo_vectors = vectors 989 else: 990 if np.prod(np.shape(vectors)) % 3 != 0: 991 print('\n please provide vector(s) from R³ \n ') 992 raise MTError(' !! ') 993 994 if np.shape(vectors)[0] == 3: 995 for ii in range(np.shape(vectors)[1]): 996 lo_vectors.append(vectors[:, ii]) 997 else: 998 for ii in range(np.shape(vectors)[0]): 999 lo_vectors.append(vectors[:, ii].transpose()) 1000 1001 lo_vecs_to_show = [] 1002 1003 for vec in lo_vectors: 1004 if system.upper() == 'NED': 1005 if fancy: 1006 lo_vecs_to_show.append(fancy_vector(vec)) 1007 else: 1008 lo_vecs_to_show.append(vec) 1009 elif system.upper() == 'USE': 1010 if fancy: 1011 lo_vecs_to_show.append(fancy_vector(NED2USE(vec))) 1012 else: 1013 lo_vecs_to_show.append(NED2USE(vec)) 1014 elif system.upper() == 'XYZ': 1015 if fancy: 1016 lo_vecs_to_show.append(fancy_vector(NED2XYZ(vec))) 1017 else: 1018 lo_vecs_to_show.append(NED2XYZ(vec)) 1019 elif system.upper() == 'NWU': 1020 if fancy: 1021 lo_vecs_to_show.append(fancy_vector(NED2NWU(vec))) 1022 else: 1023 lo_vecs_to_show.append(NED2NWU(vec)) 1024 1025 if len(lo_vecs_to_show) == 1: 1026 return lo_vecs_to_show[0] 1027 else: 1028 if fancy: 1029 return ''.join(lo_vecs_to_show) 1030 else: 1031 return lo_vecs_to_show 1032 1033 def get_M(self, system='NED', style='n'): 1034 """ 1035 Returns the moment tensor in matrix representation. 1036 1037 Call with arguments to set ouput in other basis system or in fancy 1038 style (to be viewed with 'print') 1039 """ 1040 if style == 'f': 1041 print('\n Full moment tensor in %s-coordinates:' % (system)) 1042 return self._matrix_w_style_and_system(self._M, system, style) 1043 else: 1044 return self._matrix_w_style_and_system(self._M, system, style) 1045 1046 def get_decomposition(self, in_system='NED', out_system='NED', style='n'): 1047 """ 1048 Returns a tuple of the decomposition results. 1049 1050 Order: 1051 - 1 - basis of the provided input (string) 1052 - 2 - basis of the representation (string) 1053 - 3 - chosen decomposition type (integer) 1054 1055 - 4 - full moment tensor (matrix) 1056 1057 - 5 - isotropic part (matrix) 1058 - 6 - isotropic percentage (float) 1059 - 7 - deviatoric part (matrix) 1060 - 8 - deviatoric percentage (float) 1061 1062 - 9 - DC part (matrix) 1063 -10 - DC percentage (float) 1064 -11 - DC2 part (matrix) 1065 -12 - DC2 percentage (float) 1066 -13 - DC3 part (matrix) 1067 -14 - DC3 percentage (float) 1068 1069 -15 - CLVD part (matrix) 1070 -16 - CLVD percentage (matrix) 1071 1072 -17 - seismic moment (float) 1073 -18 - moment magnitude (float) 1074 1075 -19 - eigenvectors (3-array) 1076 -20 - eigenvalues (list) 1077 -21 - p-axis (3-array) 1078 -22 - neutral axis (3-array) 1079 -23 - t-axis (3-array) 1080 -24 - faultplanes (list of two 3-arrays) 1081 """ 1082 return [in_system, out_system, self.get_decomp_type(), 1083 self.get_M(system=out_system), 1084 self.get_iso(system=out_system), self.get_iso_percentage(), 1085 self.get_devi(system=out_system), self.get_devi_percentage(), 1086 self.get_DC(system=out_system), self.get_DC_percentage(), 1087 self.get_DC2(system=out_system), self.get_DC2_percentage(), 1088 self.get_DC3(system=out_system), self.get_DC3_percentage(), 1089 self.get_CLVD(system=out_system), self.get_CLVD_percentage(), 1090 self.get_moment(), self.get_mag(), 1091 self.get_eigvecs(system=out_system), 1092 self.get_eigvals(system=out_system), 1093 self.get_p_axis(system=out_system), 1094 self.get_null_axis(system=out_system), 1095 self.get_t_axis(system=out_system), 1096 self.get_fps()] 1097 1098 def get_full_decomposition(self): 1099 """ 1100 Nice compilation of decomposition result to be viewed in the shell 1101 (call with 'print'). 1102 """ 1103 mexp = pow(10, np.ceil(np.log10(np.max(np.abs(self._M))))) 1104 m = self._M / mexp 1105 s = '\nScalar Moment: M0 = %g Nm (Mw = %3.1f)\n' 1106 s += 'Moment Tensor: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,\n' 1107 s += ' Mne = %6.3f, Mnd = %6.3f, Med = %6.3f ' 1108 s += '[ x %g ]\n\n' 1109 s = s % (self._seismic_moment, self._moment_magnitude, m[0, 0], 1110 m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2], mexp) 1111 s += self._fault_planes_as_str() 1112 return s 1113 1114 def _fault_planes_as_str(self): 1115 """ 1116 Internal setup of a nice string, containing information about the fault 1117 planes. 1118 """ 1119 s = '\n' 1120 for i, sdr in enumerate(self.get_fps()): 1121 s += 'Fault plane %i: ' % (i + 1) 1122 s += 'strike = %3.0f°, dip = %3.0f°, slip-rake = %4.0f°\n' % \ 1123 (sdr[0], sdr[1], sdr[2]) 1124 return s 1125 1126 def get_input_system(self, style='n', **kwargs): 1127 """ 1128 Returns the basis system of the input. 1129 """ 1130 if style == 'f': 1131 print('\n Basis system of the input:\n ') 1132 return self._input_basis 1133 1134 def get_output_system(self, style='n', **kwargs): 1135 """ 1136 Returns the basis system of the input. 1137 """ 1138 if style == 'f': 1139 print('\n Basis system of the output: \n ') 1140 return self._output_basis 1141 1142 def get_decomp_type(self, style='n', **kwargs): 1143 """ 1144 Returns the decomposition type. 1145 """ 1146 decomp_dict = dict(zip(('20', '21', '31'), 1147 ('ISO + DC + CLVD', 1148 'ISO + major DC + minor DC', 1149 'ISO + DC1 + DC2 + DC3'))) 1150 if style == 'f': 1151 print('\n Decomposition type: \n ') 1152 return decomp_dict[str(self._decomposition_key)] 1153 1154 return self._decomposition_key 1155 1156 def get_iso(self, system='NED', style='n'): 1157 """ 1158 Returns the isotropic part of the moment tensor in matrix 1159 representation. 1160 1161 Call with arguments to set ouput in other basis system or in fancy 1162 style (to be viewed with 'print') 1163 """ 1164 if style == 'f': 1165 print('\n Isotropic part in %s-coordinates: ' % (system)) 1166 return self._matrix_w_style_and_system(self._isotropic, system, style) 1167 1168 def get_devi(self, system='NED', style='n'): 1169 """ 1170 Returns the deviatoric part of the moment tensor in matrix 1171 representation. 1172 1173 Call with arguments to set ouput in other basis system or in fancy 1174 style (to be viewed with 'print') 1175 """ 1176 if style == 'f': 1177 print('\n Deviatoric part in %s-coordinates: ' % (system)) 1178 return self._matrix_w_style_and_system(self._deviatoric, system, style) 1179 1180 def get_DC(self, system='NED', style='n'): 1181 """ 1182 Returns the Double Couple part of the moment tensor in matrix 1183 representation. 1184 1185 Call with arguments to set ouput in other basis system or in fancy 1186 style (to be viewed with 'print') 1187 """ 1188 if style == 'f': 1189 print('\n Double Couple part in %s-coordinates:' % (system)) 1190 return self._matrix_w_style_and_system(self._DC, system, style) 1191 1192 def get_DC2(self, system='NED', style='n'): 1193 """ 1194 Returns the second Double Couple part of the moment tensor in matrix 1195 representation. 1196 1197 Call with arguments to set ouput in other basis system or in fancy 1198 style (to be viewed with 'print') 1199 """ 1200 if style == 'f': 1201 print('\n second Double Couple part in %s-coordinates:' % (system)) 1202 if self._DC2 is None: 1203 if style == 'f': 1204 print(' not available in this decomposition type ') 1205 return '' 1206 1207 return self._matrix_w_style_and_system(self._DC2, system, style) 1208 1209 def get_DC3(self, system='NED', style='n'): 1210 """ 1211 Returns the third Double Couple part of the moment tensor in matrix 1212 representation. 1213 1214 Call with arguments to set ouput in other basis system or in fancy 1215 style (to be viewed with 'print') 1216 """ 1217 if style == 'f': 1218 print('\n third Double Couple part in %s-coordinates:' % (system)) 1219 1220 if self._DC3 is None: 1221 if style == 'f': 1222 print(' not available in this decomposition type ') 1223 return '' 1224 return self._matrix_w_style_and_system(self._DC3, system, style) 1225 1226 def get_CLVD(self, system='NED', style='n'): 1227 """ 1228 Returns the CLVD part of the moment tensor in matrix representation. 1229 1230 Call with arguments to set ouput in other basis system or in fancy 1231 style (to be viewed with 'print') 1232 """ 1233 if style == 'f': 1234 print('\n CLVD part in %s-coordinates: \n' % (system)) 1235 if self._CLVD is None: 1236 if style == 'f': 1237 print(' not available in this decomposition type ') 1238 return '' 1239 1240 return self._matrix_w_style_and_system(self._CLVD, system, style) 1241 1242 def get_DC_percentage(self, system='NED', style='n'): 1243 """ 1244 Returns the percentage of the DC part of the moment tensor in matrix 1245 representation. 1246 """ 1247 if style == 'f': 1248 print('\n Double Couple percentage: \n') 1249 return self._DC_percentage 1250 1251 def get_CLVD_percentage(self, system='NED', style='n'): 1252 """ 1253 Returns the percentage of the DC part of the moment tensor in matrix 1254 representation. 1255 """ 1256 if style == 'f': 1257 print('\n CLVD percentage: \n') 1258 if self._CLVD is None: 1259 if style == 'f': 1260 print(' not available in this decomposition type ') 1261 return '' 1262 return int(100 - self._DC_percentage - self._iso_percentage) 1263 1264 def get_DC2_percentage(self, system='NED', style='n'): 1265 """ 1266 Returns the percentage of the second DC part of the moment tensor in 1267 matrix representation. 1268 """ 1269 if style == 'f': 1270 print("\n second Double Couple's percentage: \n") 1271 if self._DC2 is None: 1272 if style == 'f': 1273 print(' not available in this decomposition type ') 1274 return '' 1275 return self._DC2_percentage 1276 1277 def get_DC3_percentage(self, system='NED', style='n'): 1278 """ 1279 Returns the percentage of the third DC part of the moment tensor in 1280 matrix representation. 1281 """ 1282 if style == 'f': 1283 print("\n third Double Couple percentage: \n") 1284 if self._DC3 is None: 1285 if style == 'f': 1286 print(' not available in this decomposition type ') 1287 return '' 1288 return int(100 - self._DC2_percentage - self._DC_percentage) 1289 1290 def get_iso_percentage(self, system='NED', style='n'): 1291 """ 1292 Returns the percentage of the isotropic part of the moment tensor in 1293 matrix representation. 1294 """ 1295 if style == 'f': 1296 print('\n Isotropic percentage: \n') 1297 return self._iso_percentage 1298 1299 def get_devi_percentage(self, system='NED', style='n'): 1300 """ 1301 Returns the percentage of the deviatoric part of the moment tensor in 1302 matrix representation. 1303 """ 1304 if style == 'f': 1305 print('\n Deviatoric percentage: \n') 1306 return int(100 - self._iso_percentage) 1307 1308 def get_moment(self, system='NED', style='n'): 1309 """ 1310 Returns the seismic moment (in Nm) of the moment tensor. 1311 """ 1312 if style == 'f': 1313 print('\n Seismic moment (in Nm) : \n ') 1314 return self._seismic_moment 1315 1316 def get_mag(self, system='NED', style='n'): 1317 """ 1318 Returns the moment magnitude M_w of the moment tensor. 1319 """ 1320 if style == 'f': 1321 print('\n Moment magnitude Mw: \n ') 1322 return self._moment_magnitude 1323 1324 def get_decomposition_key(self, system='NED', style='n'): 1325 """ 1326 10 = standard decomposition (Jost & Herrmann) 1327 """ 1328 if style == 'f': 1329 print('\n Decomposition key (standard = 10): \n ') 1330 return self._decomposition_key 1331 1332 def get_eigvals(self, system='NED', style='n', **kwargs): 1333 """ 1334 Returns a list of the eigenvalues of the moment tensor. 1335 """ 1336 if style == 'f': 1337 if self._plot_clr_order < 0: 1338 print('\n Eigenvalues T N P :\n') 1339 else: 1340 print('\n Eigenvalues P N T :\n') 1341 # in the order HNS: 1342 return self._eigenvalues 1343 1344 def get_eigvecs(self, system='NED', style='n'): 1345 """ 1346 Returns the eigenvectors of the moment tensor. 1347 1348 Call with arguments to set ouput in other basis system or in fancy 1349 style (to be viewed with 'print') 1350 """ 1351 if style == 'f': 1352 1353 if self._plot_clr_order < 0: 1354 print('\n Eigenvectors T N P (in basis system %s): ' % 1355 (system)) 1356 else: 1357 print('\n Eigenvectors P N T (in basis system %s): ' % 1358 (system)) 1359 1360 return self._vector_w_style_and_system(self._eigenvectors, system, 1361 style) 1362 1363 def get_null_axis(self, system='NED', style='n'): 1364 """ 1365 Returns the neutral axis of the moment tensor. 1366 1367 Call with arguments to set ouput in other basis system or in fancy 1368 style (to be viewed with 'print') 1369 """ 1370 if style == 'f': 1371 print('\n Null-axis in %s -coordinates: ' % (system)) 1372 return self._vector_w_style_and_system(self._null_axis, system, style) 1373 1374 def get_t_axis(self, system='NED', style='n'): 1375 """ 1376 Returns the tension axis of the moment tensor. 1377 1378 Call with arguments to set ouput in other basis system or in fancy 1379 style (to be viewed with 'print') 1380 """ 1381 if style == 'f': 1382 print('\n Tension-axis in %s -coordinates: ' % (system)) 1383 return self._vector_w_style_and_system(self._t_axis, system, style) 1384 1385 def get_p_axis(self, system='NED', style='n'): 1386 """ 1387 Returns the pressure axis of the moment tensor. 1388 1389 Call with arguments to set ouput in other basis system or in fancy 1390 style (to be viewed with 'print') 1391 """ 1392 if style == 'f': 1393 print('\n Pressure-axis in %s -coordinates: ' % (system)) 1394 return self._vector_w_style_and_system(self._p_axis, system, style) 1395 1396 def get_transform_matrix(self, system='NED', style='n'): 1397 """ 1398 Returns the transformation matrix (input system to principal axis 1399 system. 1400 1401 Call with arguments to set ouput in other basis system or in fancy 1402 style (to be viewed with 'print') 1403 """ 1404 if style == 'f': 1405 print('\n rotation matrix in %s -coordinates: ' % (system)) 1406 return self._matrix_w_style_and_system(self._rotation_matrix, system, 1407 style) 1408 1409 def get_fps(self, **kwargs): 1410 """ 1411 Returns a list of the two faultplane 3-tuples, each showing strike, 1412 dip, slip-rake. 1413 """ 1414 fancy_key = kwargs.get('style', '0') 1415 if fancy_key[0].lower() == 'f': 1416 return self._fault_planes_as_str() 1417 else: 1418 return self._faultplanes 1419 1420 def get_colour_order(self, **kwargs): 1421 """ 1422 Returns the value of the plotting order (only important in BeachBall 1423 instances). 1424 """ 1425 style = kwargs.get('style', '0')[0].lower() 1426 if style == 'f': 1427 print('\n Colour order key: ') 1428 return self._plot_clr_order 1429 1430 1431# --------------------------------------------------------------- 1432# 1433# external functions: 1434# 1435# --------------------------------------------------------------- 1436 1437def _puzzle_basis_transformation(mat_tup_arr_vec, in_basis, out_basis): 1438 lo_bases = ['NED', 'USE', 'XYZ', 'NWU'] 1439 if (in_basis not in lo_bases) and (out_basis in lo_bases): 1440 sys.exit('wrong basis chosen') 1441 1442 if in_basis == out_basis: 1443 transformed_in = mat_tup_arr_vec 1444 elif in_basis == 'NED': 1445 if out_basis == 'USE': 1446 transformed_in = NED2USE(mat_tup_arr_vec) 1447 if out_basis == 'XYZ': 1448 transformed_in = NED2XYZ(mat_tup_arr_vec) 1449 if out_basis == 'NWU': 1450 transformed_in = NED2NWU(mat_tup_arr_vec) 1451 elif in_basis == 'USE': 1452 if out_basis == 'NED': 1453 transformed_in = USE2NED(mat_tup_arr_vec) 1454 if out_basis == 'XYZ': 1455 transformed_in = USE2XYZ(mat_tup_arr_vec) 1456 if out_basis == 'NWU': 1457 transformed_in = USE2NWU(mat_tup_arr_vec) 1458 elif in_basis == 'XYZ': 1459 if out_basis == 'NED': 1460 transformed_in = XYZ2NED(mat_tup_arr_vec) 1461 if out_basis == 'USE': 1462 transformed_in = XYZ2USE(mat_tup_arr_vec) 1463 if out_basis == 'NWU': 1464 transformed_in = XYZ2NWU(mat_tup_arr_vec) 1465 elif in_basis == 'NWU': 1466 if out_basis == 'NED': 1467 transformed_in = NWU2NED(mat_tup_arr_vec) 1468 if out_basis == 'USE': 1469 transformed_in = NWU2USE(mat_tup_arr_vec) 1470 if out_basis == 'XYZ': 1471 transformed_in = NWU2XYZ(mat_tup_arr_vec) 1472 1473 if len(mat_tup_arr_vec) == 3 and np.prod(np.shape(mat_tup_arr_vec)) != 9: 1474 tmp_array = np.array([0, 0, 0]) 1475 tmp_array[:] = transformed_in 1476 return tmp_array 1477 else: 1478 return transformed_in 1479 1480 1481def _return_matrix_vector_array(ma_ve_ar, basis_change_matrix): 1482 """ 1483 Generates the output for the functions, yielding matrices, vectors, and 1484 arrays in new basis systems. 1485 1486 Allowed input are 3x3 matrices, 3-vectors, 3-vector collections, 1487 3-arrays, and 6-tuples. Matrices are transformed directly, 1488 3-vectors the same. 1489 1490 6-arrays are interpreted as 6 independent components of a moment 1491 tensor, so they are brought into symmetric 3x3 matrix form. This 1492 is transformed, and the 6 standard components 11,22,33,12,13,23 1493 are returned. 1494 """ 1495 if (not np.prod(np.shape(ma_ve_ar)) in [3, 6, 9]) or \ 1496 (not len(np.shape(ma_ve_ar)) in [1, 2]): 1497 print('\n wrong input - ', end=' ') 1498 print('provide either 3x3 matrix or 3-element vector \n') 1499 raise MTError(' !! ') 1500 1501 if np.prod(np.shape(ma_ve_ar)) == 9: 1502 return np.dot(basis_change_matrix, 1503 np.dot(ma_ve_ar, basis_change_matrix.T)) 1504 elif np.prod(np.shape(ma_ve_ar)) == 6: 1505 m_in = ma_ve_ar 1506 orig_matrix = np.matrix([[m_in[0], m_in[3], m_in[4]], 1507 [m_in[3], m_in[1], m_in[5]], 1508 [m_in[4], m_in[5], m_in[2]]], dtype=np.float) 1509 m_out_mat = np.dot(basis_change_matrix, 1510 np.dot(orig_matrix, basis_change_matrix.T)) 1511 1512 return m_out_mat[0, 0], m_out_mat[1, 1], m_out_mat[2, 2], \ 1513 m_out_mat[0, 1], m_out_mat[0, 2], m_out_mat[1, 2] 1514 else: 1515 if np.shape(ma_ve_ar)[0] == 1: 1516 return np.dot(basis_change_matrix, ma_ve_ar.transpose()) 1517 else: 1518 return np.dot(basis_change_matrix, ma_ve_ar) 1519 1520 1521def USE2NED(some_matrix_or_vector): 1522 """ 1523 Function for basis transform from basis USE to NED. 1524 1525 Input: 1526 3x3 matrix or 3-element vector or 6-element array in USE basis 1527 representation 1528 1529 Output: 1530 3x3 matrix or 3-element vector or 6-element array in NED basis 1531 representation 1532 """ 1533 basis_change_matrix = np.matrix([[0., -1., 0.], 1534 [0., 0., 1.], 1535 [-1., 0., 0.]], dtype=np.float) 1536 return _return_matrix_vector_array(some_matrix_or_vector, 1537 basis_change_matrix) 1538 1539 1540def XYZ2NED(some_matrix_or_vector): 1541 """ 1542 Function for basis transform from basis XYZ to NED. 1543 1544 Input: 1545 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1546 representation 1547 1548 Output: 1549 3x3 matrix or 3-element vector or 6-element array in NED basis 1550 representation 1551 """ 1552 basis_change_matrix = np.matrix([[0., 1., 0.], 1553 [1., 0., 0.], 1554 [0., 0., -1.]], dtype=np.float) 1555 return _return_matrix_vector_array(some_matrix_or_vector, 1556 basis_change_matrix) 1557 1558 1559def NWU2NED(some_matrix_or_vector): 1560 """ 1561 Function for basis transform from basis NWU to NED. 1562 1563 Input: 1564 3x3 matrix or 3-element vector or 6-element array in NWU basis 1565 representation 1566 1567 Output: 1568 3x3 matrix or 3-element vector or 6-element array in NED basis 1569 representation 1570 """ 1571 basis_change_matrix = np.matrix([[1., 0., 0.], 1572 [0., -1., 0.], 1573 [0., 0., -1.]], dtype=np.float) 1574 return _return_matrix_vector_array(some_matrix_or_vector, 1575 basis_change_matrix) 1576 1577 1578def NED2USE(some_matrix_or_vector): 1579 """ 1580 Function for basis transform from basis NED to USE. 1581 1582 Input: 1583 3x3 matrix or 3-element vector or 6-element array in NED basis 1584 representation 1585 1586 Output: 1587 3x3 matrix or 3-element vector or 6-element array in USE basis 1588 representation 1589 """ 1590 basis_change_matrix = np.matrix([[0., -1., 0.], 1591 [0., 0., 1.], 1592 [-1., 0., 0.]], dtype=np.float).I 1593 return _return_matrix_vector_array(some_matrix_or_vector, 1594 basis_change_matrix) 1595 1596 1597def XYZ2USE(some_matrix_or_vector): 1598 """ 1599 Function for basis transform from basis XYZ to USE. 1600 1601 Input: 1602 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1603 representation 1604 1605 Output: 1606 3x3 matrix or 3-element vector or 6-element array in USE basis 1607 representation 1608 """ 1609 basis_change_matrix = np.matrix([[0., 0., 1.], 1610 [0., -1., 0.], 1611 [1., 0., 0.]], dtype=np.float) 1612 return _return_matrix_vector_array(some_matrix_or_vector, 1613 basis_change_matrix) 1614 1615 1616def NED2XYZ(some_matrix_or_vector): 1617 """ 1618 Function for basis transform from basis NED to XYZ. 1619 1620 Input: 1621 3x3 matrix or 3-element vector or 6-element array in NED basis 1622 representation 1623 1624 Output: 1625 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1626 representation 1627 """ 1628 basis_change_matrix = np.matrix([[0., 1., 0.], 1629 [1., 0., 0.], 1630 [0., 0., -1.]], dtype=np.float).I 1631 return _return_matrix_vector_array(some_matrix_or_vector, 1632 basis_change_matrix) 1633 1634 1635def NED2NWU(some_matrix_or_vector): 1636 """ 1637 Function for basis transform from basis NED to NWU. 1638 1639 Input: 1640 3x3 matrix or 3-element vector or 6-element array in NED basis 1641 representation 1642 1643 Output: 1644 3x3 matrix or 3-element vector or 6-element array in NWU basis 1645 representation 1646 """ 1647 basis_change_matrix = np.matrix([[1., 0., 0.], 1648 [0., -1., 0.], 1649 [0., 0., -1.]], dtype=np.float).I 1650 return _return_matrix_vector_array(some_matrix_or_vector, 1651 basis_change_matrix) 1652 1653 1654def USE2XYZ(some_matrix_or_vector): 1655 1656 """ 1657 Function for basis transform from basis USE to XYZ. 1658 1659 Input: 1660 3x3 matrix or 3-element vector or 6-element array in USE basis 1661 representation 1662 1663 Output: 1664 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1665 representation 1666 """ 1667 basis_change_matrix = np.matrix([[0., 0., 1.], 1668 [0., -1., 0.], 1669 [1., 0., 0.]], dtype=np.float).I 1670 return _return_matrix_vector_array(some_matrix_or_vector, 1671 basis_change_matrix) 1672 1673 1674def NWU2XYZ(some_matrix_or_vector): 1675 1676 """ 1677 Function for basis transform from basis USE to XYZ. 1678 1679 Input: 1680 3x3 matrix or 3-element vector or 6-element array in USE basis 1681 representation 1682 1683 Output: 1684 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1685 representation 1686 """ 1687 basis_change_matrix = np.matrix([[0., -1., 0.], 1688 [1., 0., 0.], 1689 [0., 0., 1.]], dtype=np.float) 1690 return _return_matrix_vector_array(some_matrix_or_vector, 1691 basis_change_matrix) 1692 1693 1694def NWU2USE(some_matrix_or_vector): 1695 1696 """ 1697 Function for basis transform from basis USE to XYZ. 1698 1699 Input: 1700 3x3 matrix or 3-element vector or 6-element array in USE basis 1701 representation 1702 1703 Output: 1704 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1705 representation 1706 """ 1707 basis_change_matrix = np.matrix([[0., 0., 1.], 1708 [-1., 0., 0.], 1709 [0., -1., 0.]], dtype=np.float) 1710 return _return_matrix_vector_array(some_matrix_or_vector, 1711 basis_change_matrix) 1712 1713 1714def XYZ2NWU(some_matrix_or_vector): 1715 """ 1716 Function for basis transform from basis USE to XYZ. 1717 1718 Input: 1719 3x3 matrix or 3-element vector or 6-element array in USE basis 1720 representation 1721 1722 Output: 1723 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1724 representation 1725 """ 1726 basis_change_matrix = np.matrix([[0., -1., 0.], 1727 [1., 0., 0.], 1728 [0., 0., 1.]], dtype=np.float).I 1729 return _return_matrix_vector_array(some_matrix_or_vector, 1730 basis_change_matrix) 1731 1732 1733def USE2NWU(some_matrix_or_vector): 1734 """ 1735 Function for basis transform from basis USE to XYZ. 1736 1737 Input: 1738 3x3 matrix or 3-element vector or 6-element array in USE basis 1739 representation 1740 1741 Output: 1742 3x3 matrix or 3-element vector or 6-element array in XYZ basis 1743 representation 1744 """ 1745 basis_change_matrix = np.matrix([[0., 0., 1.], 1746 [-1., 0., 0.], 1747 [0., -1., 0.]], dtype=np.float).I 1748 return _return_matrix_vector_array(some_matrix_or_vector, 1749 basis_change_matrix) 1750 1751 1752def strikediprake_2_moments(strike, dip, rake): 1753 """ 1754 angles are defined as in Jost&Herman (given in degrees) 1755 1756 strike 1757 angle clockwise between north and plane ( in [0,360] ) 1758 dip 1759 angle between surface and dipping plane ( in [0,90] ), 0 = horizontal, 1760 90 = vertical 1761 rake 1762 angle on the rupture plane between strike vector and actual movement 1763 (defined mathematically positive: ccw rotation is positive) 1764 1765 basis for output is NED (= X,Y,Z) 1766 1767 output: 1768 M = M_nn, M_ee, M_dd, M_ne, M_nd, M_ed 1769 """ 1770 S_rad = strike / rad2deg 1771 D_rad = dip / rad2deg 1772 R_rad = rake / rad2deg 1773 1774 for ang in S_rad, D_rad, R_rad: 1775 if abs(ang) < epsilon: 1776 ang = 0. 1777 1778 M1 = -(np.sin(D_rad) * np.cos(R_rad) * np.sin(2 * S_rad) + 1779 np.sin(2 * D_rad) * np.sin(R_rad) * np.sin(S_rad) ** 2) 1780 M2 = (np.sin(D_rad) * np.cos(R_rad) * np.sin(2 * S_rad) - 1781 np.sin(2 * D_rad) * np.sin(R_rad) * np.cos(S_rad) ** 2) 1782 M3 = (np.sin(2 * D_rad) * np.sin(R_rad)) 1783 M4 = (np.sin(D_rad) * np.cos(R_rad) * np.cos(2 * S_rad) + 1784 0.5 * np.sin(2 * D_rad) * np.sin(R_rad) * np.sin(2 * S_rad)) 1785 M5 = -(np.cos(D_rad) * np.cos(R_rad) * np.cos(S_rad) + 1786 np.cos(2 * D_rad) * np.sin(R_rad) * np.sin(S_rad)) 1787 M6 = -(np.cos(D_rad) * np.cos(R_rad) * np.sin(S_rad) - 1788 np.cos(2 * D_rad) * np.sin(R_rad) * np.cos(S_rad)) 1789 1790 Moments = [M1, M2, M3, M4, M5, M6] 1791 1792 return tuple(Moments) 1793 1794 1795def fancy_matrix(m_in): 1796 """ 1797 Returns a given 3x3 matrix or array in a cute way on the shell, if you 1798 use 'print' on the return value. 1799 """ 1800 m = m_in.copy() 1801 1802 norm_factor = round(np.max(np.abs(m.flatten())), 5) 1803 1804 try: 1805 if (norm_factor < 0.1) or (norm_factor >= 10): 1806 if not abs(norm_factor) == 0: 1807 m = m / norm_factor 1808 out = "\n / %5.2F %5.2F %5.2F \\\n" % \ 1809 (m[0, 0], m[0, 1], m[0, 2]) 1810 out += " | %5.2F %5.2F %5.2F | x %F\n" % \ 1811 (m[1, 0], m[1, 1], m[1, 2], norm_factor) 1812 out += " \\ %5.2F %5.2F %5.2F /\n" % \ 1813 (m[2, 0], m[2, 1], m[2, 2]) 1814 return out 1815 except Exception: 1816 pass 1817 1818 return "\n / %5.2F %5.2F %5.2F \\\n" % (m[0, 0], m[0, 1], m[0, 2]) + \ 1819 " | %5.2F %5.2F %5.2F |\n" % (m[1, 0], m[1, 1], m[1, 2]) + \ 1820 " \\ %5.2F %5.2F %5.2F /\n" % (m[2, 0], m[2, 1], m[2, 2]) 1821 1822 1823def fancy_vector(v): 1824 """ 1825 Returns a given 3-vector or array in a cute way on the shell, if you 1826 use 'print' on the return value. 1827 """ 1828 return "\n / %5.2F \\\n" % (v[0]) + \ 1829 " | %5.2F |\n" % (v[1]) + \ 1830 " \\ %5.2F /\n" % (v[2]) 1831 1832 1833# --------------------------------------------------------------- 1834# 1835# Class for plotting: 1836# 1837# --------------------------------------------------------------- 1838 1839class BeachBall: 1840 """ 1841 Class for generating a beachball projection for a provided moment tensor 1842 object. 1843 1844 Input: a MomentTensor object 1845 1846 Output can be plots of 1847 - the eigensystem 1848 - the complete sphere 1849 - the projection to a unit sphere 1850 .. either lower (standard) or upper half 1851 1852 Beside the plots, the unit sphere projection may be saved in a given file. 1853 1854 Alternatively, only the file can be provided without showing anything 1855 directly. 1856 """ 1857 def __init__(self, MT=MomentTensor, kwargs_dict={}, npoints=360): 1858 self.MT = MT 1859 self._M = MT._M 1860 self._set_standard_attributes() 1861 self._update_attributes(kwargs_dict) 1862 self._plot_n_points = npoints 1863 self._nodallines_in_NED_system() 1864 self.arange_1 = np.arange(3 * npoints) - 1 1865 # self._identify_faultplanes() 1866 1867 def ploBB(self, kwargs, ax=None): 1868 """ 1869 Plots the projection of the beachball onto a unit sphere. 1870 """ 1871 self._update_attributes(kwargs) 1872 self._setup_BB() 1873 self._plot_US(ax=ax) 1874 1875 def save_BB(self, kwargs): 1876 """ 1877 Saves the 2D projection of the beachball without plotting. 1878 1879 :param outfile: name of outfile, addressing w.r.t. current directory 1880 :param format: if no implicit valid format is provided within the 1881 filename, add file format 1882 """ 1883 self._update_attributes(kwargs) 1884 self._setup_BB() 1885 self._just_save_bb() 1886 1887 def _just_save_bb(self): 1888 """ 1889 Saves the beachball unit sphere plot into a given file. 1890 """ 1891 import matplotlib 1892 1893 if self._plot_outfile_format == 'svg': 1894 try: 1895 matplotlib.use('SVG') 1896 except Exception: 1897 matplotlib.use('Agg') 1898 elif self._plot_outfile_format == 'pdf': 1899 try: 1900 matplotlib.use('PDF') 1901 except Exception: 1902 matplotlib.use('Agg') 1903 pass 1904 elif self._plot_outfile_format == 'ps': 1905 try: 1906 matplotlib.use('PS') 1907 except Exception: 1908 matplotlib.use('Agg') 1909 pass 1910 elif self._plot_outfile_format == 'eps': 1911 try: 1912 matplotlib.use('Agg') 1913 except Exception: 1914 matplotlib.use('PS') 1915 pass 1916 elif self._plot_outfile_format == 'png': 1917 try: 1918 matplotlib.use('AGG') 1919 except Exception: 1920 mp_out = matplotlib.use('GTKCairo') 1921 if mp_out: 1922 mp_out2 = matplotlib.use('Cairo') 1923 if mp_out2: 1924 matplotlib.use('GDK') 1925 1926 import matplotlib.pyplot as plt 1927 1928 plotfig = self._setup_plot_US(plt) 1929 1930 outfile_format = self._plot_outfile_format 1931 outfile_name = self._plot_outfile 1932 1933 outfile_abs_name = os.path.realpath( 1934 os.path.abspath(os.path.join(os.curdir, outfile_name))) 1935 1936 try: 1937 plotfig.savefig(outfile_abs_name, dpi=self._plot_dpi, 1938 transparent=True, facecolor='k', 1939 format=outfile_format) 1940 except Exception: 1941 print('ERROR!! -- Saving of plot not possible') 1942 return 1943 plt.close(667) 1944 del plt 1945 1946 def get_psxy(self, kwargs): 1947 """ 1948 Returns one string, to be piped into psxy of GMT. 1949 1950 :param GMT_type: fill/lines/EVs (select type of string), 1951 default is 'fill' 1952 :param GMT_scaling: scale the beachball - default radius is 1.0 1953 :param GMT_tension_colour: tension area of BB - colour flag for -Z in 1954 psxy, default is 1 1955 :param GMT_pressure_colour: pressure area of BB - colour flag for -Z in 1956 psxy, default is 0 1957 :param GMT_show_2FPs: flag, if both faultplanes are to be shown, 1958 default is 0 1959 :param GMT_show_1FP: flag, if one faultplane is to be shown, default 1960 is 1 1961 :param GMT_FP_index: 1 or 2, default is 2 1962 """ 1963 self._GMT_type = 'fill' 1964 self._GMT_2fps = False 1965 self._GMT_1fp = 0 1966 1967 self._GMT_psxy_fill = None 1968 self._GMT_psxy_nodals = None 1969 self._GMT_psxy_EVs = None 1970 self._GMT_scaling = 1. 1971 1972 self._GMT_tension_colour = 1 1973 self._GMT_pressure_colour = 0 1974 1975 self._update_attributes(kwargs) 1976 1977 self._setup_BB() 1978 1979 self._set_GMT_attributes() 1980 1981 if self._GMT_type == 'fill': 1982 self._GMT_psxy_fill.seek(0) 1983 GMT_string = self._GMT_psxy_fill.getvalue() 1984 elif self._GMT_type == 'lines': 1985 self._GMT_psxy_nodals.seek(0) 1986 GMT_string = self._GMT_psxy_nodals.getvalue() 1987 else: 1988 GMT_string = self._GMT_psxy_EVs.getvalue() 1989 1990 return GMT_string 1991 1992 def _add_2_GMT_string(self, FH_string, curve, colour): 1993 """ 1994 Writes coordinate pair list of given curve as string into temporal 1995 file handler. 1996 """ 1997 colour_Z = colour 1998 wstring = bytes('> -Z%i\n' % (colour_Z), encoding='utf-8') 1999 FH_string.write(wstring) 2000 np.savetxt(FH_string, self._GMT_scaling * curve.transpose()) 2001 2002 def _set_GMT_attributes(self): 2003 """ 2004 Set the beachball lines and nodals as strings into a file handler. 2005 """ 2006 neg_nodalline = self._nodalline_negative_final_US 2007 pos_nodalline = self._nodalline_positive_final_US 2008 FP1_2_plot = self._FP1_final_US 2009 FP2_2_plot = self._FP2_final_US 2010 EV_2_plot = self._all_EV_2D_US[:, :2].transpose() 2011 US = self._unit_sphere 2012 2013 tension_colour = self._GMT_tension_colour 2014 pressure_colour = self._GMT_pressure_colour 2015 2016 # build strings for possible GMT-output, used by 'psxy' 2017 GMT_string_FH = io.BytesIO() 2018 GMT_linestring_FH = io.BytesIO() 2019 GMT_EVs_FH = io.BytesIO() 2020 2021 self._add_2_GMT_string(GMT_EVs_FH, EV_2_plot, tension_colour) 2022 GMT_EVs_FH.flush() 2023 2024 if self._plot_clr_order > 0: 2025 self._add_2_GMT_string(GMT_string_FH, US, pressure_colour) 2026 self._add_2_GMT_string(GMT_string_FH, neg_nodalline, 2027 tension_colour) 2028 self._add_2_GMT_string(GMT_string_FH, pos_nodalline, 2029 tension_colour) 2030 GMT_string_FH.flush() 2031 2032 if self._plot_curve_in_curve != 0: 2033 self._add_2_GMT_string(GMT_string_FH, US, tension_colour) 2034 2035 if self._plot_curve_in_curve < 1: 2036 self._add_2_GMT_string(GMT_string_FH, neg_nodalline, 2037 pressure_colour) 2038 self._add_2_GMT_string(GMT_string_FH, pos_nodalline, 2039 tension_colour) 2040 GMT_string_FH.flush() 2041 else: 2042 self._add_2_GMT_string(GMT_string_FH, pos_nodalline, 2043 pressure_colour) 2044 self._add_2_GMT_string(GMT_string_FH, neg_nodalline, 2045 tension_colour) 2046 GMT_string_FH.flush() 2047 else: 2048 self._add_2_GMT_string(GMT_string_FH, US, tension_colour) 2049 self._add_2_GMT_string(GMT_string_FH, neg_nodalline, 2050 pressure_colour) 2051 self._add_2_GMT_string(GMT_string_FH, pos_nodalline, 2052 pressure_colour) 2053 GMT_string_FH.flush() 2054 2055 if self._plot_curve_in_curve != 0: 2056 self._add_2_GMT_string(GMT_string_FH, US, pressure_colour) 2057 if self._plot_curve_in_curve < 1: 2058 self._add_2_GMT_string(GMT_string_FH, neg_nodalline, 2059 tension_colour) 2060 self._add_2_GMT_string(GMT_string_FH, pos_nodalline, 2061 pressure_colour) 2062 GMT_string_FH.flush() 2063 else: 2064 self._add_2_GMT_string(GMT_string_FH, pos_nodalline, 2065 tension_colour) 2066 self._add_2_GMT_string(GMT_string_FH, neg_nodalline, 2067 pressure_colour) 2068 2069 GMT_string_FH.flush() 2070 2071 # set all nodallines and faultplanes for plotting: 2072 self._add_2_GMT_string(GMT_linestring_FH, neg_nodalline, 2073 tension_colour) 2074 self._add_2_GMT_string(GMT_linestring_FH, pos_nodalline, 2075 tension_colour) 2076 2077 if self._GMT_2fps: 2078 self._add_2_GMT_string(GMT_linestring_FH, FP1_2_plot, 2079 tension_colour) 2080 self._add_2_GMT_string(GMT_linestring_FH, FP2_2_plot, 2081 tension_colour) 2082 2083 elif self._GMT_1fp: 2084 if not int(self._GMT_1fp) in [1, 2]: 2085 print('no fault plane specified for being plotted...continue', 2086 end=' ') 2087 print('without fault plane(s)') 2088 pass 2089 else: 2090 if int(self._GMT_1fp) == 1: 2091 self._add_2_GMT_string(GMT_linestring_FH, FP1_2_plot, 2092 tension_colour) 2093 else: 2094 self._add_2_GMT_string(GMT_linestring_FH, FP2_2_plot, 2095 tension_colour) 2096 2097 self._add_2_GMT_string(GMT_linestring_FH, US, tension_colour) 2098 2099 GMT_linestring_FH.flush() 2100 2101 setattr(self, '_GMT_psxy_nodals', GMT_linestring_FH) 2102 setattr(self, '_GMT_psxy_fill', GMT_string_FH) 2103 setattr(self, '_GMT_psxy_EVs', GMT_EVs_FH) 2104 2105 def get_MT(self): 2106 """ 2107 Returns the original moment tensor object, handed over to the class at 2108 generating this instance. 2109 """ 2110 return self.MT 2111 2112 def full_sphere_plot(self, kwargs): 2113 """ 2114 Plot of the full beachball, projected on a circle with a radius 2. 2115 """ 2116 self._update_attributes(kwargs) 2117 self._setup_BB() 2118 self._aux_plot() 2119 2120 def _aux_plot(self): 2121 """ 2122 Generates the final plot of the total sphere (according to the chosen 2123 2D-projection. 2124 """ 2125 import matplotlib.pyplot as plt 2126 2127 plt.close('all') 2128 plotfig = plt.figure(665, figsize=(self._plot_aux_plot_size, 2129 self._plot_aux_plot_size)) 2130 2131 plotfig.subplots_adjust(left=0, bottom=0, right=1, top=1) 2132 ax = plotfig.add_subplot(111, aspect='equal') 2133 # plt.axis([-1.1, 1.1, -1.1, 1.1], 'equal') 2134 ax.axison = False 2135 2136 EV_2_plot = getattr(self, '_all_EV' + '_final') 2137 BV_2_plot = getattr(self, '_all_BV' + '_final').transpose() 2138 curve_pos_2_plot = getattr(self, '_nodalline_positive' + '_final') 2139 curve_neg_2_plot = getattr(self, '_nodalline_negative' + '_final') 2140 FP1_2_plot = getattr(self, '_FP1' + '_final') 2141 FP2_2_plot = getattr(self, '_FP2' + '_final') 2142 2143 tension_colour = self._plot_tension_colour 2144 pressure_colour = self._plot_pressure_colour 2145 2146 if self._plot_clr_order > 0: 2147 if self._plot_fill_flag: 2148 2149 alpha = self._plot_fill_alpha * self._plot_total_alpha 2150 ax.fill(self._outer_circle[0, :], self._outer_circle[1, :], 2151 fc=pressure_colour, alpha=alpha) 2152 ax.fill(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 2153 fc=tension_colour, alpha=alpha) 2154 ax.fill(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 2155 fc=tension_colour, alpha=alpha) 2156 2157 if self._plot_curve_in_curve != 0: 2158 ax.fill(self._outer_circle[0, :], self._outer_circle[1, :], 2159 fc=tension_colour, alpha=alpha) 2160 if self._plot_curve_in_curve < 1: 2161 ax.fill(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 2162 fc=pressure_colour, alpha=alpha) 2163 ax.fill(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 2164 fc=tension_colour, alpha=alpha) 2165 else: 2166 ax.fill(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 2167 fc=pressure_colour, alpha=alpha) 2168 ax.fill(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 2169 fc=tension_colour, alpha=alpha) 2170 2171 if self._plot_show_princ_axes: 2172 alpha = self._plot_princ_axes_alpha * self._plot_total_alpha 2173 ax.plot([EV_2_plot[0, 0]], [EV_2_plot[1, 0]], 'm^', 2174 ms=self._plot_princ_axes_symsize, 2175 lw=self._plot_princ_axes_lw, alpha=alpha) 2176 ax.plot([EV_2_plot[0, 3]], [EV_2_plot[1, 3]], 'mv', 2177 ms=self._plot_princ_axes_symsize, 2178 lw=self._plot_princ_axes_lw, alpha=alpha) 2179 ax.plot([EV_2_plot[0, 1]], [EV_2_plot[1, 1]], 'b^', 2180 ms=self._plot_princ_axes_symsize, 2181 lw=self._plot_princ_axes_lw, alpha=alpha) 2182 ax.plot([EV_2_plot[0, 4]], [EV_2_plot[1, 4]], 'bv', 2183 ms=self._plot_princ_axes_symsize, 2184 lw=self._plot_princ_axes_lw, alpha=alpha) 2185 ax.plot([EV_2_plot[0, 2]], [EV_2_plot[1, 2]], 'g^', 2186 ms=self._plot_princ_axes_symsize, 2187 lw=self._plot_princ_axes_lw, alpha=alpha) 2188 ax.plot([EV_2_plot[0, 5]], [EV_2_plot[1, 5]], 'gv', 2189 ms=self._plot_princ_axes_symsize, 2190 lw=self._plot_princ_axes_lw, alpha=alpha) 2191 else: 2192 if self._plot_fill_flag: 2193 alpha = self._plot_fill_alpha * self._plot_total_alpha 2194 ax.fill(self._outer_circle[0, :], self._outer_circle[1, :], 2195 fc=tension_colour, alpha=alpha) 2196 ax.fill(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 2197 fc=pressure_colour, alpha=alpha) 2198 ax.fill(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 2199 fc=pressure_colour, alpha=alpha) 2200 2201 if self._plot_curve_in_curve != 0: 2202 ax.fill(self._outer_circle[0, :], self._outer_circle[1, :], 2203 fc=pressure_colour, alpha=alpha) 2204 if self._plot_curve_in_curve < 0: 2205 ax.fill(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 2206 fc=tension_colour, alpha=alpha) 2207 ax.fill(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 2208 fc=pressure_colour, alpha=alpha) 2209 pass 2210 else: 2211 ax.fill(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 2212 fc=tension_colour, alpha=alpha) 2213 ax.fill(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 2214 fc=pressure_colour, alpha=alpha) 2215 pass 2216 2217 if self._plot_show_princ_axes: 2218 alpha = self._plot_princ_axes_alpha * self._plot_total_alpha 2219 ax.plot([EV_2_plot[0, 0]], [EV_2_plot[1, 0]], 'g^', 2220 ms=self._plot_princ_axes_symsize, 2221 lw=self._plot_princ_axes_lw, alpha=alpha) 2222 ax.plot([EV_2_plot[0, 3]], [EV_2_plot[1, 3]], 'gv', 2223 ms=self._plot_princ_axes_symsize, 2224 lw=self._plot_princ_axes_lw, alpha=alpha) 2225 ax.plot([EV_2_plot[0, 1]], [EV_2_plot[1, 1]], 'b^', 2226 ms=self._plot_princ_axes_symsize, 2227 lw=self._plot_princ_axes_lw, alpha=alpha) 2228 ax.plot([EV_2_plot[0, 4]], [EV_2_plot[1, 4]], 'bv', 2229 ms=self._plot_princ_axes_symsize, 2230 lw=self._plot_princ_axes_lw, alpha=alpha) 2231 ax.plot([EV_2_plot[0, 2]], [EV_2_plot[1, 2]], 'm^', 2232 ms=self._plot_princ_axes_symsize, 2233 lw=self._plot_princ_axes_lw, alpha=alpha) 2234 ax.plot([EV_2_plot[0, 5]], [EV_2_plot[1, 5]], 'mv', 2235 ms=self._plot_princ_axes_symsize, 2236 lw=self._plot_princ_axes_lw, alpha=alpha) 2237 2238 self._plot_nodalline_colour = 'y' 2239 2240 ax.plot(curve_neg_2_plot[0, :], curve_neg_2_plot[1, :], 'o', 2241 c=self._plot_nodalline_colour, lw=self._plot_nodalline_width, 2242 alpha=self._plot_nodalline_alpha * self._plot_total_alpha, 2243 ms=3) 2244 2245 self._plot_nodalline_colour = 'b' 2246 2247 ax.plot(curve_pos_2_plot[0, :], curve_pos_2_plot[1, :], 'D', 2248 c=self._plot_nodalline_colour, lw=self._plot_nodalline_width, 2249 alpha=self._plot_nodalline_alpha * self._plot_total_alpha, 2250 ms=3) 2251 2252 if self._plot_show_1faultplane: 2253 if self._plot_show_FP_index == 1: 2254 ax.plot(FP1_2_plot[0, :], FP1_2_plot[1, :], '+', 2255 c=self._plot_faultplane_colour, 2256 lw=self._plot_faultplane_width, 2257 alpha=self._plot_faultplane_alpha * 2258 self._plot_total_alpha, ms=5) 2259 elif self._plot_show_FP_index == 2: 2260 ax.plot(FP2_2_plot[0, :], FP2_2_plot[1, :], '+', 2261 c=self._plot_faultplane_colour, 2262 lw=self._plot_faultplane_width, 2263 alpha=self._plot_faultplane_alpha * 2264 self._plot_total_alpha, ms=5) 2265 2266 elif self._plot_show_faultplanes: 2267 ax.plot(FP1_2_plot[0, :], FP1_2_plot[1, :], '+', 2268 c=self._plot_faultplane_colour, 2269 lw=self._plot_faultplane_width, 2270 alpha=self._plot_faultplane_alpha * self._plot_total_alpha, 2271 ms=4) 2272 ax.plot(FP2_2_plot[0, :], FP2_2_plot[1, :], '+', 2273 c=self._plot_faultplane_colour, 2274 lw=self._plot_faultplane_width, 2275 alpha=self._plot_faultplane_alpha * self._plot_total_alpha, 2276 ms=4) 2277 else: 2278 pass 2279 2280 # if isotropic part shall be displayed, fill the circle completely with 2281 # the appropriate colour 2282 if self._pure_isotropic: 2283 if abs(np.trace(self._M)) > epsilon: 2284 if self._plot_clr_order < 0: 2285 ax.fill(self._outer_circle[0, :], self._outer_circle[1, :], 2286 fc=tension_colour, alpha=1, zorder=100) 2287 else: 2288 ax.fill(self._outer_circle[0, :], self._outer_circle[1, :], 2289 fc=pressure_colour, alpha=1, zorder=100) 2290 2291 # plot NED basis vectors 2292 if self._plot_show_basis_axes: 2293 plot_size_in_points = self._plot_size * 2.54 * 72 2294 points_per_unit = plot_size_in_points / 2. 2295 2296 fontsize = plot_size_in_points / 66. 2297 symsize = plot_size_in_points / 77. 2298 2299 direction_letters = 'NSEWDU' 2300 for idx, val in enumerate(BV_2_plot): 2301 x_coord = val[0] 2302 y_coord = val[1] 2303 np_letter = direction_letters[idx] 2304 2305 rot_angle = -np.arctan2(y_coord, x_coord) + pi / 2. 2306 original_rho = np.sqrt(x_coord ** 2 + y_coord ** 2) 2307 2308 marker_x = (original_rho - (3 * symsize / points_per_unit)) * \ 2309 np.sin(rot_angle) 2310 marker_y = (original_rho - (3 * symsize / points_per_unit)) * \ 2311 np.cos(rot_angle) 2312 annot_x = (original_rho - (8.5 * fontsize / points_per_unit)) \ 2313 * np.sin(rot_angle) 2314 annot_y = (original_rho - (8.5 * fontsize / points_per_unit)) \ 2315 * np.cos(rot_angle) 2316 2317 ax.text(annot_x, annot_y, np_letter, 2318 horizontalalignment='center', size=fontsize, 2319 weight='bold', verticalalignment='center', 2320 bbox=dict(edgecolor='white', facecolor='white', 2321 alpha=1)) 2322 2323 if original_rho > epsilon: 2324 ax.scatter([marker_x], [marker_y], 2325 marker=(3, 0, rot_angle), s=symsize ** 2, c='k', 2326 facecolor='k', zorder=300) 2327 else: 2328 ax.scatter([x_coord], [y_coord], marker=(4, 1, rot_angle), 2329 s=symsize ** 2, c='k', facecolor='k', 2330 zorder=300) 2331 2332 # plot both circle lines (radius 1 and 2) 2333 ax.plot(self._unit_sphere[0, :], self._unit_sphere[1, :], 2334 c=self._plot_outerline_colour, lw=self._plot_outerline_width, 2335 alpha=self._plot_outerline_alpha * self._plot_total_alpha) 2336 ax.plot(self._outer_circle[0, :], self._outer_circle[1, :], 2337 c=self._plot_outerline_colour, lw=self._plot_outerline_width, 2338 alpha=self._plot_outerline_alpha * self._plot_total_alpha) 2339 2340 # dummy points for setting plot plot size more accurately 2341 ax.plot([0, 2.1, 0, -2.1], [2.1, 0, -2.1, 0], ',', alpha=0.) 2342 2343 ax.autoscale_view(tight=True, scalex=True, scaley=True) 2344 2345 if self._plot_save_plot: 2346 try: 2347 plotfig.savefig(self._plot_outfile + '.' + 2348 self._plot_outfile_format, dpi=self._plot_dpi, 2349 transparent=True, 2350 format=self._plot_outfile_format) 2351 except Exception: 2352 print('saving of plot not possible') 2353 2354 plt.show() 2355 2356 def pa_plot(self, kwargs): 2357 """ 2358 Plot of the solution in the principal axes system. 2359 """ 2360 import matplotlib.pyplot as plt 2361 2362 self._update_attributes(kwargs) 2363 2364 r_hor = self._r_hor_for_pa_plot 2365 r_hor_FP = self._r_hor_FP_for_pa_plot 2366 2367 plt.rc('grid', color='#316931', linewidth=0.5, linestyle='-.') 2368 plt.rc('xtick', labelsize=12) 2369 plt.rc('ytick', labelsize=10) 2370 2371 width, height = plt.rcParams['figure.figsize'] 2372 size = min(width, height) 2373 2374 fig = plt.figure(34, figsize=(size, size)) 2375 # TODO remove once minimum required matplotlib version reaches 2.0 2376 # see matplotlib/matplotlib#5501 2377 if MATPLOTLIB_VERSION < [2, 0]: 2378 axis_facecolor_kwargs = dict(axisbg='#d5de9c') 2379 else: 2380 axis_facecolor_kwargs = dict(facecolor='#d5de9c') 2381 ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], polar=True, 2382 **axis_facecolor_kwargs) 2383 2384 r_steps = [0.000001] 2385 for i in (np.arange(4) + 1) * 0.2: 2386 r_steps.append(i) 2387 r_labels = ['S'] 2388 for ii in range(len(r_steps)): 2389 if (ii + 1) % 2 == 0: 2390 r_labels.append(str(r_steps[ii])) 2391 else: 2392 r_labels.append(' ') 2393 2394 t_angles = np.arange(0., 360., 90) 2395 t_labels = [' N ', ' H ', ' - N', ' - H'] 2396 2397 plt.thetagrids(t_angles, labels=t_labels) 2398 2399 ax.plot(self._phi_curve, r_hor, color='r', lw=3) 2400 ax.plot(self._phi_curve, r_hor_FP, color='b', lw=1.5) 2401 ax.set_rmax(1.0) 2402 plt.grid(True) 2403 2404 plt.rgrids((r_steps), labels=r_labels) 2405 2406 ax.set_title("beachball in eigenvector system", fontsize=15) 2407 2408 if self._plot_save_plot: 2409 try: 2410 fig.savefig(self._plot_outfile + '.' + 2411 self._plot_outfile_format, dpi=self._plot_dpi, 2412 transparent=True, 2413 format=self._plot_outfile_format) 2414 except Exception: 2415 print('saving of plot not possible') 2416 plt.show() 2417 2418 def _set_standard_attributes(self): 2419 """ 2420 Sets default values of mandatory arguments. 2421 """ 2422 # plot basis system and view point: 2423 self._plot_basis = 'NED' 2424 self._plot_projection = 'stereo' 2425 self._plot_viewpoint = [0., 0., 0.] 2426 self._plot_basis_change = None 2427 2428 # flag, if upper hemisphere is seen instead 2429 self._plot_show_upper_hemis = False 2430 2431 # flag, if isotropic part shall be considered 2432 self._plot_isotropic_part = False 2433 self._pure_isotropic = False 2434 2435 # number of minimum points per line and full circle (number/360 is 2436 # minimum of points per degree at rounded lines) 2437 self._plot_n_points = 360 2438 2439 # nodal line of pressure and tension regimes: 2440 self._plot_nodalline_width = 2 2441 self._plot_nodalline_colour = 'k' 2442 self._plot_nodalline_alpha = 1. 2443 2444 # outer circle line 2445 self._plot_outerline_width = 2 2446 self._plot_outerline_colour = 'k' 2447 self._plot_outerline_alpha = 1. 2448 2449 # faultplane(s) 2450 self._plot_faultplane_width = 4 2451 self._plot_faultplane_colour = 'b' 2452 self._plot_faultplane_alpha = 1. 2453 2454 self._plot_show_faultplanes = False 2455 self._plot_show_1faultplane = False 2456 self._plot_show_FP_index = 1 2457 2458 # principal axes: 2459 self._plot_show_princ_axes = False 2460 self._plot_princ_axes_symsize = 10 2461 self._plot_princ_axes_lw = 3 2462 self._plot_princ_axes_alpha = 0.5 2463 2464 # NED basis: 2465 self._plot_show_basis_axes = False 2466 2467 # filling of the area: 2468 self._plot_clr_order = self.MT.get_colour_order() 2469 self._plot_curve_in_curve = 0 2470 self._plot_fill_flag = True 2471 self._plot_tension_colour = 'r' 2472 self._plot_pressure_colour = 'w' 2473 self._plot_fill_alpha = 1. 2474 2475 # general plot options 2476 self._plot_size = 5 2477 self._plot_aux_plot_size = 5 2478 self._plot_dpi = 200 2479 2480 self._plot_total_alpha = 1. 2481 2482 # possibility to add external data (e.g. measured polarizations) 2483 self._plot_external_data = False 2484 self._external_data = None 2485 2486 # if, howto, whereto save the plot 2487 self._plot_save_plot = False 2488 self._plot_outfile = './BB_plot_example' 2489 self._plot_outfile_format = 'svg' 2490 2491 def _update_attributes(self, kwargs): 2492 """ 2493 Makes an internal update of the object's attributes with the 2494 provided list of keyword arguments. 2495 2496 If the keyword (extended by a leading _ ) is in the dict of 2497 the object, the value is updated. Otherwise, the keyword is 2498 ignored. 2499 """ 2500 for key in kwargs.keys(): 2501 if key[0] == '_': 2502 kw = key[1:] 2503 else: 2504 kw = key 2505 if '_' + kw in dir(self): 2506 setattr(self, '_' + kw, kwargs[key]) 2507 if kwargs.get('plot_only_lines', False): 2508 setattr(self, '_plot_fill_flag', False) 2509 2510 def _setup_BB(self, unit_circle=True): 2511 """ 2512 Setup of the beachball, when a plotting method is evoked. 2513 2514 Contains all the technical stuff for generating the final view of the 2515 beachball: 2516 2517 - Finding a rotation matrix, describing the given viewpoint onto the 2518 beachball projection 2519 - Rotating all elements (lines, points) w.r.t. the given viewpoint 2520 - Projecting the 3D sphere into the 2D plane 2521 - Building circle lines in radius r=1 and r=2 2522 - Correct the order of line points, yielding a consecutive set of 2523 points for drawing lines 2524 - Smoothing of all curves, avoiding nasty sectioning connection lines 2525 - Checking, if the two nodalline curves are laying completely within 2526 each other ( cahnges plotting order of overlay plot construction) 2527 - Projection of final smooth solution onto the standard unit sphere 2528 """ 2529 self._find_basis_change_2_new_viewpoint() 2530 self._rotate_all_objects_2_new_view() 2531 self._vertical_2D_projection() 2532 2533 if unit_circle: 2534 self._build_circles() 2535 2536 if not self.MT._iso_percentage == 100: 2537 self._correct_curves() 2538 self._smooth_curves() 2539 self._check_curve_in_curve() 2540 2541 self._projection_2_unit_sphere() 2542 2543 if self.MT._iso_percentage == 100: 2544 if np.trace(self.MT.get_M()) < 0: 2545 self._plot_clr_order = 1 2546 else: 2547 self._plot_clr_order = -1 2548 2549 def _correct_curves(self): 2550 """ 2551 Correcting potentially wrong curves. 2552 2553 Checks, if the order of the given coordinates of the lines must be 2554 re-arranged, allowing for an automatical line plotting. 2555 """ 2556 list_of_curves_2_correct = ['nodalline_negative', 'nodalline_positive', 2557 'FP1', 'FP2'] 2558 n_curve_points = self._plot_n_points 2559 2560 for obj in list_of_curves_2_correct: 2561 obj2cor_name = '_' + obj + '_2D' 2562 obj2cor = getattr(self, obj2cor_name) 2563 2564 obj2cor_in_right_order = self._sort_curve_points(obj2cor) 2565 2566 # logger.debug( 'curve: ', str(obj)) 2567 # check, if curve closed !!!!!! 2568 start_r = np.sqrt(obj2cor_in_right_order[0, 0] ** 2 + 2569 obj2cor_in_right_order[1, 0] ** 2) 2570 r_last_point = np.sqrt(obj2cor_in_right_order[0, -1] ** 2 + 2571 obj2cor_in_right_order[1, -1] ** 2) 2572 dist_last_first_point = \ 2573 np.sqrt((obj2cor_in_right_order[0, -1] - 2574 obj2cor_in_right_order[0, 0]) ** 2 + 2575 (obj2cor_in_right_order[1, -1] - 2576 obj2cor_in_right_order[1, 0]) ** 2) 2577 2578 # check, if distance between last and first point is smaller than 2579 # the distance between last point and the edge (at radius=2) 2580 if dist_last_first_point > (2 - r_last_point): 2581 # add points on edge to polygon, if it is an open curve 2582 # logger.debug( str(obj)+' not closed - closing over edge... ') 2583 phi_end = np.arctan2(obj2cor_in_right_order[0, -1], 2584 obj2cor_in_right_order[1, -1]) % (2 * pi) 2585 R_end = r_last_point 2586 phi_start = np.arctan2(obj2cor_in_right_order[0, 0], 2587 obj2cor_in_right_order[1, 0]) % (2 * pi) 2588 R_start = start_r 2589 2590 # add one point on the edge every fraction of degree given by 2591 # input parameter, increase the radius linearly 2592 phi_end_larger = np.sign(phi_end - phi_start) 2593 angle_smaller_pi = np.sign(pi - np.abs(phi_end - phi_start)) 2594 2595 if phi_end_larger * angle_smaller_pi > 0: 2596 go_ccw = True 2597 openangle = (phi_end - phi_start) % (2 * pi) 2598 else: 2599 go_ccw = False 2600 openangle = (phi_start - phi_end) % (2 * pi) 2601 2602 radius_interval = R_start - R_end # closing from end to start 2603 2604 n_edgepoints = int(openangle * rad2deg * 2605 n_curve_points / 360.) - 1 2606 # logger.debug( 'open angle %.2f degrees - filling with %i 2607 # points on the edge\n'%(openangle/pi*180,n_edgepoints)) 2608 if go_ccw: 2609 obj2cor_in_right_order = \ 2610 list(obj2cor_in_right_order.transpose()) 2611 for kk in range(n_edgepoints + 1): 2612 current_phi = phi_end - kk * openangle / \ 2613 (n_edgepoints + 1) 2614 current_radius = R_end + kk * radius_interval / \ 2615 (n_edgepoints + 1) 2616 temp = [current_radius * math.sin(current_phi), 2617 current_radius * np.cos(current_phi)] 2618 obj2cor_in_right_order.append(temp) 2619 obj2cor_in_right_order = \ 2620 np.array(obj2cor_in_right_order).transpose() 2621 else: 2622 obj2cor_in_right_order = \ 2623 list(obj2cor_in_right_order.transpose()) 2624 for kk in range(n_edgepoints + 1): 2625 current_phi = phi_end + kk * openangle / \ 2626 (n_edgepoints + 1) 2627 current_radius = R_end + kk * radius_interval / \ 2628 (n_edgepoints + 1) 2629 temp = [current_radius * math.sin(current_phi), 2630 current_radius * np.cos(current_phi)] 2631 obj2cor_in_right_order.append(temp) 2632 obj2cor_in_right_order = \ 2633 np.array(obj2cor_in_right_order).transpose() 2634 setattr(self, '_' + obj + '_in_order', obj2cor_in_right_order) 2635 return 1 2636 2637 def _nodallines_in_NED_system(self): 2638 """ 2639 The two nodal lines between the areas on a beachball are given by the 2640 points, where tan²(alpha) = (-EWs/(EWN*cos(phi)**2 + EWh*sin(phi)**2)) 2641 is fulfilled. 2642 2643 This solution is gained in the principal axes system and then expressed 2644 in terms of the NED basis system 2645 2646 output: 2647 - set of points, building the first nodal line, coordinates in the 2648 input basis system (standard NED) 2649 - set of points, building the second nodal line, coordinates in the 2650 input basis system (standard NED) 2651 - array with 6 points, describing positive and negative part of 3 2652 principal axes 2653 - array with partition of full circle (angle values in degrees) 2654 fraction is given by parameter n_curve_points 2655 """ 2656 # build the nodallines of positive/negative areas in the principal axes 2657 # system 2658 n_curve_points = self._plot_n_points 2659 2660 # phi is the angle between neutral axis and horizontal projection 2661 # of the curve point to the surface, spanned by H- and 2662 # N-axis. Running mathematically negative (clockwise) around the 2663 # SIGMA-axis. Stepsize is given by the parameter for number of 2664 # curve points 2665 phi = (np.arange(n_curve_points) / float(n_curve_points) + 2666 1. / n_curve_points) * 2 * pi 2667 self._phi_curve = phi 2668 2669 # analytical/geometrical solution for separatrix curve - alpha is 2670 # opening angle between principal axis SIGMA and point of curve. (alpha 2671 # is 0, if curve lies directly on the SIGMA axis) 2672 2673 # CASE: including isotropic part 2674 # sigma axis flips, if EWn flips sign 2675 2676 EWh_devi = self.MT.get_eigvals()[0] - 1. / 3 * np.trace(self._M) 2677 EWn_devi = self.MT.get_eigvals()[1] - 1. / 3 * np.trace(self._M) 2678 EWs_devi = self.MT.get_eigvals()[2] - 1. / 3 * np.trace(self._M) 2679 2680 if not self._plot_isotropic_part: 2681 EWh = EWh_devi 2682 EWn = EWn_devi 2683 EWs = EWs_devi 2684 else: 2685 EWh_tmp = self.MT.get_eigvals()[0] 2686 EWn_tmp = self.MT.get_eigvals()[1] 2687 EWs_tmp = self.MT.get_eigvals()[2] 2688 2689 trace_m = np.sum(self.MT.get_eigvals()) 2690 EWh = EWh_tmp.copy() 2691 EWs = EWs_tmp.copy() 2692 2693 if trace_m != 0: 2694 if (self._plot_clr_order > 0 and EWn_tmp >= 0 and 2695 abs(EWs_tmp) > abs(EWh_tmp)) or \ 2696 (self._plot_clr_order < 0 and 2697 EWn_tmp <= 0 and abs(EWs_tmp) > abs(EWh_tmp)): 2698 2699 EWs = EWh_tmp.copy() 2700 EWh = EWs_tmp.copy() 2701 print('changed order!!\n') 2702 EVs_tmp = self.MT._rotation_matrix[:, 2].copy() 2703 EVh_tmp = self.MT._rotation_matrix[:, 0].copy() 2704 2705 self.MT._rotation_matrix[:, 0] = EVs_tmp 2706 self.MT._rotation_matrix[:, 2] = EVh_tmp 2707 self._plot_clr_order *= -1 2708 2709 EWn = EWn_tmp.copy() 2710 2711 if abs(EWn) < epsilon: 2712 EWn = 0 2713 norm_factor = max(np.abs([EWh, EWn, EWs])) 2714 2715 # norm_factor is be zero in some cases 2716 with np.errstate(all='ignore'): 2717 [EWh, EWn, EWs] = [xx / norm_factor for xx in [EWh, EWn, EWs]] 2718 2719 RHS = -EWs / (EWn * np.cos(phi) ** 2 + EWh * np.sin(phi) ** 2) 2720 2721 if np.all([np.sign(xx) >= 0 for xx in RHS]): 2722 alpha = np.arctan(np.sqrt(RHS)) * rad2deg 2723 else: 2724 alpha = phi.copy() 2725 alpha[:] = 90 2726 self._pure_isotropic = 1 2727 2728 # fault planes: 2729 RHS_FP = 1. / (np.sin(phi) ** 2) 2730 alpha_FP = np.arctan(np.sqrt(RHS_FP)) * rad2deg 2731 2732 # horizontal coordinates of curves 2733 r_hor = np.sin(alpha / rad2deg) 2734 r_hor_FP = np.sin(alpha_FP / rad2deg) 2735 2736 self._r_hor_for_pa_plot = r_hor 2737 self._r_hor_FP_for_pa_plot = r_hor_FP 2738 2739 H_values = np.sin(phi) * r_hor 2740 N_values = np.cos(phi) * r_hor 2741 H_values_FP = np.sin(phi) * r_hor_FP 2742 N_values_FP = np.cos(phi) * r_hor_FP 2743 2744 # set vertical value of curve point coordinates - two symmetric curves 2745 # exist 2746 S_values_positive = np.cos(alpha / rad2deg) 2747 S_values_negative = -np.cos(alpha / rad2deg) 2748 S_values_positive_FP = np.cos(alpha_FP / rad2deg) 2749 S_values_negative_FP = -np.cos(alpha_FP / rad2deg) 2750 2751 ############# 2752 # change basis back to original input reference system 2753 ######### 2754 2755 chng_basis = self.MT._rotation_matrix 2756 2757 line_tuple_pos = np.zeros((3, n_curve_points)) 2758 line_tuple_neg = np.zeros((3, n_curve_points)) 2759 2760 for ii in range(n_curve_points): 2761 pos_vec_in_EV_basis = np.array([H_values[ii], N_values[ii], 2762 S_values_positive[ii]]).transpose() 2763 neg_vec_in_EV_basis = np.array([H_values[ii], N_values[ii], 2764 S_values_negative[ii]]).transpose() 2765 line_tuple_pos[:, ii] = np.dot(chng_basis, pos_vec_in_EV_basis) 2766 line_tuple_neg[:, ii] = np.dot(chng_basis, neg_vec_in_EV_basis) 2767 2768 EVh = self.MT.get_eigvecs()[0] 2769 EVn = self.MT.get_eigvecs()[1] 2770 EVs = self.MT.get_eigvecs()[2] 2771 2772 all_EV = np.zeros((3, 6)) 2773 2774 all_EV[:, 0] = EVh.transpose() 2775 all_EV[:, 1] = EVn.transpose() 2776 all_EV[:, 2] = EVs.transpose() 2777 all_EV[:, 3] = -EVh.transpose() 2778 all_EV[:, 4] = -EVn.transpose() 2779 all_EV[:, 5] = -EVs.transpose() 2780 2781 # basis vectors: 2782 all_BV = np.zeros((3, 6)) 2783 all_BV[:, 0] = np.array((1, 0, 0)) 2784 all_BV[:, 1] = np.array((-1, 0, 0)) 2785 all_BV[:, 2] = np.array((0, 1, 0)) 2786 all_BV[:, 3] = np.array((0, -1, 0)) 2787 all_BV[:, 4] = np.array((0, 0, 1)) 2788 all_BV[:, 5] = np.array((0, 0, -1)) 2789 2790 # re-sort the two 90 degree nodal lines to 2 fault planes - cut each at 2791 # halves and merge pairs 2792 # additionally change basis system to NED reference system 2793 2794 midpoint_idx = int(n_curve_points / 2.) 2795 2796 FP1 = np.zeros((3, n_curve_points)) 2797 FP2 = np.zeros((3, n_curve_points)) 2798 2799 for ii in range(midpoint_idx): 2800 FP1_vec = np.array([H_values_FP[ii], N_values_FP[ii], 2801 S_values_positive_FP[ii]]).transpose() 2802 FP2_vec = np.array([H_values_FP[ii], N_values_FP[ii], 2803 S_values_negative_FP[ii]]).transpose() 2804 FP1[:, ii] = np.dot(chng_basis, FP1_vec) 2805 FP2[:, ii] = np.dot(chng_basis, FP2_vec) 2806 2807 for jj in range(midpoint_idx): 2808 ii = n_curve_points - jj - 1 2809 2810 FP1_vec = np.array([H_values_FP[ii], N_values_FP[ii], 2811 S_values_negative_FP[ii]]).transpose() 2812 FP2_vec = np.array([H_values_FP[ii], N_values_FP[ii], 2813 S_values_positive_FP[ii]]).transpose() 2814 FP1[:, ii] = np.dot(chng_basis, FP1_vec) 2815 FP2[:, ii] = np.dot(chng_basis, FP2_vec) 2816 2817 # identify with faultplane index, gotten from 'get_fps': 2818 self._FP1 = FP1 2819 self._FP2 = FP2 2820 2821 self._all_EV = all_EV 2822 self._all_BV = all_BV 2823 self._nodalline_negative = line_tuple_neg 2824 self._nodalline_positive = line_tuple_pos 2825 2826 def _identify_faultplanes(self): 2827 """ 2828 See, if the 2 faultplanes, given as attribute of the moment 2829 tensor object, handed to this instance, are consistent with 2830 the faultplane lines, obtained from the basis solution. If 2831 not, interchange the indices of the newly found ones. 2832 """ 2833 # TODO !!!!!! 2834 pass 2835 2836 def _find_basis_change_2_new_viewpoint(self): 2837 """ 2838 Finding the Eulerian angles, if you want to rotate an object. 2839 2840 Your original view point is the position (0,0,0). Input are the 2841 coordinates of the new point of view, equivalent to geographical 2842 coordinates. 2843 2844 Example: 2845 2846 Original view onto the Earth is from right above lat=0, lon=0 with 2847 north=upper edge, south=lower edge. Now you want to see the Earth 2848 from a position somewhere near Baku. So lat=45, 2849 lon=45, azimuth=0. 2850 2851 The Earth must be rotated around some axis, not to be determined. 2852 The rotation matrixx is the matrix for the change of basis to the 2853 new local orthonormal system. 2854 2855 input: 2856 - latitude in degrees from -90 (south) to 90 (north) 2857 - longitude in degrees from -180 (west) to 180 (east) 2858 - azimuth in degrees from 0 (heading north) to 360 (north again) 2859 """ 2860 new_latitude = self._plot_viewpoint[0] 2861 new_longitude = self._plot_viewpoint[1] 2862 new_azimuth = self._plot_viewpoint[2] 2863 2864 s_lat = np.sin(new_latitude / rad2deg) 2865 if abs(s_lat) < epsilon: 2866 s_lat = 0 2867 c_lat = np.cos(new_latitude / rad2deg) 2868 if abs(c_lat) < epsilon: 2869 c_lat = 0 2870 s_lon = np.sin(new_longitude / rad2deg) 2871 if abs(s_lon) < epsilon: 2872 s_lon = 0 2873 c_lon = np.cos(new_longitude / rad2deg) 2874 if abs(c_lon) < epsilon: 2875 c_lon = 0 2876 # assume input basis as NED!!! 2877 2878 # original point of view therein is (0,0,-1) 2879 # new point at lat=latitude, lon=longitude, az=0, given in old 2880 # NED-coordinates: 2881 # (cos(latitude), sin(latitude)*sin(longitude), 2882 # sin(latitude)*cos(longitude) ) 2883 # 2884 # new " down' " is given by the negative position vector, so pointing 2885 # inwards to the centre point 2886 # down_prime = - ( np.array( ( s_lat, c_lat*c_lon, -c_lat*s_lon ) ) ) 2887 down_prime = -(np.array((s_lat, c_lat * s_lon, -c_lat * c_lon))) 2888 2889 # normalise: 2890 down_prime /= np.sqrt(np.dot(down_prime, down_prime)) 2891 2892 # get second local basis vector " north' " by orthogonalising 2893 # (Gram-Schmidt method) the original north w.r.t. the new " down' " 2894 north_prime_not_normalised = np.array((1., 0., 0.)) - \ 2895 (np.dot(down_prime, np.array((1., 0., 0.))) / 2896 (np.dot(down_prime, down_prime)) * down_prime) 2897 2898 len_north_prime_not_normalised = \ 2899 np.sqrt(np.dot(north_prime_not_normalised, 2900 north_prime_not_normalised)) 2901 # check for poles: 2902 if np.abs(len_north_prime_not_normalised) < epsilon: 2903 # case: north pole 2904 if s_lat > 0: 2905 north_prime = np.array((0., 0., 1.)) 2906 # case: south pole 2907 else: 2908 north_prime = np.array((0., 0., -1.)) 2909 else: 2910 north_prime = \ 2911 north_prime_not_normalised / len_north_prime_not_normalised 2912 2913 # third basis vector is obtained by a cross product of the first two 2914 east_prime = np.cross(down_prime, north_prime) 2915 2916 # normalise: 2917 east_prime /= np.sqrt(np.dot(east_prime, east_prime)) 2918 2919 rotmat_pos_raw = np.zeros((3, 3)) 2920 rotmat_pos_raw[:, 0] = north_prime 2921 rotmat_pos_raw[:, 1] = east_prime 2922 rotmat_pos_raw[:, 2] = down_prime 2923 2924 rotmat_pos = np.asmatrix(rotmat_pos_raw).T 2925 # this matrix gives the coordinates of a given point in the old 2926 # coordinates w.r.t. the new system 2927 2928 # up to here, only the position has changed, the angle of view 2929 # (azimuth) has to be added by an additional rotation around the 2930 # down'-axis (in the frame of the new coordinates) 2931 # set up the local rotation around the new down'-axis by the given 2932 # angle 'azimuth'. Positive values turn view counterclockwise from the 2933 # new north' 2934 only_rotation = np.zeros((3, 3)) 2935 s_az = np.sin(new_azimuth / rad2deg) 2936 if abs(s_az) < epsilon: 2937 s_az = 0. 2938 c_az = np.cos(new_azimuth / rad2deg) 2939 if abs(c_az) < epsilon: 2940 c_az = 0. 2941 2942 only_rotation[2, 2] = 1 2943 only_rotation[0, 0] = c_az 2944 only_rotation[1, 1] = c_az 2945 only_rotation[0, 1] = -s_az 2946 only_rotation[1, 0] = s_az 2947 2948 local_rotation = np.asmatrix(only_rotation) 2949 2950 # apply rotation from left!! 2951 total_rotation_matrix = np.dot(local_rotation, rotmat_pos) 2952 2953 # yields the complete matrix for representing the old coordinates in 2954 # the new (rotated) frame: 2955 self._plot_basis_change = total_rotation_matrix 2956 2957 def _rotate_all_objects_2_new_view(self): 2958 """ 2959 Rotate all relevant parts of the solution - namely the 2960 eigenvector-projections, the 2 nodallines, and the faultplanes 2961 - so that they are seen from the new viewpoint. 2962 """ 2963 objects_2_rotate = ['all_EV', 'all_BV', 'nodalline_negative', 2964 'nodalline_positive', 'FP1', 'FP2'] 2965 2966 for obj in objects_2_rotate: 2967 object2rotate = getattr(self, '_' + obj).transpose() 2968 2969 rotated_thing = object2rotate.copy() 2970 for i in range(len(object2rotate)): 2971 rotated_thing[i] = np.dot(self._plot_basis_change, 2972 object2rotate[i]) 2973 2974 rotated_object = rotated_thing.copy() 2975 setattr(self, '_' + obj + '_rotated', rotated_object.transpose()) 2976 2977 # --------------------------------------------------------------- 2978 2979 def _vertical_2D_projection(self): 2980 """ 2981 Start the vertical projection of the 3D beachball onto the 2D plane. 2982 2983 The projection is chosen according to the attribute '_plot_projection' 2984 """ 2985 list_of_possible_projections = ['stereo', 'ortho', 'lambert', 'gnom'] 2986 2987 if self._plot_projection not in list_of_possible_projections: 2988 print('desired projection not possible - choose from:\n ', end=' ') 2989 print(list_of_possible_projections) 2990 raise MTError(' !! ') 2991 2992 if self._plot_projection == 'stereo': 2993 if not self._stereo_vertical(): 2994 print('ERROR in stereo_vertical') 2995 raise MTError(' !! ') 2996 elif self._plot_projection == 'ortho': 2997 if not self._orthographic_vertical(): 2998 print('ERROR in stereo_vertical') 2999 raise MTError(' !! ') 3000 elif self._plot_projection == 'lambert': 3001 if not self._lambert_vertical(): 3002 print('ERROR in stereo_vertical') 3003 raise MTError(' !! ') 3004 elif self._plot_projection == 'gnom': 3005 if not self._gnomonic_vertical(): 3006 print('ERROR in stereo_vertical') 3007 raise MTError(' !! ') 3008 3009 def _stereo_vertical(self): 3010 """ 3011 Stereographic/azimuthal conformal 2D projection onto a plane, tangent 3012 to the lowest point (0,0,1). 3013 3014 Keeps the angles constant! 3015 3016 The parts in the lower hemisphere are projected to the unit 3017 sphere, the upper half to an annular region between radii r=1 3018 and r=2. If the attribute '_show_upper_hemis' is set, the 3019 projection is reversed. 3020 """ 3021 objects_2_project = ['all_EV', 'all_BV', 'nodalline_negative', 3022 'nodalline_positive', 'FP1', 'FP2'] 3023 3024 available_coord_systems = ['NED'] 3025 3026 if self._plot_basis not in available_coord_systems: 3027 print('desired plotting projection not possible - choose from :\n', 3028 end=' ') 3029 print(available_coord_systems) 3030 raise MTError(' !! ') 3031 3032 plot_upper_hem = self._plot_show_upper_hemis 3033 3034 for obj in objects_2_project: 3035 obj_name = '_' + obj + '_rotated' 3036 o2proj = getattr(self, obj_name) 3037 coords = o2proj.copy() 3038 3039 n_points = len(o2proj[0, :]) 3040 stereo_coords = np.zeros((2, n_points)) 3041 3042 for ll in range(n_points): 3043 # second component is EAST 3044 co_x = coords[1, ll] 3045 # first component is NORTH 3046 co_y = coords[0, ll] 3047 # z given in DOWN 3048 co_z = -coords[2, ll] 3049 3050 rho_hor = np.sqrt(co_x ** 2 + co_y ** 2) 3051 3052 if rho_hor == 0: 3053 new_y = 0 3054 new_x = 0 3055 if plot_upper_hem: 3056 if co_z < 0: 3057 new_x = 2 3058 else: 3059 if co_z > 0: 3060 new_x = 2 3061 else: 3062 if co_z < 0: 3063 new_rho = rho_hor / (1. - co_z) 3064 if plot_upper_hem: 3065 new_rho = 2 - (rho_hor / (1. - co_z)) 3066 3067 new_x = co_x / rho_hor * new_rho 3068 new_y = co_y / rho_hor * new_rho 3069 else: 3070 new_rho = 2 - (rho_hor / (1. + co_z)) 3071 if plot_upper_hem: 3072 new_rho = rho_hor / (1. + co_z) 3073 3074 new_x = co_x / rho_hor * new_rho 3075 new_y = co_y / rho_hor * new_rho 3076 3077 stereo_coords[0, ll] = new_x 3078 stereo_coords[1, ll] = new_y 3079 3080 setattr(self, '_' + obj + '_2D', stereo_coords) 3081 setattr(self, '_' + obj + '_final', stereo_coords) 3082 3083 return 1 3084 3085 def _orthographic_vertical(self): 3086 """ 3087 Orthographic 2D projection onto a plane, tangent to the lowest 3088 point (0,0,1). 3089 3090 Shows the natural view on a 2D sphere from large distances (assuming 3091 parallel projection) 3092 3093 The parts in the lower hemisphere are projected to the unit 3094 sphere, the upper half to an annular region between radii r=1 3095 and r=2. If the attribute '_show_upper_hemis' is set, the 3096 projection is reversed. 3097 """ 3098 3099 objects_2_project = ['all_EV', 'all_BV', 'nodalline_negative', 3100 'nodalline_positive', 'FP1', 'FP2'] 3101 3102 available_coord_systems = ['NED'] 3103 3104 if self._plot_basis not in available_coord_systems: 3105 print('desired plotting projection not possible - choose from :\n', 3106 end=' ') 3107 print(available_coord_systems) 3108 raise MTError(' !! ') 3109 3110 plot_upper_hem = self._plot_show_upper_hemis 3111 3112 for obj in objects_2_project: 3113 obj_name = '_' + obj + '_rotated' 3114 o2proj = getattr(self, obj_name) 3115 coords = o2proj.copy() 3116 3117 n_points = len(o2proj[0, :]) 3118 coords2D = np.zeros((2, n_points)) 3119 3120 for ll in range(n_points): 3121 # second component is EAST 3122 co_x = coords[1, ll] 3123 # first component is NORTH 3124 co_y = coords[0, ll] 3125 # z given in DOWN 3126 co_z = -coords[2, ll] 3127 3128 rho_hor = np.sqrt(co_x ** 2 + co_y ** 2) 3129 3130 if rho_hor == 0: 3131 new_y = 0 3132 new_x = 0 3133 if plot_upper_hem: 3134 if co_z < 0: 3135 new_x = 2 3136 else: 3137 if co_z > 0: 3138 new_x = 2 3139 else: 3140 if co_z < 0: 3141 new_rho = rho_hor 3142 if plot_upper_hem: 3143 new_rho = 2 - rho_hor 3144 3145 new_x = co_x / rho_hor * new_rho 3146 new_y = co_y / rho_hor * new_rho 3147 else: 3148 new_rho = 2 - rho_hor 3149 if plot_upper_hem: 3150 new_rho = rho_hor 3151 3152 new_x = co_x / rho_hor * new_rho 3153 new_y = co_y / rho_hor * new_rho 3154 3155 coords2D[0, ll] = new_x 3156 coords2D[1, ll] = new_y 3157 3158 setattr(self, '_' + obj + '_2D', coords2D) 3159 setattr(self, '_' + obj + '_final', coords2D) 3160 3161 return 1 3162 3163 def _lambert_vertical(self): 3164 """ 3165 Lambert azimuthal equal-area 2D projection onto a plane, tangent to the 3166 lowest point (0,0,1). 3167 3168 Keeps the area constant! 3169 3170 The parts in the lower hemisphere are projected to the unit 3171 sphere (only here the area is kept constant), the upper half to an 3172 annular region between radii r=1 and r=2. If the attribute 3173 '_show_upper_hemis' is set, the projection is reversed. 3174 """ 3175 objects_2_project = ['all_EV', 'all_BV', 'nodalline_negative', 3176 'nodalline_positive', 'FP1', 'FP2'] 3177 3178 available_coord_systems = ['NED'] 3179 3180 if self._plot_basis not in available_coord_systems: 3181 print('desired plotting projection not possible - choose from :\n', 3182 end=' ') 3183 print(available_coord_systems) 3184 raise MTError(' !! ') 3185 3186 plot_upper_hem = self._plot_show_upper_hemis 3187 3188 for obj in objects_2_project: 3189 obj_name = '_' + obj + '_rotated' 3190 o2proj = getattr(self, obj_name) 3191 coords = o2proj.copy() 3192 3193 n_points = len(o2proj[0, :]) 3194 coords2D = np.zeros((2, n_points)) 3195 3196 for ll in range(n_points): 3197 # second component is EAST 3198 co_x = coords[1, ll] 3199 # first component is NORTH 3200 co_y = coords[0, ll] 3201 # z given in DOWN 3202 co_z = -coords[2, ll] 3203 3204 rho_hor = np.sqrt(co_x ** 2 + co_y ** 2) 3205 3206 if rho_hor == 0: 3207 new_y = 0 3208 new_x = 0 3209 if plot_upper_hem: 3210 if co_z < 0: 3211 new_x = 2 3212 else: 3213 if co_z > 0: 3214 new_x = 2 3215 else: 3216 if co_z < 0: 3217 new_rho = rho_hor / np.sqrt(1. - co_z) 3218 3219 if plot_upper_hem: 3220 new_rho = 2 - (rho_hor / np.sqrt(1. - co_z)) 3221 3222 new_x = co_x / rho_hor * new_rho 3223 new_y = co_y / rho_hor * new_rho 3224 3225 else: 3226 new_rho = 2 - (rho_hor / np.sqrt(1. + co_z)) 3227 3228 if plot_upper_hem: 3229 new_rho = rho_hor / np.sqrt(1. + co_z) 3230 3231 new_x = co_x / rho_hor * new_rho 3232 new_y = co_y / rho_hor * new_rho 3233 3234 coords2D[0, ll] = new_x 3235 coords2D[1, ll] = new_y 3236 3237 setattr(self, '_' + obj + '_2D', coords2D) 3238 setattr(self, '_' + obj + '_final', coords2D) 3239 3240 return 1 3241 3242 def _gnomonic_vertical(self): 3243 """ 3244 Gnomonic 2D projection onto a plane, tangent to the lowest 3245 point (0,0,1). 3246 3247 Keeps the great circles as straight lines (geodetics constant) ! 3248 3249 The parts in the lower hemisphere are projected to the unit 3250 sphere, the upper half to an annular region between radii r=1 3251 and r=2. If the attribute '_show_upper_hemis' is set, the 3252 projection is reversed. 3253 """ 3254 3255 objects_2_project = ['all_EV', 'all_BV', 'nodalline_negative', 3256 'nodalline_positive', 'FP1', 'FP2'] 3257 3258 available_coord_systems = ['NED'] 3259 3260 if self._plot_basis not in available_coord_systems: 3261 print('desired plotting projection not possible - choose from :\n', 3262 end=' ') 3263 print(available_coord_systems) 3264 raise MTError(' !! ') 3265 3266 plot_upper_hem = self._plot_show_upper_hemis 3267 3268 for obj in objects_2_project: 3269 obj_name = '_' + obj + '_rotated' 3270 o2proj = getattr(self, obj_name) 3271 coords = o2proj.copy() 3272 3273 n_points = len(o2proj[0, :]) 3274 coords2D = np.zeros((2, n_points)) 3275 3276 for ll in range(n_points): 3277 # second component is EAST 3278 co_x = coords[1, ll] 3279 # first component is NORTH 3280 co_y = coords[0, ll] 3281 # z given in DOWN 3282 co_z = -coords[2, ll] 3283 3284 rho_hor = np.sqrt(co_x ** 2 + co_y ** 2) 3285 3286 if rho_hor == 0: 3287 new_y = 0 3288 new_x = 0 3289 if co_z > 0: 3290 new_x = 2 3291 if plot_upper_hem: 3292 new_x = 0 3293 else: 3294 if co_z < 0: 3295 new_rho = np.cos(np.arcsin(rho_hor)) * \ 3296 np.tan(np.arcsin(rho_hor)) 3297 3298 if plot_upper_hem: 3299 new_rho = 2 - (np.cos(np.arcsin(rho_hor)) * 3300 np.tan(np.arcsin(rho_hor))) 3301 3302 new_x = co_x / rho_hor * new_rho 3303 new_y = co_y / rho_hor * new_rho 3304 3305 else: 3306 new_rho = 2 - (np.cos(np.arcsin(rho_hor)) * 3307 np.tan(np.arcsin(rho_hor))) 3308 3309 if plot_upper_hem: 3310 new_rho = np.cos(np.arcsin(rho_hor)) * \ 3311 np.tan(np.arcsin(rho_hor)) 3312 3313 new_x = co_x / rho_hor * new_rho 3314 new_y = co_y / rho_hor * new_rho 3315 3316 coords2D[0, ll] = new_x 3317 coords2D[1, ll] = new_y 3318 3319 setattr(self, '_' + obj + '_2D', coords2D) 3320 setattr(self, '_' + obj + '_final', coords2D) 3321 3322 return 1 3323 3324 def _build_circles(self): 3325 """ 3326 Sets two sets of points, describing the unit sphere and the outer 3327 circle with r=2. 3328 3329 Added as attributes '_unit_sphere' and '_outer_circle'. 3330 """ 3331 phi = self._phi_curve 3332 3333 UnitSphere = np.zeros((2, len(phi))) 3334 UnitSphere[0, :] = np.cos(phi) 3335 UnitSphere[1, :] = np.sin(phi) 3336 3337 # outer circle ( radius for stereographic projection is set to 2 ) 3338 outer_circle_points = 2 * UnitSphere 3339 3340 self._unit_sphere = UnitSphere 3341 self._outer_circle = outer_circle_points 3342 3343 def _sort_curve_points(self, curve): 3344 """ 3345 Checks, if curve points are in right order for line plotting. 3346 3347 If not, a re-arranging is carried out. 3348 """ 3349 sorted_curve = np.zeros((2, len(curve[0, :]))) 3350 # in polar coordinates 3351 r_phi_curve = np.zeros((len(curve[0, :]), 2)) 3352 for ii in range(curve.shape[1]): 3353 r_phi_curve[ii, 0] = \ 3354 math.sqrt(curve[0, ii] ** 2 + curve[1, ii] ** 2) 3355 r_phi_curve[ii, 1] = \ 3356 math.atan2(curve[0, ii], curve[1, ii]) % (2 * pi) 3357 # find index with highest r 3358 largest_r_idx = np.argmax(r_phi_curve[:, 0]) 3359 3360 # check, if perhaps more values with same r - if so, take point with 3361 # lowest phi 3362 other_idces = \ 3363 np.where(r_phi_curve[:, 0] == r_phi_curve[largest_r_idx, 0]) 3364 if len(other_idces) > 1: 3365 best_idx = np.argmin(r_phi_curve[other_idces, 1]) 3366 start_idx_curve = other_idces[best_idx] 3367 else: 3368 start_idx_curve = largest_r_idx 3369 3370 if not start_idx_curve == 0: 3371 pass 3372 3373 # check orientation - want to go inwards 3374 start_r = r_phi_curve[start_idx_curve, 0] 3375 next_idx = (start_idx_curve + 1) % len(r_phi_curve[:, 0]) 3376 prep_idx = (start_idx_curve - 1) % len(r_phi_curve[:, 0]) 3377 next_r = r_phi_curve[next_idx, 0] 3378 3379 keep_direction = True 3380 if next_r <= start_r: 3381 # check, if next R is on other side of area - look at total 3382 # distance - if yes, reverse direction 3383 dist_first_next = \ 3384 (curve[0, next_idx] - curve[0, start_idx_curve]) ** 2 + \ 3385 (curve[1, next_idx] - curve[1, start_idx_curve]) ** 2 3386 dist_first_other = \ 3387 (curve[0, prep_idx] - curve[0, start_idx_curve]) ** 2 + \ 3388 (curve[1, prep_idx] - curve[1, start_idx_curve]) ** 2 3389 3390 if dist_first_next > dist_first_other: 3391 keep_direction = False 3392 3393 if keep_direction: 3394 # direction is kept 3395 for jj in range(curve.shape[1]): 3396 running_idx = (start_idx_curve + jj) % len(curve[0, :]) 3397 sorted_curve[0, jj] = curve[0, running_idx] 3398 sorted_curve[1, jj] = curve[1, running_idx] 3399 else: 3400 # direction is reversed 3401 for jj in range(curve.shape[1]): 3402 running_idx = (start_idx_curve - jj) % len(curve[0, :]) 3403 sorted_curve[0, jj] = curve[0, running_idx] 3404 sorted_curve[1, jj] = curve[1, running_idx] 3405 3406 # check if step of first to second point does not have large angle 3407 # step (problem caused by projection of point (pole) onto whole 3408 # edge - if this first angle step is larger than the one between 3409 # points 2 and three, correct position of first point: keep R, but 3410 # take angle with same difference as point 2 to point 3 3411 3412 angle_point_1 = (math.atan2(sorted_curve[0, 0], 3413 sorted_curve[1, 0]) % (2 * pi)) 3414 angle_point_2 = (math.atan2(sorted_curve[0, 1], 3415 sorted_curve[1, 1]) % (2 * pi)) 3416 angle_point_3 = (math.atan2(sorted_curve[0, 2], 3417 sorted_curve[1, 2]) % (2 * pi)) 3418 3419 angle_diff_23 = (angle_point_3 - angle_point_2) 3420 if angle_diff_23 > pi: 3421 angle_diff_23 = (-angle_diff_23) % (2 * pi) 3422 3423 angle_diff_12 = (angle_point_2 - angle_point_1) 3424 if angle_diff_12 > pi: 3425 angle_diff_12 = (-angle_diff_12) % (2 * pi) 3426 3427 if abs(angle_diff_12) > abs(angle_diff_23): 3428 r_old = \ 3429 math.sqrt(sorted_curve[0, 0] ** 2 + sorted_curve[1, 0] ** 2) 3430 new_angle = (angle_point_2 - angle_diff_23) % (2 * pi) 3431 sorted_curve[0, 0] = r_old * math.sin(new_angle) 3432 sorted_curve[1, 0] = r_old * math.cos(new_angle) 3433 3434 return sorted_curve 3435 3436 def _smooth_curves(self): 3437 """ 3438 Corrects curves for potential large gaps, resulting in strange 3439 intersection lines on nodals of round and irreagularly shaped 3440 areas. 3441 3442 At least one coordinte point on each degree on the circle is assured. 3443 """ 3444 list_of_curves_2_smooth = ['nodalline_negative', 'nodalline_positive', 3445 'FP1', 'FP2'] 3446 3447 points_per_degree = self._plot_n_points / 360. 3448 3449 for curve2smooth in list_of_curves_2_smooth: 3450 obj_name = curve2smooth + '_in_order' 3451 obj = getattr(self, '_' + obj_name).transpose() 3452 3453 smoothed_array = np.zeros((1, 2)) 3454 smoothed_array[0, :] = obj[0] 3455 smoothed_list = [smoothed_array] 3456 3457 # now in shape (n_points,2) 3458 for idx, val in enumerate(obj[:-1]): 3459 r1 = math.sqrt(val[0] ** 2 + val[1] ** 2) 3460 r2 = math.sqrt(obj[idx + 1][0] ** 2 + obj[idx + 1][1] ** 2) 3461 phi1 = math.atan2(val[0], val[1]) 3462 phi2 = math.atan2(obj[idx + 1][0], obj[idx + 1][1]) 3463 3464 phi2_larger = np.sign(phi2 - phi1) 3465 angle_smaller_pi = np.sign(pi - abs(phi2 - phi1)) 3466 3467 if phi2_larger * angle_smaller_pi > 0: 3468 go_cw = True 3469 openangle = (phi2 - phi1) % (2 * pi) 3470 else: 3471 go_cw = False 3472 openangle = (phi1 - phi2) % (2 * pi) 3473 3474 openangle_deg = openangle * rad2deg 3475 radius_diff = r2 - r1 3476 3477 if openangle_deg > 1. / points_per_degree: 3478 3479 n_fillpoints = int(openangle_deg * points_per_degree) 3480 fill_array = np.zeros((n_fillpoints, 2)) 3481 if go_cw: 3482 angles = ((np.arange(n_fillpoints) + 1) * openangle / 3483 (n_fillpoints + 1) + phi1) % (2 * pi) 3484 else: 3485 angles = (phi1 - (np.arange(n_fillpoints) + 1) * 3486 openangle / (n_fillpoints + 1)) % (2 * pi) 3487 3488 radii = (np.arange(n_fillpoints) + 1) * \ 3489 radius_diff / (n_fillpoints + 1) + r1 3490 3491 fill_array[:, 0] = radii * np.sin(angles) 3492 fill_array[:, 1] = radii * np.cos(angles) 3493 3494 smoothed_list.append(fill_array) 3495 3496 smoothed_list.append([obj[idx + 1]]) 3497 3498 smoothed_array = np.vstack(smoothed_list) 3499 setattr(self, '_' + curve2smooth + '_final', 3500 smoothed_array.transpose()) 3501 3502 def _check_curve_in_curve(self): 3503 """ 3504 Checks, if one of the two nodallines contains the other one 3505 completely. If so, the order of colours is re-adapted, 3506 assuring the correct order when doing the overlay plotting. 3507 """ 3508 lo_points_in_pos_curve = \ 3509 list(self._nodalline_positive_final.transpose()) 3510 lo_points_in_pos_curve_array = \ 3511 self._nodalline_positive_final.transpose() 3512 lo_points_in_neg_curve = \ 3513 list(self._nodalline_negative_final.transpose()) 3514 lo_points_in_neg_curve_array = \ 3515 self._nodalline_negative_final.transpose() 3516 3517 # check, if negative curve completely within positive curve 3518 mask_neg_in_pos = 0 3519 for neg_point in lo_points_in_neg_curve: 3520 mask_neg_in_pos += self._pnpoly(lo_points_in_pos_curve_array, 3521 neg_point[:2]) 3522 if mask_neg_in_pos > len(lo_points_in_neg_curve) - 3: 3523 self._plot_curve_in_curve = 1 3524 3525 # check, if positive curve completely within negative curve 3526 mask_pos_in_neg = 0 3527 for pos_point in lo_points_in_pos_curve: 3528 mask_pos_in_neg += self._pnpoly(lo_points_in_neg_curve_array, 3529 pos_point[:2]) 3530 if mask_pos_in_neg > len(lo_points_in_pos_curve) - 3: 3531 self._plot_curve_in_curve = -1 3532 3533 # correct for ONE special case: double couple with its 3534 # eigensystem = NED basis system: 3535 testarray = [1., 0, 0, 0, 1, 0, 0, 0, 1] 3536 if np.prod(self.MT._rotation_matrix.A1 == testarray) and \ 3537 (self.MT._eigenvalues[1] == 0): 3538 self._plot_curve_in_curve = -1 3539 self._plot_clr_order = 1 3540 3541 def _point_inside_polygon(self, x, y, poly): 3542 """ 3543 Determine if a point is inside a given polygon or not. 3544 3545 Polygon is a list of (x,y) pairs. 3546 """ 3547 n = len(poly) 3548 inside = False 3549 3550 p1x, p1y = poly[0] 3551 for i in range(n + 1): 3552 p2x, p2y = poly[i % n] 3553 if y > min(p1y, p2y): 3554 if y <= max(p1y, p2y): 3555 if x <= max(p1x, p2x): 3556 if p1y != p2y: 3557 xinters = \ 3558 (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x 3559 if p1x == p2x or x <= xinters: 3560 inside = not inside 3561 p1x, p1y = p2x, p2y 3562 3563 return inside 3564 3565 def _pnpoly(self, verts, point): 3566 """ 3567 Check whether point is in the polygon defined by verts. 3568 3569 verts - 2xN array 3570 point - (2,) array 3571 3572 See 3573 https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html 3574 """ 3575 # using take instead of getitem, about ten times faster, see 3576 # http://wesmckinney.com/blog/?p=215 3577 verts = np.require(verts, dtype=np.float64) 3578 x, y = point 3579 3580 xpi = verts[:, 0] 3581 ypi = verts[:, 1] 3582 # shift 3583 xpj = xpi.take(self.arange_1[:xpi.size]) 3584 ypj = ypi.take(self.arange_1[:ypi.size]) 3585 3586 possible_crossings = \ 3587 ((ypi <= y) & (y < ypj)) | ((ypj <= y) & (y < ypi)) 3588 3589 xpi = xpi[possible_crossings] 3590 ypi = ypi[possible_crossings] 3591 xpj = xpj[possible_crossings] 3592 ypj = ypj[possible_crossings] 3593 3594 crossings = x < (xpj - xpi) * (y - ypi) / (ypj - ypi) + xpi 3595 3596 return crossings.sum() % 2 3597 3598 def _projection_2_unit_sphere(self): 3599 """ 3600 Brings the complete solution (from stereographic projection) 3601 onto the unit sphere by just shrinking the maximum radius of 3602 all points to 1. 3603 3604 This keeps the area definitions, so the colouring is not affected. 3605 """ 3606 list_of_objects_2_project = ['nodalline_positive_final', 3607 'nodalline_negative_final'] 3608 lo_fps = ['FP1_final', 'FP2_final'] 3609 3610 for obj2proj in list_of_objects_2_project: 3611 obj = getattr(self, '_' + obj2proj).transpose().copy() 3612 for idx, val in enumerate(obj): 3613 old_radius = np.sqrt(val[0] ** 2 + val[1] ** 2) 3614 if old_radius > 1: 3615 obj[idx, 0] = val[0] / old_radius 3616 obj[idx, 1] = val[1] / old_radius 3617 3618 setattr(self, '_' + obj2proj + '_US', obj.transpose()) 3619 3620 for fp in lo_fps: 3621 obj = getattr(self, '_' + fp).transpose().copy() 3622 3623 tmp_obj = [] 3624 for idx, val in enumerate(obj): 3625 old_radius = np.sqrt(val[0] ** 2 + val[1] ** 2) 3626 if old_radius <= 1 + epsilon: 3627 tmp_obj.append(val) 3628 tmp_obj2 = np.array(tmp_obj).transpose() 3629 tmp_obj3 = self._sort_curve_points(tmp_obj2) 3630 3631 setattr(self, '_' + fp + '_US', tmp_obj3) 3632 3633 lo_visible_EV = [] 3634 3635 for idx, val in enumerate(self._all_EV_2D.transpose()): 3636 r_ev = np.sqrt(val[0] ** 2 + val[1] ** 2) 3637 if r_ev <= 1: 3638 lo_visible_EV.append([val[0], val[1], idx]) 3639 visible_EVs = np.array(lo_visible_EV) 3640 3641 self._all_EV_2D_US = visible_EVs 3642 3643 lo_visible_BV = [] 3644 dummy_list1 = [] 3645 direction_letters = 'NSEWDU' 3646 3647 for idx, val in enumerate(self._all_BV_2D.transpose()): 3648 r_bv = math.sqrt(val[0] ** 2 + val[1] ** 2) 3649 if r_bv <= 1: 3650 if idx == 1 and 'N' in dummy_list1: 3651 continue 3652 elif idx == 3 and 'E' in dummy_list1: 3653 continue 3654 elif idx == 5 and 'D' in dummy_list1: 3655 continue 3656 else: 3657 lo_visible_BV.append([val[0], val[1], idx]) 3658 dummy_list1.append(direction_letters[idx]) 3659 3660 visible_BVs = np.array(lo_visible_BV) 3661 3662 self._all_BV_2D_US = visible_BVs 3663 3664 def _plot_US(self, ax=None): 3665 """ 3666 Generates the final plot of the beachball projection on the unit 3667 sphere. 3668 3669 Additionally, the plot can be saved in a file on the fly. 3670 """ 3671 import matplotlib.pyplot as plt 3672 3673 plotfig = self._setup_plot_US(plt, ax=ax) 3674 3675 if self._plot_save_plot: 3676 try: 3677 plotfig.savefig(self._plot_outfile + '.' + 3678 self._plot_outfile_format, dpi=self._plot_dpi, 3679 transparent=True, 3680 format=self._plot_outfile_format) 3681 except Exception: 3682 print('saving of plot not possible') 3683 plt.show() 3684 plt.close('all') 3685 3686 def _setup_plot_US(self, plt, ax=None): 3687 """ 3688 Setting up the figure with the final plot of the unit sphere. 3689 3690 Either called by _plot_US or by _just_save_bb 3691 """ 3692 plt.close(667) 3693 if ax is None: 3694 plotfig = plt.figure(667, 3695 figsize=(self._plot_size, self._plot_size)) 3696 plotfig.subplots_adjust(left=0, bottom=0, right=1, top=1) 3697 ax = plotfig.add_subplot(111, aspect='equal') 3698 3699 ax.axison = False 3700 3701 neg_nodalline = self._nodalline_negative_final_US 3702 pos_nodalline = self._nodalline_positive_final_US 3703 FP1_2_plot = self._FP1_final_US 3704 FP2_2_plot = self._FP2_final_US 3705 3706 US = self._unit_sphere 3707 3708 tension_colour = self._plot_tension_colour 3709 pressure_colour = self._plot_pressure_colour 3710 3711 if self._plot_fill_flag: 3712 if self._plot_clr_order > 0: 3713 alpha = self._plot_fill_alpha * self._plot_total_alpha 3714 3715 ax.fill(US[0, :], US[1, :], fc=pressure_colour, alpha=alpha) 3716 ax.fill(neg_nodalline[0, :], neg_nodalline[1, :], 3717 fc=tension_colour, alpha=alpha) 3718 ax.fill(pos_nodalline[0, :], pos_nodalline[1, :], 3719 fc=tension_colour, alpha=alpha) 3720 3721 if self._plot_curve_in_curve != 0: 3722 ax.fill(US[0, :], US[1, :], fc=tension_colour, 3723 alpha=alpha) 3724 3725 if self._plot_curve_in_curve < 1: 3726 ax.fill(neg_nodalline[0, :], neg_nodalline[1, :], 3727 fc=pressure_colour, alpha=alpha) 3728 ax.fill(pos_nodalline[0, :], pos_nodalline[1, :], 3729 fc=tension_colour, alpha=alpha) 3730 pass 3731 else: 3732 ax.fill(pos_nodalline[0, :], pos_nodalline[1, :], 3733 fc=pressure_colour, alpha=alpha) 3734 ax.fill(neg_nodalline[0, :], neg_nodalline[1, :], 3735 fc=tension_colour, alpha=alpha) 3736 pass 3737 3738 EV_sym = ['m^', 'b^', 'g^', 'mv', 'bv', 'gv'] 3739 3740 if self._plot_show_princ_axes: 3741 alpha = \ 3742 self._plot_princ_axes_alpha * self._plot_total_alpha 3743 3744 for val in self._all_EV_2D_US: 3745 ax.plot([val[0]], [val[1]], EV_sym[int(val[2])], 3746 ms=self._plot_princ_axes_symsize, 3747 lw=self._plot_princ_axes_lw, alpha=alpha) 3748 3749 else: 3750 alpha = self._plot_fill_alpha * self._plot_total_alpha 3751 3752 ax.fill(US[0, :], US[1, :], fc=tension_colour, alpha=alpha) 3753 ax.fill(neg_nodalline[0, :], neg_nodalline[1, :], 3754 fc=pressure_colour, alpha=alpha) 3755 ax.fill(pos_nodalline[0, :], pos_nodalline[1, :], 3756 fc=pressure_colour, alpha=alpha) 3757 3758 if self._plot_curve_in_curve != 0: 3759 ax.fill(US[0, :], US[1, :], fc=pressure_colour, 3760 alpha=alpha) 3761 3762 if self._plot_curve_in_curve < 1: 3763 ax.fill(neg_nodalline[0, :], neg_nodalline[1, :], 3764 fc=tension_colour, alpha=alpha) 3765 ax.fill(pos_nodalline[0, :], pos_nodalline[1, :], 3766 fc=pressure_colour, alpha=alpha) 3767 pass 3768 else: 3769 ax.fill(pos_nodalline[0, :], pos_nodalline[1, :], 3770 fc=tension_colour, alpha=alpha) 3771 ax.fill(neg_nodalline[0, :], neg_nodalline[1, :], 3772 fc=pressure_colour, alpha=alpha) 3773 pass 3774 3775 EV_sym = ['g^', 'b^', 'm^', 'gv', 'bv', 'mv'] 3776 if self._plot_show_princ_axes: 3777 alpha = self._plot_princ_axes_alpha * self._plot_total_alpha 3778 3779 for val in self._all_EV_2D_US: 3780 ax.plot([val[0]], [val[1]], EV_sym[int(val[2])], 3781 ms=self._plot_princ_axes_symsize, 3782 lw=self._plot_princ_axes_lw, alpha=alpha) 3783 3784 # 3785 # set all nodallines and faultplanes for plotting: 3786 # 3787 3788 ax.plot(neg_nodalline[0, :], neg_nodalline[1, :], 3789 c=self._plot_nodalline_colour, ls='-', 3790 lw=self._plot_nodalline_width, 3791 alpha=self._plot_nodalline_alpha * self._plot_total_alpha) 3792 # ax.plot( neg_nodalline[0,:] ,neg_nodalline[1,:],'go') 3793 3794 ax.plot(pos_nodalline[0, :], pos_nodalline[1, :], 3795 c=self._plot_nodalline_colour, ls='-', 3796 lw=self._plot_nodalline_width, 3797 alpha=self._plot_nodalline_alpha * self._plot_total_alpha) 3798 3799 if self._plot_show_faultplanes: 3800 ax.plot(FP1_2_plot[0, :], FP1_2_plot[1, :], 3801 c=self._plot_faultplane_colour, ls='-', 3802 lw=self._plot_faultplane_width, 3803 alpha=self._plot_faultplane_alpha * self._plot_total_alpha) 3804 ax.plot(FP2_2_plot[0, :], FP2_2_plot[1, :], 3805 c=self._plot_faultplane_colour, ls='-', 3806 lw=self._plot_faultplane_width, 3807 alpha=self._plot_faultplane_alpha * self._plot_total_alpha) 3808 3809 elif self._plot_show_1faultplane: 3810 if self._plot_show_FP_index not in [1, 2]: 3811 print('no fault plane specified for being plotted... ', 3812 end=' ') 3813 print('continue without faultplane') 3814 pass 3815 else: 3816 alpha = self._plot_faultplane_alpha * self._plot_total_alpha 3817 if self._plot_show_FP_index == 1: 3818 ax.plot(FP1_2_plot[0, :], FP1_2_plot[1, :], 3819 c=self._plot_faultplane_colour, ls='-', 3820 lw=self._plot_faultplane_width, alpha=alpha) 3821 else: 3822 ax.plot(FP2_2_plot[0, :], FP2_2_plot[1, :], 3823 c=self._plot_faultplane_colour, ls='-', 3824 lw=self._plot_faultplane_width, alpha=alpha) 3825 3826 # if isotropic part shall be displayed, fill the circle completely with 3827 # the appropriate colour 3828 3829 if self._pure_isotropic: 3830 # f abs( np.trace( self._M )) > epsilon: 3831 if self._plot_clr_order < 0: 3832 ax.fill(US[0, :], US[1, :], fc=tension_colour, alpha=1, 3833 zorder=100) 3834 else: 3835 ax.fill(US[0, :], US[1, :], fc=pressure_colour, alpha=1, 3836 zorder=100) 3837 3838 # plot outer circle line of US 3839 ax.plot(US[0, :], US[1, :], c=self._plot_outerline_colour, ls='-', 3840 lw=self._plot_outerline_width, 3841 alpha=self._plot_outerline_alpha * self._plot_total_alpha) 3842 3843 # plot NED basis vectors 3844 if self._plot_show_basis_axes: 3845 plot_size_in_points = self._plot_size * 2.54 * 72 3846 points_per_unit = plot_size_in_points / 2. 3847 3848 fontsize = plot_size_in_points / 40. 3849 symsize = plot_size_in_points / 61. 3850 3851 direction_letters = 'NSEWDU' 3852 3853 for val in self._all_BV_2D_US: 3854 x_coord = val[0] 3855 y_coord = val[1] 3856 np_letter = direction_letters[int(val[2])] 3857 3858 rot_angle = -np.arctan2(y_coord, x_coord) + pi / 2. 3859 original_rho = np.sqrt(x_coord ** 2 + y_coord ** 2) 3860 3861 marker_x = (original_rho - (1.5 * symsize / points_per_unit)) \ 3862 * np.sin(rot_angle) 3863 marker_y = (original_rho - (1.5 * symsize / points_per_unit)) \ 3864 * np.cos(rot_angle) 3865 annot_x = (original_rho - (4.5 * fontsize / points_per_unit)) \ 3866 * np.sin(rot_angle) 3867 annot_y = (original_rho - (4.5 * fontsize / points_per_unit)) \ 3868 * np.cos(rot_angle) 3869 3870 ax.text(annot_x, annot_y, np_letter, 3871 horizontalalignment='center', size=fontsize, 3872 weight='bold', verticalalignment='center', 3873 bbox=dict(edgecolor='white', facecolor='white', 3874 alpha=1)) 3875 3876 if original_rho > epsilon: 3877 ax.scatter([marker_x], [marker_y], 3878 marker=(3, 0, rot_angle), s=symsize ** 2, 3879 c='k', facecolor='k', zorder=300) 3880 else: 3881 ax.scatter([x_coord], [y_coord], marker=(4, 1, rot_angle), 3882 s=symsize ** 2, c='k', facecolor='k', 3883 zorder=300) 3884 3885 # plot 4 fake points, guaranteeing full visibility of the sphere 3886 ax.plot([0, 1.05, 0, -1.05], [1.05, 0, -1.05, 0], ',', alpha=0.) 3887 # scaling behavior 3888 ax.autoscale_view(tight=True, scalex=True, scaley=True) 3889 3890 return plotfig 3891 3892 3893# ------------------------------------------------------------------- 3894# 3895# input and call management 3896# 3897# ------------------------------------------------------------------- 3898 3899def main(argv=None): 3900 """ 3901 Usage 3902 3903 obspy-mopad plot M 3904 obspy-mopad save M 3905 obspy-mopad gmt M 3906 obspy-mopad convert M 3907 obspy-mopad decompose M 3908 """ 3909 def _handle_input(M_in, args): 3910 """ 3911 take the original method and its arguments, the source mechanism, 3912 and the dictionary with proper parsers for each call, 3913 """ 3914 # construct a dict with consistent keyword args suited for the current 3915 # call 3916 kwargs = args.build(args) 3917 # set the fitting input basis system 3918 in_system = kwargs.get('in_system', 'NED') 3919 # build the moment tensor object 3920 mt = MomentTensor(M=M_in, system=in_system) 3921 # call the main routine to handle the moment tensor 3922 return args.call(mt, kwargs) 3923 3924 def _call_plot(MT, kwargs_dict): 3925 """ 3926 """ 3927 bb2plot = BeachBall(MT, kwargs_dict) 3928 3929 if kwargs_dict['plot_save_plot']: 3930 bb2plot.save_BB(kwargs_dict) 3931 return 3932 3933 if kwargs_dict['plot_pa_plot']: 3934 bb2plot.pa_plot(kwargs_dict) 3935 return 3936 3937 if kwargs_dict['plot_full_sphere']: 3938 bb2plot.full_sphere_plot(kwargs_dict) 3939 return 3940 3941 bb2plot.ploBB(kwargs_dict) 3942 3943 return 3944 3945 def _call_convert(MT, kwargs_dict): 3946 """ 3947 """ 3948 if kwargs_dict['type_conversion']: 3949 3950 if kwargs_dict['type_conversion'] == 'SDR': 3951 if kwargs_dict['fancy_conversion']: 3952 return MT.get_fps(style='f') 3953 else: 3954 return MT.get_fps(style='n') 3955 elif kwargs_dict['type_conversion'] == 'T': 3956 if kwargs_dict['basis_conversion']: 3957 out_system = kwargs_dict['out_system'] 3958 if kwargs_dict['fancy_conversion']: 3959 return MT.get_M(system=out_system, style='f') 3960 else: 3961 return MT.get_M(system=out_system, style='n') 3962 else: 3963 if kwargs_dict['fancy_conversion']: 3964 return MT.get_M(style='f') 3965 else: 3966 return MT.get_M(style='n') 3967 3968 if kwargs_dict['basis_conversion']: 3969 if len(MT._original_M) in [6, 7]: 3970 if len(MT._original_M) == 6: 3971 M_converted = _puzzle_basis_transformation( 3972 MT.get_M(), 3973 'NED', kwargs_dict['out_system']) 3974 if kwargs_dict['fancy_conversion']: 3975 print('\n Moment tensor in basis %s:\n ' % 3976 (kwargs_dict['in_system'])) 3977 print(fancy_matrix(MT.get_M( 3978 system=kwargs_dict['in_system']))) 3979 print() 3980 print('\n Moment tensor in basis %s:\n ' % 3981 (kwargs_dict['out_system'])) 3982 return fancy_matrix(M_converted) 3983 else: 3984 return M_converted[0, 0], M_converted[1, 1], \ 3985 M_converted[2, 2], M_converted[0, 1], \ 3986 M_converted[0, 2], M_converted[1, 2] 3987 else: 3988 M_converted = _puzzle_basis_transformation( 3989 MT.get_M(), 'NED', kwargs_dict['out_system']) 3990 if kwargs_dict['fancy_conversion']: 3991 print('\n Moment tensor in basis %s:\n ' % 3992 (kwargs_dict['in_system'])) 3993 print(fancy_matrix(MT.get_M( 3994 system=kwargs_dict['in_system']))) 3995 print() 3996 print('\n Moment tensor in basis %s:\n ' % 3997 (kwargs_dict['out_system'])) 3998 return fancy_matrix(M_converted) 3999 else: 4000 return M_converted[0, 0], M_converted[1, 1], \ 4001 M_converted[2, 2], M_converted[0, 1], \ 4002 M_converted[0, 2], M_converted[1, 2], \ 4003 MT._original_M[6] 4004 elif len(MT._original_M) == 9: 4005 new_M = np.asarray(MT._original_M).reshape(3, 3).copy() 4006 if kwargs_dict['fancy_conversion']: 4007 return fancy_matrix(_puzzle_basis_transformation( 4008 new_M, kwargs_dict['in_system'], 4009 kwargs_dict['out_system'])) 4010 else: 4011 return _puzzle_basis_transformation( 4012 new_M, kwargs_dict['in_system'], 4013 kwargs_dict['out_system']) 4014 elif len(MT._original_M) == 3: 4015 M_converted = _puzzle_basis_transformation( 4016 MT.get_M(), 'NED', kwargs_dict['out_system']) 4017 if kwargs_dict['fancy_conversion']: 4018 print('\n Moment tensor in basis %s: ' % 4019 (kwargs_dict['out_system'])) 4020 return fancy_matrix(M_converted) 4021 else: 4022 return M_converted[0, 0], M_converted[1, 1], \ 4023 M_converted[2, 2], M_converted[0, 1], \ 4024 M_converted[0, 2], M_converted[1, 2] 4025 elif len(MT._original_M) == 4: 4026 M_converted = MT._original_M[3] * \ 4027 _puzzle_basis_transformation(MT.get_M(), 'NED', 4028 kwargs_dict['out_system']) 4029 if kwargs_dict['fancy_conversion']: 4030 print('\n Momemnt tensor in basis %s: ' % 4031 (kwargs_dict['out_system'])) 4032 return fancy_matrix(M_converted) 4033 else: 4034 return M_converted[0, 0], M_converted[1, 1], \ 4035 M_converted[2, 2], M_converted[0, 1], \ 4036 M_converted[0, 2], M_converted[1, 2] 4037 else: 4038 print('this try is meaningless - read the possible', end=' ') 4039 print('choices!\n(perhaps you want option "-v"(convert a', 4040 end=' ') 4041 print('vector) or "-t"(convert strike, dip, rake to a matrix', 4042 end=' ') 4043 print('and show THAT matrix in another basis system)', end=' ') 4044 print('instead!?!?)\n') 4045 sys.exit(-1) 4046 4047 if kwargs_dict['vector_conversion']: 4048 if kwargs_dict['fancy_conversion']: 4049 print('\n Vector in basis %s:\n ' % 4050 (kwargs_dict['vector_in_system'])) 4051 print(fancy_vector(MT._original_M)) 4052 print() 4053 print('\n Vector in basis %s:\n ' % 4054 (kwargs_dict['vector_out_system'])) 4055 return fancy_vector(_puzzle_basis_transformation( 4056 MT._original_M, kwargs_dict['vector_in_system'], 4057 kwargs_dict['vector_out_system'])) 4058 else: 4059 return _puzzle_basis_transformation( 4060 MT._original_M, 4061 kwargs_dict['vector_in_system'], 4062 kwargs_dict['vector_out_system']) 4063 else: 4064 msg = 'provide either option -t with one argument, or -b with ' 4065 msg += 'two arguments, or -v with 2 arguments' 4066 sys.exit(msg) 4067 4068 def _call_gmt(MT, kwargs_dict): 4069 """ 4070 """ 4071 bb = BeachBall(MT, kwargs_dict) 4072 return bb.get_psxy(kwargs_dict).decode('utf-8') 4073 4074 def _call_decompose(MT, kwargs_dict): 4075 """ 4076 """ 4077 MT._isotropic = None 4078 MT._deviatoric = None 4079 MT._DC = None 4080 MT._iso_percentage = None 4081 MT._DC_percentage = None 4082 MT._DC2 = None 4083 MT._DC3 = None 4084 MT._DC2_percentage = None 4085 MT._CLVD = None 4086 MT._seismic_moment = None 4087 MT._moment_magnitude = None 4088 4089 out_system = kwargs_dict['decomp_out_system'] 4090 MT._output_basis = out_system 4091 MT._decomposition_key = kwargs_dict['decomposition_key'] 4092 4093 MT._decompose_M() 4094 4095 # if total decomposition: 4096 if kwargs_dict['decomp_out_complete']: 4097 if kwargs_dict['decomp_out_fancy']: 4098 decomp = MT.get_full_decomposition() 4099 try: 4100 print(decomp) 4101 except Exception: 4102 print(decomp.replace('°', ' deg')) 4103 return 4104 else: 4105 return MT.get_decomposition(in_system=kwargs_dict['in_system'], 4106 out_system=out_system) 4107 # otherwise: 4108 else: 4109 # argument dictionary - setting the appropriate calls 4110 do_calls = dict(zip(('in', 'out', 4111 'type', 4112 'full', 'm', 4113 'iso', 'iso_perc', 4114 'dev', 'devi', 'devi_perc', 4115 'dc', 'dc_perc', 4116 'dc2', 'dc2_perc', 4117 'dc3', 'dc3_perc', 4118 'clvd', 'clvd_perc', 4119 'mom', 'mag', 4120 'eigvals', 'eigvecs', 4121 't', 'n', 'p', 4122 'fps', 'faultplanes', 'fp', 4123 'decomp_key'), 4124 ('input_system', 'output_system', 4125 'decomp_type', 4126 'M', 'M', 4127 'iso', 'iso_percentage', 4128 'devi', 'devi', 'devi_percentage', 4129 'DC', 'DC_percentage', 4130 'DC2', 'DC2_percentage', 4131 'DC3', 'DC3_percentage', 4132 'CLVD', 'CLVD_percentage', 4133 'moment', 'mag', 4134 'eigvals', 'eigvecs', 4135 't_axis', 'null_axis', 'p_axis', 4136 'fps', 'fps', 'fps', 4137 'decomp_type') 4138 )) 4139 4140 # build argument for local call within MT object: 4141 lo_args = kwargs_dict['decomp_out_part'] 4142 4143 # for single element output: 4144 if len(lo_args) == 1: 4145 if kwargs_dict['decomp_out_fancy']: 4146 print(getattr(MT, 'get_' + do_calls[lo_args[0]])( 4147 style='f', system=out_system)) 4148 return 4149 else: 4150 return getattr(MT, 'get_' + do_calls[lo_args[0]])( 4151 style='n', system=out_system) 4152 # for list of elements: 4153 else: 4154 outlist = [] 4155 for arg in lo_args: 4156 if kwargs_dict['decomp_out_fancy']: 4157 print(getattr(MT, 'get_' + do_calls[arg])( 4158 style='f', system=out_system)) 4159 else: 4160 outlist.append(getattr(MT, 'get_' + do_calls[arg])( 4161 style='n', system=out_system)) 4162 if kwargs_dict['decomp_out_fancy']: 4163 return 4164 else: 4165 return outlist 4166 4167 def _build_gmt_dict(args): 4168 """ 4169 """ 4170 consistent_kwargs_dict = {} 4171 4172 if args.GMT_show_1FP: 4173 consistent_kwargs_dict['_GMT_1fp'] = args.GMT_show_1FP 4174 4175 if args.GMT_show_2FP2: 4176 args.GMT_show_1FP = 0 4177 4178 consistent_kwargs_dict['_GMT_2fps'] = True 4179 consistent_kwargs_dict['_GMT_1fp'] = 0 4180 4181 if args.GMT_string_type[0] not in ['F', 'L', 'E']: 4182 print('type of desired string not known - taking "fill" instead') 4183 consistent_kwargs_dict['_GMT_type'] = 'fill' 4184 4185 else: 4186 if args.GMT_string_type[0] == 'F': 4187 consistent_kwargs_dict['_GMT_type'] = 'fill' 4188 elif args.GMT_string_type[0] == 'L': 4189 consistent_kwargs_dict['_GMT_type'] = 'lines' 4190 else: 4191 consistent_kwargs_dict['_GMT_type'] = 'EVs' 4192 4193 if args.GMT_scaling < epsilon: 4194 print('GMT scaling factor must be a factor larger than') 4195 print('%f - set to 1, due to obviously stupid input value' % 4196 (epsilon)) 4197 args.GMT_scaling = 1.0 4198 4199 if args.plot_viewpoint: 4200 try: 4201 vp = args.plot_viewpoint.split(',') 4202 if not len(vp) == 3: 4203 raise 4204 if not -90 <= float(vp[0]) <= 90: 4205 raise 4206 if not -180 <= float(vp[1]) <= 180: 4207 raise 4208 if not 0 <= float(vp[2]) % 360 <= 360: 4209 raise 4210 consistent_kwargs_dict['plot_viewpoint'] = \ 4211 [float(vp[0]), float(vp[1]), float(vp[2])] 4212 except Exception: 4213 pass 4214 4215 if args.GMT_projection: 4216 lo_allowed_projections = ['STEREO', 'ORTHO', 'LAMBERT'] # ,'GNOM'] 4217 do_allowed_projections = {x[0]: x.lower() 4218 for x in lo_allowed_projections} 4219 try: 4220 gmtp = args.GMT_projection 4221 if gmtp in lo_allowed_projections: 4222 consistent_kwargs_dict['plot_projection'] = gmtp.lower() 4223 elif gmtp in do_allowed_projections.keys(): 4224 consistent_kwargs_dict['plot_projection'] = \ 4225 do_allowed_projections[gmtp] 4226 else: 4227 consistent_kwargs_dict['plot_projection'] = 'stereo' 4228 except Exception: 4229 pass 4230 4231 consistent_kwargs_dict['_GMT_scaling'] = args.GMT_scaling 4232 consistent_kwargs_dict['_GMT_tension_colour'] = args.GMT_tension_colour 4233 consistent_kwargs_dict['_GMT_pressure_colour'] = \ 4234 args.GMT_pressure_colour 4235 consistent_kwargs_dict['_plot_isotropic_part'] = \ 4236 args.GMT_plot_isotropic_part 4237 4238 return consistent_kwargs_dict 4239 4240 def _build_decompose_dict(args): 4241 """ 4242 """ 4243 consistent_kwargs_dict = {} 4244 4245 consistent_kwargs_dict['in_system'] = args.decomp_in_system 4246 4247 consistent_kwargs_dict['decomp_out_system'] = args.decomp_out_system 4248 4249 consistent_kwargs_dict['decomposition_key'] = args.decomp_key 4250 4251 # if option 'partial' is not chosen, take complete decomposition: 4252 if not args.decomp_out_part: 4253 consistent_kwargs_dict['decomp_out_complete'] = 1 4254 consistent_kwargs_dict['decomp_out_part'] = 0 4255 4256 if args.decomp_out_fancy: 4257 consistent_kwargs_dict['decomp_out_fancy'] = 1 4258 else: 4259 consistent_kwargs_dict['decomp_out_fancy'] = 0 4260 # otherwise take only partial decomposition: 4261 else: 4262 lo_allowed_attribs = ['in', 'out', 'type', 4263 'full', 'm', 4264 'iso', 'iso_perc', 4265 'dev', 'devi', 'devi_perc', 4266 'dc', 'dc_perc', 4267 'dc2', 'dc2_perc', 4268 'dc3', 'dc3_perc', 4269 'clvd', 'clvd_perc', 4270 'mom', 'mag', 4271 'eigvals', 'eigvecs', 4272 't', 'n', 'p', 4273 'fps', 'faultplanes', 'fp', 4274 'decomp_key'] 4275 lo_desired_attribs = args.decomp_out_part.split(',') 4276 lo_correct_attribs = [] 4277 4278 # check for allowed parts of decomposition: 4279 for da in lo_desired_attribs: 4280 if da.lower() in lo_allowed_attribs: 4281 lo_correct_attribs.append(da.lower()) 4282 4283 # if only wrong or no arguments are handed over, change to complete 4284 # decomposition: 4285 if len(lo_correct_attribs) == 0: 4286 print(' no correct attributes for partial decomposition - ' 4287 'returning complete decomposition') 4288 consistent_kwargs_dict['decomp_out_complete'] = 1 4289 consistent_kwargs_dict['decomp_out_part'] = 0 4290 if args.decomp_out_fancy: 4291 consistent_kwargs_dict['decomp_out_fancy'] = 1 4292 else: 4293 consistent_kwargs_dict['decomp_out_fancy'] = 0 4294 4295 # if only one part is desired to be shown, fancy style is possible 4296 elif len(lo_correct_attribs) == 1: 4297 consistent_kwargs_dict['decomp_out_complete'] = 0 4298 consistent_kwargs_dict['decomp_out_part'] = lo_correct_attribs 4299 if args.decomp_out_fancy: 4300 consistent_kwargs_dict['decomp_out_fancy'] = 1 4301 else: 4302 consistent_kwargs_dict['decomp_out_fancy'] = 0 4303 4304 # if several parts are desired to be shown, fancy style is 4305 # NOT possible: 4306 else: 4307 consistent_kwargs_dict['decomp_out_complete'] = 0 4308 consistent_kwargs_dict['decomp_out_part'] = lo_correct_attribs 4309 if args.decomp_out_fancy: 4310 consistent_kwargs_dict['decomp_out_fancy'] = 1 4311 else: 4312 consistent_kwargs_dict['decomp_out_fancy'] = 0 4313 4314 consistent_kwargs_dict['style'] = 'n' 4315 if consistent_kwargs_dict['decomp_out_fancy']: 4316 consistent_kwargs_dict['style'] = 'f' 4317 4318 return consistent_kwargs_dict 4319 4320 def _build_convert_dict(args): 4321 """ 4322 """ 4323 consistent_kwargs_dict = {} 4324 lo_allowed_options = ['type_conversion', 'basis_conversion', 4325 'vector_conversion', 'fancy_conversion'] 4326 # check for allowed options: 4327 for ao in lo_allowed_options: 4328 if hasattr(args, ao): 4329 consistent_kwargs_dict[ao] = getattr(args, ao) 4330 4331 consistent_kwargs_dict['in_system'] = 'NED' 4332 if 'out_system' not in consistent_kwargs_dict: 4333 consistent_kwargs_dict['out_system'] = 'NED' 4334 4335 if args.type_conversion and args.vector_conversion: 4336 print('decide for ONE option of "-t" OR "-v"') 4337 sys.exit(-1) 4338 4339 if args.basis_conversion: 4340 consistent_kwargs_dict['in_system'] = args.basis_conversion[0] 4341 consistent_kwargs_dict['out_system'] = args.basis_conversion[1] 4342 4343 if args.type_conversion and args.type_conversion == 'SDR': 4344 if args.basis_conversion: 4345 if args.basis_conversion[1] != 'NED': 4346 print('output "sdr" from type conversion cannot be ' 4347 'displayed in another basis system!') 4348 consistent_kwargs_dict['out_system'] = 'NED' 4349 4350 if args.vector_conversion: 4351 consistent_kwargs_dict['vector_in_system'] = \ 4352 args.vector_conversion[0] 4353 consistent_kwargs_dict['vector_out_system'] = \ 4354 args.vector_conversion[1] 4355 4356 return consistent_kwargs_dict 4357 4358 def _build_plot_dict(args): 4359 """ 4360 """ 4361 consistent_kwargs_dict = {} 4362 4363 consistent_kwargs_dict['plot_save_plot'] = False 4364 if args.plot_outfile: 4365 consistent_kwargs_dict['plot_save_plot'] = True 4366 lo_possible_formats = ['svg', 'png', 'eps', 'pdf', 'ps'] 4367 4368 try: 4369 (filepath, filename) = os.path.split(args.plot_outfile) 4370 if not filename: 4371 filename = 'dummy_filename.svg' 4372 (shortname, extension) = os.path.splitext(filename) 4373 if not shortname: 4374 shortname = 'dummy_shortname' 4375 4376 if extension[1:].lower() in lo_possible_formats: 4377 consistent_kwargs_dict['plot_outfile_format'] = \ 4378 extension[1:].lower() 4379 4380 if shortname.endswith('.'): 4381 consistent_kwargs_dict['plot_outfile'] = \ 4382 os.path.realpath(os.path.abspath(os.path.join( 4383 os.curdir, filepath, 4384 shortname + extension[1:].lower()))) 4385 else: 4386 consistent_kwargs_dict['plot_outfile'] = \ 4387 os.path.realpath(os.path.abspath(os.path.join( 4388 os.curdir, filepath, shortname + '.' + 4389 extension[1:].lower()))) 4390 else: 4391 if filename.endswith('.'): 4392 consistent_kwargs_dict['plot_outfile'] = \ 4393 os.path.realpath(os.path.abspath(os.path.join( 4394 os.curdir, filepath, 4395 filename + lo_possible_formats[0]))) 4396 else: 4397 consistent_kwargs_dict['plot_outfile'] = \ 4398 os.path.realpath(os.path.abspath(os.path.join( 4399 os.curdir, filepath, filename + '.' + 4400 lo_possible_formats[0]))) 4401 consistent_kwargs_dict['plot_outfile_format'] = \ 4402 lo_possible_formats[0] 4403 4404 except Exception: 4405 msg = 'please provide valid filename: <name>.<format> !!\n' 4406 msg += ' <format> must be svg, png, eps, pdf, or ps ' 4407 exit(msg) 4408 4409 if args.plot_pa_plot: 4410 consistent_kwargs_dict['plot_pa_plot'] = True 4411 else: 4412 consistent_kwargs_dict['plot_pa_plot'] = False 4413 4414 if args.plot_full_sphere: 4415 consistent_kwargs_dict['plot_full_sphere'] = True 4416 consistent_kwargs_dict['plot_pa_plot'] = False 4417 else: 4418 consistent_kwargs_dict['plot_full_sphere'] = False 4419 4420 if args.plot_viewpoint: 4421 try: 4422 vp = args.plot_viewpoint.split(',') 4423 if not len(vp) == 3: 4424 raise 4425 if not -90 <= float(vp[0]) <= 90: 4426 raise 4427 if not -180 <= float(vp[1]) <= 180: 4428 raise 4429 if not 0 <= float(vp[2]) % 360 <= 360: 4430 raise 4431 consistent_kwargs_dict['plot_viewpoint'] = \ 4432 [float(vp[0]), float(vp[1]), float(vp[2])] 4433 except Exception: 4434 pass 4435 4436 if args.plot_projection: 4437 lo_allowed_projections = ['STEREO', 'ORTHO', 'LAMBERT'] # ,'GNOM'] 4438 do_allowed_projections = {x[0]: x.lower() 4439 for x in lo_allowed_projections} 4440 try: 4441 ppl = args.plot_projection 4442 if ppl in lo_allowed_projections: 4443 consistent_kwargs_dict['plot_projection'] = ppl.lower() 4444 elif ppl in do_allowed_projections.keys(): 4445 consistent_kwargs_dict['plot_projection'] = \ 4446 do_allowed_projections[ppl] 4447 else: 4448 consistent_kwargs_dict['plot_projection'] = 'stereo' 4449 except Exception: 4450 pass 4451 4452 if args.plot_show_upper_hemis: 4453 consistent_kwargs_dict['plot_show_upper_hemis'] = True 4454 4455 if args.plot_n_points and args.plot_n_points > 360: 4456 consistent_kwargs_dict['plot_n_points'] = args.plot_n_points 4457 4458 if args.plot_size: 4459 try: 4460 if 0.01 < args.plot_size <= 1: 4461 consistent_kwargs_dict['plot_size'] = \ 4462 args.plot_size * 10 / 2.54 4463 elif 1 < args.plot_size < 45: 4464 consistent_kwargs_dict['plot_size'] = \ 4465 args.plot_size / 2.54 4466 else: 4467 consistent_kwargs_dict['plot_size'] = 5 4468 consistent_kwargs_dict['plot_aux_plot_size'] = \ 4469 consistent_kwargs_dict['plot_size'] 4470 except Exception: 4471 pass 4472 4473 if args.plot_pressure_colour: 4474 try: 4475 sec_colour_raw = args.plot_pressure_colour.split(',') 4476 if len(sec_colour_raw) == 1: 4477 if sec_colour_raw[0].lower()[0] in 'bgrcmykw': 4478 consistent_kwargs_dict['plot_pressure_colour'] = \ 4479 sec_colour_raw[0].lower()[0] 4480 else: 4481 raise 4482 elif len(sec_colour_raw) == 3: 4483 for sc in sec_colour_raw: 4484 if not 0 <= (int(sc)) <= 255: 4485 raise 4486 consistent_kwargs_dict['plot_pressure_colour'] = \ 4487 (float(sec_colour_raw[0]) / 255., 4488 float(sec_colour_raw[1]) / 255., 4489 float(sec_colour_raw[2]) / 255.) 4490 else: 4491 raise 4492 except Exception: 4493 pass 4494 4495 if args.plot_tension_colour: 4496 try: 4497 sec_colour_raw = args.plot_tension_colour.split(',') 4498 if len(sec_colour_raw) == 1: 4499 if sec_colour_raw[0].lower()[0] in 'bgrcmykw': 4500 consistent_kwargs_dict['plot_tension_colour'] = \ 4501 sec_colour_raw[0].lower()[0] 4502 else: 4503 raise 4504 elif len(sec_colour_raw) == 3: 4505 for sc in sec_colour_raw: 4506 if not 0 <= (int(float(sc))) <= 255: 4507 raise 4508 consistent_kwargs_dict['plot_tension_colour'] = \ 4509 (float(sec_colour_raw[0]) / 255., 4510 float(sec_colour_raw[1]) / 255., 4511 float(sec_colour_raw[2]) / 255.) 4512 else: 4513 raise 4514 except Exception: 4515 pass 4516 4517 if args.plot_total_alpha: 4518 if not 0 <= args.plot_total_alpha <= 1: 4519 consistent_kwargs_dict['plot_total_alpha'] = 1 4520 else: 4521 consistent_kwargs_dict['plot_total_alpha'] = \ 4522 args.plot_total_alpha 4523 4524 if args.plot_show_1faultplane: 4525 consistent_kwargs_dict['plot_show_1faultplane'] = True 4526 try: 4527 fp_args = args.plot_show_1faultplane 4528 4529 if not int(fp_args[0]) in [1, 2]: 4530 consistent_kwargs_dict['plot_show_FP_index'] = 1 4531 else: 4532 consistent_kwargs_dict['plot_show_FP_index'] = \ 4533 int(fp_args[0]) 4534 4535 if not 0 < float(fp_args[1]) <= 20: 4536 consistent_kwargs_dict['plot_faultplane_width'] = 2 4537 else: 4538 consistent_kwargs_dict['plot_faultplane_width'] = \ 4539 float(fp_args[1]) 4540 4541 try: 4542 sec_colour_raw = fp_args[2].split(',') 4543 if len(sec_colour_raw) == 1: 4544 sc = sec_colour_raw[0].lower()[0] 4545 if sc not in 'bgrcmykw': 4546 raise 4547 consistent_kwargs_dict['plot_faultplane_colour'] = \ 4548 sec_colour_raw[0].lower()[0] 4549 elif len(sec_colour_raw) == 3: 4550 for sc in sec_colour_raw: 4551 if not 0 <= (int(sc)) <= 255: 4552 raise 4553 consistent_kwargs_dict['plot_faultplane_colour'] = \ 4554 (float(sec_colour_raw[0]) / 255., 4555 float(sec_colour_raw[1]) / 255., 4556 float(sec_colour_raw[2]) / 255.) 4557 else: 4558 raise 4559 except Exception: 4560 consistent_kwargs_dict['plot_faultplane_colour'] = 'k' 4561 4562 try: 4563 if 0 <= float(fp_args[3]) <= 1: 4564 consistent_kwargs_dict['plot_faultplane_alpha'] = \ 4565 float(fp_args[3]) 4566 except Exception: 4567 consistent_kwargs_dict['plot_faultplane_alpha'] = 1 4568 except Exception: 4569 pass 4570 4571 if args.plot_show_faultplanes: 4572 consistent_kwargs_dict['plot_show_faultplanes'] = True 4573 consistent_kwargs_dict['plot_show_1faultplane'] = False 4574 4575 if args.plot_dpi: 4576 if 200 <= args.plot_dpi <= 2000: 4577 consistent_kwargs_dict['plot_dpi'] = args.plot_dpi 4578 4579 if args.plot_only_lines: 4580 consistent_kwargs_dict['plot_fill_flag'] = False 4581 4582 if args.plot_outerline: 4583 consistent_kwargs_dict['plot_outerline'] = True 4584 try: 4585 fp_args = args.plot_outerline 4586 if not 0 < float(fp_args[0]) <= 20: 4587 consistent_kwargs_dict['plot_outerline_width'] = 2 4588 else: 4589 consistent_kwargs_dict['plot_outerline_width'] = \ 4590 float(fp_args[0]) 4591 try: 4592 sec_colour_raw = fp_args[1].split(',') 4593 if len(sec_colour_raw) == 1: 4594 if sec_colour_raw[0].lower()[0] in 'bgrcmykw': 4595 consistent_kwargs_dict['plot_outerline_colour'] = \ 4596 sec_colour_raw[0].lower()[0] 4597 else: 4598 raise 4599 elif len(sec_colour_raw) == 3: 4600 for sc in sec_colour_raw: 4601 if not 0 <= (int(sc)) <= 255: 4602 raise 4603 consistent_kwargs_dict['plot_outerline_colour'] = \ 4604 (float(sec_colour_raw[0]) / 255., 4605 float(sec_colour_raw[1]) / 255., 4606 float(sec_colour_raw[2]) / 255.) 4607 else: 4608 raise 4609 except Exception: 4610 consistent_kwargs_dict['plot_outerline_colour'] = 'k' 4611 4612 try: 4613 if 0 <= float(fp_args[2]) <= 1: 4614 consistent_kwargs_dict['plot_outerline_alpha'] = \ 4615 float(fp_args[2]) 4616 except Exception: 4617 consistent_kwargs_dict['plot_outerline_alpha'] = 1 4618 except Exception: 4619 pass 4620 4621 if args.plot_nodalline: 4622 consistent_kwargs_dict['plot_nodalline'] = True 4623 try: 4624 fp_args = args.plot_nodalline 4625 4626 if not 0 < float(fp_args[0]) <= 20: 4627 consistent_kwargs_dict['plot_nodalline_width'] = 2 4628 else: 4629 consistent_kwargs_dict['plot_nodalline_width'] = \ 4630 float(fp_args[0]) 4631 try: 4632 sec_colour_raw = fp_args[1].split(',') 4633 if len(sec_colour_raw) == 1: 4634 if sec_colour_raw[0].lower()[0] in 'bgrcmykw': 4635 consistent_kwargs_dict['plot_nodalline_colour'] = \ 4636 sec_colour_raw[0].lower()[0] 4637 else: 4638 raise 4639 elif len(sec_colour_raw) == 3: 4640 for sc in sec_colour_raw: 4641 if not 0 <= (int(sc)) <= 255: 4642 raise 4643 consistent_kwargs_dict['plot_nodalline_colour'] = \ 4644 (float(sec_colour_raw[0]) / 255., 4645 float(sec_colour_raw[1]) / 255., 4646 float(sec_colour_raw[2]) / 255.) 4647 else: 4648 raise 4649 except Exception: 4650 consistent_kwargs_dict['plot_nodalline_colour'] = 'k' 4651 try: 4652 if 0 <= float(fp_args[2]) <= 1: 4653 consistent_kwargs_dict['plot_nodalline_alpha'] = \ 4654 float(fp_args[2]) 4655 except Exception: 4656 consistent_kwargs_dict['plot_nodalline_alpha'] = 1 4657 except Exception: 4658 pass 4659 4660 if args.plot_show_princ_axes: 4661 consistent_kwargs_dict['plot_show_princ_axes'] = True 4662 try: 4663 fp_args = args.plot_show_princ_axes 4664 4665 if not 0 < float(fp_args[0]) <= 40: 4666 consistent_kwargs_dict['plot_princ_axes_symsize'] = 10 4667 else: 4668 consistent_kwargs_dict['plot_princ_axes_symsize'] = \ 4669 float(fp_args[0]) 4670 4671 if not 0 < float(fp_args[1]) <= 20: 4672 consistent_kwargs_dict['plot_princ_axes_lw '] = 3 4673 else: 4674 consistent_kwargs_dict['plot_princ_axes_lw '] = \ 4675 float(fp_args[1]) 4676 try: 4677 if 0 <= float(fp_args[2]) <= 1: 4678 consistent_kwargs_dict['plot_princ_axes_alpha'] = \ 4679 float(fp_args[2]) 4680 except Exception: 4681 consistent_kwargs_dict['plot_princ_axes_alpha'] = 1 4682 except Exception: 4683 pass 4684 4685 if args.plot_show_basis_axes: 4686 consistent_kwargs_dict['plot_show_basis_axes'] = True 4687 4688 consistent_kwargs_dict['in_system'] = args.plot_input_system 4689 4690 if args.plot_isotropic_part: 4691 consistent_kwargs_dict['plot_isotropic_part'] = \ 4692 args.plot_isotropic_part 4693 4694 return consistent_kwargs_dict 4695 4696 def _build_parsers(): 4697 """ 4698 build dictionary with 4 (5 incl. 'save') sets of options, belonging to 4699 the 4 (5) possible calls 4700 """ 4701 from argparse import (ArgumentParser, 4702 RawDescriptionHelpFormatter, 4703 RawTextHelpFormatter, 4704 SUPPRESS) 4705 4706 parser = ArgumentParser(prog='obspy-mopad', 4707 formatter_class=RawDescriptionHelpFormatter, 4708 description=""" 4709############################################################################### 4710################################ MoPaD ################################ 4711################ Moment tensor Plotting and Decomposition tool ################ 4712############################################################################### 4713 4714Multi method tool for: 4715 4716- Plotting and saving of focal sphere diagrams ('Beachballs'). 4717 4718- Decomposition and Conversion of seismic moment tensors. 4719 4720- Generating coordinates, describing a focal sphere diagram, to be 4721 piped into GMT's psxy (Useful where psmeca or pscoupe fail.) 4722 4723For more help, please run ``%(prog)s {command} --help''. 4724 4725------------------------------------------------------------------------------- 4726 4727Example: 4728 4729To generate a beachball for a normal faulting mechanism (a snake's eye type): 4730 4731 %(prog)s plot 0,45,-90 or %(prog)s plot p 0,1,-1,0,0,0 4732 """) 4733 4734 parser.add_argument('-V', '--version', action='version', 4735 version='%(prog)s ' + __version__) 4736 4737 mechanism = ArgumentParser(add_help=False, 4738 formatter_class=RawTextHelpFormatter) 4739 mechanism.add_argument('mechanism', metavar='source-mechanism', 4740 help=""" 4741The 'source mechanism' as a comma-separated list of length: 4742 47433: 4744 strike, dip, rake; 47454: 4746 strike, dip, rake, moment; 47476: 4748 M11, M22, M33, M12, M13, M23; 47497: 4750 M11, M22, M33, M12, M13, M23, moment; 47519: 4752 full moment tensor 4753 4754(With all angles to be given in degrees) 4755 """) 4756 4757 subparsers = parser.add_subparsers(title='commands') 4758 4759 # Case-insensitive typing 4760 class caps(str): 4761 def __new__(self, content): 4762 return str.__new__(self, content.upper()) 4763 4764 # Possible basis systems 4765 ALLOWED_BASES = ['NED', 'USE', 'XYZ', 'NWU'] 4766 4767 # gmt 4768 help = "return the beachball as a string, to be piped into GMT's psxy" 4769 desc = """Tool providing strings to be piped into the 'psxy' from GMT. 4770 4771 Either a string describing the fillable area (to be used with option 4772 '-L' within psxy) or the nodallines or the coordinates of the principle 4773 axes are given. 4774 """ 4775 4776 parser_gmt = subparsers.add_parser('gmt', help=help, description=desc, 4777 parents=[mechanism]) 4778 4779 group_type = parser_gmt.add_argument_group('Output') 4780 group_show = parser_gmt.add_argument_group('Appearance') 4781 group_geo = parser_gmt.add_argument_group('Geometry') 4782 4783 group_type.add_argument( 4784 '-t', '--type', dest='GMT_string_type', metavar='<type>', 4785 type=caps, default='FILL', 4786 help='choice of psxy data: area to fill (fill), nodal lines ' 4787 '(lines), or eigenvector positions (ev)') 4788 group_show.add_argument( 4789 '-s', '--scaling', dest='GMT_scaling', metavar='<scaling factor>', 4790 type=float, default=1.0, 4791 help='spatial scaling of the beachball') 4792 group_show.add_argument( 4793 '-r', '--color1', '--colour1', dest='GMT_tension_colour', 4794 metavar='<tension colour>', type=int, default=1, 4795 help="-Z option's key for the tension colour of the beachball - " 4796 'type: integer') 4797 group_show.add_argument( 4798 '-w', '--color2', '--colour2', dest='GMT_pressure_colour', 4799 metavar='<pressure colour>', type=int, default=0, 4800 help="-Z option's key for the pressure colour of the beachball - " 4801 'type: integer') 4802 group_show.add_argument( 4803 '-D', '--faultplanes', dest='GMT_show_2FP2', action='store_true', 4804 help='boolean key, if 2 faultplanes shall be shown') 4805 group_show.add_argument( 4806 '-d', '--show-1fp', dest='GMT_show_1FP', metavar='<FP index>', 4807 type=int, choices=[1, 2], default=False, 4808 help='integer key (1,2), what faultplane shall be shown ' 4809 '[%(default)s]') 4810 group_geo.add_argument( 4811 '-V', '--viewpoint', dest='plot_viewpoint', 4812 metavar='<lat,lon,azi>', 4813 help='coordinates of the viewpoint - 3-tuple of angles in degree') 4814 group_geo.add_argument( 4815 '-p', '--projection', dest='GMT_projection', 4816 metavar='<projection>', type=caps, 4817 help='projection of the sphere') 4818 group_show.add_argument( 4819 '-I', '--show-isotropic-part', dest='GMT_plot_isotropic_part', 4820 action='store_true', 4821 help='if isotropic part shall be considered for plotting ' 4822 '[%(default)s]') 4823 4824 parser_gmt.set_defaults(call=_call_gmt, build=_build_gmt_dict) 4825 4826 # convert 4827 help = 'convert a mechanism to/in (strike,dip,slip-rake) from/to ' \ 4828 'matrix form *or* convert a matrix, vector, tuple into ' \ 4829 'different basis representations' 4830 desc = """Tool providing converted input. 4831 4832 Choose between the conversion from/to matrix-moment-tensor 4833 form (-t), the change of the output basis system for a given 4834 moment tensor (-b), or the change of basis for a 3D vector 4835 (-v). 4836 """ 4837 parser_convert = subparsers.add_parser('convert', help=help, 4838 description=desc, 4839 parents=[mechanism]) 4840 4841 group_type = parser_convert.add_argument_group('Type conversion') 4842 group_basis = parser_convert.add_argument_group('M conversion') 4843 group_vector = parser_convert.add_argument_group('Vector conversion') 4844# group_show = parser_convert.add_argument_group('Appearance') 4845 4846 group_type.add_argument( 4847 '-t', '--type', dest='type_conversion', 4848 type=caps, choices=['SDR', 'T'], 4849 help='type conversion - convert to: strike,dip,rake (sdr) or ' 4850 'Tensor (T)') 4851 group_basis.add_argument( 4852 '-b', '--basis', dest='basis_conversion', 4853 type=caps, choices=ALLOWED_BASES, nargs=2, 4854 help='basis conversion for M - provide 2 arguments: input and ' 4855 'output bases') 4856 group_vector.add_argument( 4857 '-v', '--vector', dest='vector_conversion', 4858 type=caps, choices=ALLOWED_BASES, nargs=2, 4859 help='basis conversion for a vector - provide M as a 3Dvector ' 4860 'and 2 option-arguments of -v: input and output bases') 4861 parser_convert.add_argument( 4862 '-y', '--fancy', dest='fancy_conversion', action='store_true', 4863 help='output in a stylish way') 4864 4865 parser_convert.set_defaults(call=_call_convert, 4866 build=_build_convert_dict) 4867 4868 # plot 4869 help = 'plot a beachball projection of the provided mechanism' 4870 desc = """Plots a beachball diagram of the provided mechanism. 4871 4872 Several styles and configurations are available. Also saving 4873 on the fly can be enabled. 4874 """ 4875 parser_plot = subparsers.add_parser('plot', help=help, 4876 description=desc, 4877 parents=[mechanism]) 4878 4879 group_save = parser_plot.add_argument_group('Saving') 4880 group_type = parser_plot.add_argument_group('Type of plot') 4881 group_quality = parser_plot.add_argument_group('Quality') 4882 group_colours = parser_plot.add_argument_group('Colours') 4883 group_misc = parser_plot.add_argument_group('Miscellaneous') 4884 group_dc = parser_plot.add_argument_group('Fault planes') 4885 group_geo = parser_plot.add_argument_group('Geometry') 4886 group_app = parser_plot.add_argument_group('Appearance') 4887 4888 group_save.add_argument( 4889 '-f', '--output-file', dest='plot_outfile', metavar='<filename>', 4890 help='filename for saving') 4891 group_type.add_argument( 4892 '-P', '--pa-system', dest='plot_pa_plot', action='store_true', 4893 help='if principal axis system shall be plotted instead') 4894 group_type.add_argument( 4895 '-O', '--full-sphere', dest='plot_full_sphere', 4896 action='store_true', 4897 help='if full sphere shall be plotted instead') 4898 group_geo.add_argument( 4899 '-V', '--viewpoint', dest='plot_viewpoint', 4900 metavar='<lat,lon,azi>', 4901 help='coordinates of the viewpoint - 3-tuple') 4902 group_geo.add_argument( 4903 '-p', '--projection', dest='plot_projection', 4904 metavar='<projection>', type=caps, 4905 help='projection of the sphere') 4906 group_type.add_argument( 4907 '-U', '--upper', dest='plot_show_upper_hemis', action='store_true', 4908 help='if upper hemisphere shall be shown ') 4909 group_quality.add_argument( 4910 '-N', '--points', dest='plot_n_points', metavar='<no. of points>', 4911 type=int, 4912 help='minimum number of points, used for nodallines') 4913 group_app.add_argument( 4914 '-s', '--size', dest='plot_size', metavar='<size in cm>', 4915 type=float, 4916 help='size of plot in cm') 4917 group_colours.add_argument( 4918 '-w', '--pressure-color', '--pressure-colour', 4919 dest='plot_pressure_colour', metavar='<colour>', 4920 help='colour of the tension area') 4921 group_colours.add_argument( 4922 '-r', '--tension-color', '--tension-colour', 4923 dest='plot_tension_colour', metavar='<colour>', 4924 help='colour of the pressure area ') 4925 group_app.add_argument( 4926 '-a', '--alpha', dest='plot_total_alpha', metavar='<alpha>', 4927 type=float, 4928 help='alpha value for the plot - float from 1=opaque to ' 4929 '0=transparent') 4930 group_dc.add_argument( 4931 '-D', '--dc', dest='plot_show_faultplanes', action='store_true', 4932 help='if double couple faultplanes shall be plotted ') 4933 group_dc.add_argument( 4934 '-d', '--show-1fp', dest='plot_show_1faultplane', nargs=4, 4935 metavar=('<index>', '<linewidth>', '<colour>', '<alpha>'), 4936 help='plot 1 faultplane - arguments are: index [1,2] of the ' 4937 'resp. FP, linewidth (float), line colour (string or ' 4938 'rgb-tuple), and alpha value (float between 0 and 1)') 4939 group_misc.add_argument( 4940 '-e', '--eigenvectors', dest='plot_show_princ_axes', 4941 metavar=('<size>', '<linewidth>', '<alpha>'), nargs=3, 4942 help='show eigenvectors - if used, provide 3 arguments: symbol ' 4943 'size, symbol linewidth, and symbol alpha value') 4944 group_misc.add_argument( 4945 '-b', '--basis-vectors', dest='plot_show_basis_axes', 4946 action='store_true', 4947 help='show NED basis in plot') 4948 group_app.add_argument( 4949 '-l', '--lines', dest='plot_outerline', 4950 metavar=('<linewidth>', '<colour>', '<alpha>'), nargs=3, 4951 help='gives the style of the outer line - 3 arguments needed: ' 4952 'linewidth (float), line colour (string or RGB-tuple), and ' 4953 'alpha value (float between 0 and 1)') 4954 group_app.add_argument( 4955 '-n', '--nodals', dest='plot_nodalline', 4956 metavar=('<linewidth>', '<colour>', '<alpha>'), nargs=3, 4957 help='gives the style of the nodal lines - 3 arguments needed: ' 4958 'linewidth (float), line colour (string or RGB-tuple), and ' 4959 'alpha value (float between 0 and 1)') 4960 group_quality.add_argument( 4961 '-q', '--quality', dest='plot_dpi', metavar='<dpi>', type=int, 4962 help='changes the quality for the plot in terms of dpi ' 4963 '(minimum=200)') 4964 group_type.add_argument( 4965 '-L', '--lines-only', dest='plot_only_lines', action='store_true', 4966 help='if only lines are shown (no fill - so overwrites ' 4967 '"fill"-related options)') 4968 group_misc.add_argument( 4969 '-i', '--input-system', dest='plot_input_system', 4970 type=caps, choices=ALLOWED_BASES, default='NED', 4971 help='if source mechanism is given as tensor in a system other ' 4972 'than NED') 4973 group_type.add_argument( 4974 '-I', '--show-isotropic-part', dest='plot_isotropic_part', 4975 action='store_true', 4976 help='if isotropic part shall be considered for plotting ' 4977 '[%(default)s]') 4978 4979 parser_plot.set_defaults(call=_call_plot, build=_build_plot_dict) 4980 4981 # decompose 4982 help = 'decompose a given mechanism into its components' 4983 desc = """Returns a decomposition of the input moment tensor.\n 4984 4985 Different decompositions are available (following Jost & 4986 Herrmann). Either the complete decomposition or only parts are 4987 returned; in- and output basis systema can be chosen. The 4988 'fancy' option is available for better human reading. 4989 """ 4990 4991 parser_decompose = subparsers.add_parser('decompose', help=help, 4992 description=desc, 4993 parents=[mechanism]) 4994 4995 group_type = parser_decompose.add_argument_group( 4996 'Type of decomposition') 4997 group_part = parser_decompose.add_argument_group( 4998 'Partial decomposition') 4999 group_system = parser_decompose.add_argument_group( 5000 'Basis systems') 5001 5002 helpstring11 = """ 5003 Returns a list of the decomposition results. 5004 5005 Order: 5006 1: 5007 basis of the provided input (string); 5008 2: 5009 basis of the representation (string); 5010 3: 5011 chosen decomposition type (integer); 5012 5013 4: 5014 full moment tensor (matrix); 5015 5016 5: 5017 isotropic part (matrix); 5018 6: 5019 isotropic percentage (float); 5020 7: 5021 deviatoric part (matrix); 5022 8: 5023 deviatoric percentage (float); 5024 5025 9: 5026 DC part (matrix); 5027 10: 5028 DC percentage (float); 5029 11: 5030 DC2 part (matrix); 5031 12: 5032 DC2 percentage (float); 5033 13: 5034 DC3 part (matrix); 5035 14: 5036 DC3 percentage (float); 5037 5038 15: 5039 CLVD part (matrix); 5040 16: 5041 CLVD percentage (matrix); 5042 5043 17: 5044 seismic moment (float); 5045 18: 5046 moment magnitude (float); 5047 5048 19: 5049 eigenvectors (3-array); 5050 20: 5051 eigenvalues (list); 5052 21: 5053 p-axis (3-array); 5054 22: 5055 neutral axis (3-array); 5056 23: 5057 t-axis (3-array); 5058 24: 5059 faultplanes (list of two 3-arrays). 5060 5061 5062 If option 'fancy' is set, only a small overview about geometry and 5063 strength is provided instead. 5064 """ 5065 group_part.add_argument( 5066 '-c', '--complete', dest='decomp_out_complete', 5067 action='store_true', help=helpstring11) 5068 parser_decompose.add_argument( 5069 '-y', '--fancy', dest='decomp_out_fancy', action='store_true', 5070 help='key for a stylish output') 5071 group_part.add_argument( 5072 '-p', '--partial', dest='decomp_out_part', 5073 metavar='<part1,part2,... >', 5074 help='provide an argument, what part(s) shall be displayed (if ' 5075 'multiple, separate by commas): in, out, type, full, iso, ' 5076 'iso_perc, devi, devi_perc, dc, dc_perc, dc2, dc2_perc, dc3, ' 5077 'dc3_perc, clvd, clvd_perc, mom, mag, eigvals, eigvecs, t, ' 5078 'n, p, faultplanes') 5079 group_system.add_argument( 5080 '-i', '--input-system', dest='decomp_in_system', 5081 type=caps, choices=ALLOWED_BASES, default='NED', 5082 help='set to provide input in a system other than NED') 5083 group_system.add_argument( 5084 '-o', '--output-system', dest='decomp_out_system', 5085 type=caps, choices=ALLOWED_BASES, default='NED', 5086 help='set to return output in a system other than NED') 5087 group_type.add_argument( 5088 '-t', '--type', dest='decomp_key', metavar='<decomposition key>', 5089 type=int, choices=[20, 21, 31], default=20, 5090 help='integer key to choose the type of decomposition - 20: ' 5091 'ISO+DC+CLVD ; 21: ISO+major DC+ minor DC ; 31: ISO + 3 DCs') 5092 5093 parser_decompose.set_defaults(call=_call_decompose, 5094 build=_build_decompose_dict) 5095 5096 return parser 5097 5098 parser = _build_parsers() 5099 args = parser.parse_args(argv) 5100 5101 try: 5102 M_raw = [float(xx) for xx in args.mechanism.split(',')] 5103 except Exception: 5104 parser.error('invalid source mechanism') 5105 5106 if not len(M_raw) in [3, 4, 6, 7, 9]: 5107 parser.error('invalid source mechanism') 5108 if len(M_raw) in [4, 6, 7, 9] and len(np.array(M_raw).nonzero()[0]) == 0: 5109 parser.error('invalid source mechanism') 5110 5111 aa = _handle_input(M_raw, args) 5112 if aa is not None: 5113 try: 5114 print(aa) 5115 except Exception: 5116 print(str(aa).replace('°', ' deg')) 5117 5118 5119if __name__ == '__main__': 5120 main() 5121