pygribでgrbsファイルを開いてみた

grbsを使ってみた。
最初sqliteみたいな感じなのかなと思っていた。

import pygrib
grbs = pygrib.open("dataset/gfs_4_20200921_0000_000.grb2")
g = grbs.select(name="V component of wind")[0]
g

9:V component of wind:m s**-1 (instant):regular_ll:unknown:level 0:fcst time 0 hrs:from 202009210000

g.latlons()で、grbsはクラスでfunctionを持っているんだと気づく。

lats, lons = g.latlons()
lats, lons

(array([[ 90. ,  90. ,  90. , ...,  90. ,  90. ,  90. ],
        [ 89.5,  89.5,  89.5, ...,  89.5,  89.5,  89.5],
        [ 89. ,  89. ,  89. , ...,  89. ,  89. ,  89. ],
        ...,
        [-89. , -89. , -89. , ..., -89. , -89. , -89. ],
        [-89.5, -89.5, -89.5, ..., -89.5, -89.5, -89.5],
        [-90. , -90. , -90. , ..., -90. , -90. , -90. ]]),
 array([[  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
        [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
        [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
        ...,
        [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
        [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
        [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5]]))

Gribとは

そもそもGribって何ってなった。
GRIB は、グリッド化されたデータを配布するための世界気象機関の標準規格で、pygribがECMWFのGRIB API C ライブラリへのpythonインターフェースらしい。

githubで調べたら多分これかなというのがあった。
latlons()のメソッドのURL

インターフェースざっくり見てもどう使うのかピンとこなかった。
こっちのExampleのが分かりやすいかも
https://jswhit.github.io/pygrib/docs/

keyをたくさん持っているみたいなので、keyを調べてみた。
めっちゃ多いから結果は省略した

g.keys()

['globalDomain',
 'GRIBEditionNumber',
 'tablesVersionLatest',
 'grib2divider',
 'is_efas',
 'angularPrecision',

 ...(省略)

 'md5Section7',
 'section8Length',
 'analDate',
 'validDate']

numpy arrayとして扱えるらしい

  maxt = grb.values
  maxt.shape, maxt.min(), maxt.max()

(94, 192) 223.7 319.9

(94, 192)は、何キロかごとに区切ったのグリッドだと推測している。

Levelを並べ替える

levels = np.array([g['level'] for g in grb])
indexes = np.argsort(levels)[::-1]
indexes

array([22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10,  9,  8,  7,  6,
        5,  4,  3,  2,  1,  0])

グリッドからLevelを加えた3次元のデータを作る

x = grb[0].values.shape[0]
y = grb[1].values.shape[1]
cube = np.zeros([len(grb), x, y])
for i in range(len(grb)):
    cube[i,:,:] = grb[indexes[i]].values
cube_dict = {
    'data': cube,
    'units': grb[0]['units'],
    'levels': levels[indexes]
}    
cube.shape

(23, 361, 720)

とりあえず欲しい値は取れた。
BaseMapに描画してみたり、等圧線とか書きたい

この辺りを参考にしたら出来そう
https://github.com/scollis/HRRR/blob/master/notebooks/HRRR%20Grib.ipynb