psd一日改訂版

使い方

$ python psd_1day.py 2001/200109_N.MWDH_Uhp01.s

コード

psd_1week.py

import obspy
from obspy import read
from matplotlib import pylab as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
import sys
import matplotlib.ticker as tick


# コマンドラインから引数
file = sys.argv[1]

# サンプリング周波数
samplerate = 100

# sacファイル読み込み
sac1 = read(file, debug_headers=True)#読み込みたいファイルをここに入力
# 振幅データ（縦軸）
y = sac1[0].data

# データの頭の時刻
starttime = (sac1[0].stats.starttime)
# 観測点名
station = (sac1[0].stats.station)
# グラフタイトル用
title = str(starttime) +'  ' + str(station)

# 時間（横軸）．
time = np.arange(0, len(y))/samplerate

# グラフとグラフの隙間
plt.subplots_adjust(hspace=0.1)

# 波形
ax1=plt.subplot(2,1,1)
plt.title(title)
plt.plot(time, y, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(-2e-06, 2e-06)
plt.ylabel("Ground velocity [m/s]")

# スペクトログラム
ax2=plt.subplot(2,1,2, sharex=ax1)
plt.specgram(y, Fs=samplerate,scale='dB', vmin = -200, vmax = -150, cmap = 'jet',NFFT = 1024, noverlap = 512)
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)
# plt.gca().xaxis.set_minor_locator(tick.MultipleLocator(3600))
plt.gca().yaxis.set_minor_locator(tick.MultipleLocator(2))
plt.ylabel("Frequency [Hz]")
plt.xlabel("Time [hour]")
plt.xlim(0, 86400)
# カラーバー
ax=plt.colorbar(orientation='horizontal', aspect=10, shrink=0.2)
ax.set_label('Intensity [dB]')

plt.show()