24. Visualizing Probabilistic Power Spectral Densities

The following code example shows how to use the PPSD class defined in obspy.signal. The routine is useful for interpretation of e.g. noise measurements for site quality control checks. For more information on the topic see [McNamara2004].

>>> from obspy import read
>>> from obspy.io.xseed import Parser
>>> from obspy.signal import PPSD

Read data and select a trace with the desired station/channel combination:

>>> st = read("https://examples.obspy.org/BW.KW1..EHZ.D.2011.037")
>>> tr = st.select(id="BW.KW1..EHZ")[0]

Metadata can be provided as an Inventory (e.g. from a StationXML file or from a request to a FDSN web service), a Parser (e.g. from a dataless SEED file), a filename of a local RESP file (or a legacy poles and zeros dictionary). Then we initialize a new PPSD instance. The ppsd object will then make sure that only appropriate data go into the probabilistic psd statistics.

>>> parser = Parser("https://examples.obspy.org/dataless.seed.BW_KW1")
>>> ppsd = PPSD(tr.stats, metadata=parser)

Now we can add data (either trace or stream objects) to the ppsd estimate. This step may take a while. The return value True indicates that the data was successfully added to the ppsd estimate.

>>> ppsd.add(st)

We can check what time ranges are represented in the ppsd estimate. ppsd.times contains a sorted list of start times of the one hour long slices that the psds are computed from (here only the first two are printed).

>>> print(ppsd.times[:2])
[UTCDateTime(2011, 2, 6, 0, 0, 0, 935000), UTCDateTime(2011, 2, 6, 0, 30, 0, 935000)]
>>> print("number of psd segments:", len(ppsd.times))
number of psd segments: 47

Adding the same stream again will do nothing (return value False), the ppsd object makes sure that no overlapping data segments go into the ppsd estimate.

>>> ppsd.add(st)
>>> print("number of psd segments:", len(ppsd.times))
number of psd segments: 47

Additional information from other files/sources can be added step by step.

>>> st = read("https://examples.obspy.org/BW.KW1..EHZ.D.2011.038")
>>> ppsd.add(st)

The graphical representation of the ppsd can be displayed in a matplotlib window..

>>> ppsd.plot()

..or saved to an image file:

>>> ppsd.plot("/tmp/ppsd.png")  
>>> ppsd.plot("/tmp/ppsd.pdf")  

(Source code, png, hires.png)


A (for each frequency bin) cumulative version of the histogram can also be visualized:

>>> ppsd.plot(cumulative=True)

(Source code, png, hires.png)


To use the colormap used by PQLX / [McNamara2004] you can import and use that colormap from obspy.imaging.cm:

>>> from obspy.imaging.cm import pqlx
>>> ppsd.plot(cmap=pqlx)

(Source code, png, hires.png)


Below the actual PPSD (for a detailed discussion see [McNamara2004]) is a visualization of the data basis for the PPSD (can also be switched off during plotting). The top row shows data fed into the PPSD, green patches represent available data, red patches represent gaps in streams that were added to the PPSD. The bottom row in blue shows the single psd measurements that go into the histogram. The default processing method fills gaps with zeros, these data segments then show up as single outlying psd lines.


Providing metadata from e.g. a Dataless SEED or StationXML volume is safer than specifying static poles and zeros information (see PPSD).