Hi all,
When I use ppsd() (PPSD) for calculating the probabilistic power spectral densities of seismic record, ObsPy raise a warning,
/home/junliu/anaconda3/envs/obspy/lib/python3.10/site-packages/obspy/signal/spectral_estimation.py:886: UserWarning: Already covered time spans detected (e.g. 2024-04-16T00:00:10.000000Z), skipping these slices.
warnings.warn(msg)
and there is no output generated. Here is the code snippet used for processing:
from obspy.core import read
from obspy.signal import PPSD
st = read("test.sac")
paz = {'gain': 60077000.0,
'poles': [-0.037004+0.037016j, -0.037004-0.037016j,
-251.33+0j, -131.04-467.29j, -131.04+467.29j],
'sensitivity': 2516778400.0,
'zeros': [0j, 0j]}
ppsd = PPSD(st[0].stats, paz, ppsd_length=10)
# we use 10 second to crash it quickly
ppsd.add(st[0])
Here is the demo SAC file.
test.sac (98.3 KB)
Based on my observation, this issue seems to stem from potential rounding errors when reading the SAC file. st[0].stats indicates a sampling rate of 249.99998474121094. Upon inspecting obspy/signal/spectral_estimation.py, line 733 returns True in the condition:
if overlap_seconds / self.ppsd_length > self.overlap:
In this case, overlap_seconds / self.ppsd_length evaluates to 0.500000512, which is slightly greater than 0.5. Resampling the trace resolves the issue:
st[0].resample(250.00)
ppsd = PPSD(st[0].stats, paz, ppsd_length=10)
ppsd.add(st[0])
ppsd() works pretty well.
I found that the sampling rate issue has been fixed SAC: Potential problem with floating point accuracy in sampling rate / sample spacing · Issue #3408 · obspy/obspy (github.com), and I just wondered if it is needed for modifing line 733 (such as raise the right side) to avoid this problem.