Visualizing power spectral density using Obspy [Python]

Utpal Kumar     1 minute read

For details, visit Visualizing Probabilistic Power Spectral Densities.

Contents

  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 obspy.io.xseed 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 obspy.imaging.cm import pqlx
import warnings
warnings.filterwarnings('ignore')

Download stream using Obspy

## Downloading inventory
net = 'PB' 
sta = 'B075' 
loc='*'
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')
else:
    st = read(filename_prefix+".mseed")
    inventory = read_inventory(filename_prefix+'_stations.xml')

Add data to the ppsd estimate

tr = st.select(channel="EHZ")[0]
print(st)
st.plot(outfile=filename_prefix+"traces.png",show=False)
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)
    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] 
    plt.close('all')
    ppsd.plot(filename_prefix+"-ppsd_cumulative.png",cumulative=True,cmap=pqlx) #cumulative version of the histogram
    plt.close('all')
    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
    plt.close('all')
    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


References

  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