地上天気図解析とプロット図作成コード
地上天気図解析の概要と、この解析のために地上観測をプロットした地図を作成するPythonのコードについて解説します。
この最後に、実行可能なサンプルコードを添付しています。説明用に示すコードは、わかりやすさ重視でそのままでは動作しません。ご留意ください。
地上天気図解析について
地上天気図には、気象台や空港で観測された気圧や気温、風、雲量などの情報がプロットされています。これら観測データから等圧線や等温度線を解析し、雲量や天気、風の場なども把握して、低気圧や高気圧、前線、気圧の谷などを解析します。天気図解析を３ないし６時間ごとに行うことで、天気の変化や前線などの擾乱の動きなどを把握することができ、知識や経験があれば大きな場における天気変化を予想することも可能でしょう。
解析する際に知っておくべき知識や留意点があります。参考資料として、気象庁の量的予報技術資料（予報技術研修テキスト）があります。「現業作業における総観場の把握と局地気象解析について」や「アジア太平洋地上天気図の標準的な解析手法」などが参考になるでしょう。
特に気をつけてもらいたいことは、観測値に忠実に等圧線や等温度線を解析すると不適切な解析となる場合があります。観測値に忠実に解析することより、解析する天気図の範囲に見合ったスケールの気象擾乱を適切に解析することが重要です。観測値には、様々なスケールの現象が重なった結果の値であり、その位置や高度に伴う局地性や、観測誤差も一定程度含まれているためです。近隣の複数の観測地点のデータとその変化も確認しつつ、前線などの気象擾乱の時空間的な連続性を考慮して解析してください。
地上天気図には、地上気象観測の結果がプロットされています。このプロットの形式には国際式と日本式があります。ここでは、国際式のプロットした図を作成していきます。国際式の天気図の説明は、この気象庁の解説資料をご確認ください。
2021年10月11日9時、東北地方周辺の地上天気図の解析例を示します（下の図）。等圧線を黒実線、等温度線を赤破線で解析しました。寒冷前線が東北地方を通過中です。前線の北側では等温度線が集中し、温度変化が大きく、その南端に寒冷前線を解析しています。前線に伴う気圧の谷も明瞭で、気圧の傾きも比較的大きくなっていて、北日本の太平洋側の海上を中心に風が強くなっていることが想定できます。
地上気象観測のデータ読み込み
気象台などで観測された地上気象観測(SYNOP)と空港で行われている航空気象観測(METAR)、両方の観測データを読み込むコードを下に示します。unidataのサイトにあるデータセット(最新1ヶ月分のみ)をダウンロードし、SYNOPはnetCDF形式、MEATRはtext形式の観測データを読み込んでいきます。METARのデコードは、MetPyの関数metar.parse_metar_file()を利用しています。
from datetime import datetime,timedelta
from siphon.catalog import TDSCatalog
import netCDF4
from metpy.io import metar
### データ読み込み
# 表示する時刻(UTC) 読み込むデータの日時(最近1ヶ月しかない)
tg_date = datetime(year, month, day, hour, 0, 0)
## Metar & SYNOPデータの読み込みカタログ
# metar:text形式
metar_catalog_url = 'https://thredds-test.unidata.ucar.edu/thredds/catalog/noaaport/text/metar/catalog.xml'
# synop:netCDF形式
synop_catalog_url = 'https://thredds.ucar.edu/thredds/catalog/nws/synoptic/ncdecoded/files/catalog.xml'
#
## カタログ読み込み
s_cat = TDSCatalog(synop_catalog_url)
m_cat = TDSCatalog(metar_catalog_url)
# 読み込むデータの日時を指定し、urlなど取得
s_ds = s_cat.datasets.filter_time_nearest(tg_date)
m_ds = m_cat.datasets.filter_time_nearest(tg_date)
# データのダウンロード
s_ds.download()
m_ds.download()
# SYNOP netCDF読み込み
nc=netCDF4.Dataset(s_ds.name,'r')
#
### METAR Parse
m_df = metar.parse_metar_file(m_ds.name)MetPyの描画関数を利用するために、次の手順で、netCDF形式のSYNOPのデータをpandasのデータセット形式に変換します。
1. netCDFから必要なデータを配列として取り出す
2. データ量を減らすために、必要な地域のデータを抽出
3. pandas データセットの作成
下のコードは1にあたり、データ種類ごとに配列を作成しています。
### Synop pandasデータの作成
# netCDFからpandasへ変換するために、要素毎の配列作成
# 1) WMO Identification number
wmoId = nc.variables['wmoId'][:]
station_id = [str(n) for n in wmoId]
# 2) Station Name
stnName = nc.variables['stnName'][:]
# 3) time of Observation seconds since 1970-01-01 00 UTC
# => Timestamp('2018-09-21 03:04:50')
time_obs = [ pd.Timestamp(n,unit='s') for n in nc.variables['time_obs'][:]]
# 4) time nominal seconds since 1970-01-01 00 UTC
# => Timestamp('2018-09-21 03:04:50')
time_nominal = [pd.Timestamp(n,unit='s') for n in nc.variables['time_nominal'][:]]
# 5) Station
latitude = nc.variables['Lat'][:]
longitude=nc.variables['Lon'][:]
elevation=nc.variables['elev'][:] # meters
stnType = nc.variables['stnType'][:] # Table 1860
# 6) 視程
vis=nc.variables['VIS'][:] # Table 4377
# 7) 風速
wind_direction=nc.variables['DIR'][:] # degree_true 360
wind_speed=nc.variables['SPD'][:] # m/s
# 8) 気温、露点温度
air_temperature=nc.variables['T'][:]
dew_point_temperature=nc.variables['TD'][:]
# 9) 気圧 hectopascals
air_pressure_at_sea_level=nc.variables['SLP'][:]
air_pressure_at_station_level=nc.variables['PRES'][:]
# Change in Pressure in past 3 hours
ptend=nc.variables['Ptend'][:]
# Character Pressure tendency(0|1|2|3|4|5|6|7|8 Table 12-5)
cptend=nc.variables['char_Ptend'][:]
# 10) precipitation
# 雨量(mm) Table 3590
precip_amt=nc.variables['PRECIP_amt'][:]
# 雨量の期間 (hours) Table 4019
precip_period=nc.variables['PRECIP_period'][:]
# 11) weather
# Manned stn Table 4677, autoTable 4680
present_weather=nc.variables['WXpresent'][:]
# Manned stn Table4561, autoTable4531
past_weather=nc.variables['WXpast'][:]
# 12) cloud
# Total cloud cover in oktas Table 2700
cloud_coverage=nc.variables['cloudCover'][:]
# Lower level clouds Table 0513
low_cloud_type=nc.variables['cloudLow'][:]
# Middle level clouds Table 0515
medium_cloud_type=nc.variables['cloudMiddle'][:]
# Higher level clouds Table 0509
highest_cloud_type=nc.variables['cloudHigh'][:]下のコードではデータ量を減らすために、緯度・経度の範囲を指定し、その範囲内の観測地点を抽出します。
## データカット領域の指定
latMax = 90.0
latMin = 0.0
lonMax = 180.0
lonMin = 90.0
## 配列切り出し用の関数
def pickup_array(ary, ok_a):
tary=[]
for i in ok_a:
tary.append(ary[i])
return tary
## 配列のサイズを小さくするために、エリア内にあるデータに絞る
ok_a=[]
for i in range(len(station_id)):
lat = latitude[i]
lon = longitude[i]
if (lat < latMax and lat > latMin and lon > lonMin and
lon < lonMax and tg_date == time_obs[i]):
ok_a.append(i)
#
station_id = pickup_array(station_id,ok_a)
latitude = pickup_array(latitude,ok_a)
longitude = pickup_array(longitude,ok_a)
elevation = pickup_array(elevation,ok_a)
time_obs = pickup_array(time_obs,ok_a)
air_temperature = pickup_array(air_temperature,ok_a)
dew_point_temperature = pickup_array(dew_point_temperature,ok_a)
cloud_coverage = pickup_array(cloud_coverage,ok_a)
air_pressure_at_sea_level = pickup_array(air_pressure_at_sea_level,ok_a)
present_weather= pickup_array(present_weather,ok_a)
wind_direction = pickup_array(wind_direction,ok_a)
wind_speed = pickup_array(wind_speed,ok_a)
ptend = pickup_array(ptend,ok_a)
cptend = pickup_array(cptend,ok_a)
vis=pickup_array(vis,ok_a)
low_cloud_type = pickup_array(low_cloud_type,ok_a)
medium_cloud_type = pickup_array(medium_cloud_type,ok_a)
highest_cloud_type = pickup_array(highest_cloud_type,ok_a)次に、pandas Datasetを作成します。ただし、利用する描画関数の仕様上、風のデータを東成分と北成分として与える必要があるため変換しています。
# 風の西・北風成分へ変換
# Maskedデータを描画させるとエラーが出るので、とりあえず風速0とする
eastward_wind = []
northward_wind = []
for i in range(len(wind_direction)):
if (type(wind_direction[i]) != np.ma.core.MaskedConstant and
type(wind_speed[i]) != np.ma.core.MaskedConstant):
dr = wind_direction[i]
sp = wind_speed[i] * 1.94384 # m/s => kt
ew = -1.0 * sp * math.sin(dr * math.pi / 180.0)
nw = -1.0 * sp * math.cos(dr * math.pi / 180.0)
eastward_wind.append(ew)
northward_wind.append(nw)
else:
eastward_wind.append(0.0)
northward_wind.append(0.0)
#
## SYNOP pandas DataFrameの作成
dic_arr = {'station_id':station_id,
'latitude': latitude, 'longitude': longitude,
'elevation': elevation, 'date_time': time_obs,
'air_temperature':air_temperature,
'dew_point_temperature':dew_point_temperature,
'cloud_coverage':cloud_coverage,
'air_pressure_at_sea_level':air_pressure_at_sea_level,
'eastward_wind':eastward_wind,'northward_wind':northward_wind,
'present_weather':present_weather, 'ptend':ptend,
'cptend':cptend, 'vis':vis,
'low_cloud_type':low_cloud_type,'medium_cloud_type':medium_cloud_type,
'highest_cloud_type':highest_cloud_type}
data = pd.core.frame.DataFrame(dic_arr, index=station_id)観測データプロットの描画指定
PlotObsクラスを利用して、METARとSYNOP別々に、これまでで読み込んだ観測データを入力し、観測プロットではどのデータをどの位置に表示するかを指定する必要があります。位置を示すlocationsの指定は'NW', 'SW', 'C' などとしています。この'NW'は、観測地点の位置からみてなら北西側（左上）に表示するという意味です。'C'を指定すると、観測地点の位置を中心に表示されます。下のコードでは気温を観測ポイントの北西側に赤色で表示するよう指定しています。また、reduce_pointsに0.0より大きな値を入れると、近接した観測地点が表示対象外となり、間引かれて表示されることになります。
## METAR：Plotするための準備 PlotObsインスタンス作成
m_obs = PlotObs()
m_obs.data = m_df
m_obs.time = tg_date
m_obs.level = None
#
# plotするデータの指定、表示位置、文字の色、Format、矢羽
m_obs.fields = ['air_temperature', 'dew_point_temperature',
'cloud_coverage','air_pressure_at_sea_level']
m_obs.locations = ['NW', 'SW', 'C','NE']
m_obs.colors = ['tab:red', 'tab:green', 'black', 'blue']
m_obs.formats = [None, None, 'sky_cover',lambda v: format(10 * v, '.0f')[-3:], ]
m_obs.vector_field = ['eastward_wind', 'northward_wind']
# 値を大きくすると描画するデータが間引かれます
m_obs.reduce_points = 0.0
## SYNOP：Plotするための準備 PlotObsインスタンス作成
obs = PlotObs()
obs.data = data
obs.time = tg_date
obs.level = None
#
# plotするデータの指定、表示位置、文字の色、Format、矢羽
obs.fields = ['air_temperature', 'dew_point_temperature', 'cloud_coverage',
'air_pressure_at_sea_level',
'present_weather','ptend','cptend','vis',
'low_cloud_type', 'medium_cloud_type', 'highest_cloud_type']
obs.locations = ['NW', 'SW', 'C', [1.8,1.0], [-1.3,0], [1.5,0], [2.5,0],
[-2.3,0], [0, -1.0], [0, 1.0],[0, 2.0]]
obs.colors = ['tab:red', 'tab:green', 'black', 'blue','blue','blue','black',
'black','blue','black','black']
obs.formats = [None, None, 'sky_cover', lambda v: format(10 * v, '.0f')[-3:],
'current_weather',
lambda v: format(10 * v, '.0f')[-2:], 'pressure_tendency',None,
'low_clouds', 'mid_clouds', 'high_clouds']
obs.vector_field = ['eastward_wind', 'northward_wind']
# 値を大きくすると描画するデータが間引かれます
obs.reduce_points = 0.0地上天気図解析のためのプロット
ここまでで描画準備は完了しています。後は、表示する地図の図法や、図の大きさ、描画する海岸線、タイトルなどを指定すれば、プロット図を作成できます。panel.plotsで、SYNOPとMETARの両方を指定することで、両方の観測データをプロットできます。
## 図のsize指定 インチ
fig_size = (13,13)
## 地図の表示範囲指定
fig_area = [136,145,35,43] # 東北付近
## MapPanelクラスを利用し、地図のエリア、投影、描画する地図要素を指定
panel = MapPanel()
panel.layout = (1, 1, 1)
panel.area = fig_area
#
## 投影する図法
crs = ccrs.Mercator(central_longitude=140.0)
panel.projection = crs
#
## 表示指定
panel.layers = ['coastline', 'borders', 'states']
panel.plots = [m_obs, obs] # METAR & SYNOP
#
## タイトル
panel.title = 'SYNOP & METAR PLOT ' + obs.time.strftime("%Y/%m/%d %H:%M:%SZ")
#
## 一般的なパネルクラスを利用して、図の大きさ、表示するパネルを指定
pc = PanelContainer()
pc.size = fig_size
pc.panels = [panel]
#
pc.show()サンプルコード
下にサンプルコードを掲載します。このコードを利用し地上観測プロット図を作成し、気象庁ホームページにある高層天気図を参照しながら、地上天気図解析してみてください。残念なことに、紹介したサイトからではSYNOPは毎日日本時間９時のデータしかありません。METATは毎時ある様です。ご留意ください。
なお、この記事では、MetPy Vesionは1.1を利用しています。