GRIB2をPythonで読み込んで気象データを地図に描く方法

はじめに
使用ライブラリ
サンプルデータ
GRIB2を触ってみる
- まずは読み込む
  - ①データを開く
  - ②データを取得する
- 地図上にプロットしてみる
参考

はじめに

気象データを利用しようとするとテキストからバイナリまでさまざまなデータ形式に遭遇しますが、なかでも数値計算の出力に使われるGRIB2は気象分野独特のデータ形式の一つだと思います。

GRIB2はWMOが策定したグリッドデータの標準フォーマットです。形式としてはバイナリに分類されるため、人間が中身を読むには専用のツールやライブラリを使ってデータを抽出・変換する必要があります。代表的な変換ツールとしてはNOAAのwgrib2などが挙げられますが、いつのまにかPythonでも簡単に処理できるようになっていて便利だったので、備忘も兼ねて記事にすることにしました。

使用ライブラリ

pygrib
GRIB2を読み出す（デコードする）ライブラリです。 ECMWFのGRIBデータ用CライブラリをPythonで扱えるようにしたものです。
Cartopy
地図への描画ライブラリで、Matplotlibの拡張ライブラリBasemapの後継版にあたります。もともとの開発元は英国気象局です。
NumPy、Matplotlib
説明がいらないくらいどこででもよく使われるライブラリたちです。

サンプルデータ

GSMデータ
気象庁の全球数値予報モデルで計算されたグリッドデータで、日々GRIB2形式で配信されています。
今回は日本域の地表面データを使用しました（GSMデータ自体の詳細については気象庁HPの各種技術資料をご参照ください）。

GRIB2を触ってみる

まずは読み込む

①データを開く

GRIB2データを開きます。read()でデータの情報をリストで確認できます。

import pygrib

grbs = pygrib.open("Z__C_RJTD_YYYYMMDDhh0000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0000-0100_grib2.bin")
grbs.read()

②データを取得する

必要な物理パラメータを指定し、データを読み込みます。
pygribのdata()を用いることで、data, lats, lonsに全グリッドの物理パラメータと緯度経度をそれぞれ取得できます。

#t=0のある物理パラメータを読み込む
grb = grbs.select(parameterName='Low cloud cover')[0]
grb.data()
data, lats, lons = grb.data()

これであとは二次元配列データとして扱えます。

ちなみにパラメータ名を指定するのに使える引数はいくつかあるのですが（e.g. nameやshortName）、GSMはなぜか一部に「unknown」が入っていたりして処理に困ったので、parameterNameを使うのがよさそうでした。

各パラメータがどんな名前で登録されているかはリストで書き出したりすると確認できると思います。

check = list()
for ck in grbs:
    check.append(ck.parameterName)
print(check)

地図上にプロットしてみる

②で取得したデータを、Cartopyを使って試しに作図してみます。
操作は基本的にMatplotlibと同じです。サブプロットで呼び出して使います。

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

parameter_name = 'Low cloud cover'

#作図領域の指定とCartopyの呼び出し
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})

#データを描画
contour = ax.contourf(lons, lats, data, levels=np.linspace(0, 100, 11), cmap='Blues_r')
plt.colorbar(contour, ax=ax, label=f'{parameter_name}', orientation='vertical')

#海岸線を描画
ax.coastlines()

#軸関係を諸々指定
gl = ax.gridlines(draw_labels=True, linestyle='-')
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.top_labels = False
gl.right_labels = False

ax.set_title(f'{parameter_name}')
plt.show()