4
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-06-20

目的

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

#風
気象庁のサイトからまず生情報をダウンロードする

image.png

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

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

image.png

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

image.png

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

image.png

image.png

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

image.png

なるほど、元データに何も入っていないわけか、いずれ日別の平均を取るので、そのデータは捨ててもいい。
合わせて、'静穏'の場合、風速が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)

ここまできたら、こんな感じとなる
image.png

風には風速以外、風向もあるので、そのまま風速だけ平均を取るのはだめだよね。
高校時代に習ったベクトルの知識を思い出して、まず北と東の二つ垂直の方向の分解をしよう。
幸い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)

こうなる
image.png

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

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

image.png

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

#潮
潮の情報も気象庁のサイトからダウンロードできる
親切な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)

image.png

データの無いところは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))

image.png

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

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

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

幸いピーク検出はかなりのアルゴリズムがある。ここ参照
その中に簡単なものを一個選んで試そう

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

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

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

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

image.png

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

4
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?