- 追加された行はこのように表示されます。
- 削除された行は
このように表示されます。
!! 防災研の Hi-Net のデータを Python で
(どっかにあるか?)
http://www.hinet.bosai.go.jp/faq/?LANG=ja#Q09
!ダウンロード(イベント(地震)データ
Login -> ベント波形データダウンロード -> 日時入力 -> 地震を選択
! SAMPLE
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import subprocess
import numpy as np
import scipy as sp
from scipy import fftpack
import pylab as pl
# win32tool をインストールしてね〜〜
# cmds = ['./dewin_32', '-c','センサ名', 'データフアイル.evt']
# cmds = ['./dewin_32', '-c','54a0', './U20180219000073_20.evt']
cmds = ['./dewin_32', '-c','5452', './U20180219000073_20.evt']
time_step=0.01 # 100Hz データから撮りたい!!
# win32tool を実行
try:
res = subprocess.check_output(cmds)
except:
print ( "Error." )
##
# nlen = len(res)
rdat = res.split()
nlen = len(rdat)
#
# Data の配列化
# data = float(Scale) / float(Factor) * np.array(rdat, dtype=np.float64)
data = np.array(rdat,dtype=np.float64)
##
## FFT
sample_freq = fftpack.fftfreq(data.size, d=time_step)
sig_fft = fftpack.fft(data)
pidxs = np.where(sample_freq > 0)
(freqs, power) = sample_freq[pidxs], np.abs(sig_fft)[pidxs]
freq = freqs[power.argmax()]
#実行結果
pl.figure()
pl.subplot(1,2,1),pl.plot(freqs, power)
pl.ylabel('plower')
pl.xlabel('Frequency [Hz]')
#
pl.subplot(1, 2, 2), pl.plot(data)
#
pl.show()
pl.savefig('EQ.png')
#