一週間分の波形

コマンド

$ python waveform_nm_3.py 190107_N.MMOH_U_bp2-30_sqr.s 190108_N.MMOH_U_bp2-30_sqr.s 190109_N.MMOH_U_bp2-30_sqr.s 190110_N.MMOH_U_bp2-30_sqr.s 190111_N.MMOH_U_bp2-30_sqr.s 190112_N.MMOH_U_bp2-30_sqr.s 190113_N.MMOH_U_bp2-30_sqr.s

コード

waveform_nm_3.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

# コマンドラインから引数
file1 = sys.argv[1]
file2 = sys.argv[2]
file3 = sys.argv[3]
file4 = sys.argv[4]
file5 = sys.argv[5]
file6 = sys.argv[6]
file7 = sys.argv[7]

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

# sacファイル読み込み
sac1 = read(file1, debug_headers=True)#読み込みたいファイルをここに入力
sac2 = read(file2, debug_headers=True)#読み込みたいファイルをここに入力
sac3 = read(file3, debug_headers=True)#読み込みたいファイルをここに入力
sac4 = read(file4, debug_headers=True)#読み込みたいファイルをここに入力
sac5 = read(file5, debug_headers=True)#読み込みたいファイルをここに入力
sac6 = read(file6, debug_headers=True)#読み込みたいファイルをここに入力
sac7 = read(file7, debug_headers=True)#読み込みたいファイルをここに入力

# 振幅データ（縦軸）
y1 = sac1[0].data
y1 = y1*1e-18
y2 = sac2[0].data
y2 = y2*1e-18
y3 = sac3[0].data
y3 = y3*1e-18
y4 = sac4[0].data
y4 = y4*1e-18
y5 = sac5[0].data
y5 = y5*1e-18
y6 = sac6[0].data
y6 = y6*1e-18
y7 = sac7[0].data
y7 = y7*1e-18


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

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

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

# 波形
ax1=plt.subplot(7,1,1)
# plt.title(title)
plt.plot(time, y1, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

ax2=plt.subplot(7,1,2, sharex=ax1)
plt.plot(time, y2, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

ax3=plt.subplot(7,1,3, sharex=ax1)
plt.plot(time, y3, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

ax4=plt.subplot(7,1,4, sharex=ax1)
plt.plot(time, y4, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

ax5=plt.subplot(7,1,5, sharex=ax1)
plt.plot(time, y5, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

ax6=plt.subplot(7,1,6, sharex=ax1)
plt.plot(time, y6, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

ax7=plt.subplot(7,1,7, sharex=ax1)
plt.plot(time, y7, color='black', linestyle='solid', linewidth = 0.3)
plt.tick_params(labelbottom=False)
plt.xlim(0, 86400)
plt.ylim(0, 1e-13)
# plt.ylabel("Ground velocity [m^2/s^2]")
# plt.xlabel("Time [hour]")

plt.xticks(np.arange(0, 86400, step=3600))
plt.tick_params(labelbottom=False)

plt.show()

出力結果

一週間分の地震波形