やりたいこと
PythonのMatplotlibで風配図を描きます．Pythonでの風配図の描き方をまとめている記事があんまりなかったので書きました．ぜひ参考にしてみてください．
風配図とは
サクッと風配図についての説明をします．風配図（Wind rose）はある期間における各方位（大体8，16方位）の風向，風速の頻度を円形にプロットした図です．これを見ればその地点の卓越風や風環境の特徴がよくわかるので，建築分野や気象シミュレーションの精度評価などに使われます．今回は簡単のため，風向のみをプロットした風配図を描画してみようと思います．
使用データ
今回は気象庁の過去の気象データから風配図を作ります．データは以下のURLから地点，期間等を指定してダウンロードできます．この記事では高知市地方気象台の2020年7月1ヶ月ぶんの1時間ごとのデータを用います．これをテキストデータとして書き出して使います．
ここで使用するテキストデータをDataFrameとして表示するとこんな感じになります．
気象データは以下のリンクからダウンロードできます．
プログラム全体
説明の前にプログラム全体を載せておきます
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# データの読み込み
file = './input/Kochi_202007_1hur.txt'
df = pd.read_csv(file, delimiter=' ')
# データの配列化， 風向の抽出
data = df.values
wdir = data[:, 12]
# 16方位をリストアップ
dir_type = ['N','NNE','NE','ENE',
'E','ESE','SE','SSE',
'S','SSW','SW','WSW',
'W','WNW','NW','NNW']
# 各風向ごとの出現頻度を集計する
val_arr = np.zeros((16,))
count = 0
for i in dir_type:
index = np.where(wdir == dir_type[count])
num = index[0].size
per = num / len(wdir)
val_arr[count] = per
count += 1
# 風配図を閉じた円形にするために， 風向の軸と， 出現頻度の配列の終点と視点を一致させる
angles = np.linspace(0.0, 2 * np.pi, len(dir_type)+1, endpoint=True)
val = np.concatenate((val_arr, [val_arr[0]]))
# 風配図を描画する
plt.rcParams["font.size"] = 15 # フォントサイズを15に統一する
fig = plt.figure(figsize=(8, 8))
# 円形のグラフを作成する
ax = fig.add_subplot(111, projection='polar')
# 時計回りに指定する
ax.set_theta_direction(-1)
# 北（N）を始点にする
ax.set_theta_zero_location('N')
# プロットする
ax.plot(angles, val, color='b')
# 図形を塗りつぶす
ax.fill(angles, val, alpha=0.25, color='b')
# 円形のグリッドを表示させる
ax.set_thetagrids(angles[:-1]*180 / np.pi, dir_type)
# グラフの径を出現頻度の最大値に合わせる
ax.set_rlim(0, np.max(val))
# グラフタイトル
plt.title('Wind Rose \n Kochi City, 2020, Jul.')
# 表示
plt.show()
# 保存
fig.savefig('./output/windrose.png')
説明
ライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
データの読み込み
まずpandasで気象データを読み込みます．
# データの読み込み
file = './input/Kochi_202007_1hur.txt'
df = pd.read_csv(file, delimiter=' ')
# データの配列化， 風向の抽出
data = df.values
wdir = data[:, 12]
pandasで読み込んだテキストデータのDataFrameを配列化し，スライスで風向部分を抽出しました．
風配図作成のための準備
風配図作成には16方位の風向の出現頻度を算出する必要があります．まずはこの配列を作ります．
# 16方位をリストアップ
dir_type = ['N','NNE','NE','ENE',
'E','ESE','SE','SSE',
'S','SSW','SW','WSW',
'W','WNW','NW','NNW']
# 各風向ごとの出現頻度を集計する
val_arr = np.zeros((16,))
count = 0
for i in dir_type:
index = np.where(wdir == dir_type[count])
num = index[0].size
per = num / len(wdir)
val_arr[count] = per
count += 1
あらかじめ16方位のリストを作成します．次に16方位ごとの出現頻度を算出します．
for文で16方位ごとに風向の出現回数を数え上げ，全数で割ることで頻度を計算し，16要素の空の配列(val_arr)に格納していきます．
これで16方位ごとの出現頻度を配列化することができました．
データが揃ったのでグラフ化といきたいところですが，今回描きたいグラフは円形なので，入力するデータもこれに対応させないといけません．つまり，データの始点と終点を一致させる作業が必要です．
# 風配図を閉じた円形にするために， 風向の軸と， 出現頻度の配列の終点と視点を一致させる
angles = np.linspace(0.0, 2 * np.pi, len(dir_type)+1, endpoint=True)
val = np.concatenate((val_arr, [val_arr[0]]))
これで風配図を描くためのデータが完成しました．
風配図の描画
グラフを描いていきましょう．add_subplotのprojectionを'polar'に指定することで円形のグラフが描けます．
# 風配図を描画する
plt.rcParams["font.size"] = 15 # フォントサイズを15に統一する
fig = plt.figure(figsize=(8, 8))
# 円形のグラフを作成する
ax = fig.add_subplot(111, projection='polar')
風配図は基本的に北(N)から北北西(NNE)，北西(NE)の順で時計回りにプロットされます．デフォルトでは反時計回りで始点が90°となっているので設定をいじってやります．ちなみにこの状態でプロットするとこんな感じになります．
# 時計回りに指定する
ax.set_theta_direction(-1)
# 北（N）を始点にする
ax.set_theta_zero_location('N')
出現頻度をプロットしていきます．
# プロットする
ax.plot(angles, val, color='b')
# 図形を塗りつぶす
ax.fill(angles, val, alpha=0.25, color='b')
グリッドやタイトル，軸の最大値を設定して，完成です．
# 円形のグリッドを表示させる
ax.set_thetagrids(angles[:-1]*180 / np.pi, dir_type)
# グラフの径を出現頻度の最大値に合わせる
ax.set_rlim(0, np.max(val))
# グラフタイトル
plt.title('Wind Rose \n Kochi City, 2020, Jul.')
# 表示
plt.show()
# 保存
fig.savefig('./output/windrose.png')
完成
最後に
実際の気象データの読み込みから風配図作成まで解説してみました．参考にしていただけると幸いです．