ベクトル形式の天気図から前線記号などを抜き出す（アメリカ気象局天気図）

0.はじめに

　訳あってしばらくアメリカで生活をすることになったので、アメリカの気象データについても色々と見ていきたいと思います。アメリカの気象機関はなんといってもまずはNOAA、National Oceanography and Atomspheric Administry（アメリカ海洋大気庁）です。NOAAはアメリカ商務省に属する機関です。ちなみに日本の気象庁は昔の運輸省、現在は国土交通省の傘下ですね。

NOAAの下には日本の気象庁に相当するNational Weather Service（国立気象局）など、いくつかの組織があります。観測や予報のデータを広く公開していて、ありがたいことに"free and available to the public without restriction on use"であるとの記載が見られます。

今回は、ベクトル形式天気図から前線記号などを抜き出す（1）気象庁SVG形式天気図編に引き続き、アメリカの天気図から前線記号を抜き出す話です。ベクトル形式天気図と書きましたが、正確にはベクトル形式天気図ではありません。ですが、前線記号などとその位置を取り出す方法であることには変わりありません。

1. アメリカの地上解析天気図

気象予報センター(Weather Prediction Center)が出している天気図のアーカイブページ(地上気象解析, Surface Analysis)がこちらです。

北米大陸全体・アメリカ合衆国部分など範囲を変えたり、地上観測点のプロット・衛星・レーダーなどレイヤーをオーバーラップさせたりと、色々な形式で天気図を見ることができます。選択肢が多く無料で見られるて素晴らしい。
私の目的は「気象データから前線を自動描画する」というものですので、この目的に沿った天気図を使用します。

下記は、2014年3月5日0時（世界協定時）のものです。

各種の前線（温暖前線、寒冷前線、停滞前線と閉塞前線）や、高気圧（H）および低気圧（L）が解析されています。
また、ところどころにオレンジ色の点線が記載されています。これはトラフ、つまり「気圧の谷」ラインです。

しかし前線解析の流儀は国によって違います。例えば西海岸ワシントン州の沖合には停滞前線と、閉塞前線＋温暖前線が並行に解析されています。このようなパターンの描画は、日本ではあまり見かけないと思います。

それにしても国土が広い。広範囲な高気圧・低気圧と前線パターンがすっぽりと国土の上に載っています。
余談ですが、大陸の国アメリカでは起きる気象現象も結構違っています。私が住むマサチューセッツ州付近は海に近いとはいえ、日本に比べると乾燥しています。癖毛の家族もサラサラで過ごしています。その分肌の乾燥も強めですが。

2.気象要素の切り抜き

2.1 画像形式データからの切り抜き

この天気図データは画像形式（ラスタ形式）です。ピクセルの値が赤・青・紫・オレンジのものを抜き出すと、以下のように気象要素のみを抜き出すことができます。

左：元の天気図(GIF形式）右：色を元に抜き出したもの

色で切り抜くコードは下に折り畳んで示します。
ポイントは、図の大きさが日によって微妙に異なったりしますのでその対応処理が入っていること（笑）。
アメリカのデータにはこういうケースがままあり、多少データフォーマットに対するレジリエンシーが必要になります。

ピクセルの色での切り抜きコード

アメリカ天気図切り抜き

import copy
import random
import sys 
import matplotlib.pyplot as plt 
import matplotlib.cm as cm
from matplotlib.animation import ArtistAnimation
import numpy as np
import math
import csv 
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from PIL import Image,ImageFilter,ImageOps,ImageChops

#
    # 天気図ファイル t_imgfileをオープン
    t_img = Image.open(t_imgfile)

    img_size = t_img.size
    
    #RGB形式に変換する
    m_img = t_img.convert('RGB')

    #最後に出力するファイルをRGB形式で作成する
    mask1 =  Image.new('RGB', img_size)

    #Pixelごとに値を読み取る
	for x in range(0,img_size[0]):
		for y in range(0,img_size[1]):
    		r0,g0,b0 = m_img.getpixel((x,y))

			"""
			 Typical Front: Red(243/21/2) Blue(24/0/255) Purple(125/0/234)
			 Other: Number(118/10/-0) Line(245/108/56) Border(255/255/255) 
			"""
			r,g,b=255,255,255

			if r0 > 250 and g0 <  230 and b0 < 230   : # Set to BLOWN
    			r,g,b = 255,130,71
			if r0 > 240 and g0 <  30 and b0 < 30   : # Set to RED
				r,g,b = 255, 0, 0
			if r0 < 230 and g0 < 200 and b0 > 220 : # Set to Purple
				r,g,b = 125, 0, 234
			if r0 < 220 and g0 < 220 and b0 > 250 : # Set to BLUE
				r,g,b = 0,   0, 255

			# NOAAのロゴなどのエリアを白で塗りつぶす
			if x > 0 and x < 60 and y > 424 and y < 470:
				r,g,b = 255, 255, 255

				mask1.putpixel((x,y),(r,g,b))

    # 大きさが日によって違うことがあるので、調整する。
	tmp_img_to_array = np.asarray(mask1)

	if tmp_img_to_array.shape[1] == 750:
		add_img_to_array = tmp_img_to_array[:,1:749,:]
	else:
		add_img_to_array = tmp_img_to_array
			

    # イメージ形式に戻す
	m_s_img = Image.fromarray(np.uint8(add_img_to_array),mode='RGB')
	m_s_img.save(o_filelist[c_sfile])

しかしながら、ピクセルの色では気象要素毎に抜き出すことは難しく、温暖前線・寒冷前線と停滞前線や高低気圧記号が同じ分類となってしまいます。

2.2 Coded Surface Bulletinからの切り抜き

様々な天気図が掲載されているページに、「Coded Surface Bulletin」という情報が掲載されています。「符号化地表解析通報」とでも訳す感じですかね。以下は単に通報と称します。

日本でも、SYNOPとかTEMPといった符号化された気象通報式がありますが、同じように前線や高気圧低気圧の位置を一定の規約に従ってテキスト化したものです。

例を示します。

CODED SUS


CODED SURFACE FRONTAL POSITIONS
NWS WEATHER PREDICTION CENTER COLLEGE PARK MD
922 PM EDT SUN OCT 29 2023
 
VALID 103000Z
HIGHS 1031 53130 1022 4671 1023 4873 1017 3383 1029 64116 1021 68152 1029 
62145 1030 64137 1033 57128 1035 52120 1037 51117 1038 46113 1036 44108 1034 
41104 1033 42114 1027 57115 1024 39130 1019 3279 
LOWS 1007 58161 1007 56157 1000 25103 1015 4083 1014 3877 1014 3589 1013 3293 
1010 28100 1019 38107 1011 34107 1011 34109 1008 5877 1013 6596 1010 6977 
1011 56147 1024 61138 1012 54145 1015 4079 961 5652 966 5856 1002 7571 
TROF 1795 1893 1991 2090 2189 
COLD WK 28100 28101 28103 28104 28105 29106 30106 31107 
COLD WK 33109 32108 31109 30110 29112 

　（中略）

OCFNT WK 5856 5956 6054 6051 5844 5539 5136 
TROF 7572 7384 7298 
TROF 8171 7974 7671 
COLD WK 3962 3864 3966 
 
$$

この通報の復号方法が、

というところに解説されています。

比較的簡単な複合で読み解けます。低解像度版（緯度経度を度単位）と高解像度版（緯度経度を小数点1位の度単位）がありますが、低解像度版はHIGHS 1031 53130なら　1031hPaの高気圧が、北緯53度、西経130度にある。という意味です。高精度版はHIGHS 1031 5291303ならば、北緯53.9度西経130.3度です。
前線については、通過地点の緯度経度が同様に記載されています。

そして、この通報を解読して天気図を描画するプログラムがMetPyの1.5版に登場しています。

2.3 MetPyを用いた気象テキストの解読

MetPyは気象に関する計算や図の描画に便利なライブラリです。ライセンス条件はここに記載があります。

この中で、metpy.io.parse_wpc_surface_bulletinというライブラリがあって、通報のテキストファイルを読んでpandas形式のデータに変換して返します。

parse_wpc_surface_bulletin


>>> import os
>>> from metpy.io import parse_wpc_surface_bulletin
>>> import pandas
>>> df = parse_wpc_surface_bulletin('file_name_of_bulletin')
>>> type(df)
<class 'pandas.core.frame.DataFrame'>
>>> df
                 valid feature  strength                                       geometry
0  2023-10-30 12:00:00    HIGH    1025.0                            POINT (-117.7 35.7)
1  2023-10-30 12:00:00    HIGH    1030.0                            POINT (-119.3 38.6)
2  2023-10-30 12:00:00    HIGH    1037.0                            POINT (-119.1 51.7)
3  2023-10-30 12:00:00    HIGH    1034.0                            POINT (-123.6 53.5)
4  2023-10-30 12:00:00    HIGH    1038.0                            POINT (-113.6 45.9)
..                 ...     ...       ...                                                ...
72 2023-10-30 12:00:00    COLD       NaN            LINESTRING (-134.4 40, -134.5 36.9)
73 2023-10-30 12:00:00    TROF       NaN              LINESTRING (-174 61.6, -170 60.5)
74 2023-10-30 12:00:00    TROF       NaN              LINESTRING (-176 62, -179.7 62.8)
75 2023-10-30 12:00:00   STNRY       NaN  LINESTRING (-72.5 40.5, -73.2 40.1, -74.099999...
76 2023-10-30 12:00:00    TROF       NaN  LINESTRING (-117.6 37.8, -118.4 39.4, -119.8 4...

[77 rows x 4 columns]

pandasのDataFrame形式でデータが返ってきています。featureに応じて、geometryのPOINTやLINESTRINGを処理すれば良いのですが、これらを描画するためのスタイルもMetPyに用意されました。Basemapやmatplotlibを用いて天気図を作成することができるようになっており、MetPyのチュートリアルには以下のように描画例とサンプルコードが掲載されています。

Tutorial of Plotting Fronts -MetPy 1.5から

つまり、この通報が一種のベクトル形式の天気図として利用できるということになります。

2.4 MetPyを用いた前線記号の抜き出し

さてこれらの便利なライブラリを利用すると、簡単に記号を抜き出せます。
抜き出す、というより単に、地図を描画せずに記号を描画する、というだけ。

ソースコードはMetPyのサンプルをベースに、Basemapを用いた描画に変更しています。

通報を読み込んで、parse_wpc_surface_bulletinによってpandasのDataFrame形式にしたデータをmatplotlibを用いて描画していく、という処理です。

plot with MetPy

import os
import sys
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator

from metpy.cbook import get_test_data
from metpy.io import parse_wpc_surface_bulletin
from metpy.plots import (\
　add_metpy_logo, ColdFront, OccludedFront, StationaryFront,\
                             StationPlot, WarmFront)

def plot_bulletin(ax, data, plt):
  # 天気図要素の描画
  size = 2.5
  fontsize = 12
  complete_style =　\
   {'HIGH': {'color': 'blue', 'fontsize': fontsize}, \
    'LOW':  {'color': 'red', 'fontsize': fontsize}, \
    'WARM': {'linewidth': 0.7, 'path_effects': [WarmFront(size=size)]}, \
    'COLD': {'linewidth': 0.7, 'path_effects': [ColdFront(size=size)]}, \
    'OCFNT':{'linewidth': 0.7, 'path_effects': [OccludedFront(size=size)]}, \
    'STNRY':{'linewidth': 0.7, 'path_effects': [StationaryFront(size=size)]}, \
    'TROF': {'linewidth': 0.7, 'linestyle': 'dashed' , 'color': 'darkorange'}}

  # 高気圧・低気圧のプロット 
  for field in ('HIGH', 'LOW'):
      rows = data[data.feature == field]
      if rows.shape[0] > 0 :
        x, y = zip(*((pt.x, pt.y) for pt in rows.geometry))
        px, py = ax(x, y)
        for ppx, ppy in zip(*(px, py)) :
          if field[0] == "H":
            plt.text(ppx, ppy, field[0] , \
            fontsize=12,fontweight='bold', \
            ha='left',va='bottom',color='k') 
          if field[0] == "L":
            plt.text(ppx, ppy, field[0] , 
            fontsize=12,fontweight='bold', \
            ha='left',va='bottom',color='c') 

  # 前線やトラフの線を描画
  for field in ('WARM', 'COLD', 'STNRY', 'OCFNT', 'TROF'):
      rows = data[data.feature == field]
      if rows.shape[0] > 0 :
        for pt in rows.geometry :
          xx, yy = pt.xy
          pxx, pyy = ax(xx, yy)
          ax.plot(pxx, pyy, **complete_style[field])



# Basemapを利用して地図を作成
m = Basemap(projection='stere',\
  llcrnrlat=18, urcrnrlat=51, llcrnrlon=-126, urcrnrlon=-54, \
  lat_0=40,  lon_0=-100, resolution='i' )

fig = plt.subplots(figsize=(5.8,4.2))
plt.subplots_adjust(left=0.01, right=0.98, top=0.99, bottom=0.01)

# 地図装飾の描画
m.drawparallels(np.arange(-80.,81.,10.))   #緯度線
m.drawmeridians(np.arange(-180.,181.,10.)) #経度線
m.drawcoastlines(linewidth=0.5)            #海岸線
m.drawstates(linewidth=0.3)                #州境界
m.drawcountries(linewidth=0.3)             #国境

# 通報テキスト読み込み
df = parse_wpc_surface_bulletin(codfilename)

# 描画ルーチンの呼び出し
plot_bulletin(m, df, plt)


plt.savefig(outfilename)
plt.close()

このようにして抜き出した前線要素をBasemapを用いて地図に描画したものを以下に示します。

上段：NOAA/NWSによる解析の天気図下段：同じ時刻の描画したもの

またも余談ですが、よく見ると上下の図では地図の投影が異なっていることに気づかれるかと思います。特に右下のフロリダ半島など四隅に近い方を見るとよく分かると思います。私が描画した下段はステレオ画法なのですが、上の地図がどういう投影方法なのか、色々検索してみましたが結局見つかりませんでした。
アメリカの国土は広く、州の大きさも大小さまざまなので、単一の投影方式では目的にそぐわない場合もあることでしょう。そのため、地図の部分によって投影方法が異なるような地図の方式もあったりしてなかなか奥の深い世界です。

さて、セマンティックセグメンテーションのマスク画像としては、下記のように地図の装飾を行わずに書き出します。

左：寒冷前線　右：温暖前線

左：停滞前線　右：閉塞前線

トラフ

左：高気圧記号　右：低気圧記号

3. まとめ

　アメリカの気象機関であるNOAA/NWSが発行している天気図（地表解析）について、Coded Surface Bulletinを用いた天気図要素の抜き出しを行いました。

Coded Surface Bulletinは、MetPyが提供しているライブラリを使用することで、非常に簡単に復号して描画することができます。

私は日本の気象データ可視化画像をもとに前線を自動検出するといったことを行なってきました。
気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（5）
アメリカのデータでも、NWSのような天気図を自動描画できるか、実施中です。