GSMの500hPa面天気図作成 1

この連載は、予報士など気象に関する知識があり、様々な天気図を作成したい人のために、python上で動作するMetPyというツールを紹介しつつ、日々利用したい天気図について考えていきます。

これからの天気がどうなるかをざっくり知りたいとき、私は気象庁作成の数値予報天気図を確認します。

気象庁作成の500hPa面天気図

この数値予報天気図はよく考えられていて、通年を通して使える基本の天気図です。ただし、昔から変わっていません。白黒のため、すっと頭に入ってこないなど改良の余地があり、図をカラー化しただけでも見やすくなります。

ここ数十年で数値予報の精度も気象学も、とても進歩しています。これらを上手く活用して、新しい基盤となる数値予報天気図を作れないでしょうか？

総観気象学の教科書や論文で掲載されている天気図をお手本に、様々な数値予報天気図をMetPyを使って作成していきます。予想される気象擾乱や現象を素早く把握できて、その構造やステージなども確認できる数値予報天気図を作成していくことを目指します。

入力データは気象庁GSMのGRIB2形式を使います。この形式のデータは最新のものでなければ、インターネットからダウンロードできるサイトがいくつかあります。

この第1回ではMetPyのインストールから、GSMの予測データの読み込み、物理量の計算までを説明します。最初に作成する天気図は500hPa面天気図で、高度と渦度の分布図のカラー化を行います。

 なお、第3回にJupyter Notebook用のコードを公開しています。実際にコードを実行しながら試してみたい方は、先にこれをダウンロードした方がよいでしょう。

500hPa面の高度・渦度の天気図

500hPa面は対流圏のおよそ中間の高度にあたります。この高度や渦度を示した天気図から、天気に大きな影響を与える偏西風の蛇行やリッジ、トラフを解析できます。この解析結果と地上にある低気圧や高気圧の位置関係から、これら擾乱の盛衰を大まかに予想できることができるなど、500hPa面の天気図は地上天気図と同様に、とても基本的な天気図の一つです。

MetPyのインストールについて

MetPyを利用するには、Anaconda3をインストールするのがスムーズです。このインストールの詳細は、他の方々の記事に譲ります。

私は、met_envという環境を作成して、MetPyをインストールしています。参考のために、その際のコマンドを示します。

## conda のインストール
$ conda config --set auto_activate_base false
$ conda update anaconda
$ conda update conda
$ conda update --\all
#
## metpy, netCDFなどをインストールする、met_env環境の作成
$ conda create -n met_env python=3.8
$ conda info -e
$ conda activate met_env
$ conda install netCDF4 
$ conda install cartopy
$ conda install -c conda-forge metpy
$ conda deactivate

GRIB2ファイルの読み込み

MetPyを利用した天気図作成コードでは、GRIB2ファイルの読み込み、xarray というデータ形式などを利用するために、次のimportを行なっています。

import math
import pygrib
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import matplotlib.path as mpath
import cartopy.crs as ccrs
import datetime
import sys
#                                                                                                                                   
import metpy.calc as mpcalc
from metpy.units import units
from scipy.ndimage.filters import maximum_filter, minimum_filter
#                                                                                                                                   
import argparse

次に、初期時刻・予想時間・等圧面の指定と、データの読み込みのコードを示します。

ポイントは、下から4から6行目にあるデータ取得する際の、shortNameの指定でしょう。GSMの指定気圧面には5つの物理量が格納されていて、高度(gpm)のshortNameは"gh"、風速(西風成分 m/s)は"u"、風速(南風成分 m/s)は"v"、気温(K)は"t"、鉛直速度(pa/s)は"W"、相対湿度(%)は"r"となっています。ただし、相対湿度は地上から300hPa面までの高度にしか格納されていないことに、留意してください。

なお、このサンプルコードを実行する際は、フォルダー./data/gsm/ を作成し、その中に、ファイル
　　Z__C_RJTD_20210823000000_GSM_GPV_Rgl_FD0000_grib2.bin
をあらかじめダウンロードしてください。

# GSMの読み込む初期値の年月日時をUTCで与えます。
i_year =2021
i_month = 8
i_day = 23
i_hourZ = 0
#
# 予想時間を与える。この値は注意が必要です。下3桁目が日数、下２桁で時間で与えます。
#  初期値なら0、18時間後なら18、24時間後なら100、36時間後なら112となります。
i_ft = 0
#
# 読み込む気圧面の気圧を与えます。ここでは500hPaと指定します。
tagHp = 500
#
# データの格納先フォルダー名
data_fld="./data/gsm/"
####################################################
#
# 読み込むGRIB2形式GSMのファイル名
gsm_fn_t="Z__C_RJTD_{0:4d}{1:02d}{2:02d}{3:02d}0000_GSM_GPV_Rgl_FD{4:04d}_grib2.bin"
gr_fn= gsm_fn_t.format(i_year,i_month,i_day,i_hourZ,i_ft)
#
# データOpen
grbs = pygrib.open(data_fld + gr_fn)
#
# データ取得                                                                                                                       
grbHt = grbs(shortName="gh",typeOfLevel='isobaricInhPa',level=tagHp)[0]
grbWu = grbs(shortName="u",typeOfLevel='isobaricInhPa',level=tagHp)[0]
grbWv = grbs(shortName="v",typeOfLevel='isobaricInhPa',level=tagHp)[0]
#
# データClose
grbs.close()

読み込んだ格子点データは全球分あります。表示に必要な領域のみを切り出し保存することで、スクリプトの実行時間を短縮できます。データの切り出しのコードは次の通りです。

## GPVの切り出し領域の指定：(lonW,latS)-(lonE,latN)の矩形                                                                                                      
latS=-20
latN=80
lonW=70
lonE=190
## データ切り出し                                                                                                                   
valHt, latHt, lonHt = grbHt.data(lat1=latS,lat2=latN,lon1=lonW,lon2=lonE)
valWu, latWu, lonWu = grbWu.data(lat1=latS,lat2=latN,lon1=lonW,lon2=lonE)
valWv, latWv, lonWv = grbWv.data(lat1=latS,lat2=latN,lon1=lonW,lon2=lonE)

渦度などの物理量を計算するためのデータ変換

MetPyは、様々な物理量を計算する関数がそろっています。これら関数を利用するために、格子点データをxarrayというデータセットの形式に変換する必要があります。この変換の際に、データの単位を指定することも忘れないでください。単位を指定することで、例えばm/sからktへの変換も簡単に行うことができ、関数の引数に誤った単位の変数を指定するとアラートが出るなど便利です。

## 渦度などの算出のためにxarrayデータセットを作成                                                                                   
ds = xr.Dataset(
   {
       "Geopotential_height": (["lat", "lon"], valHt),                                                                        
       "u_wind": (["lat", "lon"], valWu),
       "v_wind": (["lat", "lon"], valWv),
   },
   coords={
       "level": [tagHp],
       "lat": latHt[:,0],
       "lon": lonHt[0,:],
       "time": [grbHt.validDate],
   },
)
# 単位も入力する
ds['Geopotential_height'].attrs['units'] = 'm'
ds['u_wind'].attrs['units']='m/s'
ds['v_wind'].attrs['units']='m/s'
ds['level'].attrs['units'] = 'hPa'
ds['lat'].attrs['units'] = 'degrees_north'
ds['lon'].attrs['units'] = 'degrees_east'

渦度の計算

相対渦度の算出コードを示します。とてもシンプルです。xarray形式の恩恵を実感できます。

## 算出                                                                                                                             
# 相対渦度                                                                                                                          
ds['vorticity'] = mpcalc.vorticity(ds['u_wind'],ds['v_wind'])

次回は作図部分のコードをご紹介します。