トップ 差分 一覧 ソース 検索 ヘルプ RSS ログイン

PRG-py-hinet

  防災研の 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')
#