SACPZ format output is not in radians per second convention

[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

[megies]

Sorry this got overlooked somehow.. nice catch and this should be an easy fix. I'm just confused right now with A0 in the example you posted. The StationXML has A0=100.295 and the correct SACPZ has A0=3.959488e+03. From memory, with 4 poles and 3 zeros shouldn't the difference also just be a factor of (2 * Pi) **(4-3)?

[chad-earthscope]

I believe it's because the conversion to radians/sec is done prior to converting from velocity to displacement (by adding a zero). The original response in m/s has 2 zeros and 4 poles, making the calculation:

In [1]: 100.295 * (2*3.14159)**(4-2)
Out[1]: 3959.4812047591577

Now, I do not readily have a good answer why that order of operations is correct. For what it's worth, reviewing the rdseed source code, it's done in that order, and likely why the service is following suit.

[megies]

Ah, I overlooked that it's also changed to displacement.. 👍
Hmm yeah, the calculation seems OK now, A0 seems correct like this. Does SACPZ always use meters / displacement? Just if you know right away, I guess I could research it. (If velocity is OK too, maybe that might it might be simpler to just leave the physical units alone and only do the radians/Hertz conversion?)

[chad-earthscope]

Does SACPZ always use meters / displacement?

Well the sac convention is that data are in nanometers / displacement, although I think it's common for folks to work in meters also and just skip converting to nanometers.

A quick search of the source code and it does not appear that sac reads the INPUT UNIT annotation. The help for the transfer command says:

The INPUT UNIT for any polezero produced by sacpz or evalresp will be "M"

but really, following the use of any SACPZ file for deconvolution, the data are in whatever units the SACPZ describe.

At the DMC we convert the time dimension to displacement and leave the length dimension as it is, but it is very commonly meters and so that becomes presumed unit.

[megies]

Should be fixed now, hope the simple fix covers all scenarios