Subscribe to our weekly newsletter

Visualizing power spectral density using Obspy in Python (codes included)

Utpal Kumar   1 minute read      

Short demonstration of the ppsd class defined in Obspy using 3 days of data for station PB-B075

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

Similar posts

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.

Disclaimer of liability

The information provided by the Earth Inversion is made available for educational purposes only.

Whilst we endeavor to keep the information up-to-date and correct. Earth Inversion makes no representations or warranties of any kind, express or implied about the completeness, accuracy, reliability, suitability or availability with respect to the website or the information, products, services or related graphics content on the website for any purpose.

UNDER NO CIRCUMSTANCE SHALL WE HAVE ANY LIABILITY TO YOU FOR ANY LOSS OR DAMAGE OF ANY KIND INCURRED AS A RESULT OF THE USE OF THE SITE OR RELIANCE ON ANY INFORMATION PROVIDED ON THE SITE. ANY RELIANCE YOU PLACED ON SUCH MATERIAL IS THEREFORE STRICTLY AT YOUR OWN RISK.

Leave a comment