PPSD errors

Hi All,

I am new to obspy and am trying to plot a PPSD of my miniseed data.

I keep getting an error saying no data accumulated.

I was able to plot the trace no problem.

The code I am currently running can be seen below ( I have added some unnecessary lines for trouble shooting).

from obspy import read
from obspy.signal import PPSD
st = read(“C:\Users\21352643\Desktop\S- Seismic Data\S- Garulp Sensor at South ITM\2017_04_07\01h00m_200_e.msd”)
tr = st.select(id=".CMG6…HHE")[0]
tr.plot()
paz = {‘gain’: 1983000, ‘poles’: [-0.02365+0.02365j, -0.02365-0.02365j, -393.011, -7.4904, -53.5979-21.7494j, -53.5979+21.7494j], ‘sensitivity’: 1259.73, ‘zeros’: [-5.03207, 0, 0]}
ppsd= PPSD(tr.stats, paz)
print(paz)
print(len(tr))
print(tr.stats)
ppsd.add(st)
#st = read(“C:\Users\21352643\Desktop\01h00m_200_e.msd”)
#print(ppsd.add(st))
print(ppsd.times_processed[:10])
print(ppsd.id)

#print(ppsd.id)
ppsd.plot()

The output of this is:


{‘zeros’: [-5.03207, 0, 0], ‘sensitivity’: 1259.73, ‘poles’: [(-0.02365+0.02365j), (-0.02365-0.02365j), -393.011, -7.4904, (-53.5979-21.7494j), (-53.5979+21.7494j)], ‘gain’: 1983000}
720000
network:
station: CMG6
location:
channel: HHE
starttime: 2017-04-07T01:00:00.000000Z
endtime: 2017-04-07T01:59:59.995000Z
sampling_rate: 200.0
delta: 0.005
npts: 720000
calib: 1.0
_format: MSEED
mseed: AttribDict({‘record_length’: 4096, ‘encoding’: u’STEIM1’, ‘filesize’: 1568768, u’dataquality’: u’D’, ‘number_of_records’: 383, ‘byteorder’: u’>’})
[]
.CMG6…HHE
Traceback (most recent call last):

File “”, line 1, in
runfile(‘C:/Users/21352643/untitled0.py’, wdir=‘C:/Users/21352643’)

File “C:\Users\21352643\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py”, line 866, in runfile
execfile(filename, namespace)

File “C:\Users\21352643\AppData\Local\Continuum\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py”, line 87, in execfile
exec(compile(scripttext, filename, ‘exec’), glob, loc)

File “C:/Users/21352643/untitled0.py”, line 18, in
ppsd.plot()

File “C:\Users\21352643\AppData\Local\Continuum\Anaconda2\envs\obspy\lib\site-packages\obspy\signal\spectral_estimation.py”, line 1541, in plot
self.__check_histogram()

File “C:\Users\21352643\AppData\Local\Continuum\Anaconda2\envs\obspy\lib\site-packages\obspy\signal\spectral_estimation.py”, line 758, in __check_histogram
raise Exception(msg)

Exception: No data accumulated

I have the calibration data for the sensor I am using (see attached).
Am I entering this correctly into the paz function?

Thanks all

Hi Joshua,

Please select data with large time windows. It takes at least one hour for a valid PPSD calculation.

Best,

Youngman

Hi Youngman,

Thanks for that. So if I add more data using ppsd.add until I have at least an hours worth I should be able to plot the PPSD?

Thanks

Joshua

Hi Joshua,

Yes. PSD is calculated once for 1-hour and continuous data. You could either import a large amount of data to calculate PPSD, or use several times of ppsd.add(st) to import discontinuous data (each data set must > 1-hour).

Best,

Youngman

Just a few more comments, most vital info was already given by
Youngman.. (thanks!)

The length of one psd piece can be controlled by setting "ppsd_length"
in PPSD() initialization
(http://docs.obspy.org/master/packages/autogen/obspy.signal.spectral_estimation.PPSD.html?highlight=ppsd_length),
in case you want to operate with shorter time slices (obviously this
will limit the frequency range of the PPSD at the long period end). The
default is one hour, as is used by PQLX and the paper by McNamara et al
and in general by most people.

In general, supplying one day streams to PPSD.add() is a good idea (as
opposed to using e.g. 1-hour streams, as PPSD is also using an overlap
when selecting time windows and the checks for which data is already
covered and will be skipped by PPSD are not ideal and might omit valid
data when using shorter time slices, see
https://github.com/obspy/obspy/pull/896)

T

Thanks for the additional information Tobias. I have decided that I will extend the data file capture length of the guralp sensor to 1 day. It is currently set to one hour but is giving me 59 minutes and 59 seconds which I believe is the cause of the error (since it is just under an hour). Will report back once I have a full days data to test ppsd with.

Thanks

Joshua

You can also load multiple 1-hour files into one stream object and then
feed the assembled (e.g.) 1-day stream object to PPSD.add(). If you're
missing data in between hours like you indicated be aware that gaps are
filled with zeros in processing by default (controlled by skip_on_gaps
parameter), same as with PQLX.

T

Hi All,

So I was able to plot the ppsd by ensuring the length was correct.

Unfortunately the ppsd does not come out correctly (I get a red line across the top).

I think the paz section of my code must be the cause.

I am entering

paz = {‘gain’: 1983000, ‘poles’: [-0.02365+0.02365j, -0.02365-0.02365j, -393.011, -7.4904, -53.5979-21.7494j, -53.5979+21.7494j], ‘sensitivity’: 1, ‘zeros’: [-5.03207, 0, 0]}

Then using it in

ppsd = PPSD(tr.stats, paz, ppsd_length=win_len, overlap=ovr_lap)

The values entered in the paz come from the following:

CMG-6TD CALIBRATION



WORKS ORDER:

|

2591

|

DATE:

|

22 October 2004

|

  • | - | - | - |




    |



    |



    |



    |


    SERIAL NUMBER:

    |

    T6294

    |

    TESTED BY:

    |

    ENS

    |




    |



    |



    |



    |

Velocity

Digitiser

Digital

Response

Output

Output

V/m/s

uV/count

m/s/count

VERTICAL

1039.49

0.2575

2.477E-10

NORTH/SOUTH

1104.32

0.2594

2.349E-10

EAST/WEST

1259.73

0.2544

2.019E-10

Mass Position

Mass Position

Mass Position

Coil Constant

Response

Digital Output

Digital Output

A/m/s2

V/m/s2

uV/count

m/s2/count

VERTICAL

0.0026

7.80

0.5846

7.498E-08

NORTH/SOUTH

0.0028

8.28

0.5832

7.041E-08

EAST/WEST

0.0031

9.45

0.5845

6.187E-08

Current Consumption: @12V ( ±10 mA)

79mA

This sensor operates from 10 to 36 Volts only.

POLES AND ZERO TABLE

WORKS ORDER NUMBER: 2591

SENSOR SERIAL NO: T6294

Velocity

Digitiser

Digital

Response

Output

Output

V/m/s

uV/count

m/s/count

VERTICAL

1039.49

0.2575

2.477E-10

NORTH/SOUTH

1104.32

0.2594

2.349E-10

EAST/WEST

1259.73

0.2544

2.019E-10

Mass Position

Mass Position

Mass Position

Coil Constant

Response

Digital Output

Digital Output

A/m/s2

V/m/s2

uV/count

m/s2/count

VERTICAL

0.0026

7.80

0.5846

7.498E-08

NORTH/SOUTH

0.0028

8.28

0.5832

7.041E-08

EAST/WEST

0.0031

9.45

0.5845

6.187E-08

Velocity response output:



POLES (HZ)

|

ZEROS HZ

|

  • | - |




    |



    |


    -23.65e-3 + 23.65e-3j

    |

    -5.03207

    |


    -23.65e-3 - 23.65e-3j

    |

    0

    |


    -393.011

    |

    0

    |


    -7.4904

    |



    |


    -53.5979 - 21.7494j

    |



    |


    -53.5979 + 21.7494j

    |



    |




    |



    |


    Normalizing factor at 1 Hz: A =1.983*10^6

    |



    |

Sensor Sensitivity: See Calibration Sheet.

Any assistance would be greatly appreciated.

Thanks

Joshua

Hi Joshua,

you need to set the sensitivity in the paz dict to the inverse of the digital output in your calibration sheet. units of sensitivity should be counts / (m/s).

best
Christian