気象庁レーダー毎極座標GPVのGRIBをnetCDF(CF/Radial 規約)に変換する

TL;DR

変換プログラムをGitHubに置いておくので自己責任でどうぞ。
二重偏波気象レーダー極座標データの変換についてはこちらの記事を参照ください。

はじめに

降水現象を調べるときに、気象庁(JMA)のレーダーデータを使いたいことがあります。通常は、1kmメッシュ全国合成レーダーエコー強度GPVや解析雨量）などの複数のレーダーからの情報を加工して作られた降水量の空間2次元データを使うことが多いでしょう。しかし降水の3次元の空間分布を調べるといった凝ったことをする場合、レーダーを中心とする仰角別の極座標格子に反射強度が格納されたレーダー毎極座標データを扱う必要があります。

このデータは他の気象庁のデータと同じくGRIB2形式で格納されており、テンプレート3.50120や4.51022といったJMA独自のテンプレートが使われていることが特徴としてあげられます（配信資料に関する仕様 No.13702を参照）。そのため、NCLやGrADSといった描画ソフトでそのままの形で読み出して描画することはできません。

描画する一つのやり方として、NICTの気象庁極座標レーダーデータのダウンロードサイトにあるサンプルプログラムを利用する方法があります。このサンプルプログラムには、wgrib2でJMA独自のテンプレートを読めるようにするパッチ、パッチを当てたwgrib2を使って変換したプレーンバイナリからレーダーサイトを中心とした直交座標に変換するCのプログラム、変換後のデータを描画するGrADSスクリプトが含まれています。

一方で、最近ではPyARTやwradlibといったPythonでのレーダーデータの解析ツールが充実してきています。これらのツールは、当然ながら独自テンプレートのGRIBを読んでくれるはずはなく、CF/Radial規約のNetCDF形式などにデータを変換して読み込ませる必要があります。このような変換プログラムを公開している人はいないようなので、作成しました。

CF/Radial規約とは何か

CF/Radial規約の説明をする前に、CF規約について知っておく必要があります。CF規約については、詳細は豊田様の netCDFとCF規約について の記事が参考になります。ここではNetCDF CF 規約 日本語訳のサイトの説明が簡潔なので引用させていただきます。

NetCDF は（ちょうど XML のように）データ形式の枠組みを決めるものです。 データセットには属性と呼ばれる名前-値ペアを与えることができますが、 利用者の合意にしたがってどのような情報でも書き込めるので、 どのように使うのかの約束が必要で、それを NetCDF 規約といいます。 その規約の中でデファクトスタンダード化しつつある CF 規約を日本語訳したものです。

CF/Radial規約は正式名称が"CF-compliant netCDF Format for Moments Data for RADAR and LIDAR in Radial Coordinates"であることから分かるように、CF規約では対応が不十分であった、極座標で表現されるようなレーダーやライダーの観測データを表現する規約です。

CF/Radial規約にはバージョン1.x系(CfRadial1)と2.x系(CfRadial2)の2種類が存在します。CfRadial1はNetCDFのclassic data modelに基づくものです。2011年に最初のバージョン1.1が策定され、2016年に策定されたバージョン1.4が最新版になります。一方のCfRadial2は、NetCDF4とgroupを用いるもので、CfRadial1との後方互換性はありません。現在の最新版のバージョン2.0は2017年に策定され、2019年に一部更新が加えられています。またドラフトの段階であるバージョン2.1が2019年に公開されています。

CfRadial2は最近出来たため、ツールによって対応状況が異なるようです。PyARTに関しては、Issue Reading CF/Radial 2.0 Files #858を読む限り、公式版ではCfRadial2に一部問題があるようで、有志で互換性のあるリーダーが作成されている模様です。wradlibについては、wradlib.io.xarray.CfRadialを見る限り、CfRadial2に対応しているようです。（いずれも実際に確認したわけではないです）

変換プログラムの作成

基本的に配信資料に関する仕様 No.13702を読んで変換するコードを書くだけです。今回は利用できる環境が多く構造が簡単なことから、CF/Radial Version 1.4に変換することにしました。NetCDFの作成には、PythonのnetCDF4パッケージを使いました。いくつかポイントを記しておきます。

• GRIBの第1節（識別節）・第3節（格子系定義節）・第4節（プロダクト定義節）はPythonでバイナリを読み出しました。第5節（資料表現節）と第7節（資料節）のランレングス圧縮されたデータについては、Pythonで解凍するコードを書くのが面倒なので、wgrib2でプレーンバイナリに変換してそれを読み込んでいます。プレーンバイナリへの出力はJMA独自テンプレートは関係しないため、通常のwgrib2でもエラーになりません。もちろん外部プログラムへの依存性が気になる人は、ランレングスの解凍をPythonで書くこともできると思います。

• CfRadial Data File Formatの2.3節を読めば分かる通り、データの格納方法には Regular 2-D storage と Staggered 2-D storage の2種類があります。前者は時間によらずレンジ方向のビン数が一定な場合で、反射強度などの変数はrangeとtimeの2次元配列として格納されます。後者はレンジ方向のビン数が時間によって可変な場合で、反射強度などの変数は1次元配列して格納され、各時刻のレンジ方向のビン数と位置についてはray_n_gatesとray_start_indexに格納されます。レーダー毎極座標データの場合は仰角によってレンジ方向のビン数が可変なため、はじめは Staggered 2-D storage を使おうと思いました。しかし出力したnetCDFがPyARTで上手く読み込めなかったので、Regular 2-D storage にして値がないところは_FillValueで埋めました。

• 2020年から気象庁のレーダーが順次更新されていて、執筆時点では東京レーダーが更新されました。この更新の目玉は二重偏波化ですが、隠れた改善点としてボリュームスキャン（上から下まで全ての仰角で3次元的にスキャンすること）にかかる時間が10分から5分になることが挙げられます。GRIB2のデータは引き続き10分で1つのため、1ファイルに2つのボリュームスキャンが含まれることになります。CF/Radial規約を読むと1ファイルに1ボリュームのようなので（本当かは自信がない）、1つのGRIBから2つのnetCDFを生成したい場合がありそうです。そこで、GRIBの格納されているデータの番号（何番目の仰角かに対応）を指定して切り出せるようにしました。

作成した変換プログラムはGitHubにMITライセンスで置いておきます。

変換したデータの描画

本題から外れますが、出力されたデータをPyARTで図を描くとこんな感じになります。実質4行で最低限の図が描けるのは便利ですね。

draw.py

import matplotlib.pyplot as plt
import pyart

radar = pyart.io.read_cfradial('./work/RS47695.nc')
display = pyart.graph.RadarDisplay(radar)
display.plot('DBZ', sweep=29)
plt.show()

謝辞

• 山崎一哉さん（東京大学理学系研究科）には、PythonでのGRIBの読み出し方法を参考にさせて頂きました。
• 松岸修平さん（東京大学大気海洋研究所）には、Pythonのレーダー解析のパッケージの情報を教えて頂きました。
• 気象庁レーダー毎極座標GPVは気象業務支援センター経由で取得しました。