気象庁のGRIB2形式のファイルを読む

# はじめに
juyter-notebookでpygribを使用し，[気象庁のＬＦＭ格子点データ](https://www.data.jma.go.jp/add/suishin/cgi-bin/catalogue/make_product_page.cgi?id=KyoModel)を読んだ．
## ファイルをダウンロード
[京都大学生存圏研究所のデータベース](http://database.rish.kyoto-u.ac.jp/)からデータをダウンロードする．
ダウンロードが成功すれば下記の「Out」のように出力される．
```python:In
import urllib.request
from bs4 import BeautifulSoup

url = 'http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2020/08/01/'

# urlの情報を取得
response = urllib.request.urlopen(url)
soup = BeautifulSoup(response)
response.close()

## a href にファイル名が入っていそうなので，１行ずつリストに入れる
file_list = []
for row in soup.findAll('a'):
    file_list.append(row.get('href'))

## binファイルのみ抽出
file_list = [file for file in file_list if file.endswith('.bin')==True]

# ダウンロードしたファイルを置く場所のPATH
DL_PATH = 'F:/'

# file_list中の始めのファイルのみダウンロード
urllib.request.urlretrieve(url+fname, DL_PATH+file_list[0])
```
```python:Out
('F:/Z__C_RJTD_20200801000000_CWM_GPV_Rjp_Gll0p05deg_FD0000-0300_grib2.bin',
 <http.client.HTTPMessage at 0x1d09ef42d30>)
```


## ファイルを開く
ファイルを開き，`read()`することでgribmessageのリストとして中身を確認できる．
```python:In
import pygrib

grbs = pygrib.open('F:/Z__C_RJTD_20171205010000_LFM_GPV_Rjp_L-pall_FH0100_grib2.bin')
grbs.read()
```
```python:Out
[1:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 1 mins:from 201712050100,
 2:U component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 1 mins:from 201712050100,
 3:V component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 1 mins:from 201712050100,
...
 90:V component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 10000.0 Pa:fcst time 1 mins:from 201712050100,
 91:Temperature:K (instant):regular_ll:isobaricInhPa:level 10000.0 Pa:fcst time 1 mins:from 201712050100,
 92:Vertical velocity:Pa s**-1 (instant):regular_ll:isobaricInhPa:level 10000.0 Pa:fcst time 1 mins:from 201712050100]

```
### selectメソッド
引数の条件でgribmessageを検索できる．
| 引数 | 入力型 |
|:-:|:-:|
| name  | str  |
|  shortName | str  |
|  units |  str |
|  forecastTime | int  |
```python:In
grbs.select(name='Geopotential height')
grbs(name='Geopotential height')          # selectメソッドと同じことができる
```
```python:Out
[1:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 1 mins:from 201712050100,
 7:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 97500.0 Pa:fcst time 1 mins:from 201712050100,
 13:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 95000.0 Pa:fcst time 1 mins:from 201712050100,
...
 78:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 20000.0 Pa:fcst time 1 mins:from 201712050100,
 83:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 15000.0 Pa:fcst time 1 mins:from 201712050100,
 88:Geopotential height:gpm (instant):regular_ll:isobaricInhPa:level 10000.0 Pa:fcst time 1 mins:from 201712050100]
```
## 変数の読み込み
`select`メソッドで選んだ変数には様々な情報が格納されている．
```python:In
grb = grbs.select(name='Geopotential height')[0]

# keyを取得し，その中に格納されている値を読む
for grb_key in grb.keys():
    if grb.valid_key(grb_key)==True:       # keyの欠落がなく，値が読み取れるものを出力
        print(grb_key, ' : ',grb[grb_key])
    else:                                  # 読めなかったkeyは赤字で出力
        print('\033[31m'+grb_key+' : Can not read'+'\033[0m')
```
```python:Out
globalDomain  :  g
GRIBEditionNumber  :  2
tablesVersionLatestOfficial  :  29
tablesVersionLatest  :  29
...
productType  :  unknown
md5Section7  :  4caa43e9fc2f9bd2ba637aed3af7dc72
section8Length  :  4
```
### valuesメソッド
数値データを含むmasked_arrayを返す．
```python:In
grb = grbs(name='Temperature')[0]
grb.values
```
```python:Out
masked_array(
  data=[[--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        [--, --, --, ..., --, --, --],
        ...,
        [291.39064025878906, 291.29689025878906, 291.23439025878906, ...,
         --, --, --],
        [291.48439025878906, 291.40626525878906, 291.34376525878906, ...,
         --, --, --],
        [291.59376525878906, 291.51564025878906, 291.48439025878906, ...,
         --, --, --]],
  mask=[[ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        [ True,  True,  True, ...,  True,  True,  True],
        ...,
        [False, False, False, ...,  True,  True,  True],
        [False, False, False, ...,  True,  True,  True],
        [False, False, False, ...,  True,  True,  True]],
  fill_value=9999.0)
```
### lationsメソッド
緯度，経度データを返す．
```python:In
grb = grbs(name='Temperature')[0]
grb.latlons()
```
```python:Out
(array([[47.6 , 47.6 , 47.6 , ..., 47.6 , 47.6 , 47.6 ],
        [47.56, 47.56, 47.56, ..., 47.56, 47.56, 47.56],
        [47.52, 47.52, 47.52, ..., 47.52, 47.52, 47.52],
        ...,
        [22.48, 22.48, 22.48, ..., 22.48, 22.48, 22.48],
        [22.44, 22.44, 22.44, ..., 22.44, 22.44, 22.44],
        [22.4 , 22.4 , 22.4 , ..., 22.4 , 22.4 , 22.4 ]]),
 array([[120.  , 120.05, 120.1 , ..., 149.9 , 149.95, 150.  ],
        [120.  , 120.05, 120.1 , ..., 149.9 , 149.95, 150.  ],
        [120.  , 120.05, 120.1 , ..., 149.9 , 149.95, 150.  ],
        ...,
        [120.  , 120.05, 120.1 , ..., 149.9 , 149.95, 150.  ],
        [120.  , 120.05, 120.1 , ..., 149.9 , 149.95, 150.  ],
        [120.  , 120.05, 120.1 , ..., 149.9 , 149.95, 150.  ]]))
```
### dataメソッド
数値データ，緯度・経度データを返す．
出力結果は`values`メソッドと`lations`メソッドの結果を合わせたような感じ．
```python:In
grb = grbs(name='Temperature')[0]
grb.data()
```
### その他のメソッド
keyに入っている情報をメソッドとして取り出せるのでよく使いそうなものをまとめた．
| メソッド | 戻り値の内容 |
|:-:|:-:|
| .name | 変数名  |
| .shortName  | 変数名の省略形  |
| .units  | 単位  |
| .maximum  |  最大値 |
| .minimum  |  最小値 |
| .average  | 平均値  |
| .validDate  |  予報時刻(UTC) |
| .analDate  |  解析基準時刻(UTC) |
| .forecastTime  | 予報時間  |
| .missingValue  |欠損値|
## おまけ
ファイルの中にどのような変数があるのかを簡単に見るために，気になる要素を抽出して出力するようにした．
```python:In
import pygrib
import pandas as pd

grbs = pygrib.open('F:/Z__C_RJTD_20171205010000_LFM_GPV_Rjp_L-pall_FH0100_grib2.bin')


name = []
s_name  = []
unit = []
level = []
level_type = []

# 各リストに要素を追加
for grb in grbs:
    name.append(grb.name)
    s_name.append(grb.shortName)
    unit.append(grb.units)
    level.append(grb.level)
    level_type.append(grb.typeOfLevel)


# listをDataFrameに変換
df = pd.DataFrame([s_name, name, unit, level, level_type])
df = df.T
df.columns=['shortName', 'name', 'units', 'level', 'typeOfLevel']

# 行と列をすべて表示させる
pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)

print(df)

```
```python:Out
   shortName                 name     units level    typeOfLevel
0         gh  Geopotential height       gpm  1000  isobaricInhPa
1          u  U component of wind   m s**-1  1000  isobaricInhPa
2          v  V component of wind   m s**-1  1000  isobaricInhPa
3          t          Temperature         K  1000  isobaricInhPa
4          w    Vertical velocity  Pa s**-1  1000  isobaricInhPa
...
87        gh  Geopotential height       gpm   100  isobaricInhPa
88         u  U component of wind   m s**-1   100  isobaricInhPa
89         v  V component of wind   m s**-1   100  isobaricInhPa
90         t          Temperature         K   100  isobaricInhPa
91         w    Vertical velocity  Pa s**-1   100  isobaricInhPa
```
# 参考文献
- [pygrib documentation >> API Example usage](https://jswhit.github.io/pygrib/api.html#example-usage)
- [pygribで気象庁の数値予報GPVデータを読み込む](https://predora005.hatenablog.com/entry/2020/10/31/000000)
- [pygribで気象庁のGRIB2ファイルを読む](https://qiita.com/kurukuruz/items/6fc0be9efa34a2fd6741)