Boom Event Report

Developers
#1

I thought I’d share a Python module I’ve written for the Raspberry (Shake and) Boom.

Though I have programming experience from pre-windows days (Fortran, Basic, Pascal, and machine code) and some Visual Basic experience 15 to 20 years ago, I’ve only played with Python since getting the RS&B. So Python and Obspy is all relatively new to me. This is the first real result since starting to learn Python and how to handle the instrument response.

I hope it helps someone to learn or get started with Python and Obspy. I am also happy for suggestions and feedback about any errors. For example, y axis range on the PSD doesn’t seem to agree with SWARM or Dataview, so I may be missing something there.

image
image1335×1725 248 KB

Here’s what the output looks like:

R21C0 HDF 20220411 092000 UTC
R21C0 HDF 20220411 092000 UTC1200×800 116 KB

1 Like
#2

Thank you for the code share, sheeny!

I am sure that other users will find this very helpful for their first steps in Python/Obspy programming, and the fact that the output provides so much information it’s a welcome bonus.

#3

One comment:

Presenting the code as an image prevents people copying it to try out themselves.
If that was the intent, much better to include it as text.

Philip

#4

Hi Philip,

I tried posting the code as text between quotes [quote][/quote] but it turned the #comments into headers and did some strange things, so I went back to the screenshot.

OK, after having a look at some page source code here’s the code:


from obspy.clients.fdsn import Client
from obspy import UTCDateTime
import matplotlib.pyplot as plt
from obspy.signal import filter

# define start and end times
start = UTCDateTime(2022, 4, 11, 9, 20, 00) # (YYYY, m, d, H, M, S)
end = start + 10                # start plus plot duration in seconds (recommend a minimum of 10s)

# define fdsn client to get data from
client = Client('RASPISHAKE')

# get data from the FSDN server and detrend it
STATION = "R21C0"               # station name
st = client.get_waveforms("AM", STATION, "00", "HDF", starttime=start, endtime=end, attach_response=True)
st.merge(method=0, fill_value='latest')
st.detrend(type="demean")

# get Instrument Response
inv = client.get_stations(network="AM", station=STATION, level="RESP")

# set-up figure and subplots
fig = plt.figure(figsize=(12,8))    #12 x 8 inches
ax1 = fig.add_subplot(6,2,(1,3))    #left top 2 high
ax2 = fig.add_subplot(6,2,(5,7))    #left middle 2 high
ax3 = fig.add_subplot(6,2,(2,6))    #right top 3 high
ax4 = fig.add_subplot(6,2,(8,12))   #right bottom 3 high
ax5 = fig.add_subplot(6,2,(9,11))   #left bottom 2 high
fig.suptitle("AM."+STATION+'.00.HDF Channel '+start.strftime('%d/%m/%Y %H:%M:%S UTC'))


# plot the raw curve
ax1.plot(st[0].times(reftime=start), st[0].data, lw=1)

#plot the PSD
ax2.psd(x=st[0], NFFT=512, noverlap=2, Fs=100, color='r', lw=1)
ax2.set_xscale('log')

#plot spectrogram
ax5.specgram(x=st[0], NFFT=64, noverlap=16, Fs=100, cmap='plasma')
ax5.set_yscale('log')
ax5.set_ylim(0.5,50)

# set up filter
filt = [0.3, 0.5, 10, 12]   #edit filter values to suit

# remove instrument response
resp_removed = st.remove_response(inventory=inv,pre_filt=filt,output="DEF", water_level=60)

# plot the filtered and corrected, and dB curves
ax3.plot(st[0].times(reftime=start), resp_removed[0].data, lw=1, color='g')
ax4.plot(st[0].times(reftime=start), 93.979400087*abs(resp_removed[0].data), lw=1, color='b')

# set-up some plot details
ax1.set_ylabel("Counts")
ax2.set_ylabel("Power Spectral Density")         
ax3.set_ylabel("Pascals (Pa) ("+str(filt[1])+" to "+str(filt[2])+"Hz)")
ax4.set_ylabel("Decibels [dB] ("+str(filt[1])+" to "+str(filt[2])+"Hz)")
ax5.set_ylabel("Frequency, Hz")
ax1.set_xlabel("Seconds after start")
ax3.set_xlabel("Seconds after start")
ax4.set_xlabel("Seconds after start")
ax5.set_xlabel("Seconds after start")

ax1.grid(color='dimgray', ls = '-.', lw = 0.33)
ax2.grid(color='dimgray', ls = '-.', lw = 0.33)
ax3.grid(color='dimgray', ls = '-.', lw = 0.33)
ax4.grid(color='dimgray', ls = '-.', lw = 0.33)
ax5.grid(color='dimgray', ls = '-.', lw = 0.33)

#adjust subplots for readability
plt.subplots_adjust(hspace=1)

# save the final figure
#plt.savefig(STATION+' HDF '+start.strftime('%Y%m%d %H%M%S UTC'))  #comment this line out till figure is final

# show the final figure
plt.show()

Enjoy.

Al.

3 Likes