Visualizing power spectral density using Obspy [Python]

Utpal Kumar     1 minute read

For details, visit Visualizing Probabilistic Power Spectral Densities.


  1. Import necessary libraries
  2. Download stream using Obspy
  3. Add data to the ppsd estimate
  4. Visualization using Obspy
  5. Output Figures
  6. References

Location of the event (yellow star) and station (green triangle) [Screenshot from IRIS]

Downloaded stream

Import necessary libraries

from import Parser
from obspy.signal import PPSD
from obspy.clients.fdsn import Client
from obspy import UTCDateTime, read_inventory, read
import os, glob
import matplotlib.pyplot as plt
from import pqlx
import warnings

Download stream using Obspy

## Downloading inventory
net = 'PB'
sta = 'B075'
chan = 'EH*'
filename_prefix = f"{net}_{sta}"
mseedfiles = glob.glob(filename_prefix+".mseed")
xmlfiles = glob.glob(filename_prefix+'_stations.xml')

if not len(mseedfiles) or not len(xmlfiles):
    print("--> Missing mseed / station xml file, downloading...")
    time = UTCDateTime('2008-02-19T13:30:00')
    wf_starttime = time - 60*60
    wf_endtime = time + 3 * 24 * 60 * 60 #3 days of data (requires atleast 1 hour)

    client = Client('IRIS')

    st = client.get_waveforms(net, sta, loc, chan, wf_starttime, wf_endtime)
    st.write(filename_prefix+".mseed", format="MSEED")
    inventory = client.get_stations(starttime=wf_starttime, endtime=wf_endtime,network=net, station=sta, channel=chan, level='response', location=loc)
    inventory.write(filename_prefix+'_stations.xml', 'STATIONXML')
    st = read(filename_prefix+".mseed")
    inventory = read_inventory(filename_prefix+'_stations.xml')

Add data to the ppsd estimate

tr ="EHZ")[0]
ppsd = PPSD(tr.stats, metadata=inventory)
add_status = ppsd.add(st) #add data (either trace or stream objects) to the ppsd estimate

Visualization using Obspy

if add_status:
    print(ppsd.times_processed[:2]) #check what time ranges are represented in the ppsd estimate
    print("Number of psd segments:", len(ppsd.times_processed))
    ppsd.plot(filename_prefix+"-ppsd.png",cmap=pqlx) #colormap used by PQLX / [McNamara2004]
    ppsd.plot(filename_prefix+"-ppsd_cumulative.png",cumulative=True,cmap=pqlx) #cumulative version of the histogram
    ppsd.plot_temporal(period=[0.1, 1.0, 10], filename=filename_prefix+"-ppsd_temporal_plot.png") #The central period closest to the specified period is selected
    ppsd.plot_spectrogram(filename=filename_prefix+"-spectrogram.png", show=False)

Output Figures

Probabilistic Power Spectral Densities with colormap used by [McNamara2004]

Cumulative version of the histogram

Plot the evolution of PSD value of one (or more) period bins over time.

Spectrogram of the estimate


  1. McNamara, D. E., & Buland, R. P. (2004). Ambient Noise Levels in the Continental United States. Bulletin of the Seismological Society of America, 94(4), 1517–1527.
  2. McNamara, D. E., & Boaz, R. I. (2006). Seismic noise analysis system using power spectral density probability density functions: A stand-alone software package. Citeseer.

Leave a comment


Plotting seismograms with increasing epicentral distance

How to plot the distance vs seismic waveforms?


Plotting track and trajectory of hurricanes on a topographic map

How to plot the track or trajectory of a hurricane on a map?