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
>>> 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")
Metadata can be provided as an
Inventory (e.g. from a StationXML file
or from a request to a FDSN web service – be sure to use level=’response’),
Parser (mostly legacy, dataless SEED files
can be read into Inventory objects using
read_inventory()), a filename of a
local RESP file (also legacy, RESP files can be parsed into Inventory
objects; 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.
>>> inv = read_inventory("https://examples.obspy.org/BW_KW1.xml") >>> ppsd = PPSD(tr.stats, metadata=inv)
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) True
We can check what time ranges are represented in the ppsd estimate.
ppsd.times_processed 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). Other attributes storing information about the data fed into the
PPSD object are
>>> print(ppsd.times_processed[: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_processed)) 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) False >>> print("number of psd segments:", len(ppsd.times_processed)) 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) True
The graphical representation of the ppsd can be displayed in a matplotlib window..
..or saved to an image file:
>>> ppsd.plot("/tmp/ppsd.png") >>> ppsd.plot("/tmp/ppsd.pdf")
A (for each frequency bin) cumulative version of the histogram can also be visualized:
>>> from obspy.imaging.cm import pqlx >>> ppsd.plot(cmap=pqlx)
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
Time series of psd values can also be extracted from the PPSD by accessing the
plotted using the
plot_temporal() method (temporal
restrictions can be used in the plot, see documentation):
>>> ppsd.plot_temporal([0.1, 1, 10])
Spectrogram-like plots can be done using the