天気図っぽい気象前線の自動描画を日本から世界に拡張してみる〜領域拡張と地図投影

目次
0. はじめに
1. 日本付近の気象データから「天気図っぽい」前線を自動描画する（機械学習）
2. 学習したネットワーク・パラメタを他の地域の気象データに適用する(機械学習と地図）
3. ステレオ投影された前線画像から正投影地図への座標変換(地図）
4. pcolormeshを用いての画像の貼り付け(地図と可視化）
5. まとめ

0. はじめに

以前に数値気象予報モデルのデータから、気象庁が解析した「気象前線」を、機械学習（ニューラルネットワーク）を用いて自動的に描画する、というプログラムを作成したことがありました。

これは日本付近のデータと日本付近の天気図データをもとにしたものだったのですが、オリンピックも始まったことだし、せっかく日本で学習した結果を世界に拡張したらどうなるか？と思い立ってやってみました。

実は世界各地で行われている気象前線解析にはそれぞれの国や機関の流儀があるようで、ヨーロッパやアメリカの前線の解析の仕方は日本と異なりますので、いわば「日本流解析」の世界拡張ということになります。

日本的な気象前線解析を学習したAIがどこまで世界の気象データを日本流に解析できるでしょうか。

結論からいうと現在のところ 北半球 南北両半球の中緯度帯まで拡張してみたところ、それらしく前線が出てきました。

こんな感じです。

梅雨前線が日本周辺からさらに西側に伸びていたり、太平洋北部を進む低気圧の前線も見えています。
ヨーロパや北米大陸の低気圧にも前線が解析されました。大西洋域の高気圧の西縁の前線はやや？なところですが。

2021.7.25更新　南半球の前線も描画しました。

さて今回は、この図が出てくるまでに下記のような地図の投影で少し苦労した話を、前線学習の話とともに纏めました。

地図の投影での苦労
・投影後の画像を別形式地図に投影するための座標変換
　Basemapの座標変換機能を用います
・この座標変換により使えなくなる画像貼り付け機能の代替機能の工夫
　カラー画像をグレー変換して、自作したカラーマップを用いてpcolormeshを用います

1. 日本付近の気象データから「天気図っぽい」前線を自動描画する（機械学習）

詳細は以前ここに纏めています。

気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（5) 機械学習編 Automatic Front Detection in Weather Data

簡単に言ってしまうと、

気象データ :気象庁が日々作成している数値気象予報モデルのデータを画像化したもの。
天気図っぽい前線 :気象庁が日々作成している「速報天気図(SPAS)」という天気図の前線部分を抜き出したもの。

というデータの対を大量に用意して、それらの間の関係をニューラルネットワーク（CNN）によりセマンティックセグメンテーションとして学習させるわけです。気象データ画像のピクセル単位に寒冷前線・温暖前線・閉塞前線・それ以外という分類を行わせることで、気象前線を描画するものです。学習データは3年9ヶ月（3枚/日）程度の約5000組です。

GPVデータはいつもいつものことですが、京都大学生存圏研究所様から利用させていただいております。

2. 学習したネットワーク・パラメタを他の地域の気象データに適用する

日本付近のデータで学習したものが他地域で通じるのか？という疑問も湧きますが、北半球の同じ緯度帯ならまず大丈夫では？と考えてすすめることにしました。
その後、南半球についても少し工夫することで前線描画ができました（2021.7.25追記）。

2.1 学習時のデータの地理的な範囲について

元々学習に使用したデータは日本付近のものです。例えば以下のような画像です。

2019年1月の地上気圧・気温・風

2019年1月の850hPaでの相当温位

上の例では画像を連ねてアニメーションGIFにしています。

こういう画像を2018年8月から2021年4月までの3年半強使って、「天気図っぽい」前線を描く学習を行わせています。

(上段左より 地上気圧 850hPa相当温位 700hPa湿数、中段左より 850hPa気圧・温度・風　700hPa鉛直流 500hPa気圧・温度・風、下段左より 地上気圧・温度・風 300hPa風速）

入力するデータの種類は、以前にやったものと少し違っていて(*)、下記の7種類です。
それぞれ3チャンネルの画像データなので、ニューラルネットワークの入力としてはこれらをConcatenateした21チャネルの画像情報となります。

No	高度	項目
1	地上	気圧・気温・風ベクトル
2	850hPa	気圧(ジオポテンシャル高度）・気温・風ベクトル
3	850hPa	相当温位
4	700hPa	湿数
5	700hPa	鉛直風速
6	500hPa	気圧(ジオポテンシャル高度）・気温・風ベクトル
7	300hPa	風速分布

(*)以前との違い
増やしたのはNo.7の300hPaでの風速分布で、上図の下段右のようなものです。
以前から閉塞前線も分類として増やしました。
気象学的に閉塞点と高層ジェット気流の速度の強いところには関係がありますので、ニューラルネットワークにとってヒントとなるようにと思ったためです。

2.2 他の地域の入力データの作成

2.1 で例示したデータはBasemapのステレオ投影による記法を用いて作成しています。
入力するGPVデータは全球分ありますので、描画する対象の緯度経度を変更してあげれば簡単に作成できます。

なお、格子点形式の気象データ（GPV)をpygribを用いてデコードして、Basemapやmatplotlibを用いて可視化するまでについては、気象データをもとに「天気図っぽい前線」を機械学習で描いてみる(2)に纏めております。

Basemap.py

from mpl_toolkits.basemap import Basemap

m = Basemap(projection='stere', llcrnrlat=9, urcrnrlat=54, 
   llcrnrlon=115, urcrnrlon=178, lat_0=60,  lon_0=140, resolution='i' )

（緯度、経度）= (60, 140) を投影直下点(lat_0, lon_0)として、
左下コーナー(llcrnrlat, llcrnrlon)、
右上コーナー(urcrnrlat, urcrnrlon)
について(9,115),(54,178)という指定をしています。

こうすると例えば下記のような図が描画されます（気圧図を重ねています）。

2021年6月1日0時(UTC)の地上気圧

他の地域データを作るには、これらの投影直下点、左下コーナー、右上コーナーを変更すればOKです。

北米大陸付近（投影直下点 西経100度、北緯60度）

ヨーロッパ・北アフリカ付近（投影直下点 東経20度、北緯60度）

中央アジア付近（投影直下点　東経100度、北緯60度）

このように、範囲としては経度を40度づつずらして、北半球一周で9個の領域を対象としました。

2.3 南半球の入力データ

2021.7.25更新

南半球も同じ用に経度が40度づつずれた領域を9個作成しました。

ところで南半球の気象は、北半球を東西反転したような画像になります。
今回のニューラルネットワークは北半球で学習しているため、この反転されたデータに対しては前線描画がうまくできません。そこで、南半球の画像は左右反転させてニューラルネットワークに入力して、出力された前線画像もまた左右反転させて用いることとしました。

左右反転されたオーストラリア付近(上が南になっています）

反転前の上と同じエリア（地図としてはこれが正しい）

前線検出例を下に示します。
反転させたデータの方が、前線らしいものを検出しています。

左右反転した入力データを用いて前線検出させた例

左右反転しない入力データを用いて前線検出させた例

2.4 セマンティックセグメンテーション実行結果

北半球中緯度帯ということで気象パターンが似ているせいか、日本付近しか見たことがない日本代表AIにしては（思ったよりも）それなりな前線が描画されました。

北米領域(上段左より 地上気圧 850hPa相当温位 700hPa湿数、中段左より 850hPa気圧・温度・風　700hPa鉛直流 500hPa気圧・温度・風、下段左より 地上気圧・温度・風 300hPa風速）

大西洋領域(図の並びは北米領域と同じ）

大西洋海域にも速報天気図（SPAS）と同じような前線記号、高気圧・低気圧の記号が描画されました。

日本代表AIとしては初めて見る地図に惑わされてはいないようです。やはり相当温位の等温位線の密な部分、湿数の大きなところを前線を引くべきところと思って、日本で学んだ知識を活かして描画しているようです(笑）。

3. ステレオ投影された前線画像から正投影地図への座標変換

せっかく日本代表として描画してくれた図をもう少し見やすくしたいところ。
AIは学んだサイズで、領域別に前線画像を出力しているのですが、これらを地球全体の地図に貼り付けたほうが見やすいというものです。これはAIには教えていないのでこちらでやってあげます。

ここでひとつ考えなければならない話があります。

ニューラルネットワークが出力するのは下図左のような「ステレオ投影済みの画像」です。
この画像の中のそれぞれの点の緯度・経度が分からないと、右のような地球儀への正投影ができません。緯度線、経度線を引いてみるとわかるように、「ステレオ投影済み画像」の中で緯度経度が直交座標系にはなっていないので、この画像から元の緯度経度へ座標変換をしてやる必要があります。

3.1 ステレオ投影済画像の画像座標系から緯度経度への変換

Basemapが持っている座標変換機能を用います。

まず出力画像と同じステレオ投影形式のBasemapを作成します。

stereo_basemap.py

# 左下、右上、投影直下の経度
ll_lon, ur_lon, ss_lon  = 117 , 172 , 140

# ステレオ投影のBasemapを作成  
  m0 = Basemap(projection='stere', llcrnrlat=12, urcrnrlat=54,
               llcrnrlon=ll_lon, urcrnrlon=ur_lon, 
               lat_0=60,  lon_0=ss_lon, resolution='i')

（ll_lon, ur_lon, ss_lonを変数にしているのは別領域用に入れ替えやすくしたため）

この準備をした上で、ステレオ投影の地図に(256+1)x(256+1)のグリッドを割り当てます。

makegrid.py

lons, lats, x, y = m0.makegrid(256+1,256+1,returnxy=True)

こうすることで、ステレオ投影画像の縦横に256×256の等間隔格子を作成し、かつそれぞれの格子が該当する緯度・経度を取得する事ができます。

3.2 正投影図上のxy座標への変換

次に、投影する全球地図のBasemapを作成します。
全球のデータを一度に見えるように、東半球と西半球それぞれの投影画像を作成します。

globalmap.py

  # 1行2列の図を作成します。
  fig, ax = plt.subplots(1, 2, figsize=(12,6))
  plt.subplots_adjust(left=0.01, right=0.98, top=0.99, bottom=0.01)
  tstr = タイトル文字列を設定

  # 投影直下点を北緯10度、東経120度とした半球の正投影図
  m1 = Basemap(projection='ortho', lat_0=10, lon_0=120, 
       resolution='i', ax=ax[0])
  ax[0].set_title(tstr)                # タイトルを設定
  im1=m1.drawcoastlines(linewidth=0.5) # 海岸線を描画

  # 同様に投影直下点を北緯10度、西経60度とした半球の正投影図
  m2 = Basemap(projection='ortho', lat_0=10, lon_0=-60, 
       resolution='i', ax=ax[1])
  ax[1].set_title(tstr)
  im2=m2.drawcoastlines(linewidth=0.5)

この準備をしておくことで、3.1で得られた緯度経度lons, latsから、Basemapの機能を用いて

xy_trans.py

xx, yy = m1(lons, lats)

とすることで、m1の正投影図の上でのxy座標（xx, yy)に変換されます。

4 pcolormeshを用いての画像の貼り付け

座標変換もできたので、下記のようにimshowの位置指定機能を用いて前線画像を貼り付けよう！と考えたのですが、これはうまくいきません。

imsho.py

### うまく行かないコードです！！

# 前線画像ファイルをオープンしてX座標を反転させておく
t_img   = 
(Image.open(sampleimg_file).convert('RGB')).resize((256,256),Image.HAMMING)
t_array_r = np.array(t_img)/255
t_array = t_array_r[::-1,:]

#画像をm1のBasemapに貼り付けて、array_extentにより画像の左下と右上コーナーを指定する。
x0, y0       = m1(117, 12)
x1, y1       = m1(172, 54)
array_extent = (x0, x1, y0, y1)

im1 = m1.imshow(t_array,origin="lower")
im1.set_extent(array_extent)

なぜならば、imshowで指定できるset_extentでは、左下と左上、右上と右下の経度座標は同じでなければならないからです（同様に左上と右上、左下と右下の緯度座標も同じでなければなりません）。
つまり下の図の一番左側の矩形イメージをimshowによって投影図に貼り付けることは可能なのですが、中央の図のようになってしまいます。本来は、右側の図のようにならなければなりません。左上と右下の位置が違っていることがわかります（わかりやすいように背景に地形図を重ねています）。

そこで画像配列を緯度経度にもとづき、pcolormeshを用いて球面上に「塗る」ことにしました。
（右図は実際にそのようにして作成しています）

4.2 カラーマップ(pcolormesh)による画像描画

pcolormeshでは、格子状に並んだデータと、その格子を割り当てる座標を指定することができます。
このため4.1で示したような少し歪んだ並びについてもうまく塗ることができるわけです。
ただ、データについてはRGBのような3チャンネルではなくひとつの値となっていなければなりません。

そこで、RGB形式の出力データをグレースケールに変換して1チャンネルのデータとして、その値を表現するカラーマップを自作して元の色を再現することにしました。

このあたりの方法に関しては、下記サイトで纏めておられた先達のお知恵を拝借しています。

【Python/Pillow(PIL)】カラー,モノクロ,HSVなどの変換
不連続カラーバーの作成(Python)

まずグレースケール変換後に赤、青、ピンク、白がそれぞれどんな数値に変換されるかを調査します。

colormap.py

>>> from PIL import Image

# 4x1 という調査用のRGB画像を用意
>>> img_rgb = Image.new("RGB",(4,1))

>>> img_rgb.putpixel((0,0),(255,0,0)) # 赤色（温暖前線）をセット
>>> img_rgb.putpixel((1,0),(0,0,255)) # 青色（寒冷前線）をセット
>>> img_rgb.putpixel((2,0),(255,0,255)) # ピンク（閉塞前線）をセット
>>> img_rgb.putpixel((3,0),(255,255,255)) # 白（前線以外）をセット

# 画像をグレースケールに変換
>>> img_L = img_rgb.convert("L") 

# 変換後にそれぞれの色がどういう数値に変換されたかを調査
>>> print(img_L.getpixel((0,0)))
76 # 赤は76
>>> print(img_L.getpixel((1,0)))
29 # 青は29
>>> print(img_L.getpixel((2,0)))
105 # ピンクは105
>>> print(img_L.getpixel((3,0)))
255 # 白は255

前線画像(RGB)をグレースケール変換したデータ配列がこのような値に変換されることを見越して、pcolormeshで用いるカラーマップを作成します。
変換の際は上で調べた値を挟んで色が再現される範囲をもたせるように敷居を設定します。

custom_cmap.py

# しきい値を設定。青(29:20-70),赤(76:70-95),ピンク(105:95-240),白(255:240-255)
level_rgb = [20/255,70/255,95/255,240/255,255/255]

# 各しきいのカラーマップデータをRGBで指定。白は背景画像が透過されるように透明にする。
cmap_data = [(0,0,1),(1,0,0),(1,0,1),(1,1,1,0)]

# camp_orig という名のカラーマップを作成。
cmap_orig = cls.ListedColormap(cmap_data)
cmap_orig.set_under((1,1,1,0)) # しきい値より下の値がきたら透明白で塗る
cmap_orig.set_over((1,1,1,0))  # しきい値より上の値がきたら透明白で塗る

# camp_origとしきい値をくくりつける
norm      = cls.BoundaryNorm( level_rgb, cmap_orig.N )

このカラーマップを用いると、pcolormeshで使用できるように1チャンネルの配列化したデータから元のカラー図に戻すことができます。

前線画像をオープンするときにグレースケール化します。

pcolormesh.py

# 前線画像(samplefile)をグレースケールでオープンし、
# t_array_Lというndarrayに変換する
t_img_L
 = (Image.open(samplefile).convert('L')).resize((256,256),Image.HAMMING)

t_array_L_t   = np.asarray(t_img_L)/255
t_array_L     = t_array_L_t[::-1,:]

# 描画準備
fig, ax = plt.subplots(1, 2, figsize=(18,6))
plt.subplots_adjust(left=0.01, right=0.98, top=0.99, bottom=0.01)

# Basemapで正投影（m2）とステレオ投影(m0)を作成する。
# 前線画像はm0の形式で作成されている
m2 = Basemap(projection='ortho', lat_0=35, lon_0=135, 
     resolution='i', ax=ax[1])
m0 = Basemap(projection='stere', llcrnrlat=12, urcrnrlat=54, 
     llcrnrlon=115, urcrnrlon=178, lat_0=60,  lon_0=140, 
     resolution='i' , ax=ax[0])

# 前線画像に256x256格子を割り当て、それぞれの格子の緯度経度を取得する
lons, lats, x, y = m0.makegrid(256+1,256+1,returnxy=True)

# もとの前線画像をimshowでステレオ投影の地図に貼り付ける
# これは読み取った元画像なので緯度経度と画像座標は一致している
im0 = m0.imshow(t_array_L, origin="lower", cmap="gray")
im0 = m0.drawcoastlines(linewidth=0.4)

# 緯度経度を正投影図の画像座標に変換する
xx, yy = m2(lons, lats)

# pcolormeshで作成したカラーマップを用いてグレースケールデータを表示
im2=m2.shadedrelief()
im2=m2.pcolormesh(xx, yy, t_array_L, cmap=cmap_orig, norm=norm)
fig.colorbar(im2, ax=ax[1])

下図左のようにグレースケール化したデータを、右図のカラーマップでpcolormeshを用いて表示することで元の図の色を再現しました。白色については、地上気圧図と重ねることを想定して透明にしてあります。

4.3 地上気圧図との重ね合わせ

あとは天気図らしく地上気圧図を重ね合わせます。
ここにも気象図の可視化に関しては書いていますが、pygribというライブラリを用います。

GPVデータもそれ自体が緯度・経度ごとのデータなので、正投影図用に座標変換します。

mslp.py

import pygrib

# 等気圧線の間隔を定義
levels_prs = np.arange(930.0 , 1080.0, 4.0)

# gpvfileにセットしたGRIB形式のGPVファイルをオープン
grbs = pygrib.open(gpvfile)

# GPVデータから地上気圧(海面更正気圧）を抜き出す
letter_mslp = "Pressure reduced to MSL"
grb_prs = grbs.select(parameterName=letter_mslp , level=0 , 
forecastTime=fcsttime)

# GPVデータ自体の座標を取得し、投影図用に変換
gpvlats, gpvlons = grb_prs[0].latlons()
x1, y1 = m1(gpvlons, gpvlats)
x2, y2 = m2(gpvlons, gpvlats)

# ラベル付きの等圧線を描画。GRIBの気圧はPa単位なので100で割ってhPa単位で表示。
plt.clabel(m1.contour(x1, y1, grb_prs[0].values / 100.0 , levels_prs , 
linewidths=0.3, colors='k') , fontsize=4, inline=1, fmt='%4.0f')
# 等圧線の上から色付きのコンターを重ね合わせ
m1.contourf(x1, y1, grb_prs[0].values / 100.0 , levels_prs , 
cmap="coolwarm" )

# 別半球にも同じことを実施
plt.clabel(m2.contour(x2, y2, grb_prs[0].values / 100.0 , levels_prs , 
linewidths=0.3, colors='k') , fontsize=4, inline=1, fmt='%4.0f')
m2.contourf(x2, y2, grb_prs[0].values / 100.0 , levels_prs , 
cmap="coolwarm" )

# このあとにpcolormeshを用いた前線の描画が続きます。

4.4 その他

pcolormeshのデータをBasemapの正投影図に渡した際に、見えない側の半球のデータについても描画をしようとしてしまい、不要かつ奇妙な線や色が表示されることがあります。
この問題というか仕様の回避が下記で議論されておりました。

Problem with ortho projection and pcolormesh #470

結論としては不要な部分のデータはnp.nanにセットしてマスクすべし！ということです。
このために下記のようなコードを追加します。

AvoidingOrthoProjectionProblem.py

xx, yy = m2(lons, lats)
#　不要な部分(画像座標が下記を満たす部分）のインデックスを取得
ii = ( \
      ( xx[:-1, :-1] > 1e20) | \
      ( xx[1: , :-1] > 1e20) | \
      ( xx[:-1, 1: ] > 1e20) | \
      ( xx[1: , 1: ] > 1e20) | \
      ( yy[:-1, :-1] > 1e20) | \
      ( yy[1: , :-1] > 1e20) | \
      ( yy[:-1, 1: ] > 1e20) | \
      ( yy[1: , 1: ] > 1e20)   \
     )

# 上記条件のインデックスのデータをnp.nanでマスク
t_array_2[ii] = np.nan

# pcolormesh を実行する
im2=m2.pcolormesh(xx, yy, t_array_2, cmap=cmap, norm=norm)

こうして最終的に下図のような画像ができます。
下図では予報時間別に作った図をアニメーションGIFにしています。

2021/7/2 UTC12時を初期値とした予報値に前線検出を行ったもの

2021/7/16 UTC12時を初期値とした予報値に前線検出を行ったもの

2021/6/19 UTC12時を初期値とした予報値に前線検出を行ったもの

4. まとめ

日本付近の気象図（SPAS）を用いて前線を自動描画するよう学習させたニューラルネットワークを北半球の他の領域に適用して日本的な前線解析をした全球画像を作成しました。

全球画像作成にあたっては、
・投影後の画像を別形式地図に投影するための座標変換
　Basemapの座標変換機能を用います
・この座標変換により使えなくなる画像貼り付け機能の代替機能の工夫
　カラー画像をグレー変換して、自作したカラーマップを用いてpcolormeshを用います
という工夫を行いました。

今後、南半球の同緯度帯や、別の緯度帯にも展開していきたいです。
いまのところ、南半球の同緯度帯をうまくステレオ表示できず（南が地図の上部にくるステレオ投影ができない？）、入力データの作成ができていません。

2021.7.25修正
南半球の同緯度帯も入力画像を作ることができました。普通に左下、右上を指定すると南が上に来る地図が出力されます。

長文・乱文のところ最後まで読んでいただいてありがとうございました。