気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（5）機械学習編 Automatic Front Detection in Weather Data
2021.9.17 学習したニューラルネットワークを地球上の他の地域に適用して全球前線自動描画に挑戦したことを記載しました。
変更概要
「3.4ところで・・・」に解説記事へのリンクを追記
2020.11.15 カラー版天気図データが増加したことを受けて再学習を行った結果を踏まえて一部を更新しました。
投稿してから月日が経ち、この間に気象図と気象データを地道に収集していたことから再学習を行いました。
変更概要
・教師データとして使用した天気図の期間を変更
・「3.2 生成データ(初見データ)」の中の「どのデータが寄与しているのか？」の箇所で、結果を可視化したデータを変更してアニメを追加
・「3.４ところで・・・」でGSMの予測結果に対して前線を描画させたアニメを変更
・「4.3.1 ネットワーク構造」でネットワークの変更した図を追加掲載。
1. 「天気図っぽい前線」の学習
1.1 「天気図っぽい前線」とは？
天気図の前線は、防災的な観点からか、日本に影響のあるような部分を中心に解析されているようです。
他の国の天気図を見ると、日本の天気図では描かれないような前線が登場してたりします。
例えば下記のコラムではヨーロッパの気象図に登場する前線についての記載があります。
衛星でも見えない、「隠れた」前線の話
これらから言えるのは、前線解析は必ずしも一意に決まるものでなく、機関ごとの流儀や意図がありそうということです。ということは、私が試みている機械学習による前線描画は、気象データから日本の天気図っぽい前線を描画させるもの、と言えるでしょう。
学習精度をどんどん高めて気象庁のような日本の天気図っぽい前線解析を完璧に描画できるようになると、気象実況解析のエキスパートシステムとなるものかもしれません。
永澤義嗣氏が著された気象予報と防災ー予報官の道（中公新書）には、「天気図を三千枚書いて一人前」という記載があります。今回学習に使用しているデータはおよそ2000枚です。Deep Learning予報官として、一人前というにはまだ努力不足です。
ところで、この機械学習は天気図を生成しているものの、天気予報をしているものではありません。
天気予報は現在から将来を予測するもので、前線描画は現在から現在のデータへの変換です。
2. 機械学習の流れ
2.1 何をやるのか？
CNNによって画像から画像を生成するものです。生成画像のピクセルごとに赤・青・白どの色になるべきか、確率を計算させて確率の高かった色をピクセルにセットするという手法です。
この手法は下記を参考にしています。
Exascale Deep Learning for Climate Analytics
Thorsten Kurth, Sean Treichler, Joshua Romero, Mayur Mudigonda, Nathan Luehr, Everett Phillips, Ankur Mahesh, Michael Matheson, Jack Deslippe, Massimiliano Fatica, Prabhat, Michael Houston
arXiv:1810.01993 [cs.DC]
2018年にテキサス州ダラスで開催されたSC18というスパコン分野の国際学会において、ACM Gordon Bell Prizeを受賞した論文です。私は幸いにも、SC18へ出張していて受賞講演を聴講する機会に恵まれ、公私混同しながらとても感慨深かった思い出があります。
2.2 機械学習の流れ
ニューラルネットワークによる機械学習の流れをまとめます。
(1) 入力画像を作る
気象データをダウンロードして可視化する。
6種類のカラー画像を用意しました（第2回）
(2) 教師画像を作る
「速報天気図」から色をもとに前線要素だけを抜き出した画像を作る(第3回）
教師画像を増やすために白黒天気図のカラー化なんてのもやりました(第4回)
2020.11.15更新 正式なカラー版天気図が蓄積されたので、正式カラー版のみを用いて学習するようデータを修正しました。
(3) CNNのあれこれ
・入力データ：並べ替え
CNNの入力は、6種類の入力画像(channel数は3）をconcatenateした18channelのテンソルになります。
CNNのミニバッチの中に、近接した時間のデータばかりになると学習時に偏ってしまうので、入力データと教師画像の組を時間的にランダムに並べ替えておきます。
・入力データ：one-hotベクトル
CNNは各ピクセルが赤、青、白のそれぞれに該当する確率を計算するように作成するので、RGBである教師画像については、正解が1、その他が0となった配列に変換しておきます。つまり、白なら(1,0,0)、赤なら(0,1,0)、青なら(0,0,1)というデータを作ります。
・出力データ
例えばCNNはあるピクセルの値としては、例えば(0.1, 0.7, 0.2)を出力します。この場合はこのピクセルは赤だ、として絵を作ります。
・ニューラルネットワーク
CNNはU-Netを意識した分岐経路を有するものを作成します。
・ロス関数
categorical_cross_entropyを用います。これにより赤、青、白の確率を計算します。
(4) 学習させる！(2020.11.15更新）
いよいよCNNを学習させます。
今回、2017年1月から2019年10月まで、各日から2枚（6時UTCと18時UTCのデータ）を使用して学習させています。 2018年9月から2020年10月まで各日から3枚(6時、12時、18時UTC）を使用して学習させました。
カラー版をもとにして前線要素を切り取った教師データの方が前線としてのクオリティが高いため、それらのみを教師データをとして用いたためです。白黒天気図からカラー化したものは退場頂きました。
教師データは赤、青、白のいずれかが1になっているone-hotベクトルですから学習が進むにつれて前線要素のカラーを再現できるようになっていきます。
私のMac miniは数日の間、触れないくらいの高温になります。
Mac miniのうち、最もCPUを酷使されている個体のひとつではないかと思っており、CPU冥利に尽きるというものでしょう。
(5) 予測させる！（初見データに前線を描かせる）
学習が収束したところで、学習済ネットワークを使って、初見データの前線を描画します。今回は初見データとして2019年11月から2020年1月までのデータに前線を描画させました。
3. 何はともあれ結果
3.1 学習データ
学習データの一例です。
速報天気図の前線ととほぼ同じような位置に該当する前線要素が生成できるまで収束しているように見えます。
左上から右へ、
出力結果:生成した前線要素（地上気圧の可視化画像へ重ね合わせたもの）
教師画像:速報天気図
入力:相当温位850hPa
入力:気温・気圧・風（地上）
左下から右へ、
入力:気温・高度・風（850hPa)
入力:湿数（700hPa）
入力:鉛直流（700hPa）
入力:気温・高度・風（500hPa)
3.2 生成データ（初見データ）
3.1で、2017年1月から2019年10月のデータによって学習させたネットワークを使って、前線生成させました。
うまく生成できている例
日本の東の海上に低気圧があり、寒冷前線が南西に延びています。
ニューラルネットワークでもこの寒冷前線を生成できています。また、温暖前線も途切れながらも生成されています。
かなりダメな例
速報天気図では、千島列島近海の低気圧から南西に寒冷前線が延びています。また、北緯30度線にほぼ沿って停滞前線が解析されています。
これに対して、ニューラルネットワークではどちらの前線もほぼ生成できていません。
どのデータが寄与しているのか？（2020.11.15更新)
入力データにも生成された前線を描画して作成したアニメです。
・2017年6月10日から14日までの結果
上段左から：生成した前線付き天気図、当日の速報天気図、850hPa相当温位、地上気圧・風・気温
下段左から：850hPa高度・風・気温、700hPa湿数、700hPa上昇流、500hPa高度・風・気温
6月ということで日本の南に梅雨前線（停滞前線）が生成されている時期のものです。これを見ると、850hPa相当温位の等相当温位線の集中帯のやや南側に前線が描画されているようです。気象予報士試験問題の教科書通りですね。
・2017年2月1日から6日の結果
(2月3日に天気図が入手できていない時間があるため上段の表示が乱れる日があります)
冬に低気圧が発達していくのに伴い、寒冷前線と温暖前線が生成されているところが生成できています。ここでも850hPa相当温位の等値線の集中帯との対応が大きいようです。上昇流の対応を見ると低気圧前面で上昇流（図では青色）、後面で下降流（図では赤色）となっている対応も見て取れます。
3.3 学習の進行状況
LossおよびAccuracyの遷移です。このデータを取り直すために再度学習させた800エポックまでの状況です。
200エポックごとに学習を止めて再スタートさせているのですが、その際に入力画像のランダム化（後述）を行うことで学習対象のデータが一部入れ替わったりするので、一度LossやAccuracyが悪くなります。
先ほどのうまく生成できた例についての生成画像の変化です。
左から、200、400、600、800エポック学習時点のパラメタでの生成画像、最終的に1500エポック以上
実行したときの生成画像。
前線を描くべき場所に前線を引くこと自体は早期にできるようになり、前線記号（山型のマークなど）の精度が上がってきているように見えます。
ダメな例の方です。
ん？学習が浅い方のパラメタで生成した方がまだマシだったかもしれない？
学習時間ですが、私のMac miniで1エポックで580秒前後でした。これでマル４日くらい実行し続けでした。
GPUを持っていたりするともっと先まで実行する気になるのかなあ。
3.4 ところで・・・(2020.11.15更新)
まだ前線がくっきり間違えずにできるところまで来ていませんが、データ例を増やして計算リソースを増やせば、かなり良いところまでいけるのではないか、という感触です。
それではこれは何かの役に立つのか？
前線解析は気象庁で実施している結果を入手できることができますので、プロにお任せした方がお得です。
ひとつの利用例として数値予報結果や気候計算結果に自動的に前線を描く、ということは考えられないでしょうか？
例えば上記は2020年11月5日12時UTCを初期値としたGSMの予報計算結果に対して前線を描いたものです。予想天気図は24時間後、48時間後という間隔での発表ですので、GSMの６時間ごとの結果に対して前線を描く、ということを行っています。
2021.9.17追記
ひとつの応用先として、この学習済ニューラルネットワークを、地球の他の地域のデータにも適用してみました。全球の自動前線描画というわけです。
意外に日本周辺以外の地域でも前線を描画することができています。
「天気図っぽい」気象前線の自動描画を日本から世界に拡張してみる〜領域拡張と地図投影
4.CNNのあれこれ
以下はCNN本体の説明と、多少ハマった部分のTipsです。
4.1 入力データ
4.1.1 入力画像・教師画像の組をランダムに並べ替える
入力画像と教師画像は、読み込んだ直後には日付順に並んでいます。
このまま学習させると、ミニバッチの各回は時間的に近接した気象状況となってしまうことが想定されるので、異なる気象状況が混じったバッチとするために、リストをランダムにソートします。
# i_datalist_train 年月日時順に並んだn_datalist個の入力用気象画像データ
# t_datalist_train 年月日時順に並んだn_datalist個の教師画像データ
# 0からn_datalist-1までの数列からランダムにn_datalistのインデックスを生成して
# ループ処理する
for c_list in random.sample(range(n_datalist), n_datalist ):
t_datalist_train.append(t_datalist[c_list])
# 並び替え済教師画像リスト
i_datalist_train.append(i_datalist[c_list])
# 並び替え済入力画像リスト
4.1.2 前線画像のone-hot化
前回までで作成した教師画像用の前線要素画像は、各ピクセルが赤、青、白のいずれかとなっています。
ほとんどの面積を占めるのバックグラウンドが白色、寒冷前線と高気圧記号が青色、温暖前線と閉塞前線および低気圧記号が赤色です。
CNNにより画像を生成して教師画像との間でcategorical_cross_entropyによる誤差関数を計算しますので、この前線要素画像データをone-hotベクトル化しておきます。to_categoricalというnumpyのメソッドを利用するために、RGB配列をいったん0,1,2のいずれかをとる配列に変換します。その後にto_categoricalを用いてone-hotベクトルに変換します。そのあたりを実行しているのが下記のソースです。
# t_img 読み込んだ前線画像ファイル
# img_size t_imgの画像サイズ
t_data = np.empty((img_size[1],img_size[0]))
# t_data one-hot化するために0,1,2の3値化するための配列を用意
# ピクセル毎に赤は1, 青は2, それ以外は0 とセットする
for x in range(img_size[1]):
for y in range(img_size[0]):
r,g,b = t_img.getpixel((y,x))
# r g bにそれぞれRGB値を格納する
if(r>g+20):
if(r>b+20):
t_data[x,y]=1 # 赤色ピクセル
else:
if(b>g+20):
t_data[x,y]=2 # 青色ピクセル
else:
t_data[x,y]=0 # 白色ピクセル
else:
if(b>r+20):
if(b>g+20):
t_data[x,y]=2 # 青色ピクセル
else:
t_data[x,y]=0 # 白色ピクセル
else:
t_data[x,y]=0 # 白色ピクセル
# t_data から、3要素のone-hot vector配列 T_data に変換
T_data = np_utils.to_categorical(t_data[:,:],3)
4.2 出力データ
予測させた出力データから、最も確率が大きいものを採用して３色画像に変換します。
データを出力すると下記のように、ひとつの値が大きくなっています。その値に該当する色をピクセルに配置します。
w_array
(256, 256, 3)
[[[9.99969959e-01 1.26371087e-05 1.73822737e-05]
[1.00000000e+00 8.79307649e-09 8.33461922e-09]
[1.00000000e+00 1.22459204e-12 8.95228910e-16]
...
[9.99985695e-01 6.48013793e-06 7.86928376e-06]
[9.99960303e-01 8.51386540e-06 3.12020056e-05]
[9.99833941e-01 2.61777150e-05 1.39806682e-04]]
[[9.99999881e-01 8.55169304e-08 1.83308195e-08]
[1.00000000e+00 9.66997732e-11 1.11044485e-12]
[1.00000000e+00 4.26908814e-16 1.04265986e-22]
...
ピクセルへの色配置のためのソースです。w_arrayは、(R,G,B)の確率値ベクトルが緯度x経度に並んでいます。これをもとに、画像データmask1を作成します。
# w_array 予測させた出力データのNd_array配列の1画面分が格納されている
s_img = array_to_img(w_array[:,:,:].reshape(i_dmlat,i_dmlon,3))
# s_img 画像化したデータ
new_img_size = [w_array.shape[1], w_array.shape[0]]
mask1 = Image.new('RGB', new_img_size)
# mask1 出力するための配列
for x in range(new_img_size[1]):
for y in range(new_img_size[0]):
# w_arrayの第3要素をw1/r1/b1に格納する
w1 = w_array[x,y,0]
r1 = w_array[x,y,1]
b1 = w_array[x,y,2]
if(r1>w1):
if(r1>b1): # r1>w1, r1>b1
r,g,b=255,0,0
# r1が最大の場合はRGBとして赤色をセット
else: # r1=<w1
if(b1>w1): # r1=<w1<b1
r,g,b=0,0,255
# b1 が最大の場合はRGBとして青色をセット
else: # r1=<w1, b1=<w1
r,g,b=255,255,255
# w1 が最大の場合はRGBとして白色をセット
else: # w1>=r1
if(b1>r1):
if(b1>w1): # b1>w1>r1
r,g,b=0,0,255
# b1 が最大の場合はRGBとして青色をセット
else: # w1>=b1>r1
r,g,b=255,255,255
# w1 が最大の場合はRGBとして白色をセット
else: # w1>=r1 w>=b1
r,g,b=255,255,255
# w1 が最大の場合はRGBとして白色をセット
mask1.putpixel((y,x),(r,g,b))
# RBGの値をmask1にセット
4.3 ニューラルネットワーク(Convolutional Neural Network)
4.3.1 ネットワーク構造(2020.11.15更新)
CNNの構造は、U-Netのような分岐を持たせているもので自作です。
入力が18x256x256のテンソル、出力が3x256x256のテンソルです。
画像サイズはConv2Dのstride=(2,2)パラメタの指定により半分になっていき、最終的に120x16x16になります。
通常のU-Netでは、画像サイズを半分にする際にチャネル数を倍にするのですが、今回はメモリ量の都合で、きれいに倍にはしていません。
2020.11.15更新版
入力に128x128ピクセルサイズ、64x64ピクセルサイズの画像を追加しました。エンコードの過程で256x256ピクセルサイズのデータサイズが変換される時にこれらデータをconcatします。またデコード過程でもサイズが変わるところでこれらデータを、U-Netとしてエンコーダから渡ってくるデータと合わせてconcatします。
4.3.2 Kerasによるネットワーク定義
ネットワーク定義部分のPython-Kerasのソースは下記です。
Kerasのネットワークの定義のやり方としては、Functional APIを用いて各レイヤの出力を次レイヤまたは先の方の結合部分に繋いでいくものです。
cnv2dtra = concatenate([cnv2dtra , cnv2d14 ], axis=1)
などの部分が,cnv2dtraに分岐しきたcnv2d14を第1要素（チャンネル）で結合させている部分です。
#- Neural Network Definition
num_hidden1 = 18
num_hidden2 = 32
num_hidden3 = 48
num_hidden3a = 72
num_hidden3b = 120
num_hidden4 = 128
num_hidden5 = 16
num_hidden6 = 16
num_itter = int(_num_itter) # 学習回数
#--- seting Neural Network using Keras
in_ch = 18 # 入力画像 RGB x 6種類 のチャンネル数
cnn_input = Input(shape=(in_ch, i_dmlat, i_dmlon))
# i_dmlat, i_dmlon 画像の緯度方向、経度方向の画素数
# 1st hidden layer
cnv2d1 = Conv2D(num_hidden1, data_format='channels_first', kernel_size=(3,3), dilation_rate=(3,3), activation='relu', padding='same')(cnn_input)
cnv2d2 = Conv2D(num_hidden1, data_format='channels_first', kernel_size=(3,3), dilation_rate=(3,3), activation='relu', padding='same')(cnv2d1)
cnv2d4 = Conv2D(num_hidden1, data_format='channels_first', strides=(2,2) , kernel_size=(3,3), activation='relu', padding='same')(cnv2d2)
# 2nd hidden layer
cnv2d5 = Conv2D(num_hidden2, data_format='channels_first', kernel_size=(3,3), dilation_rate=(3,3), activation='relu', padding='same')(cnv2d4)
cnv2d6 = Conv2D(num_hidden2, data_format='channels_first', kernel_size=(3,3), dilation_rate=(3,3), activation='relu', padding='same')(cnv2d5)
cnv2d8 = Conv2D(num_hidden2, data_format='channels_first', strides=(2,2) , kernel_size=(3,3), activation='relu', padding='same')(cnv2d6)
# 3rd hidden layer
cnv2d9 = Conv2D(num_hidden3, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2d8)
cnv2d10 = Conv2D(num_hidden3, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2d9)
cnv2d12 = Conv2D(num_hidden3, data_format='channels_first', strides=(2,2) , kernel_size=(3,3), activation='relu', padding='same')(cnv2d10)
# 4th hidden layer
cnv2d13 = Conv2D(num_hidden3a, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2d12)
cnv2d14 = Conv2D(num_hidden3a, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2d13)
cnv2d16 = Conv2D(num_hidden3a, data_format='channels_first', strides=(2,2) , kernel_size=(3,3), activation='relu', padding='same')(cnv2d14)
cnv2d17 = Conv2D(num_hidden3b, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2d16)
cnv2d19 = Conv2D(num_hidden3b, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2d17)
#--- decode start
cnv2dtra = Conv2DTranspose(num_hidden3a, data_format='channels_first', kernel_size=(3,3), strides=(2,2), activation='relu', padding='same')(cnv2d19)
cnv2dtra = concatenate([cnv2dtra , cnv2d14 ], axis=1)
# 分岐結合
cnv2dtrb = Conv2D(num_hidden3a, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2dtra)
cnv2dtr1 = Conv2DTranspose(num_hidden3, data_format='channels_first', kernel_size=(3,3), strides=(2,2), activation='relu', padding='same')(cnv2dtrb)
cnv2dtr1 = concatenate([cnv2dtr1 , cnv2d10 ], axis=1)
# 分岐結合
cnv2dtr2 = Conv2D(num_hidden3, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2dtr1)
cnv2dtr4 = Conv2DTranspose(num_hidden2, data_format='channels_first', kernel_size=(3,3), strides=(2,2), activation='relu', padding='same')(cnv2dtr2)
cnv2dtr4 = concatenate([cnv2dtr4, cnv2d6], axis=1)
# 分岐結合
cnv2dtr5 = Conv2D(num_hidden2, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2dtr4)
cnv2dtr7 = Conv2DTranspose(num_hidden1, data_format='channels_first', kernel_size=(3,3), strides=(2,2), activation='relu', padding='same')(cnv2dtr5)
cnv2dtr7 = concatenate([cnv2dtr7, cnv2d2], axis=1)
# 分岐結合
cnv2dtr8 = Conv2D(num_hidden1, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2dtr7)
cnv2dtr9 = Conv2D(num_hidden1, data_format='channels_first', kernel_size=(3,3), activation='relu', padding='same')(cnv2dtr8)
cnv2dtr10 = Conv2D(3, data_format='channels_first', kernel_size=(1,1), activation=softMaxAxis , padding='same')(cnv2dtr8)
# 最終的に1枚の3値の画像にする
NN_1 = Model( input=cnn_input , output=cnv2dtr10 )
# 入力はcnn_inputで最後のcnv2dtr10を出力としてネットワークを構築
NN_1.compile(optimizer='adam', loss='categorical_crossentropy' , metrics=['accuracy'])
# ロス関数 categorical_crossentropy
4.4 最終層のsoftmax
最終層では、各ピクセルについて赤青白をsoftmaxで計算します。
標準のsoftmaxでは配列の一部を指定することができないので、softMaxAxisという関数を定義しています。
この方法は
Stack Over flow:"How to specify the axis when using the softmax activation in a Keras layer?"
を参考にしています。
def softMaxAxis(x):
return softmax(x,axis=1)
このとき、モデルを保存して再読み込みさせるときに、カスタムオブジェクトの名前を教えてあげる必要があるようです。
NN_1 = load_model( paramfile , custom_objects={'softMaxAxis': softMaxAxis })
4.5 生成画像と等圧線図の重ね合わせ
ImageMagickを利用します。
convert -compose multiply が重ね合わせを行うコマンドです。
-gravity center -geometry +0+0 により2つの画像の中心を縦横のオフセット無しで合わせることを指定します。
# m_file 生成した前線画像（PNG)
# s_file 等圧線図（PNG)
composite -gravity center -geometry +0+0 -compose multiply ${m_file} ${s_file} out.png
5.まとめ
気象データ可視化画像を認識して、「日本気象庁の解析するような前線画像」を生成するニューラルネットワークを作成し、前線の自動描画をやってみました。
2021.09
日本付近の描画だけだったものを、「日本気象庁の解析するような」なやり方を学習したネットワークとして世界の他の地域へ展開し、全球自動前線描画を試みました。
「天気図っぽい」気象前線の自動描画を日本から世界に拡張してみる〜領域拡張と地図投影
この手法は前線描画に限らず、地図上にプロットされる何らかのデータで気象と関連するものであれば有効なのではないかと思っています。
一方で、当たりハズレの評価が結構難しい気がします。機械的にピクセル値のズレを見ても違う気がしますし、前線という形を捉えていることの評価をどうするか。気象予報士試験のように、人が採点するのが案外正しいかもしれません。
現在、MSMの気象データから合成レーダの画像を生成するということに取り組んでいて、それについては正解データとの評価も検討しております。
こちらも近日公開しようかと思っています。
2020.04
こちらは次のような内容になりました。
変分オートエンコーダを用いた気象データ画像の生成 Generating Weather Data Images Using Variational Autoencoder
今回Qiita初投稿でしたが、５回にわたって長文を投稿してしまいした。
読んで頂いた皆様どうもありがとうございました。