Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SACPZ: Make sure PAZ gets converted to radians/s before writing output #3434

Merged
merged 7 commits into from
Apr 30, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
3 changes: 3 additions & 0 deletions CHANGELOG.txt
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,9 @@ Changes:
- obspy.io.gcf:
* Fixed an issue in algorithm to split and encode last few data into GCF
blocks (see #3252)
- obspy.io.sac:
* Fix writing SACPZ from poles and zeros stages with Hertz transfer function
type. SAC is expecting the SACPZ data to be in radians/s (see #3334)
- obspy.io.wav:
* Fixed reading of stereo wav files (see #3399)
- obspy.io.win:
Expand Down
44 changes: 43 additions & 1 deletion obspy/core/inventory/response.py
Original file line number Diff line number Diff line change
Expand Up @@ -319,6 +319,35 @@ def pz_transfer_function_type(self, value):
else:
raise ValueError(msg)

def to_radians_per_second(self):
"""
Convert to type 'LAPLACE (RADIANS/SECOND)'
"""
if self.pz_transfer_function_type == 'LAPLACE (RADIANS/SECOND)':
return
elif self.pz_transfer_function_type == 'LAPLACE (HERTZ)':
twopi = 2 * pi
self.normalization_factor *= twopi ** (
len(self.poles) - len(self.zeros))
self.poles = [
ComplexWithUncertainties(
x.real * twopi, x.imag * twopi,
upper_uncertainty=x.upper_uncertainty * twopi,
lower_uncertainty=x.lower_uncertainty * twopi)
for x in self.poles]
self.zeros = [
ComplexWithUncertainties(
x.real * twopi, x.imag * twopi,
upper_uncertainty=x.upper_uncertainty * twopi,
lower_uncertainty=x.lower_uncertainty * twopi)
for x in self.zeros]
self.pz_transfer_function_type = 'LAPLACE (RADIANS/SECOND)'
else:
msg = (f"Can not convert transfer function type "
f"'{self.pz_transfer_function_type}' to "
f"'LAPLACE (RADIANS/SECOND)'")
raise ValueError(msg)


class CoefficientsTypeResponseStage(ResponseStage):
"""
Expand Down Expand Up @@ -1994,17 +2023,30 @@ def from_paz(cls, zeros, poles, stage_gain,
return resp


def paz_to_sacpz_string(paz, instrument_sensitivity):
def paz_to_sacpz_string(paz, instrument_sensitivity,
to_radians_per_second=True):
"""
Returns SACPZ ASCII text representation of Response.

:type paz: :class:`PolesZerosResponseStage`
:param paz: Poles and Zeros response information
:type instrument_sensitivity: :class:`InstrumentSensitivity`
:param paz: Overall instrument sensitivity of response
:type to_radians_per_second: bool
:param to_radians_per_second: Whether to convert poles and zeros to
radians/s before returning SACPZ string. This should be done in most
cases as SAC is expecting radians/s.
:rtype: str
:returns: Textual SACPZ representation of poles and zeros response stage.
"""
if to_radians_per_second:
paz = copy.deepcopy(paz)
paz.to_radians_per_second()
if paz.pz_transfer_function_type != 'LAPLACE (RADIANS/SECOND)':
msg = ('Returning a SACPZ string from a PAZResponseStage that is not '
'radians/s. SAC is expecting SACPZ data to be in radians/s '
'(see #3334).')
warnings.warn(msg)
# assemble output string
out = []
out.append("ZEROS %i" % len(paz.zeros))
Expand Down