気象庁の風と潮情報を整形する

目的

仕事の都合で風と潮汐の影響分析も必要となったので、その事前処理の過程をメモします。

#風
[気象庁のサイト](http://www.data.jma.go.jp/gmd/risk/obsdl/index.php)からまず生情報をダウンロードする

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/84cd097e-cb9b-ada0-5f3f-3e52237e079c.png)

日別の風速と**風向**が必要ですが、直接日別でダウンロードすると、風向の情報が入らない。
仕方なく時別の情報をダウンロードして自前で何とかする。

ダウンロードしたファイルはこんな感じ

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/17a5fa39-97cc-97c3-51ce-54949832c6b3.png)

ファイルがいっぱいあるわけでもないので、全てプログラムで処理するのではなく、ここは少し手作業が入ります。
不要な情報を削除し、カラム名を英語にしてから、ファイル全体をUTF-8にする（さくらで）。
こんな感じとなります

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/8f68c905-cd95-f782-fe59-faa24638a063.png)

Excelでそのまま開くと、文字化けとなるが、それはUTF-8のせいで、気にしなくていい。

今度はPythonの出番となる

```py

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()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/f0dcd94f-0c1f-ca62-a47e-204674da86bf.png)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/2fdf84de-c7f3-2e1e-5ba0-e791aec91f3a.png)

風向は漢字のままじゃ何にも出来ないよね。おまけに’静穏’とか、nanも入ったりします。
'静穏'はいいとしても、nanは計測の欠損？調べてみよう

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/064afcb8-5c94-79e9-0017-38d99cc4c8eb.png)

なるほど、元データに何も入っていないわけか、いずれ日別の平均を取るので、そのデータは捨ててもいい。
合わせて、'静穏'の場合、風速が0.3以下と調べって分かった。影響が少ないので、一応’北’に変換する。

```py
wind.direction.replace("静穏","北",inplace=True)
wind.dropna(axis=0,how='any',inplace=True)
```

それで、風向をラジアン数字で表すようにする（座標的には北方向を0にする）。

```py
#北北東などを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)
```

ここまできたら、こんな感じとなる
![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/3e790ea7-6fc2-ca48-a93e-07cd4bb3cdba.png)

風には風速以外、風向もあるので、そのまま風速だけ平均を取るのはだめだよね。
高校時代に習ったベクトルの知識を思い出して、まず北と東の二つ垂直の方向の分解をしよう。
幸い1時間等間隔のデータなので、持続時間を掛ける必要がなくなった。

```py
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)
```

こうなる
![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/5766980c-7b97-77e6-24db-d4336ea9c5b5.png)

ここまできたら、後は日別の平均をとるだけ（北と東両方向それぞれ）

```py
wind_processed = wind.groupby(wind['date'].dt.date).mean()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/10c587df-73ca-4135-489a-3061896bf039.png)

風のプレー処理は一旦ここで終了

#潮
潮の情報も[気象庁のサイト](https://www.data.jma.go.jp/kaiyou/db/tide/suisan/index.php)からダウンロードできる
親切なCSVファイルの形であってほしかったが、固定長のファイルだった。フォーマットについては
[ここ](https://www.data.jma.go.jp/kaiyou/db/tide/suisan/readme.html)にのっている。

風情報と同じようにUtf-8に変換してから処理する

```py
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)
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/df9f390c-480f-0a55-7e2d-1f5b2b513354.png)

データの無いところは9999表記なので、変換が必要

```py
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回大潮、小潮ですよね、月毎に極値をとればいいじゃと思ったら

```py
#一日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))
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/8a9eb741-a2f7-5e9d-e18b-84cdc6b0c6b4.png)

甘かった。何だそのぎざぎざ。そうか、月の引力があってからの潮汐だよね、太陽暦はそのまま通用するはずがない。しかも地形とかの影響もあって、不規則も生じるよね。

もうカレンダーはたよりにならない、その波形だけでピーク（大潮、小潮）を探そう。
単純に上がり下がりの境目がピークだよね

```py
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月の間）
![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/053f9bd9-f34e-a21d-69b4-0ee6348c2fa3.png)

幸いピーク検出はかなりのアルゴリズムがある。[ここ参照](https://blog.ytotech.com/2015/11/01/findpeaks-in-python/)
その中に簡単なものを一個選んで試そう

```py
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')
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/b2d292ae-ec9f-72c8-b614-a62b3e171099.png)

おー、先よりずっとよかったじゃ。

最後は、そのピークのラベルを元データに貼る

```py
df2015['peak'] = 0
df2015['peak'][~df['min'].isnull()] = -1
df2015['peak'][~df['max'].isnull()] = 1
```

#おまけ

実は[ここ](https://org-technology.com/posts/smoother.html)を参考して、フーリエ変換で波形を平滑化することもためしたので、備考として書き残す

```py
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))
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/186671/f61ee3c6-7a13-ad5c-26ec-40a60f0fc48d.png)

平滑化はそれなりに出来た。