はじめに
ここでは気象庁のサンプルデータ一覧のうち、以下のファイルを対象とする。
- 全球数値予報モデルGPV (GSM)
- メソ数値予報モデルGPV (MSM)
- 局所数値予報モデルGPV (LFM)
いずれもGRIB2というバイナリ形式が採用されている。
検証環境
Anacondaの仮想環境を使用した。
- Windows 10 Home
- conda 4.8.2
- Python 3.7.7
ライブラリインストール
記事タイトルの通り、pygribを利用する。pipやcondaからインストール可能。
conda install -c conda-forge pygrib
実装例
気象庁のサンプルファイル(MSM)から各時刻(JST)における神宮球場(35.6745, 139.7169)(に一番近いメッシュ地点)の気温[℃]を取得する。
※有効数字があやしいですが気にしないでください。
from datetime import timedelta
import pygrib
import pandas as pd
time_diff = timedelta(hours=9)
gpv_file = pygrib.open("Z__C_RJTD_20171205000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin")
t_messages = gpv_file.select(name="Temperature")
df = pd.DataFrame({
"validDate": [msg.validDate + time_diff for msg in t_messages],
"temperature": [
msg.data(
lat1=35.6745-0.025,
lat2=35.6745+0.025,
lon1=139.7169-0.03125,
lon2=139.7169+0.03125,
)[0][0][0] - 273.15 for msg in t_messages
]
})
print(df)
↓出力結果
validDate temperature
0 2017-12-05 09:00:00 9.810709
1 2017-12-05 10:00:00 11.409219
2 2017-12-05 11:00:00 12.358789
3 2017-12-05 12:00:00 13.116861
4 2017-12-05 13:00:00 13.770517
5 2017-12-05 14:00:00 14.698541
6 2017-12-05 15:00:00 14.488687
7 2017-12-05 16:00:00 13.063196
8 2017-12-05 17:00:00 11.467722
9 2017-12-05 18:00:00 10.320886
10 2017-12-05 19:00:00 9.592111
11 2017-12-05 20:00:00 8.797952
12 2017-12-05 21:00:00 8.171686
13 2017-12-05 22:00:00 7.832407
14 2017-12-05 23:00:00 7.461450
15 2017-12-06 00:00:00 6.884409
実装解説
ファイル接続
pygrib.openクラスを使ってファイルとのインターフェースを作成する。
import pygrib
gpv_file = pygrib.open("Z__C_RJTD_20171205000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin")
-
pygrib.openインスタンス自身がpygrib.gribmessageインスタンスのイテレータになっている - ファイルに含まれる特定の情報を
pygrib.gribmessageインスタンスまたはそのリストとして抽出するためのインターフェースとして利用できる(再利用可)
※pygrib.gribmessageについては後述
messageメソッド
引数N(int)をとり、イテレータ上のN番目(1オリジン)のpygrib.gribmessageインスタンスを返却する。
In[]
gpv_file.message(2)
↓
Out[]
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201712050000
selectメソッド
**kwargs形式の引数をとり、その条件に合致するpygrib.gribmessageのリストを返却する。
複数の検索条件を指定した場合はAND条件として評価される。
- 1条件
In[]
gpv_file.select(name="Temperature")
↓
Out[]
[5:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 0 hrs:from 201712050000,
15:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 1 hrs:from 201712050000,
27:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 2 hrs:from 201712050000,
39:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 3 hrs:from 201712050000,
51:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 4 hrs:from 201712050000,
63:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 5 hrs:from 201712050000,
75:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 6 hrs:from 201712050000,
87:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 7 hrs:from 201712050000,
99:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 8 hrs:from 201712050000,
111:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 9 hrs:from 201712050000,
123:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 10 hrs:from 201712050000,
135:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 11 hrs:from 201712050000,
147:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 12 hrs:from 201712050000,
159:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 13 hrs:from 201712050000,
171:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 14 hrs:from 201712050000,
183:Temperature:K (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 15 hrs:from 201712050000]
- 2条件
In[]
gpv_file.select(shortName="r", forecastTime=3)
↓
Out[]
[40:Relative humidity:% (instant):regular_ll:heightAboveGround:level 1.5 m:fcst time 3 hrs:from 201712050000]
現実的に検索条件に使えそうな項目は以下の通り。
| 項目名 | 入力型 |
|---|---|
| name | str |
| shortName | str |
| forecastTime(※) | int |
| validDate | datetime.datetime |
| parameterCategory | int |
| parameterNumber(※※) | int |
※forecastTimeはGSM, MSMでの単位は[時間]だがLFMでの単位は[分]なので注意
※※parameterNumberはparameterCategoryと併用して使用すべき
データ値取得
上で抽出したpygrib.gribmessageインスタンスにはある時刻・あるデータ項目の全メッシュのデータが含まれている。
主な変数
| 変数名 | 型 | 内容 |
|---|---|---|
| validDate | datetime.datetime | 予報の日時(UTC) |
| analDate | datetime.datetime | 予報の解析基準日時(UTC) |
| forecastTime | int |
analDateからvalidDateまでの時刻差 |
| distinctLatitudes | numpy.ndarray | メッシュの緯度一覧(1次元) |
| distinctLongitudes | numpy.ndarray | メッシュの経度一覧(1次元) |
| values | numpy.ndarray | 全メッシュのデータ値一覧(2次元) |
※ファイル内の日時はすべてUTCになっている。
latlonsメソッド
メッシュの緯度・経度を全メッシュごとに羅列したものを返却する。戻り値の型はどちらもnumpy.ndarray。
In[]
message = gpv_file.message(2)
lats, lons = message.latlons()
Out[]
# lats
array([[47.6 , 47.6 , 47.6 , ..., 47.6 , 47.6 , 47.6 ],
[47.55, 47.55, 47.55, ..., 47.55, 47.55, 47.55],
[47.5 , 47.5 , 47.5 , ..., 47.5 , 47.5 , 47.5 ],
...,
[22.5 , 22.5 , 22.5 , ..., 22.5 , 22.5 , 22.5 ],
[22.45, 22.45, 22.45, ..., 22.45, 22.45, 22.45],
[22.4 , 22.4 , 22.4 , ..., 22.4 , 22.4 , 22.4 ]])
# lons
array([[120. , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150. ],
[120. , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150. ],
[120. , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150. ],
...,
[120. , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150. ],
[120. , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150. ],
[120. , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150. ]])
行列として表記した場合と地図上に表示した場合が視覚的に一致するようにこのような順番になっていると思われる。
dataメソッド
引数lat1(緯度最小)、lat2(緯度最大)、lon1(経度最小)、lon2(経度最大)の範囲内にあるメッシュのデータ値とそのメッシュの緯度・経度を取得する。戻り値はデータ値・緯度・経度それぞれのメッシュ様2次元numpy.ndarrayがタプルになっている。
In[]
message.data(lat1=22.45, lat2=22.55, lon1=120.0625, lon2=120.1875)
Out[]
(array([[101985.52856445, 101985.52856445, 101985.52856445],
[101985.52856445, 101985.52856445, 101985.52856445]]),
array([[22.5 , 22.5 , 22.5 ],
[22.45, 22.45, 22.45]]),
array([[120.0625, 120.125 , 120.1875],
[120.0625, 120.125 , 120.1875]]))