Skip to content

SACPZ format output is not in radians per second convention #3334

@chad-earthscope

Description

@chad-earthscope

Avoid duplicates

  • I searched existing issues

Bug Summary

Poles and zeros are expected to be in radians per second by the sac program, and likely other programs that read these files. Furthermore, there is no annotation or indication in the format that would indicate any other units.

A simple conversion of StationXML, which allows the units to be either rad/s or Hertz, to "SACPZ" does not convert to rad/s. An example of such a conversion using a StationXML response with a PZ (Laplace) response stage in Hertz is as follows:

G_CAN__LHZ.xml.gz

import obspy
inv=obspy.read_inventory('G_CAN__LHZ.xml')
inv.write('G_CAN__LHZ.sacpz',format='SACPZ')
... UserWarning: More than one PolesZerosResponseStage encountered. Returning first one found.

Never mind the multiple PZ error, the initial set is what we need.

The first PZs are converted to "SACPZ" format, but with the PZs, A0 and CONSTANT in Hertz. There is no indication of units. Even if there were units noted, the sac program would not recognize them. Using this "SACPZ" file with sac, and possibly other programs, will result in completely incorrect transform with the user having no indication that anything is wrong.

I suggest that when writing "SACPZ" format that the units are consistently in radians per second to avoid confusion and quiet data corruption. A correct conversion from the test StationXML to SAC PZ (by irisws-sacpz) to rad/s would look like:

* **********************************
* NETWORK   (KNETWK): G
* STATION    (KSTNM): CAN
* LOCATION   (KHOLE): 
* CHANNEL   (KCMPNM): LHZ
* CREATED           : 2023-07-25T14:40:53
* START             : 1989-06-02T00:00:00
* END               : 2006-12-10T02:00:00
* DESCRIPTION       : Canberra, Australia
* LATITUDE          : -35.318715
* LONGITUDE         : 148.996325
* ELEVATION         : 700.0
* DEPTH             : 0.0
* DIP               : 0.0
* AZIMUTH           : 0.0
* SAMPLE RATE       : 1.0
* INPUT UNIT        : M
* OUTPUT UNIT       : COUNTS
* INSTTYPE          : STRECKEISEN STS1
* INSTGAIN          : 2.252000e+03 (M/S)
* COMMENT           : 
* SENSITIVITY       : 1.844840e+09 (M/S)
* A0                : 3.959488e+03
* **********************************
ZEROS	3
	+0.000000e+00	+0.000000e+00	
	+0.000000e+00	+0.000000e+00	
	+0.000000e+00	+0.000000e+00	
POLES	4
	-1.233948e-02	+1.234319e-02	
	-1.233948e-02	-1.234319e-02	
	-3.917566e+01	+4.912339e+01	
	-3.917566e+01	-4.912339e+01	
CONSTANT	7.304622e+12

There does appear to be some code to convert to rad/s, but only when attaching responses to traces. In that case the caller must know the units a priori to know if a conversion is needed. This is the case even when attaching PZs from RESP, which is weird as the RESP specifies the units.

Code to Reproduce

No response

Error Traceback

No response

ObsPy Version?

1.4.0

Operating System?

macOS

Python Version?

3.11.4

Installation Method?

None

Metadata

Metadata

Assignees

No one assigned

    Labels

    .ioissues generally related to our read/write plugins.io.sacbugconfirmed bug

    Type

    No type

    Projects

    No projects

    Milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions