GPVデータファイルからデータを抽出する
この記事は、SDGsの目標１３「気候変動に具体的な対策を」をテーマに、気象データの利用方法の普及と気象関連ソリューション開発を促進するためのナレッジを提供することを目的として作成しています。
今回は、Pythonを使って気象庁の数値予報結果をデータ処理する基本的な方法について解説していきます。
目次
数値予報モデルとは
数値予報モデルは、大気や海洋の物理的なプロセスを流体力学方程式や熱力学方程式で表現し、これらの方程式をスーパーコンピュータで数値的に解くことで、未来の気象や気候の変化を予測するための気象予報ツールであり、世界中の気象庁や気象機関で様々なモデルが運用されています。
風、気温、湿度、気圧などの物理パラメーターを時間とともに予測することで、日々の天気予報や気候変動の解析を行うことができます。
気象庁は、全球モデル（GSM）、メソモデル（MSM）、および局地モデル（LFM）を一日数回運用して、日々の天気予報の基礎となる数値予報資料を提供しています。また、1か月予報、3か月予報や寒候期／暖候期予報を作成するため、アンサンブル予報モデルも運用しています。
これらのモデルには、地上観測、ラジオゾンデや航空機による上層観測、気象衛星ひまわりなどの衛星観測の結果などが入力されており、これにより予報精度が飛躍的に向上しています。
https://www.jma.go.jp/jma/kishou/know/whitep/1-3-1.html
数値予報モデルは、地表面と大気を格子状で表現し、それぞれのメッシュごとに風、気温、湿度、気圧などの物理パラメーターを計算していきます。
これら数値予報モデルの計算結果はGPV（Grid Point Value）と呼ばれ、データファイル形式で提供されています。
GPVデータファイルとは
GPVデータファイルは、数値予報モデルが計算した風、気温、湿度、気圧などの物理パラメーターが２次元の配列として記録されたバイナリファイルです。
空間要素の場合は、2次元データを鉛直方向に重ねて３次元で表現することができます。また、予報時間の軸を含むため、GPVデータファイルには最大で4次元の物理パラメータが記録されています。
これらのGPVデータファイルは、世界気象機関（WMO）が定めるFM92 GRIB 二進形式格子点資料気象通報式（第2版）に従ってフォーマットされています。通常、このデータはGRIB2として広く知られています。
GRIB2から任意地点のデータを抽出する
ここからは、Pythonを使ってGRIB2から任意地点のデータを抽出していきます。
今回は、気象庁メソモデル（MSM）のGRIB2を使用します。企業活動等の営利目的利用の場合、日本国内では気象業務支援センターから購入することができますが、ここは気象データ利用に関する普及活動を目的にしているため、北海道大学、京都大学及び九州大学によって共同運営されている「流体電脳倶楽部」からデータをお借りしました。
Pythonの環境を準備する
この実験で使用するLinuxサーバの環境は以下のとおりです。コンテナ上で環境を作成しています。
| OS（ホスト） | Red Hat Enterprise Linux release 8.8 (Ootpa) |
| コンテナ仮想化（ホスト） | podman version 4.4.1 |
| OS（コンテナ） | AlmaLinux release 9.3 （Shamrock Pampas Cat） |
| Python（コンテナ） | Python 3.12.1 |
はじめに、Pythonの仮想環境を作成します。以下の手順を参考にして、ご自身の環境にあわせて作成してください。
$ mkdir -p Python/venv/Pygrib
$ python3.12 -V
Python 3.12.1
$ python3 -m venv Python/venv/Pygrib
$ source Python/venv/Pygrib/bin/activate
(Pygrib) [reafnex@devsrv01 ~]$ python -V
Python 3.12.1Pythonの仮想環境が起動できたら、実験に必要なPythonのモジュールをインストールしていきます。
Pythonは、GRIB2を処理するためのPygribというモジュールが利用可能です。さらにPandasとpyarrowもあわせてインストールします。
Pythonのモジュールは、pipコマンドを使用してインストールします。
(Pygrib) [reafnex@devsrv01 ~]$ pip3 install pygrib pandas pyarrow
(Pygrib) [reafnex@devsrv01 ~]$ pip3 list
Package Version
--------------- --------
certifi 2024.2.2
numpy 1.26.3
pandas 2.2.0
pip 23.2.1
pyarrow 15.0.0
pygrib 2.1.5
pyproj 3.6.1
python-dateutil 2.8.2
pytz 2024.1
setuptools 69.0.3
six 1.16.0
tzdata 2023.4Pythonの仮想環境にインストールされたPythonモジュールは、pip3 listコマンドで確認することができます。pipインストール時に指定したモジュールの他、依存関係があるモジュールも一緒にインストールされています。
ちなみに、Windowsでpipを使ってPygribをインストールすることはできませんでした。Anacondaを利用すればWindowsでもPygribを利用できるようです。
Pygribを使ったデータアクセス
今回は、Pygribモジュールを使って、気象庁MSMのGRIB2から特定の地点の気温と相対湿度を抽出してみます。
Pythonコードの全体は、以下の通りとなります。
import pygrib
import pandas as pd
from datetime import timedelta
# UTCから日本時間に変換するための時刻差を設定する。
diff_time = timedelta(hours=9)
# GPVファイル名を指定してオープンする。
gpv_file = pygrib.open("Z__C_RJTD_20240204000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin")
# 要素データを取り出す。
data_tk = gpv_file.select(name="Temperature")
data_rh = gpv_file.select(name="Relative humidity")
# 予想時刻と抽出した要素データを格納するリストを初期化する。
valid_dates = []
tc_list = []
rh_list = []
# 東京タワーの緯度経度を設定する。
req_lat = 35.65858404079
req_lon = 139.74543164468
# MSMの格子間隔の半分を設定する。
diff_lat = 0.025
diff_lon = 0.03125
# 気温と相対湿度を抽出する。
for tk, rh in zip(data_tk, data_rh):
# 予想時刻を出力する。
valid_dates.append(tk.validDate + diff_time)
# 気温を抽出して摂氏に単位変換する。
tk_pt = tk.data(
lat1 = req_lat - diff_lat, lat2 = req_lat + diff_lat,
lon1 = req_lon - diff_lon, lon2 = req_lon + diff_lon
)[0][0][0]
tc_pt = tk_pt - 273.15
tc_list.append(tc_pt)
# 相対湿度を抽出する。
rh_pt = rh.data(
lat1 = req_lat - diff_lat, lat2 = req_lat + diff_lat,
lon1 = req_lon - diff_lon, lon2 = req_lon + diff_lon
)[0][0][0]
rh_list.append(rh_pt)
# 予想時間とWBGT計算結果をPandasのデータフレームに格納する。
df = pd.DataFrame({
"validDate(JST)": valid_dates,
"Temp(C)": tc_list,
"Humidity(%)": rh_list
})
# データフレームを表示する。
print(df)上記のコードを実行すると、以下のような結果が出力されます。
(Pygrib) [reafnex@devsrv01 Pygrib]$ python3 main.py
validDate(JST) Temp(C) Humidity(%)
0 2024-02-04 09:00:00 3.856439 86.514521
1 2024-02-04 10:00:00 3.486642 88.784033
2 2024-02-04 11:00:00 3.428369 90.453701
3 2024-02-04 12:00:00 3.837823 90.388992
4 2024-02-04 13:00:00 4.632776 83.865077
5 2024-02-04 14:00:00 5.661386 73.763103
6 2024-02-04 15:00:00 6.395868 64.078033
7 2024-02-04 16:00:00 6.676202 60.489502
8 2024-02-04 17:00:00 6.421426 62.950851
9 2024-02-04 18:00:00 5.880045 68.590605
10 2024-02-04 19:00:00 5.697641 66.246988
11 2024-02-04 20:00:00 5.329889 61.383592
12 2024-02-04 21:00:00 4.583215 55.026819
13 2024-02-04 22:00:00 4.259134 54.676161
14 2024-02-04 23:00:00 4.132852 55.833216
15 2024-02-05 00:00:00 3.894037 58.245659Pygribの使い方
まずは、モジュールをインポートしてPygribを利用できるようにします。
import pygribpygrib.openで気象庁MSMのGRIB2ファイルをオープンし、ファイルハンドラ（gpv_file）を作成します。
# GPVファイル名を指定してオープンする。
gpv_file = pygrib.open("Z__C_RJTD_20240204000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin")次に、ファイルの中から気温と相対湿度のGRIBメッセージを抽出してオブジェクト変数を作成します。
# 要素データを取り出す。
data_tk = gpv_file.select(name="Temperature")
data_rh = gpv_file.select(name="Relative humidity")気象庁MSMに含まれているGRIBメッセージは、GRIB2のファイルハンドラ（ここではgpv_file）をループでまわしてprintすることで確認できます。
for msg in gpv_file:
print(msg)以下のような結果が得られます。
1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 202402040000
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202402040000
3:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202402040000
4:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202402040000
5:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 0 hrs:from 202402040000
6:Relative humidity:% (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 0 hrs:from 202402040000
7:Low cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202402040000
8:Medium cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202402040000
9:High cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202402040000
10:Total cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202402040000
11:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 1 hrs:from 202402040000
12:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 1 hrs:from 202402040000
13:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 1 hrs:from 202402040000
14:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 1 hrs:from 202402040000
15:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 1 hrs:from 202402040000
16:Relative humidity:% (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 1 hrs:from 202402040000
17:Low cloud cover:% (instant):regular_ll:surface:level 0:fcst time 1 hrs:from 202402040000
18:Medium cloud cover:% (instant):regular_ll:surface:level 0:fcst time 1 hrs:from 202402040000
19:High cloud cover:% (instant):regular_ll:surface:level 0:fcst time 1 hrs:from 202402040000
20:Total cloud cover:% (instant):regular_ll:surface:level 0:fcst time 1 hrs:from 202402040000
21:Total precipitation:kg m-2 (accum):regular_ll:surface:level 0:fcst time 0-1 hrs (accum):from 202402040000
・・<以下省略>・・東京タワーの所在地付近のデータを抽出するため、緯度経度を設定します。また、設定した緯度経度に最も近い格子点のデータを１つだけ抽出するため、気象庁MSMの格子点間隔（緯度方向：0.05°、経度方向：0.0625°）の半分を設定しています。
# 東京タワーの緯度経度を設定する。
req_lat = 35.65858404079
req_lon = 139.74543164468
# MSMの格子間隔の半分を設定する。
diff_lat = 0.025
diff_lon = 0.03125ここから、東京タワーに最も近い格子点のデータを抽出します。気温と相対湿度のGRIBメッセージを抽出したオブジェクト変数をループに入力して処理していきます。
気温データを抽出したオブジェクト変数からGRIB2メッセージに含まれている予想時間を抽出して、リスト変数に追記しています。
気温と相対湿度は、データを抽出する範囲の緯度（lat1とlat2）及び経度（lon1とlon2）を指定して抽出しています。今回は、東京タワーに最も近い1点の格子点データを抽出しています。
気象庁MSMに含まれている気温は絶対温度（K）であるため、摂氏（℃）に変換してからリスト変数に追記しています。相対湿度の単位はパーセントで、同じくリスト変数に追記しています。
# 気温と相対湿度を抽出する。
for tk, rh in zip(data_tk, data_rh):
# 予想時刻を出力する。
valid_dates.append(tk.validDate + diff_time)
# 気温を抽出して摂氏に単位変換する。
tk_pt = tk.data(
lat1 = req_lat - diff_lat, lat2 = req_lat + diff_lat,
lon1 = req_lon - diff_lon, lon2 = req_lon + diff_lon
)[0][0][0]
tc_pt = tk_pt - 273.15
tc_list.append(tc_pt)
# 相対湿度を抽出する。
rh_pt = rh.data(
lat1 = req_lat - diff_lat, lat2 = req_lat + diff_lat,
lon1 = req_lon - diff_lon, lon2 = req_lon + diff_lon
)[0][0][0]
rh_list.append(rh_pt)最後に、予想時間、気温及び相対湿度を格納したリスト変数をPandasのデータフレームに代入します。Pandasのデータフレームを利用することで、表計算したりCSVファイルやExcelファイルとして出力することが簡単にすることができます。
# 予想時間とWBGT計算結果をPandasのデータフレームに格納する。
df = pd.DataFrame({
"validDate(JST)": valid_dates,
"Temp(C)": tc_list,
"Humidity(%)": rh_list
})
# データフレームを表示する。
print(df)まとめ
いかがでしたでしょうか。
今回は、PythonのPygribモジュールを利用して、気象庁MSMのGRIB2から特定地点の気温と相対湿度のデータを抽出するコードを紹介しました。
気象データは私たちの日常生活や多くの産業にとって重要なインテリジェンスです。気象データをニーズに合わせて適切に処理し分析するテクニックは、社会課題の解決に必要不可欠なものです。
気象データを機械学習やディープラーニングを使用して処理することで、農業や様々な産業分野で新しい知見を得られる可能性があり、これにより生産性向上やフードロスの改善に寄与するかもしれません。
そのため、気象データを機械学習やディープラーニングと相性が良いPythonで処理することについては、大きな意義があります。ただし、日本では、確定論的に扱うことができない気象データをビジネス創出などに活用できるエンジニアが不足しており、この状況を改善するため、気象庁は「気象データアナリスト育成講座の認定制度」の運用を開始して、DX（デジタルトランスフォーメーション）を加速させるための気象データの適切な利用促進と人材育成を行っています。
本サイトでは、気象データの活用方法について普及活動を行い、社会課題解決とビジネス創出に基づくSDGsの目標達成に向けた活動を推進しています。
参考になれば幸いです。
気象データ利用のお悩みについてご相談ください
気象システムの設計製造及び運用保守にかかわる現役システムエンジニアが、気象予報士の視点で気象データを活用したビジネスアイディアの創出と課題解決のお手伝いをさせていただきます。
気象や波浪などの数値予報モデルを活用した気象予報システム開発や運用についても是非ご相談ください。
