Skip to content

Instantly share code, notes, and snippets.

@jamm1985
Last active July 25, 2020 04:50
Show Gist options
  • Save jamm1985/9f159e1c909828fcfea30c4497577b30 to your computer and use it in GitHub Desktop.
Save jamm1985/9f159e1c909828fcfea30c4497577b30 to your computer and use it in GitHub Desktop.
Simple waveforms (from various earthquakes) plot with trace differentiation, normalization.
"""
File: wave_plots.py
Author: Andrey Stepnov
Email: [email protected], [email protected]
Github: https://github.com/jamm1985
Description: puts waveforms from different earthquakes in single plot using obspy and matplotlib
"""
from obspy import read
from obspy import UTCDateTime
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
tracelength = 15
textsize = 12
# https://eqalert.ru/#/events/Wgvx4NMj
# LNSK SHZ - velocity
# read waveform file
ch_1 = read("../wave_data/2013-10-14-0836-05S.IMGG__021")
ch_1 = ch_1.select(station='LNSK', channel='SHZ')
# slice time
dt = UTCDateTime("2013-10-14T08:36:25.60")
ch_1 = ch_1.slice(dt, dt + tracelength)
ch_1.differentiate()
ch_1.normalize()
#ch_1.plot()
# https://eqalert.ru/#/events/GMmGPNgm
# ARGI SHZ - velocity
ch_2 = read("../wave_data/2016-06-21-1759-25M.IMGG__018")
ch_2 = ch_2.select(station='ARGI',channel='SHZ')
dt = UTCDateTime("2016-06-21T17:59:55.87")
ch_2 = ch_2.slice(dt, dt + tracelength)
ch_2.differentiate()
ch_2.normalize()
#ch_2.plot()
# https://eqalert.ru/#/events/WgvnwNgj
# STO11 Z - acceleration
ch_3 = read("../wave_data/2017-06-08-1907-39S.ST011_003")
ch_3 = ch_3.select(station='ST011',channel='C03')
dt = UTCDateTime("2017-06-08T19:08:07.80")
ch_3 = ch_3.slice(dt, dt + tracelength)
ch_3.normalize()
# custom fix zero shift :)
ch_3[0].data = ch_3[0].data - 0.08
#ch_3.plot()
# plot
matplotlib.rcParams.update({'font.size': 15})
ax3 = plt.subplot(313)
plt.plot(ch_3[0].times(), ch_3[0].data, "g-")
plt.yticks(np.arange(-1, 2, step=1))
ax3.set_title(str(ch_3[0].stats.starttime),
fontsize = textsize,
loc = 'right',
pad = 1)
plt.xlabel('Seconds')
ax1 = plt.subplot(311, sharex=ax3)
plt.plot(ch_1[0].times(), ch_1[0].data, "b-")
plt.setp(ax1.get_xticklabels(), visible=False)
plt.yticks(np.arange(-1, 2, step=1))
ax1.set_title(str(ch_1[0].stats.starttime),
fontsize = textsize,
loc = 'right',
pad = 1)
ax2 = plt.subplot(312, sharex=ax3)
plt.plot(ch_2[0].times(), ch_2[0].data, "r-")
plt.setp(ax2.get_xticklabels(), visible=False)
plt.yticks(np.arange(-1, 2, step=1))
ax2.set_title(str(ch_2[0].stats.starttime),
fontsize = textsize,
loc = 'right',
pad = 1)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment