目的
仕事の都合で風と潮汐の影響分析も必要となったので、その事前処理の過程をメモします。
風
気象庁のサイトからまず生情報をダウンロードする
日別の風速と風向が必要ですが、直接日別でダウンロードすると、風向の情報が入らない。
仕方なく時別の情報をダウンロードして自前で何とかする。
ダウンロードしたファイルはこんな感じ
ファイルがいっぱいあるわけでもないので、全てプログラムで処理するのではなく、ここは少し手作業が入ります。
不要な情報を削除し、カラム名を英語にしてから、ファイル全体をUTF-8にする（さくらで）。
こんな感じとなります
Excelでそのまま開くと、文字化けとなるが、それはUTF-8のせいで、気にしなくていい。
今度はPythonの出番となる
from pathlib import Path
import pandas as pd
import numpy as np
filename = Path('data/wind2016.csv')
wind = pd.read_csv(filename)
wind.head()
wind.direction.unique()
風向は漢字のままじゃ何にも出来ないよね。おまけに’静穏’とか、nanも入ったりします。
'静穏'はいいとしても、nanは計測の欠損？調べてみよう
なるほど、元データに何も入っていないわけか、いずれ日別の平均を取るので、そのデータは捨ててもいい。
合わせて、'静穏'の場合、風速が0.3以下と調べって分かった。影響が少ないので、一応’北’に変換する。
wind.direction.replace("静穏","北",inplace=True)
wind.dropna(axis=0,how='any',inplace=True)
それで、風向をラジアン数字で表すようにする（座標的には北方向を0にする）。
#北北東などをangel に変換
dirs_en = ["N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW"]
dirs_jp = ["北", "北北東", "北東", "東北東", "東", "東南東", "南東", "南南東", "南", "南南西", "南西", "西南西", "西", "西北西", "北西", "北北西"]
wind['dir_dig'] = wind.direction.apply(lambda x: dirs_jp.index(x)*2*np.pi/16)
ここまできたら、こんな感じとなる
風には風速以外、風向もあるので、そのまま風速だけ平均を取るのはだめだよね。
高校時代に習ったベクトルの知識を思い出して、まず北と東の二つ垂直の方向の分解をしよう。
幸い1時間等間隔のデータなので、持続時間を掛ける必要がなくなった。
wind.date = pd.to_datetime(wind.date)
wind['speed_n'] = wind.speed * np.cos(wind.dir_dig)
wind['speed_e'] = wind.speed * np.sin(wind.dir_dig)
こうなる
ここまできたら、後は日別の平均をとるだけ（北と東両方向それぞれ）
wind_processed = wind.groupby(wind['date'].dt.date).mean()
風のプレー処理は一旦ここで終了
潮
潮の情報も気象庁のサイトからダウンロードできる
親切なCSVファイルの形であってほしかったが、固定長のファイルだった。フォーマットについては
ここにのっている。
風情報と同じようにUtf-8に変換してから処理する
from pathlib import Path
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack
from scipy.signal import argrelextrema
filename = Path('data/tidefrom2015.raw')
#using pandas with a column specification
col_specification =[(0,72),(72,78),(78,80),(80,84),(84,87),(87,91),(91,94),(94,98),(98,101),(101,105),(105,108),(108,112),(112,115),(115,119),(119,122),(122,126),(126,129),(129,133),(133,136)]
columns_name = ['timetide','date','location','hight1','highh1','hight2','highh2','hight3','highh3','hight4','highh4','lowt1','lowh1','lowt2','lowh2','lowt3','lowh3','lowt4','lowh4']
data = pd.read_fwf(filename, header=None,colspecs=col_specification)
data.columns = columns_name
# 日付、満潮、干潮時刻、潮位のみが必要で、時刻別の潮位データが不要.Location も全部同じなので、不要でも
data.drop(columns=['timetide','location','highh3','hight3','highh4','hight4','lowt3','lowh3','lowt4','lowh4'],axis=1,inplace=True)
#0 padding
data.replace('\s+', '0',regex=True,inplace=True)
データの無いところは9999表記なので、変換が必要
data.replace('9999',np.nan,inplace=True)
data.replace(9999,np.nan,inplace=True)
data.replace('999',np.nan,inplace=True)
data.replace(999,np.nan,inplace=True)
data.head()
大潮、小潮の場合の動きを分析するので、まずどれが大潮、どれが小潮を洗い出す必要だった。潮汐って、月ごとに2回大潮、小潮ですよね、月毎に極値をとればいいじゃと思ったら
#一日2回切り替えがあるので、その中の潮位差の大きい方をとる
data['diff1'] = data.highh1 - data.lowh1
data['diff2'] = data.highh2 - data.lowh1
data['diff'] = data[['diff1','diff2']].max(axis=1)
#データが多いと見づらくなるので、とりあえず2015のデータだけを対象
df2015 = data.iloc[:365].copy()
df2015[['diff']].plot(figsize=(14,6))
甘かった。何だそのぎざぎざ。そうか、月の引力があってからの潮汐だよね、太陽暦はそのまま通用するはずがない。しかも地形とかの影響もあって、不規則も生じるよね。
もうカレンダーはたよりにならない、その波形だけでピーク（大潮、小潮）を探そう。
単純に上がり下がりの境目がピークだよね
df = df2015[['diff']].copy()
# Find local peaks
df['min'] = df['diff'][(df['diff'].shift(1) > df['diff']) & (df['diff'].shift(-1) > df['diff'])]
df['max'] = df['diff'][(df['diff'].shift(1) < df['diff']) & (df['diff'].shift(-1) < df['diff'])]
# Plot results
plt.figure(figsize=(14, 6), dpi=80)
plt.scatter(df.index, df['min'], c='r')
plt.scatter(df.index, df['max'], c='g')
plt.plot(df.index, df['diff'], c='b')
やはり、ぎざぎざは邪魔だ（8月、9月辺り）、しかもイコールの場合も厄介（3,4月の間）
幸いピーク検出はかなりのアルゴリズムがある。ここ参照
その中に簡単なものを一個選んで試そう
df = df2015[['diff']].copy()
n=5 # number of points to be checked before and after
df['min'] = df.iloc[argrelextrema(df['diff'].values, np.less_equal, order=n)[0]]['diff']
df['max'] = df.iloc[argrelextrema(df['diff'].values, np.greater_equal, order=n)[0]]['diff']
# Plot results
plt.figure(figsize=(14, 6), dpi=80)
plt.scatter(df.index, df['min'], c='r')
plt.scatter(df.index, df['max'], c='g')
plt.plot(df.index, df['diff'], c='b')
おー、先よりずっとよかったじゃ。
最後は、そのピークのラベルを元データに貼る
df2015['peak'] = 0
df2015['peak'][~df['min'].isnull()] = -1
df2015['peak'][~df['max'].isnull()] = 1
おまけ
実はここを参考して、フーリエ変換で波形を平滑化することもためしたので、備考として書き残す
df = df2015[['diff']].copy()
f = df.values
N = f.shape[0]
yf = fftpack.fft(f.reshape(-1))/(N/2)
dt = 0.01 # サンプリング間隔
freq = fftpack.fftfreq(N, dt)
fs = 10 # カットオフ周波数[Hz]
yf2 = np.copy(yf)
yf2[(freq > fs)] = 0
yf2[(freq < 0)] = 0
y2 = np.real(fftpack.ifft(yf2)*N) # 実数に復元
pd.DataFrame(y2).plot(figsize=(14,6))
平滑化はそれなりに出来た。