pygribで気象庁のGRIB2ファイルを読む

# はじめに

ここでは[気象庁のサンプルデータ一覧](https://www.data.jma.go.jp/developer/gpv_sample.html)のうち、以下のファイルを対象とする。

- 全球数値予報モデル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)(に一番近いメッシュ地点)の気温[℃]を取得する。
※有効数字があやしいですが気にしないでください。

```python
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](https://jswhit.github.io/pygrib/docs/pygrib.open-class.html)クラスを使ってファイルとのインターフェースを作成する。

```python
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`インスタンスを返却する。

```python: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条件

```python: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条件

```python: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](https://jswhit.github.io/pygrib/docs/pygrib.gribmessage-class.html)インスタンスには**ある時刻・あるデータ項目の全メッシュのデータ**が含まれている。

### 主な変数

| 変数名             | 型                | 内容                                  |
| ------------------ | ----------------- | ------------------------------------- |
| 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`。

```python: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`がタプルになっている。

```python: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]]))
```



# 参考

- [pygrib](https://jswhit.github.io/pygrib/docs/)
- [ファイル形式データ - 気象業務支援センター](http://www.jmbsc.or.jp/jp/online/f-online0.html)
- [配信資料に関する技術情報 第 500 号](https://www.data.jma.go.jp/add/suishin/jyouhou/pdf/500.pdf) ※pdf
- [ECMWF | Parameter database](https://apps.ecmwf.int/codes/grib/param-db)