「風配図」 (英語では windrose) とは、風向別の頻度を、それぞれの方位にあわせて表示したものである。 ヒストグラムの極座標表示ともいえる。
Python で 観測データにもとづいて風配図をかくプログラム例をしめしているウェブページがいくつかある。 そのうちには、windrose というパッケージをつかっているものがおおい。 それをつかって、風向だけでなく風速の値についても階級わけをした図をつくっている例もある。 しかし、基本的な作図パッケージ matplotlib をつかい、 作図開始の際に極座標を (「projection='polar'」という形で) 指定して、 plot や fill で作図している例もある。 ここでは (1) として後者の例、 (2) として windrose パッケージをすなおにつかった例ををしめす。
[2024-10-10] 2024-09-30 のバージョンは、結果の方位が 10°ずれていました。2024-10-10 夜に修正しました。
実習用マシンでは 2024-09-30 現在 niigata でも hamana でも /home/masudako/adv2024/windrose/ にある。これは変更するかもしれない。
作図について、つぎの記事を手本にした。ただしその記事であつかわれた風向データは16方位である。 データの読みこみや頻度の集計は別の形をとった。
材料として、気象庁アメダスの 1地点の ある期間の 1時間ごとの風向をとりあげることにする。 ただし、具体的には、2010年から2019年までについてつくられた『気象庁観測月報』 CD-ROM のうちの アメダス10分値から、毎時の正時の10分値をとりだしたものをつかう。
増田による形式変換後のファイルは、各地点ごと、月ごとのCSVファイルになっていて、 アメダス地点番号 (5けたの数字) と 年・月は、ファイル名にふくまれている。 たとえば windhourly36476_201912.csv は、地点番号 36476 の地点の2019年12月の風速・風向データである。 各行が1時刻 (時間間隔は1時間) の観測値をあらわす。 もとデータに日時があって観測値が欠測であるものは、欠測をあらわすフラグ値つきでふくまれているが、 もとデータの日時のフィールドが正しくないものはふくめないことにした。 したがって行数はかならずしもあるべき時間数にひとしくない。 CSVの列は左から、日時、風速 [m/s で小数1けたまで]、風速のフラグ、風向 [度 で整数値]、風向のフラグ である。 日時は日本標準時によるもので、たとえば「2019-12-31 01:00」のような形ではいっている。 静穏は風向 0、北風が風向 360度とされている。
このデータセットでは、風向は「度」の整数値で表現されている。 しかし、内容をこれまでに見たかぎりでは、数値はいずれも10の倍数である。 つまり、方位は、静穏を別として、36とおりの値をもっている。 (ただし、全部の記録で確認したわけではないので、10の倍数でない値をふくむ地点・年月もあるのかもれない。) ここでは、36方位で集計することにする。
【他方、2020年以後の観測値についてわれわれがつかえるのは、気象庁ウェブサイトの 「過去の地点気象データダウンロード」で得られるCSVファイルであり、こちらは「北北東」などの16方位で 収録されている。 長期を統一して見るためには、16方位に統一すべきかもしれない。 しかし、36方位を16方位にまとめると、16のうち4つの方位に他の方位よりも多くの事例がふくまれるだろう。】
対象期間は解析者側のつごうできめればよいのだが、 「月報」は月ごとにつくられているので、 月単位で、複数の年、複数の月をひとまとめにかぞえることができるようなプログラムにする。 (これまでにつくったプログラムでは、年と月とをそれぞれ独立に、はじめとおわりを指定するようにしてあり、 年にまたがる期間についての集計には、まだ対応していない。)
観測データには欠損や異常値がありうる。 「月報」では、気象要素ごと、時刻ごとにフラグ値がつけられていて、 フラグ値が 0 または 1 のばあいが正常とされている。 ここでは、(今後、風速もあつかおうとしているので) 風向と風速の両方のフラグを読み、 どちらも 0 または 1 のばあいを有効とする。 また、相対頻度は、有効な観測値 (静穏をふくむ) の件数を分母にとることにする。
【気象庁ウェブサイトからダウンロードしたCSVファイルの「品質情報」は、ここでいうフラグと 同じ趣旨のものだが、値がまったくちがう。そちらでは品質情報 0 や 1 はデータ欠損で、 8 が正常なのだ。両方をあわせてあつかうプログラムをかくのはむずかしい。】
CSVファイルをひとつずつ、pandas の read_csv でデータフレームによみこむ。 それを pandas の concat で結合し、対象期間の風速・風向をひとつのデータフレームにする。 それから、pandas の query をつかって、風向と風速のフラグ値がいずれも 0 または 1 の日時だけにした データフレームをつくる。
区間別の頻度の集計は、自分でプログラムを書くこともむずかしくないが、 ここでは Numpy の histogram 関数を利用する。 これは、(matplotlib の hist や pandas の hist がヒストグラムを作図する機能であるのとちがって) ヒストグラムになる数値の表をつくる機能である。
histogram で「density=False」とすると、各区間の要素個数がそのまま結果の配列にはいる。 「density=True」とすると、各区間の要素個数を全要素個数でわった相対値がはいる。 「density」は「確率密度」にでてくる意味での「密度」にするかどうかの真偽値なのだ。 わたしのプログラム例は、ここの True と False を書きかえるだけで両方に対応できるようにした。
ここで区別したい要素値は、0 から 360 まで 10 きざみの値である (そのうち 0 は静穏なので作図にはつかわないが)。 そこで、区間は、-5 から 365 までの 37 区間をとるような形で指定している [この部分 2024-10-10 修正]。
[注意 (2024-10-11)] このプログラム例は、風向を「度」を単位としてあらわした数値の 1の位が 0であることを前提として つくったので、そうでないと正しく動作しないことがある。 とくに、その値が 1～4 であると、実際には北からやや東寄りの風なのに、静穏と誤認してしまう。 この点は追って修正したいと考えているが、まだとりかかっていない。
Numpy の historgram 関数が、区間の端にあたる要素値をどう処理するか、わたしはまだ確認していない。 これまで見たかぎりでは風向の値の1の位は 5 にならないので、さしつかえないと判断している。
ひとまず、ひとつだけの風配図をつくるプログラムにした。 plt.show() による画面表示をする形にしてあるが、 それをコメントにして、プログラムの最後にある plt.savefig() をふくむ「##」ではじまるコメント行の「##」をとりのぞく (空白をいれない) と、 画像ファイルを保存する形になる。
風向は北から時計まわりだが、極座標の偏角は xの正の方向 (方位でいえば東) から反時計まわりが標準である。 ここでは、(Qiita の記事で知った知識をつかって) matplotlib の set_theta_direction(-1) で時計まわり、 set_theta_zero_location('N') で北を開始点とすることをあらかじめ指定する。
グラフは、36方位のそれぞれに1点をとった多角形をして作図する。 これは、ヒストグラムを折れ線グラフで表現したものに相当する。 【長方形の柱をならべたヒストグラムに相当するものは、 36方位のそれぞれに2点をとった図形にすれば実現できるだろう。追ってやってみようと思う。】
日本語文字のあつかいはOS依存性があるので、文字はASCIIとする。用語は英語をおもにする。
風向の目盛りは、16方位で、「NNE」のようにつける。
静穏の件数は、いまのところ図には表示せず、printコマンドで標準出力に 「Calm cases: 99」のような形で書きだすだけにしている。
データが /diskb/jma_kansoku_geppo/windhourly/ にあることを想定している。 はじめのほうに書かれた変数 datapath の値で変更可能。 【hamana では /diskb/jma_kansoku_geppo/ の実体は /raid/jma_kansoku_geppo/ である。】
結果の画像ファイルを置く場所として、current directory の直下に「image」があることを 想定している。 おわりのほうに書かれた変数 imagefilepath の値で変更可能。
「windrose」パッケージを追加した。
【[ローカル作業メモ] Linux マシン niigata で、anaconda3 の pip をつかって 「pip install windrose」でインストールした。 そのまえに「conda install windrose」をこころみたが、パッケージがみつからないというエラーになった。 condaの探索対象のリポジトリをふやせばよいのだろうが、pip でできたので、そちらはこころみていない。 pip でインストールしたパッケージと conda でインストールしたパッケージが混在してしまったので、 今後、もしまずいことが起きたら、インストールしなおす必要があるかもしれない。】
データを読むところは (1) と同じにし、 月単位で複数の月の観測値をまとめ、欠損をはずして有効な観測値をまとめた 風向・風速のデータフレーム dfvalid をつくるようにした。
windrose パッケージの WindroseAxes の bar 関数に dfvalid の風向・風速をあたえると、風配図をつくることができる。 ここで風向は「度」単位であたえることになっているので、気象庁観測月報CD-ROM の値のままでよい。
結果をみると、16方位で集計されている。 オプションを指定すれば36方位で集計することもできるはずだが、まだしらべていない。 方位の目盛りは、とくに指定しないと、8方位の英語略号になるが、北東は「N-E」のようにハイフンがはいる。 これもオプションで変更できるはずだが、まだしらべていない。 (1) のとおりに ax.set_thetagrids をつかったのではうまくいかなかった。
bar 関数のオプション normed は、頻度を規格化(normalize) するかどうかであり、True にすると百分率 (%) になる。
bar 関数のオプション opening は、方位別の棒 (扇型) のあいだをあけるかどうかであり、1.0 にすると ほぼすきまなしになる。
bar 関数のオプション edgecolor は、棒のない部分の色の指定である。
棒の色の指定は、あらかじめ「plt.rcParams['image.cmap'] = 'bwr'」で カラーマップを指定している。この例の 'bwr' は、青から白をへて赤にむかうカラーマップである。 これを 'jet' にかえると、青から緑をへて赤にむかうカラーマップになる。
風速の階級わけは windrose パッケージ内で自動的にきめている。 その凡例を表示させるため set_legend を つかっている。
(作図と関係ないデータ読みこみの部分で) Pandas ライブラリのつかいかたがうまくないので、つぎの warning がでる。
SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
(1) の実行例と同じデータについて実行してみた。
windroseパッケージが静穏をどのように判断しているのか、それに対するアメダスデータの与えかたが適切か、 わからないので、「静穏でない観測値」をぬきだしたデータフレーム dfnotcalm をつくって WindroseAxes の bar 関数にあたえられるようにし、 ひとまず「静穏でない観測値」の条件を「風向が 0 でない観測値」としてあたえてみた。 まだ精密に比較していないが、図を見たところでは、静穏をふくめたままあたえたばあいと同じに見えた。
ここまでの設定では、風向ごとの頻度と図の中心からの距離との関係も、 風速と色との関係も、ひとつひとつの図ごとに、データから決めるようになっている。 地点または時期のちがう複数の観測値群についてそれぞれ風配図をつくるとき、 このふたつの関係をそれぞれ統一したいことがある。
頻度と図の中心からの距離との関係については、set_rgrids をつかえばよいことがわかった。 たとえば、規格化した頻度を表示するばあいについて、 つぎのようにすると、図の円の半径が 規格化した頻度 30 % になり、 5 % ごとに円がかかれる。
rgrids = np.arange(5, 30+5, 5) ax.set_rgrids(rgrids, labels=list(rgrids))
風速と色との関係については、bar の bins オプションをつかえばよいことがわかった。 たとえば、つぎのように bins という配列変数に値をいれておき
bins = np.arange(0,16+2, 2)
ax.bar のかっこ内に「, bins=bins」というオプションを追加する。 (ここで1つめのbinsはオプション名で変えられない。 2つめのbinsは配列変数名で変えられる。) すると、風速 2 m/s きざみで 16 m/s まで色をかえて表示する。
凡例 (legend) の位置は、風配図の情報とかさならないように調整する必要がある。 これは windrose に特有なものではなく、matplotlib の legend のオプション bbox_to_anchor, loc, borderaxespad をつかえばよい。 まだ自分の例を用意できていないが、つぎのウェブページ参照。