気象データアナリスト育成講座・試験講座

PythonによるGPVのハンドリング

PythonでGRIB2形式のGPVデータをデコードし、気象要素の変換計算などを経て、欲しいデータを取得する方法を解説します。
やり方はいくつかありますが、ここではPythonモジュールwxbcgribwxparamsを使う方法を紹介します。
またwxbcgrib内部で使われている、wgrib2ソフトウェアとnetCDF4モジュールについても紹介します。


目次

  1. Python上でwgrib2を実行する
  2. PythonモジュールnetCDF4を使ってGPVを読み込む
  3. Pythonモジュールwxbcgribを使ってGPVを読み込む
  4. Pythonモジュールwxparamsを使って気象要素を計算する

1. Python上でwgrib2を実行する

MacのターミナルやWindowsのコマンドプロンプトで実行したコマンドをPythonで実行するには、subprocessモジュールを使います。
ディレクトリやファイルの操作をするのに便利なosモジュールとともに、まずはモジュールの読み込みをしましょう。

以下のセルを実行して下さい。

subprocessを使ってコマンドラインのコマンドを実行するには、以下のようにコードを書きます。

まずwgrib2のバージョン情報を出力してみましょう。
コマンドラインで実行した時と同様に、wgrib2のバージョン情報が出力されればOKです。

では次にweatherフォルダにある、数値予報モデルMSMのデータ(2020/09/05 12zベース、FT=0〜15)の概要(目録)を表示させてみましょう。

ダブルクォーテーションの中にコマンドを記述します。コマンドラインで実行した場合と同じ出力が得られればOKです。

それでは次は、wgrib2matchオプションを使って183番目の要素を、Vオプションを使ってその詳細情報を出力してみましょう。
オプションが増えるとコマンドが長くなりますので、下記のように適宜変数を活用するとコードが見やすくなります。
ここでf文字列と呼ばれる機能を使い、変数からコマンドの文字列を構成しています。

だいぶ慣れてきましたでしょうか?

ではいよいよGRIB2形式のMSMデータを、NetCDF形式に変換しましょう。
wgrib2netcdfオプションを、以下のように与えて実行します。

wgrib2 <GRIB2ファイル> -netcdf <出力するNetCDFファイル名>

ここではGRIB2ファイル名.ncという命名則でNetCDFファイルを作るようにします。
変数を活用して、ダブルクォーテーションの中にコマンドを記述し、セルを実行して下さい。

このipynbファイルがあるフォルダを開き、NetCDFファイルが作成されていることを確認して下さい。

ここでファイルサイズを見て下さい。元々のGRIB2が約70MB、NetCDFが約180MBほどではないでしょうか。
一般にはGRIB2の方が圧縮性が高いため、ファイルサイズも小さくなっています。

以上で、Python上でwgrib2を実行し、GRIB2のGPVデータをNetCDFへ変換することができました。
目次へ戻る


2. PythonモジュールnetCDF4を使ってGPVを読み込む

NetCDF形式のファイルの読み込みには、その名もnetCDF4モジュールを利用します。
また読み込んだデータを扱う際にはnumpyも利用します。

まずはモジュールのimportから行います。以下のセルを実行して下さい。

netCDF4Datasetクラスを使ってNetCDFファイルを読み込みます。
以下のセルを実行するだけで、読み込みは完了です。

変数dsには、Datasetクラスにより作成されたオブジェクトが入っています。

NetCDFファイルから読み込んだデータの要素を見てみましょう。以下のセルを実行して下さい。

緯度経度、予想対象時刻、各種気象要素が入っているようです。
気象要素は『要素名_高度情報』という命名則になっています。

ではこの中から気温を選択してみましょう。以下のセルを実行して下さい。

何やら気温のデータに関するヘッダー情報のようなものが表示されました。
これらの情報は、以下のようにしてもう少し細かく分解できます。

気温の値そのものには、以下のようにしてアクセスできます。

見るとmasked_arrayという特殊なフォーマットで値が格納されていることがわかります。

以下のようにすると、通常のnumpy.ndarray型で値を取り出すことができます。

気温の値は、そのオーダーからも絶対温度[K]であることがわかります。

またデータサイズは3次元になっており、(FT, Latitude, Longitude)の数を表します。

これで気象要素を代表して気温の値を取り出すことができました。

では次は地上気圧のデータを取得してみましょう。

以上で、netCDF4モジュールを使ってNetCDF形式のGPVデータを読み込むことができました。
目次へ戻る


3. Pythonモジュールwxbcgribを使ってGPVを読み込む

第1章・第2章ではwgrib2を使ってGRIB2形式のGPVをNetCDFに変換し、netCDF4モジュールを使ってそれを読み込んで処理する、という流れを説明してきましたが、これらの処理を一括で行えるPythonモジュールがあったら便利ですね。

そこで気象ビジネス推進コンソーシアムでは独自にPythonモジュールwxbcgribを作成し、それらの機能をまとめて実装しました。
wxbcgribを使えば、上でやってきた処理を1行で完了させることができます。

それではまずwxbcgribのインポートから行います。
wxbcgribはこのipynbファイルと同じフォルダにwxbcgrib.pyを置くことでインポートすることができます。
このipynbファイルと同じフォルダにwxbcgrib.pyがあることを確認して、以下のセルを実行して下さい。

それでは早速、GRIB2形式のMSM-GPVファイルを、直接wxbcgribを使って読み込んでみましょう。
読み込みには、wxbcgribDS_from_grnc関数を使います。

wg.DS_from_grnc(GRIB2ファイルへのパス, 取り出したい気象要素名)

それでは以下のセルを実行して下さい。

これで気温の値を読み込むことができました。簡単ですね!

なお今回は気温の要素名がTMP_1D5mabovegroundと知っていたので指定できますが、知らなかったらどうしましょうか?
その時は少々強引ですが、要素名を指定せずに実行すると、エラーメッセージとして有効な要素名を表示してくれます。
wxbcgrib.py製作者のドキュメントより)

さて変数varには、netCDF4.Datasetクラスと同じで、オブジェクトが入っています。

まず要素名と単位を表示してみましょう。

次に予想対象時刻(validtime)を見てみましょう。

validtime、はdatetime型という日時を表す特別なフォーマットで格納されています。
例えば、一番上はFT=0で、2020/09/05 21:00(JST)であることを表します。

気温のデータそのものは、以下のようにしてアクセスできます。

データサイズは3次元になっており、(FT, Latitude, Longitude)の数を表します。

緯度経度情報もあります。以下のセルを実行して下さい。

wxbcgribでは、緯度経度を指定して時系列データを取り出すことができます。
tsメソッドを使って、たった1行のコードでできてしまいます。引数として緯度・経度の順番で与えます。

変数名.ts(緯度, 経度)

今回は北緯35度・東経140度の値を取り出してみましょう。以下のセルを実行して下さい。

こんなに簡単にできてしまうなんて、便利ですね!

これで、第1章・第2章とで行った様々な処理を、wxbcgribで簡単に実行することができました。

それでは次は、北緯26.2度・東経125.75度の地上気圧の時系列データを、wxbcgribで取得してみましょう。

以上で、wxbcgribを使ってGRIB2形式のGPVデータをハンドリングすることができました。
目次へ戻る


4. Pythonモジュールwxparamsを使って気象要素を計算する

第1章〜第3章でGPVデータをPythonで読み込むことができました。
すると次に必要になる処理として、読み込んだ気象要素の単位変換や、別の気象要素への変換があります。

これらの計算を簡単に行う方法として、Weather Data Science 加藤芳樹・史葉のオリジナルPythonモジュールwxparamsを紹介します。
wxparamsは、気象要素の変換計算や単位変換の計算を行う関数群から成り立っています。

ではwxparamsと、解説で使うnumpypandasをインポートします。以下のセルを実行して下さい。

まずはwxparamsを簡単に使ってみます。

気温20度・相対湿度75%の時の露点温度を計算してみましょう。
wxparamsはインプットとしてnumpy.ndarray型を想定していますので、Numpy配列の形にして引数を渡します。
また露点温度の計算にはRH_to_Td関数を使います。

以下のセルを実行して下さい。

戻り値もnumpy.ndarray型になっています。

またpandas.Series型でもOKです。まず以下のセルを実行して、データフレームを作ります。

引数としてTEMP列とRH列を、以下のようにRH_to_Td関数へ渡すこともできます。

簡単に概要を説明しましたので、練習用データフレームでwxparamsを使ってみましょう。
まず以下のセルで練習用データを定義します。

データの内容は、気温[K], 相対湿度[%], 気圧[Pa]です。

では露点温度を計算してみましょう。注意点として、RH_to_Td関数へは気温を摂氏温度[C]で渡す必要があります。

では次にTheta_e関数で相当温位を計算してみましょう。各引数の単位に注意して下さい。


ここまでは、テーブルデータ(pandas.DataFrame)の形で気象データを持っている場合を練習してきました。
次は、GPVから取り出したデータを使い、GPVには含まれていない気象要素を計算してみましょう。

ここでは計算する気象要素として『風向・風速』を扱ってみます。
第3章で使ったGPVデータをwxbcgribを使って読み込み、風のUV成分(東西風・南北風)を入手しましょう。

読み込んだデータのデータサイズ(shape)と、単位を確認しましょう。

読み込んだデータは3次元のnumpy.ndarray型ですが、wxparamsは多次元のNumpy配列を想定しています。
よって、多次元配列のまま引数として与えることができます。

それではUV_to_SpdDir関数を用いて風向・風速を求めてみましょう。
注意点として、戻り値は2つあり、風速・風向の順番でreturnされてきます。

また、風速のデータサイズ(shape)と、FT=0における値を出力してみましょう。

以上でwxparamsを使って気象要素の変換計算と、単位変換の計算をすることができました。
目次へ戻る