Source code for obspy.io.sac.sacpz

"""
obspy.io.sac.sacpz - SACPZ write support for ObsPy
==================================================

Module for SAC poles and zero (SACPZ) file I/O.

:copyright:
    The ObsPy Development Team (devs@obspy.org)
:license:
    GNU Lesser General Public License, Version 3
    (https://www.gnu.org/copyleft/lesser.html)
"""
import numpy as np
import warnings
from copy import deepcopy
from pathlib import Path

from obspy import UTCDateTime
from obspy.core import AttribDict
from obspy.core.inventory.response import paz_to_sacpz_string


[docs] def _write_sacpz(inventory, file_or_file_object): """ Writes an inventory object to a SACPZ file. .. note:: The channel DIP in the SACPZ comment fields is specified like defined by SEED (positive down from horizontal). .. warning:: This function should NOT be called directly, it registers via the the :meth:`~obspy.core.inventory.inventory.Inventory.write` method of an ObsPy :class:`~obspy.core.inventory.inventory.Inventory` object, call this instead. :type inventory: :class:`~obspy.core.inventory.inventory.Inventory` :param inventory: The inventory instance to be written. :param file_or_file_object: The file or file-like object to be written to. """ out = [] now = UTCDateTime() for net in inventory: for sta in net: for cha in sta: resp = cha.response sens = resp.instrument_sensitivity try: paz = resp.get_paz() except Exception: msg = ( f"{net.code}.{sta.code}." f"{cha.location_code}.{cha.code} " f"{cha.start_date} has no paz. Skipping." ) warnings.warn(msg) continue input_unit = sens.input_units.upper() # convert hertz to radians/s, see #3334 paz = deepcopy(paz) paz.to_radians_per_second() # SAC seems to expect "meters" as the unit if input_unit == "M": pass elif input_unit in ["M/S", "M/SEC"]: paz.zeros.append(0j) elif input_unit in ["M/S**2", "M/SEC**2"]: paz.zeros.extend([0j, 0j]) else: msg = ( f"{net.code}.{sta.code}." f"{cha.location_code}.{cha.code} " f"{cha.start_date} " "has unrecognized input units in " f"response: {input_unit}. Skipping" ) warnings.warn(msg) continue out.append("* " + "*" * 50) out.append(f"* NETWORK : {net.code}") out.append(f"* STATION : {sta.code}") out.append(f"* LOCATION : {cha.location_code}") out.append(f"* CHANNEL : {cha.code}") out.append(f"* CREATED : {now}") out.append(f"* START : {cha.start_date}") out.append(f"* END : {cha.end_date}") out.append(f"* DESCRIPTION : {sta.site.name}") out.append(f"* LATITUDE : {cha.latitude or sta.latitude}") out.append(f"* LONGITUDE : {cha.longitude or sta.longitude}") out.append(f"* ELEVATION : {cha.elevation or sta.elevation}") out.append(f"* DEPTH : {cha.depth}") # DIP in SACPZ served by EarthScope/IRIS SACPZ web service is # systematically different from the StationXML entries. # It is defined as an incidence angle (like done in SAC for # sensor orientation), rather than an actual dip. # Add '(SEED)' to clarify that we adhere to SEED convention out.append(f"* DIP (SEED) : {cha.dip}") out.append(f"* AZIMUTH : {cha.azimuth}") out.append(f"* SAMPLE RATE : {cha.sample_rate}") out.append("* INPUT UNIT : M") out.append(f"* OUTPUT UNIT : {sens.output_units}") inst_type = cha.sensor and cha.sensor.type or '' out.append(f"* INSTTYPE : {inst_type}") out.append(f"* INSTGAIN : {paz.stage_gain} " f"({sens.input_units})") out.append(f"* SENSITIVITY : {sens.value} " f"({sens.input_units})") out.append(f"* A0 : {paz.normalization_factor}") out.append("* " + "*" * 50) out.append(paz_to_sacpz_string(paz, sens)) out.extend(["", ""]) out = "\n".join(out) + "\n\n" try: file_or_file_object.write(out) except AttributeError: with open(file_or_file_object, "wt") as fh: fh.write(out)
# UTILITIES
[docs] def attach_paz(tr, paz_file, todisp=False, tovel=False, torad=False, tohz=False): ''' Attach tr.stats.paz AttribDict to trace from SAC paz_file This is experimental code, taken from obspy.io.gse2.libgse2.attach_paz and adapted to the SAC-pole-zero conventions. Especially the conversion from velocity to displacement and vice versa is still under construction. It works but I cannot guarantee that the values are correct. For more information on the SAC-pole-zero format see: https://ds.iris.edu/files/sac-manual/commands/transfer.html. For a useful discussion on polezero files and transfer functions in general see: http://seis-uk.le.ac.uk/equipment/downloads/data_management/\ seisuk_instrument_resp_removal.pdf Also bear in mind that according to the SAC convention for pole-zero files CONSTANT is defined as: digitizer_gain*seismometer_gain*A0. This means that it does not have explicit information on the digitizer gain and seismometer gain which we therefore set to 1.0. Attaches to a trace a paz AttribDict containing poles zeros and gain. :param tr: An ObsPy :class:`~obspy.core.trace.Trace` object :param paz_file: path to pazfile or file pointer :param todisp: change a velocity transfer function to a displacement transfer function by adding another zero :param tovel: change a displacement transfer function to a velocity transfer function by removing one 0,0j zero :param torad: change to radians :param tohz: change to Hertz >>> from obspy import Trace >>> import io >>> tr = Trace() >>> f = io.StringIO("""ZEROS 3 ... -5.032 0.0 ... POLES 6 ... -0.02365 0.02365 ... -0.02365 -0.02365 ... -39.3011 0. ... -7.74904 0. ... -53.5979 21.7494 ... -53.5979 -21.7494 ... CONSTANT 2.16e18""") >>> attach_paz(tr, f,torad=True) >>> for z in tr.stats.paz['zeros']: ... print("%.2f %.2f" % (z.real, z.imag)) -31.62 0.00 0.00 0.00 0.00 0.00 ''' poles = [] zeros = [] if isinstance(paz_file, (str, Path)): paz_file = open(paz_file, 'r') is_filename = True else: is_filename = False try: while True: line = paz_file.readline() if not line: break # lines starting with * are comments if line.startswith('*'): continue if line.find('ZEROS') != -1: a = line.split() noz = int(a[1]) for _k in range(noz): line = paz_file.readline() a = line.split() if line.find('POLES') != -1 or \ line.find('CONSTANT') != -1 or \ line.startswith('*') or not line: while len(zeros) < noz: zeros.append(complex(0, 0j)) break else: zeros.append(complex(float(a[0]), float(a[1]))) if line.find('POLES') != -1: a = line.split() nop = int(a[1]) for _k in range(nop): line = paz_file.readline() a = line.split() if line.find('CONSTANT') != -1 or \ line.find('ZEROS') != -1 or \ line.startswith('*') or not line: while len(poles) < nop: poles.append(complex(0, 0j)) break else: poles.append(complex(float(a[0]), float(a[1]))) if line.find('CONSTANT') != -1: a = line.split() # in the observatory this is the seismometer gain [muVolt/nm/s] # the A0_normalization_factor is hardcoded to 1.0 constant = float(a[1]) finally: if is_filename: paz_file.close() # To convert the velocity response to the displacement response, # multiplication with jw is used. This is equivalent to one more # zero in the pole-zero representation if todisp: zeros.append(complex(0, 0j)) # To convert the displacement response to the velocity response, # division with jw is used. This is equivalent to one less zero # in the pole-zero representation if tovel: for i, zero in enumerate(list(zeros)): if zero == complex(0, 0j): zeros.pop(i) break else: raise Exception("Could not remove (0,0j) zero to change \ displacement response to velocity response") # convert poles, zeros and gain in Hertz to radians if torad: tmp = [z * 2. * np.pi for z in zeros] zeros = tmp tmp = [p * 2. * np.pi for p in poles] poles = tmp # When extracting RESP files and SAC_PZ files # from a dataless SEED using the rdseed program # where the former is in Hz and the latter in radians, # there gains seem to be unaffected by this. # According to this document: # http://www.le.ac.uk/ # seis-uk/downloads/seisuk_instrument_resp_removal.pdf # the gain should also be converted when changing from # hertz to radians or vice versa. However, the rdseed programs # does not do this. I'm not entirely sure at this stage which one is # correct or if I have missed something. I've therefore decided # to leave it out for now, in order to stay compatible with the # rdseed program and the SAC program. # constant *= (2. * np.pi) ** 3 # convert poles, zeros and gain in radian to Hertz if tohz: for i, z in enumerate(zeros): if abs(z) > 0.0: zeros[i] /= 2 * np.pi for i, p in enumerate(poles): if abs(p) > 0.0: poles[i] /= 2 * np.pi # constant /= (2. * np.pi) ** 3 # fill up ObsPy Poles and Zeros AttribDict # In SAC pole-zero files CONSTANT is defined as: # digitizer_gain*seismometer_gain*A0 tr.stats.paz = AttribDict() tr.stats.paz.seismometer_gain = 1.0 tr.stats.paz.digitizer_gain = 1.0 tr.stats.paz.poles = poles tr.stats.paz.zeros = zeros # taken from obspy.io.gse2.paz:145 tr.stats.paz.sensitivity = tr.stats.paz.digitizer_gain * \ tr.stats.paz.seismometer_gain tr.stats.paz.gain = constant
[docs] def attach_resp(tr, resp_file, todisp=False, tovel=False, torad=False, tohz=False): """ Extract key instrument response information from a RESP file, which can be extracted from a dataless SEED volume by, for example, using the script obspy-dataless2resp or the rdseed program. At the moment, you have to determine yourself if the given response is for velocity or displacement and if the values are given in rad or Hz. This is still experimental code (see also documentation for :func:`obspy.io.sac.sacpz.attach_paz`). Attaches to a trace a paz AttribDict containing poles, zeros, and gain. :param tr: An ObsPy :class:`~obspy.core.trace.Trace` object :param resp_file: path to RESP-file or file pointer :param todisp: change a velocity transfer function to a displacement transfer function by adding another zero :param tovel: change a displacement transfer function to a velocity transfer function by removing one 0,0j zero :param torad: change to radians :param tohz: change to Hertz >>> from obspy import Trace >>> import os >>> tr = Trace() >>> respfile = os.path.join(os.path.dirname(__file__), 'tests', 'data', ... 'RESP.NZ.CRLZ.10.HHZ') >>> attach_resp(tr, respfile, torad=True, todisp=False) >>> for k in sorted(tr.stats.paz): # doctest: +NORMALIZE_WHITESPACE ... print(k) digitizer_gain gain poles seismometer_gain sensitivity t_shift zeros >>> print(tr.stats.paz.poles) # doctest: +NORMALIZE_WHITESPACE +ELLIPSIS [(-0.15931...+0.15931...j), (-0.15931...-0.15931...j), (-314.159...+202.318...j), (-314.159...-202.318...j)] """ if not hasattr(resp_file, 'write'): resp_filep = open(resp_file, 'r') else: resp_filep = resp_file zeros_pat = r'B053F10-13' poles_pat = r'B053F15-18' a0_pat = r'B053F07' sens_pat = r'B058F04' t_shift_pat = r'B057F08' t_shift = 0.0 poles = [] zeros = [] while True: line = resp_filep.readline() if not line: break if line.startswith(a0_pat): a0 = float(line.split(':')[1]) if line.startswith(sens_pat): sens = float(line.split(':')[1]) if line.startswith(poles_pat): tmp = line.split() poles.append(complex(float(tmp[2]), float(tmp[3]))) if line.startswith(zeros_pat): tmp = line.split() zeros.append(complex(float(tmp[2]), float(tmp[3]))) if line.startswith(t_shift_pat): t_shift += float(line.split(':')[1]) constant = a0 * sens if not hasattr(resp_file, 'write'): resp_filep.close() if torad: tmp = [z * 2. * np.pi for z in zeros] zeros = tmp tmp = [p * 2. * np.pi for p in poles] poles = tmp if todisp: zeros.append(complex(0, 0j)) # To convert the displacement response to the velocity response, # division with jw is used. This is equivalent to one less zero # in the pole-zero representation if tovel: for i, zero in enumerate(list(zeros)): if zero == complex(0, 0j): zeros.pop(i) break else: raise Exception("Could not remove (0,0j) zero to change \ displacement response to velocity response") # convert poles, zeros and gain in radian to Hertz if tohz: for i, z in enumerate(zeros): if abs(z) > 0.0: zeros[i] /= 2 * np.pi for i, p in enumerate(poles): if abs(p) > 0.0: poles[i] /= 2 * np.pi constant /= (2. * np.pi) ** 3 # fill up ObsPy Poles and Zeros AttribDict # In SAC pole-zero files CONSTANT is defined as: # digitizer_gain*seismometer_gain*A0 tr.stats.paz = AttribDict() tr.stats.paz.seismometer_gain = sens tr.stats.paz.digitizer_gain = 1.0 tr.stats.paz.poles = poles tr.stats.paz.zeros = zeros # taken from obspy.io.gse2.paz:145 tr.stats.paz.sensitivity = tr.stats.paz.digitizer_gain * \ tr.stats.paz.seismometer_gain tr.stats.paz.gain = constant tr.stats.paz.t_shift = t_shift
if __name__ == "__main__": import doctest doctest.testmod()