一週間分のpsd

使い方

python psd_1week.py 2001/200106_N.MWDH_U.s 2001/200107_N.MWDH_U.s 2001/200108_N.MWDH_U.s 2001/200109_N.MWDH_U.s 2001/200110_N.MWDH_U.s 2001/200111_N.MWDH_U.s 2001/200112_N.MWDH_U.s

コード

psd_1week.py

import obspy
from obspy import read
from matplotlib import pylab as plt
import numpy as np
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 #サンプリング周波数

st01 = read(file1, debug_headers=True)
data01 = st01[0].data

st02 = read(file2, debug_headers=True)
data02 = st02[0].data

st03 = read(file3, debug_headers=True)
data03 = st03[0].data

st04 = read(file4, debug_headers=True)
data04 = st04[0].data

st05 = read(file5, debug_headers=True)
data05 = st05[0].data

st06 = read(file6, debug_headers=True)
data06 = st06[0].data

st07 = read(file7, debug_headers=True)
data07 = st07[0].data



## tは時間方向での離散値．stの長さをサンプルレートで切っているので，サンプル数分の配列になってます．
# t = np.arange(0, len(data))/samplerate


# data.np.extend(dataa)
y = np.append(data01, data02)
y = np.append(y, data03)
y = np.append(y, data04)
y = np.append(y, data05)
y = np.append(y, data06)
y = np.append(y, data07)

# tは時間方向での離散値．stの長さをサンプルレートで切っているので，サンプル数分の配列になってます．
time = np.arange(0, len(y))/samplerate


# スペクトログラム
plt.specgram(y, Fs=samplerate,scale='dB', vmin = -200, vmax = -150, cmap = 'jet',NFFT = 4096, noverlap = 2048)
plt.ylabel("Frequency [Hz]")
plt.xlabel("Time [day]")

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

# plt.gca().xaxis.set_minor_locator(tick.MultipleLocator(86400))
plt.gca().yaxis.set_minor_locator(tick.MultipleLocator(2))
plt.xticks(color="None")
# カラーバー
ax=plt.colorbar()
ax.set_label('Intensity [dB]')

plt.show()

2020年1月6〜12日