1# -*- coding: utf-8 -*- 2""" 3Python interface to the Seismic Analysis Code (SAC) file format. 4 5:copyright: 6 The Los Alamos National Security, LLC, Yannik Behr, C. J. Ammon, 7 C. Satriano, L. Krischer, and J. MacCarthy 8:license: 9 GNU Lesser General Public License, Version 3 10 (https://www.gnu.org/copyleft/lesser.html) 11 12 13The SACTrace object maintains consistency between SAC headers and manages 14header values in a user-friendly way. This includes some value-checking, native 15Python logicals (True, False) and nulls (None) instead of SAC's 0, 1, or 16-12345... 17 18SAC headers are implemented as properties, with appropriate getters and 19setters. 20 21 22Features 23-------- 24 251. **Read and write SAC binary or ASCII** 26 27 - autodetect or specify expected byteorder 28 - optional file size checking and/or header consistency checks 29 - header-only reading and writing 30 - "overwrite OK" checking ('lovrok' header) 31 322. **Convenient access and manipulation of relative and absolute time 33 headers** 343. **User-friendly header printing/viewing** 354. **Fast access to header values from attributes** 36 37 - With type checking, null handling, and enumerated value checking 38 395. **Convert to/from ObsPy Traces** 40 41 - Conversion from ObsPy Trace to SAC trace retains detected previous 42 SAC header values. 43 - Conversion to ObsPy Trace retains the *complete* SAC header. 44 45 46Usage examples 47-------------- 48 49Read/write SAC files 50~~~~~~~~~~~~~~~~~~~~ 51 52.. code:: python 53 54 # read from a binary file 55 sac = SACTrace.read(filename) 56 57 # read header only 58 sac = SACTrace.read(filename, headonly=True) 59 60 # write header-only, file must exist 61 sac.write(filename, headonly=True) 62 63 # read from an ASCII file 64 sac = SACTrace.read(filename, ascii=True) 65 66 # write a binary SAC file for a Sun machine 67 sac.write(filename, byteorder='big') 68 69Build a SACTrace from a header dictionary and data array 70~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 71 72.. rubric:: Example 73 74>>> header = {'kstnm': 'ANMO', 'kcmpnm': 'BHZ', 'stla': 40.5, 'stlo': -108.23, 75... 'evla': -15.123, 'evlo': 123, 'evdp': 50, 'nzyear': 2012, 76... 'nzjday': 123, 'nzhour': 13, 'nzmin': 43, 'nzsec': 17, 77... 'nzmsec': 100, 'delta': 1.0/40} 78>>> sac = SACTrace(data=np.random.random(100), **header) 79>>> print(sac) # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS 80Reference Time = 05/02/2012 (123) 13:43:17.100000 81 iztype IB: begin time 82b = 0.0 83delta = 0.0250000... 84e = 2.4750000... 85evdp = 50.0 86evla = -15.123000... 87evlo = 123.0 88iftype = itime 89internal0 = 2.0 90iztype = ib 91kcmpnm = BHZ 92kstnm = ANMO 93lcalda = False 94leven = True 95lovrok = True 96lpspol = True 97npts = 100 98nvhdr = 6 99nzhour = 13 100nzjday = 123 101nzmin = 43 102nzmsec = 100 103nzsec = 17 104nzyear = 2012 105stla = 40.5 106stlo = -108.23000... 107 108Reference-time and relative time headers 109~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 110 111.. rubric:: Example 112 113>>> sac = SACTrace(nzyear=2000, nzjday=1, nzhour=0, nzmin=0, nzsec=0, 114... nzmsec=0, t1=23.5, data=np.arange(100)) 115>>> print(sac.reftime) 1162000-01-01T00:00:00.000000Z 117>>> sac.b, sac.e, sac.t1 118(0.0, 99.0, 23.5) 119 120Move reference time by relative seconds, relative time headers are 121preserved. 122 123.. rubric:: Example 124 125>>> sac = SACTrace(nzyear=2000, nzjday=1, nzhour=0, nzmin=0, nzsec=0, 126... nzmsec=0, t1=23.5, data=np.arange(100)) 127>>> sac.reftime -= 2.5 128>>> sac.b, sac.e, sac.t1 129(2.5, 101.5, 26.0) 130 131Set reference time to new absolute time, relative time headers are 132preserved. 133 134.. rubric:: Example 135 136>>> sac = SACTrace(nzyear=2000, nzjday=1, nzhour=0, nzmin=0, nzsec=0, 137... nzmsec=0, t1=23.5, data=np.arange(100)) 138>>> # set the reftime two minutes later 139>>> sac.reftime = UTCDateTime(2000, 1, 1, 0, 2, 0, 0) 140>>> sac.b, sac.e, sac.t1 141(-120.0, -21.0, -96.5) 142 143Quick header viewing 144~~~~~~~~~~~~~~~~~~~~ 145 146Print non-null header values. 147 148.. rubric:: Example 149 150>>> sac = SACTrace() 151>>> print(sac) # doctest: +NORMALIZE_WHITESPACE 152Reference Time = 01/01/1970 (001) 00:00:00.000000 153 iztype IB: begin time 154b = 0.0 155delta = 1.0 156e = 0.0 157iftype = itime 158internal0 = 2.0 159iztype = ib 160lcalda = False 161leven = True 162lovrok = True 163lpspol = True 164npts = 0 165nvhdr = 6 166nzhour = 0 167nzjday = 1 168nzmin = 0 169nzmsec = 0 170nzsec = 0 171nzyear = 1970 172 173Print relative time header values. 174 175.. rubric:: Example 176 177>>> sac = SACTrace() 178>>> sac.lh('picks') # doctest: +NORMALIZE_WHITESPACE 179Reference Time = 01/01/1970 (001) 00:00:00.000000 180 iztype IB: begin time 181a = None 182b = 0.0 183e = 0.0 184f = None 185o = None 186t0 = None 187t1 = None 188t2 = None 189t3 = None 190t4 = None 191t5 = None 192t6 = None 193t7 = None 194t8 = None 195t9 = None 196 197Header values as attributes 198~~~~~~~~~~~~~~~~~~~~~~~~~~~ 199 200Great for interactive use, with (ipython) tab-completion... 201 202.. code:: python 203 204 sac.<tab> 205 206:: 207 208 sac.a sac.kevnm sac.nzsec 209 sac.az sac.kf sac.nzyear 210 sac.b sac.khole sac.o 211 sac.baz sac.kinst sac.odelta 212 sac.byteorder sac.knetwk sac.read 213 sac.cmpaz sac.ko sac.reftime 214 sac.cmpinc sac.kstnm sac.scale 215 sac.copy sac.kt0 sac.stdp 216 sac.data sac.kt1 sac.stel 217 sac.delta sac.kt2 sac.stla 218 sac.depmax sac.kt3 sac.stlo 219 sac.depmen sac.kt4 sac.t0 220 sac.depmin sac.kt5 sac.t1 221 sac.dist sac.kt6 sac.t2 222 sac.e sac.kt7 sac.t3 223 sac.evdp sac.kt8 sac.t4 224 sac.evla sac.kt9 sac.t5 225 sac.evlo sac.kuser0 sac.t6 226 sac.f sac.kuser1 sac.t7 227 sac.from_obspy_trace sac.kuser2 sac.t8 228 sac.gcarc sac.lcalda sac.t9 229 sac.idep sac.leven sac.to_obspy_trace 230 sac.ievreg sac.lh sac.unused23 231 sac.ievtyp sac.listhdr sac.user0 232 sac.iftype sac.lovrok sac.user1 233 sac.iinst sac.lpspol sac.user2 234 sac.imagsrc sac.mag sac.user3 235 sac.imagtyp sac.nevid sac.user4 236 sac.internal0 sac.norid sac.user5 237 sac.iqual sac.npts sac.user6 238 sac.istreg sac.nvhdr sac.user7 239 sac.isynth sac.nwfid sac.user8 240 sac.iztype sac.nzhour sac.user9 241 sac.ka sac.nzjday sac.validate 242 sac.kcmpnm sac.nzmin sac.write 243 sac.kdatrd sac.nzmsec 244 245...and documentation (in IPython)! 246 247.. code:: python 248 249 sac.iztype? 250 251:: 252 253 Type: property 254 String form: <property object at 0x106404940> 255 Docstring: 256 I Reference time equivalence: 257 258 * IUNKN (5): Unknown 259 * IB (9): Begin time 260 * IDAY (10): Midnight of reference GMT day 261 * IO (11): Event origin time 262 * IA (12): First arrival time 263 * ITn (13-22): User defined time pick n, n=0,9 264 265Convert to/from ObsPy Traces 266~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 267 268.. rubric:: Example 269 270>>> from obspy import read 271>>> tr = read()[0] 272>>> print(tr.stats) # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE 273 network: BW 274 station: RJOB 275 location: 276 channel: EHZ 277 starttime: 2009-08-24T00:20:03.000000Z 278 endtime: 2009-08-24T00:20:32.990000Z 279 sampling_rate: 100.0 280 delta: 0.01 281 npts: 3000 282 calib: 1.0 283 back_azimuth: 100.0 284 inclination: 30.0 285 response: Channel Response 286 ... 287 288 289>>> sac = SACTrace.from_obspy_trace(tr) 290>>> print(sac) # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS 291Reference Time = 08/24/2009 (236) 00:20:03.000000 292 iztype IB: begin time 293b = 0.0 294delta = 0.009999999... 295e = 29.989999... 296iftype = itime 297iztype = ib 298kcmpnm = EHZ 299knetwk = BW 300kstnm = RJOB 301lcalda = False 302leven = True 303lovrok = True 304lpspol = True 305npts = 3000 306nvhdr = 6 307nzhour = 0 308nzjday = 236 309nzmin = 20 310nzmsec = 0 311nzsec = 3 312nzyear = 2009 313scale = 1.0 314 315>>> tr2 = sac.to_obspy_trace() 316>>> print(tr2.stats) # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS 317 network: BW 318 station: RJOB 319 location: 320 channel: EHZ 321 starttime: 2009-08-24T00:20:03.000000Z 322 endtime: 2009-08-24T00:20:32.990000Z 323 sampling_rate: 100.0 324 delta: 0.01 325 npts: 3000 326 calib: 1.0 327 sac: AttribDict(...) 328 329""" 330from __future__ import (absolute_import, division, print_function, 331 unicode_literals) 332from future.builtins import * # NOQA 333from future.utils import native_str 334 335import sys 336import warnings 337from copy import deepcopy 338from itertools import chain 339 340import numpy as np 341from obspy import Trace, UTCDateTime 342from obspy.geodetics import gps2dist_azimuth, kilometer2degrees 343 344from . import header as HD # noqa 345from .util import SacError, SacHeaderError 346from . import util as _ut 347from . import arrayio as _io 348 349 350# ------------- HEADER DESCRIPTORS -------------------------------------------- 351# 352# A descriptor is a class that manages an object attribute, using the 353# descriptor protocol. A single instance of a descriptor class (FloatHeader, 354# for example) will exist for both the host class and all instances of that 355# host class (i.e. SACTrace and all its instances). As a result, we must 356# implement logic that can tell if the methods are being called on the host 357# class or on an instance. This looks like "if instance is None" on methods. 358# 359# See: 360# https://docs.python.org/3.5/howto/descriptor.html 361# https://nbviewer.jupyter.org/urls/gist.github.com/ChrisBeaumont/ 362# 5758381/raw/descriptor_writeup.ipynb 363 364class SACHeader(object): 365 def __init__(self, name): 366 try: 367 self.__doc__ = HD.DOC[name] 368 except KeyError: 369 # header doesn't have a docstring entry in HD.DOC 370 pass 371 self.name = name 372 373 374class FloatHeader(SACHeader): 375 def __get__(self, instance, instance_type): 376 if instance is None: 377 # a FloatHeader on the owner class was requested. 378 # return the descriptor itself. 379 value = self 380 else: 381 # a FloatHeader on an instance was requested. 382 # return the descriptor value. 383 value = float(instance._hf[HD.FLOATHDRS.index(self.name)]) 384 if value == HD.FNULL: 385 value = None 386 return value 387 388 def __set__(self, instance, value): 389 if value is None: 390 value = HD.FNULL 391 instance._hf[HD.FLOATHDRS.index(self.name)] = value 392 393 394# Descriptor for setting relative time headers with either a relative 395# time float or an absolute UTCDateTime 396# used for: b, o, a, f, t0-t9 397class RelativeTimeHeader(FloatHeader): 398 def __set__(self, instance, value): 399 """ 400 Intercept the set value to make sure it is an offset from the SAC 401 reference time. 402 403 """ 404 if isinstance(value, UTCDateTime): 405 offset = value - instance.reftime 406 else: 407 offset = value 408 # reuse the normal floatheader setter. 409 super(RelativeTimeHeader, self).__set__(instance, offset) 410 411 412# Factory function for setting geographic header values 413# (evlo, evla, stalo, stalat) 414# that will check lcalda and calculate and set dist, az, baz, gcarc 415class GeographicHeader(FloatHeader): 416 def __set__(self, instance, value): 417 super(GeographicHeader, self).__set__(instance, value) 418 if instance.lcalda: 419 try: 420 instance._set_distances() 421 except SacHeaderError: 422 pass 423 424 425class IntHeader(SACHeader): 426 def __get__(self, instance, instance_type): 427 if instance is None: 428 value = self 429 else: 430 value = int(instance._hi[HD.INTHDRS.index(self.name)]) 431 if value == HD.INULL: 432 value = None 433 return value 434 435 def __set__(self, instance, value): 436 if value is None: 437 value = HD.INULL 438 if value % 1: 439 warnings.warn("Non-integers may be truncated. ({}: {})".format( 440 self.name, value)) 441 instance._hi[HD.INTHDRS.index(self.name)] = value 442 443 444class BoolHeader(IntHeader): 445 def __get__(self, instance, instance_type): 446 # value can be an int or None 447 value = super(BoolHeader, self).__get__(instance, instance_type) 448 return bool(value) if value in (0, 1) else value 449 450 def __set__(self, instance, value): 451 if value not in (True, False, 1, 0): 452 msg = "Logical header values must be {True, False, 1, 0}" 453 raise ValueError(msg) 454 # booleans are subclasses of integers. They will be set (cast) 455 # directly into an integer array as 0 or 1. 456 super(BoolHeader, self).__set__(instance, value) 457 if self.name == 'lcalda': 458 if value: 459 try: 460 instance._set_distances() 461 except SacHeaderError: 462 pass 463 464 465class EnumHeader(IntHeader): 466 def __get__(self, instance, instance_type): 467 value = super(EnumHeader, self).__get__(instance, instance_type) 468 # value is int or None 469 if value is None: 470 name = None 471 elif _ut.is_valid_enum_int(self.name, value): 472 name = HD.ENUM_NAMES[value] 473 else: 474 msg = """Unrecognized enumerated value {} for header "{}". 475 See .header for allowed values.""".format(value, 476 self.name) 477 warnings.warn(msg) 478 name = None 479 return name 480 481 def __set__(self, instance, value): 482 if value is None: 483 value = HD.INULL 484 elif _ut.is_valid_enum_str(self.name, value): 485 if self.name == 'iztype': 486 reftime = _iztype_reftime(instance, value) 487 instance.reftime = reftime 488 # this also shifts all non-null relative times (instance._allt) 489 value = HD.ENUM_VALS[value] 490 else: 491 msg = 'Unrecognized enumerated value "{}" for header "{}"' 492 raise ValueError(msg.format(value, self.name)) 493 super(EnumHeader, self).__set__(instance, value) 494 495 496class StringHeader(SACHeader): 497 def __get__(self, instance, instance_type): 498 if instance is None: 499 value = self 500 else: 501 value = instance._hs[HD.STRHDRS.index(self.name)] 502 try: 503 # value is a bytes 504 value = value.decode() 505 except AttributeError: 506 # value is a str 507 pass 508 509 if value == HD.SNULL: 510 value = None 511 512 try: 513 value = value.strip() 514 except AttributeError: 515 # it's None. no .strip method 516 pass 517 return value 518 519 def __set__(self, instance, value): 520 if value is None: 521 value = HD.SNULL 522 elif len(value) > 8: 523 msg = ("Alphanumeric headers longer than 8 characters are " 524 "right-truncated.") 525 warnings.warn(msg) 526 # values will truncate themselves, since _hs is dtype '|S8' 527 try: 528 instance._hs[HD.STRHDRS.index(self.name)] = value.encode('ascii', 529 'strict') 530 except AttributeError: 531 instance._hs[HD.STRHDRS.index(self.name)] = value 532 533 534# Headers for functions of .data (min, max, mean, len) 535class DataHeader(SACHeader): 536 def __init__(self, name, func): 537 self.func = func 538 super(DataHeader, self).__init__(name) 539 540 def __get__(self, instance, instance_type): 541 if instance is None: 542 value = self 543 else: 544 try: 545 value = self.func(instance.data) 546 # convert to native Python types 547 if self.name in HD.FLOATHDRS: 548 value = float(value) 549 elif self.name in HD.INTHDRS: 550 value = int(value) 551 except TypeError: 552 # instance.data is None, get value from header 553 if self.name in HD.FLOATHDRS: 554 value = instance._hf[HD.FLOATHDRS.index(self.name)].item() 555 value = None if value == HD.FNULL else value 556 elif self.name in HD.INTHDRS: 557 value = instance._hi[HD.INTHDRS.index(self.name)].item() 558 value = None if value == HD.INULL else value 559 return value 560 561 def __set__(self, instance, value): 562 msg = "{} is read-only".format(self.name) 563 raise AttributeError(msg) 564 565 566# OTHER GETTERS/SETTERS 567def _get_e(self): 568 try: 569 if self.npts: 570 e = self.b + (self.npts - 1) * self.delta 571 else: 572 e = self.b 573 except TypeError: 574 # b, npts, and/or delta are None/null 575 # TODO: assume "b" is 0.0? 576 e = None 577 return e 578 579 580def _iztype_reftime(sactrace, iztype): 581 """ 582 Get the new reftime for a given iztype. 583 584 Setting the iztype will shift the relative time headers, such that the 585 header that iztype points to is (near) zero, and all others are shifted 586 together by the difference. 587 588 Affected headers: b, o, a, f, t0-t9 589 590 :param sactrace: 591 :type sactrace: SACTrace 592 :param iztype: One of the following strings: 593 'iunkn' 594 'ib', begin time 595 'iday', midnight of reference GMT day 596 'io', event origin time 597 'ia', first arrival time 598 'it0'-'it9', user defined pick t0-t9. 599 :type iztype: str 600 :rtype reftime: UTCDateTime 601 :return: The new SAC reference time. 602 603 """ 604 # The Plan: 605 # 1. find the seconds needed to shift the old reftime to the new one. 606 # 2. shift reference time onto the iztype header using that shift value. 607 # 3. this triggers an _allt shift of all relative times by the same amount. 608 # 4. If all goes well, actually set the iztype in the header. 609 610 # 1. 611 if iztype == 'iunkn': 612 # no shift 613 ref_val = 0.0 614 elif iztype == 'iday': 615 # seconds since midnight of reference day 616 reftime = sactrace.reftime 617 ref_val = reftime - UTCDateTime(year=reftime.year, 618 julday=reftime.julday) 619 else: 620 # a relative time header. 621 # remove the 'i' (first character) in the iztype to get the header name 622 ref_val = getattr(sactrace, iztype[1:]) 623 if ref_val is None: 624 msg = "Reference header for iztype '{}' is not set".format(iztype) 625 raise SacError(msg) 626 627 # 2. set a new reference time, 628 # 3. which also shifts all non-null relative times (sactrace._allt). 629 # remainder microseconds may be in the reference header value, because 630 # nzmsec can't hold them. 631 new_reftime = sactrace.reftime + ref_val 632 633 return new_reftime 634 635 636# kevnm is 16 characters, split into two 8-character fields 637# intercept and handle in while getting and setting 638def _get_kevnm(self): 639 kevnm = self._hs[HD.STRHDRS.index('kevnm')] 640 kevnm2 = self._hs[HD.STRHDRS.index('kevnm2')] 641 try: 642 kevnm = kevnm.decode() 643 kevnm2 = kevnm2.decode() 644 except AttributeError: 645 # kevnm is a str 646 pass 647 648 if kevnm == HD.SNULL: 649 kevnm = '' 650 if kevnm2 == HD.SNULL: 651 kevnm2 = '' 652 653 value = (kevnm + kevnm2).strip() 654 655 if not value: 656 value = None 657 658 return value 659 660 661def _set_kevnm(self, value): 662 if value is None: 663 value = HD.SNULL + HD.SNULL 664 elif len(value) > 16: 665 msg = "kevnm over 16 characters. Truncated to {}.".format(value[:16]) 666 warnings.warn(msg) 667 kevnm = '{:<8s}'.format(value[0:8]) 668 kevnm2 = '{:<8s}'.format(value[8:16]) 669 self._hs[HD.STRHDRS.index('kevnm')] = kevnm 670 self._hs[HD.STRHDRS.index('kevnm2')] = kevnm2 671 672# TODO: move get/set reftime up here, make it a property 673 674 675# -------------------------- SAC OBJECT INTERFACE ----------------------------- 676class SACTrace(object): 677 __doc__ = """ 678 Convenient and consistent in-memory representation of Seismic Analysis Code 679 (SAC) files. 680 681 This is the human-facing interface for making a valid instance. For 682 file-based or other constructors, see class methods .read and 683 .from_obspy_trace. SACTrace instances preserve relationships between 684 header values. 685 686 :param data: Associated time-series data vector. Optional. If omitted, None 687 is set as the instance data attribute. 688 :type data: :class:`numpy.ndarray` of float32 689 690 Any valid header key/value pair is also an optional input keyword argument. 691 If not provided, minimum required headers are set to valid default values. 692 The default instance is an evenly-space trace, with a sample rate of 1.0, 693 and len(data) or 0 npts, starting at 1970-01-01T00:00:00.000000. 694 695 :var reftime: Read-only reference time. Calculated from nzyear, nzjday, 696 nzhour, nzmin, nzsec, nzmsec. 697 :var byteorder: The byte order of the underlying header/data arrays. 698 Raises :class:`SacError` if array byte orders are inconsistent, even in 699 the case where '<' is your native order and byteorders look like '<', 700 '=', '='. 701 702 Any valid header name is also an attribute. See below, :mod:`header`, 703 or individial attribution docstrings for more header information. 704 705 THE SAC HEADER 706 707 NOTE: All header names and string values are lowercase. Header value 708 access should be through instance attributes. 709 710 """ + HD.HEADER_DOCSTRING 711 712 # ------------------------------- SAC HEADERS ----------------------------- 713 # SAC header values are defined as managed attributes, either as 714 # descriptors or as properties, with getters and setters. 715 # 716 # Managed attributes are defined at the class leval, and therefore shared 717 # across all instances, not attribute data themselves. 718 # 719 # This section looks ugly, but it allows for the following: 720 # 1. Relationships/checks between header variables can be done in setters 721 # 2. The underlying header array structure is retained, for quick writing. 722 # 3. Header access looks like simple attribute-access syntax. 723 # Looks funny to read here, but natural to use. 724 # 725 # FLOATS 726 delta = FloatHeader('delta') 727 depmin = DataHeader('depmin', min) 728 depmax = DataHeader('depmax', max) 729 scale = FloatHeader('scale') 730 odelta = FloatHeader('odelta') 731 b = RelativeTimeHeader('b') 732 e = property(_get_e, doc=HD.DOC['e']) 733 o = RelativeTimeHeader('o') 734 a = RelativeTimeHeader('a') 735 internal0 = FloatHeader('internal0') 736 t0 = RelativeTimeHeader('t0') 737 t1 = RelativeTimeHeader('t1') 738 t2 = RelativeTimeHeader('t2') 739 t3 = RelativeTimeHeader('t3') 740 t4 = RelativeTimeHeader('t4') 741 t5 = RelativeTimeHeader('t5') 742 t6 = RelativeTimeHeader('t6') 743 t7 = RelativeTimeHeader('t7') 744 t8 = RelativeTimeHeader('t8') 745 t9 = RelativeTimeHeader('t9') 746 f = RelativeTimeHeader('f') 747 stla = GeographicHeader('stla') 748 stlo = GeographicHeader('stlo') 749 stel = FloatHeader('stel') 750 stdp = FloatHeader('stdp') 751 evla = GeographicHeader('evla') 752 evlo = GeographicHeader('evlo') 753 evdp = FloatHeader('evdp') 754 mag = FloatHeader('mag') 755 user0 = FloatHeader('user0') 756 user1 = FloatHeader('user1') 757 user2 = FloatHeader('user2') 758 user3 = FloatHeader('user3') 759 user4 = FloatHeader('user4') 760 user5 = FloatHeader('user5') 761 user6 = FloatHeader('user6') 762 user7 = FloatHeader('user7') 763 user8 = FloatHeader('user8') 764 user9 = FloatHeader('user9') 765 dist = FloatHeader('dist') 766 az = FloatHeader('az') 767 baz = FloatHeader('baz') 768 gcarc = FloatHeader('gcarc') 769 depmen = DataHeader('depmen', np.mean) 770 cmpaz = FloatHeader('cmpaz') 771 cmpinc = FloatHeader('cmpinc') 772 # 773 # INTS 774 nzyear = IntHeader('nzyear') 775 nzjday = IntHeader('nzjday') 776 nzhour = IntHeader('nzhour') 777 nzmin = IntHeader('nzmin') 778 nzsec = IntHeader('nzsec') 779 nzmsec = IntHeader('nzmsec') 780 nvhdr = IntHeader('nvhdr') 781 norid = IntHeader('norid') 782 nevid = IntHeader('nevid') 783 npts = DataHeader('npts', len) 784 nwfid = IntHeader('nwfid') 785 iftype = EnumHeader('iftype') 786 idep = EnumHeader('idep') 787 iztype = EnumHeader('iztype') 788 iinst = IntHeader('iinst') 789 istreg = IntHeader('istreg') 790 ievreg = IntHeader('ievreg') 791 ievtyp = EnumHeader('ievtyp') 792 iqual = IntHeader('iqual') 793 isynth = EnumHeader('isynth') 794 imagtyp = EnumHeader('imagtyp') 795 imagsrc = EnumHeader('imagsrc') 796 leven = BoolHeader('leven') 797 lpspol = BoolHeader('lpspol') 798 lovrok = BoolHeader('lovrok') 799 lcalda = BoolHeader('lcalda') 800 unused23 = IntHeader('unused23') 801 # 802 # STRINGS 803 kstnm = StringHeader('kstnm') 804 kevnm = property(_get_kevnm, _set_kevnm, doc=HD.DOC['kevnm']) 805 khole = StringHeader('khole') 806 ko = StringHeader('ko') 807 ka = StringHeader('ka') 808 kt0 = StringHeader('kt0') 809 kt1 = StringHeader('kt1') 810 kt2 = StringHeader('kt2') 811 kt3 = StringHeader('kt3') 812 kt4 = StringHeader('kt4') 813 kt5 = StringHeader('kt5') 814 kt6 = StringHeader('kt6') 815 kt7 = StringHeader('kt7') 816 kt8 = StringHeader('kt8') 817 kt9 = StringHeader('kt9') 818 kf = StringHeader('kf') 819 kuser0 = StringHeader('kuser0') 820 kuser1 = StringHeader('kuser1') 821 kuser2 = StringHeader('kuser2') 822 kcmpnm = StringHeader('kcmpnm') 823 knetwk = StringHeader('knetwk') 824 kdatrd = StringHeader('kdatrd') 825 kinst = StringHeader('kinst') 826 827 def __init__(self, leven=True, delta=1.0, b=0.0, e=0.0, iztype='ib', 828 nvhdr=6, npts=0, iftype='itime', nzyear=1970, nzjday=1, 829 nzhour=0, nzmin=0, nzsec=0, nzmsec=0, lcalda=False, 830 lpspol=True, lovrok=True, internal0=2.0, data=None, **kwargs): 831 """ 832 Initialize a SACTrace object using header key-value pairs and a 833 numpy.ndarray for the data, both optional. 834 835 ..rubric:: Example 836 837 >>> sac = SACTrace(nzyear=1995, nzmsec=50, data=np.arange(100)) 838 >>> print(sac) # doctest: +NORMALIZE_WHITESPACE 839 Reference Time = 01/01/1995 (001) 00:00:00.050000 840 iztype IB: begin time 841 b = 0.0 842 delta = 1.0 843 e = 99.0 844 iftype = itime 845 internal0 = 2.0 846 iztype = ib 847 lcalda = False 848 leven = True 849 lovrok = True 850 lpspol = True 851 npts = 100 852 nvhdr = 6 853 nzhour = 0 854 nzjday = 1 855 nzmin = 0 856 nzmsec = 50 857 nzsec = 0 858 nzyear = 1995 859 860 """ 861 # The Plan: 862 # 1. Build the default header dictionary and update with provided 863 # values. 864 # 2. Convert header dict to arrays (util.dict_to_header_arrays 865 # initializes the arrays and fills in without checking. 866 # 3. set the _h[fis] and data arrays on self. 867 868 # 1. 869 # build the required header from provided or default values 870 header = {'leven': leven, 'npts': npts, 'delta': delta, 'b': b, 'e': e, 871 'iztype': iztype, 'nvhdr': nvhdr, 'iftype': iftype, 872 'nzyear': nzyear, 'nzjday': nzjday, 'nzhour': nzhour, 873 'nzmin': nzmin, 'nzsec': nzsec, 'nzmsec': nzmsec, 874 'lcalda': lcalda, 'lpspol': lpspol, 'lovrok': lovrok, 875 'internal0': internal0} 876 877 # combine header with remaining non-required args. 878 # user can put non-SAC key:value pairs into the header, but they're 879 # ignored on write. 880 header.update(kwargs) 881 882 # -------------------------- DATA ARRAY ------------------------------- 883 if data is None: 884 # this is like "headonly=True" 885 pass 886 else: 887 if not isinstance(data, np.ndarray): 888 raise TypeError("data needs to be a numpy.ndarray") 889 else: 890 # Only copy the data if they are not of the required type 891 # XXX: why require little endian instead of native byte order? 892 # data = np.require(data, native_str('<f4')) 893 pass 894 895 # --------------------------- HEADER ARRAYS --------------------------- 896 # 2. 897 # TODO: this is done even when we're reading a file. 898 # if it's too much overhead, it may need to change 899 900 # swap enum names for integer values in the header dictionary 901 header = _ut.enum_string_to_int(header) 902 903 # XXX: will these always be little endian? 904 hf, hi, hs = _io.dict_to_header_arrays(header) 905 906 # we now have data and headers, either default or provided. 907 908 # 3. 909 # this completely sidesteps any checks provided by class properties 910 self._hf = hf 911 self._hi = hi 912 self._hs = hs 913 self.data = data 914 915 self._set_distances() 916 917 @property 918 def _header(self): 919 """ 920 Convenient read-only dictionary of non-null header array values. 921 922 Header value access should be through instance attributes. All header 923 names and string values are lowercase. Computed every time, so use 924 frugally. See class docstring for header descriptions. 925 926 """ 927 out = _io.header_arrays_to_dict(self._hf, self._hi, self._hs, 928 nulls=False) 929 return out 930 931 @property 932 def byteorder(self): 933 """ 934 The byte order of the underlying header/data arrays. 935 936 Raises SacError if array byte orders are inconsistent, even in the 937 case where '<' is your native order and byteorders look like 938 '<', '=', '='. 939 940 """ 941 try: 942 if self.data is None: 943 assert self._hf.dtype.byteorder == self._hi.dtype.byteorder 944 else: 945 assert self._hf.dtype.byteorder == self._hi.dtype.byteorder ==\ 946 self.data.dtype.byteorder 947 except AssertionError: 948 msg = 'Inconsistent header/data byteorders.' 949 raise SacError(msg) 950 951 bo = self._hf.dtype.byteorder 952 953 if bo == '=': 954 byteorder = sys.byteorder 955 elif bo == '<': 956 byteorder = 'little' 957 elif bo == '>': 958 byteorder = 'big' 959 960 return byteorder 961 962 # TODO: make a byteorder setter? 963 def _byteswap(self): 964 """ 965 Change the underlying byte order and dtype interpretation of the float, 966 int, and (if present) data arrays. 967 968 """ 969 try: 970 self._hf = self._hf.byteswap(True).newbyteorder('S') 971 self._hi = self._hi.byteswap(True).newbyteorder('S') 972 if self.data is not None: 973 self.data = self.data.byteswap(True).newbyteorder('S') 974 except Exception as e: 975 # if this fails, roll it back? 976 raise e 977 978 @property 979 def reftime(self): 980 """ 981 Get or set the SAC header reference time as a UTCDateTime instance. 982 983 reftime is not an attribute, but is constructed and dismantled each 984 time directly to/from the SAC "nz"-time fields. 985 986 Setting a new reftime shifts all non-null relative time headers 987 accordingly. It accepts a UTCDateTime object, from which time shifts 988 are calculated. 989 990 ..rubric:: notes 991 992 The reftime you supply will be robbed of its remainder microseconds, 993 which are then pushed into the relative time header shifts. This means 994 that the reftime you observe after you set it here may not exactly 995 match the reftime you supplied; it may be `remainder microseconds` 996 earlier. Nor will the iztype reference header value be exactly zero; 997 it will be equal to `remainder microseconds` (as seconds). 998 999 """ 1000 return _ut.get_sac_reftime(self._header) 1001 1002 @reftime.setter 1003 def reftime(self, new_reftime): 1004 try: 1005 old_reftime = self.reftime 1006 1007 # snap the new reftime to the most recent milliseconds 1008 # (subtract the leftover microseconds) 1009 ns = new_reftime.ns 1010 utc = UTCDateTime(ns=(ns - ns % 1000000)) 1011 1012 self.nzyear = utc.year 1013 self.nzjday = utc.julday 1014 self.nzhour = utc.hour 1015 self.nzmin = utc.minute 1016 self.nzsec = utc.second 1017 self.nzmsec = utc.microsecond / 1000 1018 1019 # get the float seconds between the old and new reftimes 1020 shift = old_reftime - utc 1021 1022 # shift the relative time headers 1023 self._allt(np.float32(shift)) 1024 1025 except AttributeError: 1026 msg = "New reference time must be an obspy.UTCDateTime instance." 1027 raise TypeError(msg) 1028 1029 # --------------------------- I/O METHODS --------------------------------- 1030 @classmethod 1031 def read(cls, source, headonly=False, ascii=False, byteorder=None, 1032 checksize=False, debug_strings=False, encoding='ASCII'): 1033 """ 1034 Construct an instance from a binary or ASCII file on disk. 1035 1036 :param source: Full path string for File-like object from a SAC binary 1037 file on disk. If it is an open File object, open 'rb'. 1038 :type source: str or file 1039 :param headonly: If headonly is True, only read the header arrays not 1040 the data array. 1041 :type headonly: bool 1042 :param ascii: If True, file is a SAC ASCII/Alphanumeric file. 1043 :type ascii: bool 1044 :param byteorder: If omitted or None, automatic byte-order checking is 1045 done, starting with native order. If byteorder is specified and 1046 incorrect, a :class:`SacIOError` is raised. Only valid for binary 1047 files. 1048 :type byteorder: str {'little', 'big'}, optional 1049 :param checksize: If True, check that the theoretical file size from 1050 the header matches the size on disk. Only valid for binary files. 1051 :type checksize: bool 1052 :param debug_strings: By default, non-ASCII and null-termination 1053 characters are removed from character header fields, and those 1054 beginning with '-12345' are considered unset. If True, they 1055 are instead passed without modification. Good for debugging. 1056 :type debug_strings: bool 1057 :param encoding: Encoding string that passes the user specified 1058 encoding scheme. 1059 :type encoding: str 1060 1061 :raises: :class:`SacIOError` if checksize failed, byteorder was wrong, 1062 or header arrays are wrong size. 1063 1064 .. rubric:: Example 1065 1066 >>> from obspy.core.util import get_example_file 1067 >>> from obspy.io.sac.util import SacInvalidContentError 1068 >>> file_ = get_example_file("test.sac") 1069 >>> sac = SACTrace.read(file_, headonly=True) 1070 >>> sac.data is None 1071 True 1072 >>> sac = SACTrace.read(file_, headonly=False) 1073 >>> sac.data # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE 1074 array([ -8.74227766e-08, -3.09016973e-01, -5.87785363e-01, 1075 -8.09017122e-01, -9.51056600e-01, -1.00000000e+00, 1076 -9.51056302e-01, -8.09016585e-01, -5.87784529e-01, 1077 ... 1078 8.09022486e-01, 9.51059461e-01, 1.00000000e+00, 1079 9.51053500e-01, 8.09011161e-01, 5.87777138e-01, 1080 3.09007347e-01], dtype=float32) 1081 1082 See also: :meth:`SACTrace.validate` 1083 1084 """ 1085 if ascii: 1086 hf, hi, hs, data = _io.read_sac_ascii(source, headonly=headonly) 1087 else: 1088 hf, hi, hs, data = _io.read_sac(source, headonly=headonly, 1089 byteorder=byteorder, 1090 checksize=checksize) 1091 if not debug_strings: 1092 for i, val in enumerate(hs): 1093 val = _ut._clean_str(val.decode(encoding, 'replace'), 1094 strip_whitespace=False) 1095 if val.startswith(native_str('-12345')): 1096 val = HD.SNULL 1097 hs[i] = val.encode(encoding, 'replace') 1098 1099 sac = cls._from_arrays(hf, hi, hs, data) 1100 if sac.dist is None: 1101 sac._set_distances() 1102 1103 return sac 1104 1105 def write(self, dest, headonly=False, ascii=False, byteorder=None, 1106 flush_headers=True): 1107 """ 1108 Write the header and (optionally) data arrays to a SAC binary file. 1109 1110 :param dest: Full path or File-like object to SAC binary file on disk. 1111 :type dest: str or file 1112 :param headonly: If headonly is True, only read the header arrays not 1113 the data array. 1114 :type headonly: bool 1115 :param ascii: If True, file is a SAC ASCII/Alphanumeric file. 1116 :type ascii: bool 1117 :param byteorder: If omitted or None, automatic byte-order checking is 1118 done, starting with native order. If byteorder is specified and 1119 incorrect, a :class:`SacIOError` is raised. Only valid for binary 1120 files. 1121 :type byteorder: str {'little', 'big'}, optional 1122 :param flush_headers: If True, update data headers like 'depmin' and 1123 'depmax' with values from the data array. 1124 :type flush_headers: bool 1125 1126 """ 1127 if headonly: 1128 data = None 1129 else: 1130 # do a check for float32 data here instead of arrayio.write_sac? 1131 data = self.data 1132 if flush_headers: 1133 self._flush_headers() 1134 1135 if ascii: 1136 _io.write_sac_ascii(dest, self._hf, self._hi, self._hs, data) 1137 else: 1138 byteorder = byteorder or self.byteorder 1139 _io.write_sac(dest, self._hf, self._hi, self._hs, data, 1140 byteorder=byteorder) 1141 1142 @classmethod 1143 def _from_arrays(cls, hf=None, hi=None, hs=None, data=None): 1144 """ 1145 Low-level array-based constructor. 1146 1147 This constructor is good for getting a "blank" SAC object, and is used 1148 in other, perhaps more useful, alternate constructors ("See Also"). 1149 No value checking is done and header values are completely overwritten 1150 with the provided arrays, which is why this is a hidden constructor. 1151 1152 :param hf: SAC float header array 1153 :type hf: :class:`numpy.ndarray` of floats 1154 :param hi: SAC int header array 1155 :type hi: :class:`numpy.ndarray` of ints 1156 :param hs: SAC string header array 1157 :type hs: :class:`numpy.ndarray` of str 1158 :param data: SAC data array, optional. 1159 1160 If omitted or None, the header arrays are intialized according to 1161 :func:`arrayio.init_header_arrays`. If data is omitted, it is 1162 simply set to None on the corresponding :class:`SACTrace`. 1163 1164 .. rubric:: Example 1165 1166 >>> sac = SACTrace._from_arrays() 1167 >>> print(sac) # doctest: +NORMALIZE_WHITESPACE +SKIP 1168 Reference Time = XX/XX/XX (XXX) XX:XX:XX.XXXXXX 1169 iztype not set 1170 lcalda = True 1171 leven = False 1172 lovrok = False 1173 lpspol = False 1174 1175 """ 1176 # use the first byteorder we find, or system byteorder if we 1177 # never find any 1178 bo = '=' 1179 for arr in (hf, hi, hs, data): 1180 try: 1181 bo = arr.dtype.byteorder 1182 break 1183 except AttributeError: 1184 # arr is None (not supplied) 1185 pass 1186 hf0, hi0, hs0 = _io.init_header_arrays(byteorder=bo) 1187 # TODO: hf0, hi0, hs0 = _io.init_header_array_values(hf0, hi0, hs0) 1188 1189 if hf is None: 1190 hf = hf0 1191 if hi is None: 1192 hi = hi0 1193 if hs is None: 1194 hs = hs0 1195 1196 # get the default instance, but completely replace the arrays 1197 # initializes arrays twice, but it beats converting empty arrays to a 1198 # dict and then passing it to __init__, i think...maybe... 1199 sac = cls() 1200 sac._hf = hf 1201 sac._hi = hi 1202 sac._hs = hs 1203 sac.data = data 1204 1205 return sac 1206 1207 # TO/FROM OBSPY TRACES 1208 @classmethod 1209 def from_obspy_trace(cls, trace, keep_sac_header=True): 1210 """ 1211 Construct an instance from an ObsPy Trace. 1212 1213 :param trace: Source Trace object 1214 :type trace: :class:`~obspy.core.Trace` instance 1215 :param keep_sac_header: If True, any old stats.sac header values are 1216 kept as is, and only a minimal set of values are updated from the 1217 stats dictionary: npts, e, and data. If an old iztype and a valid 1218 reftime are present, b and e will be properly referenced to it. If 1219 False, a new SAC header is constructed from only information found 1220 in the stats dictionary, with some other default values introduced. 1221 :type keep_sac_header: bool 1222 1223 """ 1224 header = _ut.obspy_to_sac_header(trace.stats, keep_sac_header) 1225 1226 # handle the data headers 1227 data = trace.data 1228 try: 1229 if len(data) == 0: 1230 # data is a empty numpy.array 1231 data = None 1232 except TypeError: 1233 # data is None 1234 data = None 1235 1236 try: 1237 byteorder = data.dtype.byteorder 1238 except AttributeError: 1239 # data is None 1240 byteorder = '=' 1241 1242 hf, hi, hs = _io.dict_to_header_arrays(header, byteorder=byteorder) 1243 sac = cls._from_arrays(hf, hi, hs, data) 1244 # sac._flush_headers() 1245 1246 return sac 1247 1248 def to_obspy_trace(self, debug_headers=False, encoding='ASCII'): 1249 """ 1250 Return an ObsPy Trace instance. 1251 1252 Required headers: nz-time fields, npts, delta, calib, kcmpnm, kstnm, 1253 ...? 1254 1255 :param debug_headers: Include _all_ SAC headers into the 1256 Trace.stats.sac dictionary. 1257 :type debug_headers: bool 1258 :param encoding: Encoding string that passes the user specified 1259 encoding scheme. 1260 :type encoding: str 1261 1262 .. rubric:: Example 1263 1264 >>> from obspy.core.util import get_example_file 1265 >>> file_ = get_example_file("test.sac") 1266 >>> sac = SACTrace.read(file_, headonly=True) 1267 >>> tr = sac.to_obspy_trace() 1268 >>> print(tr) # doctest: +ELLIPSIS 1269 .STA..Q | 1978-07-18T08:00:10.000000Z - ... | 1.0 Hz, 100 samples 1270 1271 """ 1272 # make the obspy test for tests/data/testxy.sac pass 1273 # ObsPy does not require a valid reftime 1274 # try: 1275 # self.validate('reftime') 1276 # except SacInvalidContentError: 1277 # if not self.nzyear: 1278 # self.nzyear = 1970 1279 # if not self.nzjday: 1280 # self.nzjday = 1 1281 # for hdr in ['nzhour', 'nzmin', 'nzsec', 'nzmsec']: 1282 # if not getattr(self, hdr): 1283 # setattr(self, hdr, 0) 1284 self.validate('delta') 1285 if self.data is None: 1286 # headonly is True 1287 # Make it something palatable to ObsPy 1288 data = np.array([], dtype=self._hf.dtype.byteorder + 'f4') 1289 else: 1290 data = self.data 1291 1292 sachdr = _io.header_arrays_to_dict(self._hf, self._hi, self._hs, 1293 nulls=debug_headers, 1294 encoding=encoding) 1295 # TODO: logic to use debug_headers for real 1296 1297 stats = _ut.sac_to_obspy_header(sachdr) 1298 1299 return Trace(data=data, header=stats) 1300 1301 # ---------------------- other properties/methods ------------------------- 1302 def validate(self, *tests): 1303 """ 1304 Check validity of loaded SAC file content, such as header/data 1305 consistency. 1306 1307 :param tests: One or more of the following validity tests: 1308 'delta' : Time step "delta" is positive. 1309 'logicals' : Logical values are 0, 1, or null 1310 'data_hdrs' : Length, min, mean, max of data array match header 1311 values. 1312 'enums' : Check validity of enumerated values. 1313 'reftime' : Reference time values in header are all set. 1314 'reltime' : Relative time values in header are absolutely 1315 referenced. 1316 'all' : Do all tests. 1317 :type tests: str 1318 1319 :raises: :class:`SacInvalidContentError` if any of the specified tests 1320 fail. :class:`ValueError` if 'data_hdrs' is specified and data is 1321 None, empty array, or no tests specified. 1322 1323 .. rubric:: Example 1324 1325 >>> from obspy.core.util import get_example_file 1326 >>> from obspy.io.sac.util import SacInvalidContentError 1327 >>> file_ = get_example_file("LMOW.BHE.SAC") 1328 >>> sac = SACTrace.read(file_) 1329 >>> # make the time step invalid, catch it, and fix it 1330 >>> sac.delta *= -1.0 1331 >>> try: 1332 ... sac.validate('delta') 1333 ... except SacInvalidContentError as e: 1334 ... sac.delta *= -1.0 1335 ... sac.validate('delta') 1336 >>> # make the data and depmin/men/max not match, catch the validation 1337 >>> # error, then fix (flush) the headers so that they validate 1338 >>> sac.data += 5.0 1339 >>> try: 1340 ... sac.validate('data_hdrs') 1341 ... except SacInvalidContentError: 1342 ... sac._flush_headers() 1343 ... sac.validate('data_hdrs') 1344 1345 """ 1346 _io.validate_sac_content(self._hf, self._hi, self._hs, self.data, 1347 *tests) 1348 1349 def _format_header_str(self, hdrlist='all'): 1350 """ 1351 Produce a print-friendly string of header values for __repr__ , 1352 .listhdr(), and .lh() 1353 1354 """ 1355 # interpret hdrlist 1356 if hdrlist == 'all': 1357 hdrlist = sorted(self._header.keys()) 1358 elif hdrlist == 'picks': 1359 hdrlist = ('a', 'b', 'e', 'f', 'o', 't0', 't1', 't2', 't3', 't4', 1360 't5', 't6', 't7', 't8', 't9') 1361 else: 1362 msg = "Unrecognized hdrlist '{}'".format(hdrlist) 1363 raise ValueError(msg) 1364 1365 # start building header string 1366 # 1367 # reference time 1368 header_str = [] 1369 try: 1370 timefmt = "Reference Time = %m/%d/%Y (%j) %H:%M:%S.%f" 1371 header_str.append(self.reftime.strftime(timefmt)) 1372 except (ValueError, SacError): 1373 msg = "Reference time information incomplete." 1374 warnings.warn(msg) 1375 notime_str = "Reference Time = XX/XX/XX (XXX) XX:XX:XX.XXXXXX" 1376 header_str.append(notime_str) 1377 # 1378 # reftime type 1379 # TODO: use enumerated value dict here? 1380 iztype = self.iztype 1381 if iztype is None: 1382 header_str.append("\tiztype not set") 1383 elif iztype == 'ib': 1384 header_str.append("\tiztype IB: begin time") 1385 elif iztype == 'io': 1386 header_str.append("\tiztype IO: origin time") 1387 elif iztype == 'ia': 1388 header_str.append("\tiztype IA: first arrival time") 1389 elif iztype[1] == 't': 1390 vals = (iztype.upper(), iztype[1:]) 1391 izfmt = "\tiztype {}: user-defined time {}" 1392 header_str.append(izfmt.format(*vals)) 1393 elif iztype == 'iunkn': 1394 header_str.append("\tiztype IUNKN (Unknown)") 1395 else: 1396 header_str.append("\tunrecognized iztype: {}".format(iztype)) 1397 # 1398 # non-null headers 1399 hdrfmt = "{:10.10s} = {}" 1400 for hdr in hdrlist: 1401 # XXX: non-null header values might have no property for getattr 1402 try: 1403 header_str.append(hdrfmt.format(hdr, getattr(self, hdr))) 1404 except AttributeError: 1405 header_str.append(hdrfmt.format(hdr, self._header[hdr])) 1406 1407 return '\n'.join(header_str) 1408 1409 def listhdr(self, hdrlist='all'): 1410 """ 1411 Print header values. 1412 1413 Default is all non-null values. 1414 1415 :param hdrlist: Which header fields to you want to list. Choose one of 1416 {'all', 'picks'} or iterable of header fields. An iterable of 1417 header fields can look like 'bea' or ('b', 'e', 'a'). 1418 1419 'all' (default) prints all non-null values. 1420 'picks' prints fields which are used to define time picks. 1421 1422 An iterable of header fields can look like 'bea' or ('b', 'e', 'a'). 1423 1424 .. rubric:: Example 1425 1426 >>> from obspy.core.util import get_example_file 1427 >>> file_ = get_example_file("LMOW.BHE.SAC") 1428 >>> sac = SACTrace.read(file_) 1429 >>> sac.lh() # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS 1430 Reference Time = 04/10/2001 (100) 00:23:00.465000 1431 iztype IB: begin time 1432 a = 0.0 1433 b = 0.0 1434 delta = 0.009999999... 1435 depmax = 0.003305610... 1436 depmen = 0.00243799... 1437 depmin = 0.00148824... 1438 e = 0.98999997... 1439 iftype = itime 1440 iztype = ib 1441 kcmpnm = BHE 1442 kevnm = None 1443 kstnm = LMOW 1444 lcalda = True 1445 leven = True 1446 lpspol = False 1447 nevid = 0 1448 norid = 0 1449 npts = 100 1450 nvhdr = 6 1451 nzhour = 0 1452 nzjday = 100 1453 nzmin = 23 1454 nzmsec = 465 1455 nzsec = 0 1456 nzyear = 2001 1457 stla = -39.409999... 1458 stlo = 175.75 1459 unused23 = 0 1460 """ 1461 # https://ds.iris.edu/files/sac-manual/commands/listhdr.html 1462 print(self._format_header_str(hdrlist)) 1463 1464 def lh(self, *args, **kwargs): 1465 """Alias of listhdr method.""" 1466 self.listhdr(*args, **kwargs) 1467 1468 def __str__(self): 1469 return self._format_header_str() 1470 1471 def __repr__(self): 1472 # XXX: run self._flush_headers first? 1473 # TODO: make this somehow more readable. 1474 h = sorted(self._header.items()) 1475 fmt = ", {}={!r}" * len(h) 1476 argstr = fmt.format(*chain.from_iterable(h))[2:] 1477 return self.__class__.__name__ + "(" + argstr + ")" 1478 1479 def copy(self): 1480 return deepcopy(self) 1481 1482 def _flush_headers(self): 1483 """ 1484 Flush to the header arrays any header property values that may not be 1485 reflected there, such as data min/max/mean, npts, e. 1486 1487 """ 1488 # XXX: do I really care which byte order it is? 1489 # self.data = np.require(self.data, native_str('<f4')) 1490 self._hi[HD.INTHDRS.index('npts')] = self.npts 1491 self._hf[HD.FLOATHDRS.index('e')] = self.e 1492 self._hf[HD.FLOATHDRS.index('depmin')] = self.depmin 1493 self._hf[HD.FLOATHDRS.index('depmax')] = self.depmax 1494 self._hf[HD.FLOATHDRS.index('depmen')] = self.depmen 1495 1496 def _allt(self, shift): 1497 """ 1498 Shift all relative time headers by some value (addition). 1499 1500 Similar to SAC's "chnhdr allt". 1501 1502 Note 1503 ---- 1504 This method is triggered by setting an instance's iztype or changing 1505 its reference time, which is the most likely use case for this 1506 functionality. If what you're trying to do is set an origin time and 1507 make a file origin-based: 1508 1509 SAC> CHNHDR O GMT 1982 123 13 37 10 103 1510 SAC> LISTHDR O 1511 O 123.103 1512 SAC> CHNHDR ALLT -123.103 IZTYPE IO 1513 1514 ...it is recommended to just make sure your target reference header is 1515 set and correct, and set the iztype: 1516 1517 >>> from obspy import UTCDateTime 1518 >>> from obspy.core.util import get_example_file 1519 >>> file_ = get_example_file("test.sac") 1520 >>> sac = SACTrace.read(file_) 1521 >>> sac.o = UTCDateTime(year=1982, julday=123, 1522 ... hour=13, minute=37, 1523 ... second=10, microsecond=103) 1524 >>> sac.iztype = 'io' 1525 1526 The iztype setter will deal with shifting the time values. 1527 1528 """ 1529 for hdr in ['b', 'o', 'a', 'f'] + ['t' + str(i) for i in range(10)]: 1530 val = getattr(self, hdr) 1531 if val is not None: 1532 setattr(self, hdr, val + shift) 1533 1534 def _set_distances(self, force=False): 1535 """ 1536 Calculate dist, az, baz, gcarc. If force=True, ignore lcalda. 1537 Raises SacHeaderError if force=True and geographic headers are unset. 1538 1539 """ 1540 if self.lcalda or force: 1541 try: 1542 m, az, baz = gps2dist_azimuth(self.evla, self.evlo, self.stla, 1543 self.stlo) 1544 dist = m / 1000.0 1545 gcarc = kilometer2degrees(dist) 1546 self._hf[HD.FLOATHDRS.index('az')] = az 1547 self._hf[HD.FLOATHDRS.index('baz')] = baz 1548 self._hf[HD.FLOATHDRS.index('dist')] = dist 1549 self._hf[HD.FLOATHDRS.index('gcarc')] = gcarc 1550 except (ValueError, TypeError): 1551 # one or more of the geographic values is None 1552 if force: 1553 msg = ("Not enough information to calculate distance, " 1554 "azimuth.") 1555 raise SacHeaderError(msg) 1556