Skip to content

Commit

Permalink
ppsd: add a switch for hydrophone data, unify with special handling of
Browse files Browse the repository at this point in the history
rotational data (from ringloaser)

ppsd: trying to be nice to ppsd pickles with is_rotational_data
attributes

PPSD: better workaround for older PPSD pickle files
  • Loading branch information
megies committed Sep 7, 2015
1 parent 80c7db1 commit 719142d
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 13 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.txt
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,11 @@ master:
- new option `PPSD.plot(..., cumulative=True)` for a cumulative plot of
the histogram, i.e. a non-exceedence percentage visualization, similar
to the `percentile` option.
- changes to special handling of rotational: now handled by kwarg
`special_handling="ringlaser"` (kwarg `is_rotational_data` is
deprecated, see #916)
- special handling option for hydrophone data (no differentiation, see
#916)

0.10.x:
- obspy.fdsn:
Expand Down
52 changes: 39 additions & 13 deletions obspy/signal/spectral_estimation.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,6 @@
'delta',
'freq',
'id',
'is_rotational_data',
'len',
'location',
'merge_method',
Expand All @@ -84,6 +83,7 @@
'ppsd_length',
'sampling_rate',
'spec_bins',
'special_handling',
'station',
]

Expand Down Expand Up @@ -270,8 +270,8 @@ class PPSD(object):
@deprecated_keywords({'paz': 'metadata', 'parser': 'metadata',
'water_level': None})
def __init__(self, stats, metadata, skip_on_gaps=False,
is_rotational_data=False, db_bins=(-200, -50, 1.),
ppsd_length=3600., overlap=0.5, **kwargs):
db_bins=(-200, -50, 1.), ppsd_length=3600., overlap=0.5,
special_handling=None, **kwargs):
"""
Initialize the PPSD object setting all fixed information on the station
that should not change afterwards to guarantee consistent spectral
Expand All @@ -292,10 +292,8 @@ def __init__(self, stats, metadata, skip_on_gaps=False,
is changing over the timespans that are added to the PPSD.
Use with caution!
:note: When using `is_rotational_data=True` the applied processing
steps are changed (and it is assumed that a dictionary is
provided as `metadata`).
Differentiation of data (converting velocity
:note: When using `special_handling="ringlaser"` the applied processing
steps are changed. Differentiation of data (converting velocity
to acceleration data) will be omitted and a flat instrument
response is assumed, leaving away response removal and only
dividing by `metadata['sensitivity']` specified in the provided
Expand All @@ -317,9 +315,6 @@ def __init__(self, stats, metadata, skip_on_gaps=False,
`skip_on_gaps=True` for not filling gaps with zeros which might
result in some data segments shorter than `ppsd_length` not
used in the PPSD.
:type is_rotational_data: bool, optional
:param is_rotational_data: If set to True adapt processing of data to
rotational data. See note for details.
:type db_bins: tuple of three ints/floats
:param db_bins: Specify the lower and upper boundary and the width of
the db bins. The bin width might get adjusted to fit a number
Expand All @@ -334,15 +329,32 @@ def __init__(self, stats, metadata, skip_on_gaps=False,
values between 0 and 1 and is given as fraction of the length
of one segment, e.g. `ppsd_length=3600` and `overlap=0.5`
result in an overlap of 1800s of the segments.
:type special_handling: str, optional
:param special_handling: Switches on customized handling for
data other than seismometer recordings. Can be one of: 'ringlaser'
(no instrument correction, just division by
`metadata["sensitivity"]` of provided metadata dictionary),
'hydrophone' (no differentiation after instrument correction).
"""
# remove after release of 0.11.0
if kwargs.pop("is_rotational_data", None) is True:
msg = ("Keyword 'is_rotational_data' is deprecated. Please use "
"'special_handling=\"ringlaser\"' instead.")
warnings.warn(msg, ObsPyDeprecationWarning)
special_handling = "ringlaser"

self.id = "%(network)s.%(station)s.%(location)s.%(channel)s" % stats
self.network = stats.network
self.station = stats.station
self.location = stats.location
self.channel = stats.channel
self.sampling_rate = stats.sampling_rate
self.delta = 1.0 / self.sampling_rate
self.is_rotational_data = is_rotational_data
self.special_handling = special_handling and special_handling.lower()
if self.special_handling not in (None, "ringlaser", "hydrophone"):
msg = "Unsupported value for 'special_handling' parameter: %s"
msg = msg % self.special_handling
raise ValueError(msg)
self.ppsd_length = ppsd_length
self.overlap = overlap
# trace length for one segment
Expand Down Expand Up @@ -621,9 +633,11 @@ def __process(self, tr):
# Here we remove the response using the same conventions
# since the power is squared we want to square the sensitivity
# we can also convert to acceleration if we have non-rotational data
if self.is_rotational_data:
if self.special_handling == "ringlaser":
# in case of rotational data just remove sensitivity
spec /= self.metadata['sensitivity'] ** 2
# special_handling "hydrophone" does instrument correction same as
# "normal" data
else:
# determine instrument response from metadata
try:
Expand All @@ -644,7 +658,11 @@ def __process(self, tr):
w = 2.0 * math.pi * _freq[1:]
w = w[::-1]
# Here we do the response removal
spec = (w ** 2) * spec / respamp
# Do not differentiate when `special_handling="hydrophone"`
if self.special_handling == "hydrophone":
spec = spec / respamp
else:
spec = (w ** 2) * spec / respamp
# avoid calculating log of zero
idx = spec < dtiny
spec[idx] = dtiny
Expand Down Expand Up @@ -850,6 +868,14 @@ def load(filename):
with open(filename, 'rb') as file_:
ppsd = pickle.load(file_)

# some workarounds for older PPSD pickle files
if hasattr(ppsd, "is_rotational_data"):
if ppsd.is_rotational_data is True:
ppsd.special_handling = "ringlaser"
delattr(ppsd, "is_rotational_data")
if not hasattr(ppsd, "special_handling"):
ppsd.special_handling = None

return ppsd

def save_npz(self, filename):
Expand Down

0 comments on commit 719142d

Please sign in to comment.