pygribでオープンできないGRIB2形式の気象データをnetCDFに変換して可視化する

# pygribでオープンできないGRIB2形式の気象データをnetCDFに変換して可視化する

### 2024.5　一部修正

この記事で使用しているデータ(1kmメッシュ全国合成レーダーGPV)について、気象庁からの配信が終了しました。
詳細は気象庁資料をご覧ください。

[～ 10分毎１kmメッシュ全国合成レーダーGPVの配信終了 ～
（配信資料に関する技術情報第 568 号関連）](https://www.data.jma.go.jp/add/suishin/oshirase/pdf/20240124.pdf)

これは技術の進歩によって、データ配信間隔が10分から5分に変更となったり、メッシュ間隔が1kmから250mになったりしたためです。
気象庁によるとしばらくの間、250m間隔メッシュを低解像度化した1kmメッシュが提供されるそうです。
これに伴い、元の投稿での以下が変更になっています。

・ファイル名
 ファイル名途中の `_Prr10lv_`が、`_Prr05lv_`に変更
・パラメータ名
 Netcdfファイルでの指定パラメータが`var0_1_201_surface`から`var0_1_203_surface`に変更



## はじめに

面的立体的な広がりを持った気象データのフォーマットは「GRIB2」形式が業界標準です。
GRIB2形式については@e_toyodaさんによる詳しい紹介があります。
[WMO（世界気象機関）の格子データ形式GRIB2について](https://qiita.com/e_toyoda/items/ce7497e1a633b16f1ff1)

GRIB2形式のファイルからデータを読み出すための様々なツールがあります。
pythonならばpygribというモジュールを利用できます。

ところが、全国合成レーダーや推計気象分布などは全国1km間隔でのデータであることからデータ量が大きく、GSMやMSMとは違った、日本気象庁が独自に指定している圧縮形式が用いられています。
これがpygribでサポートされていないようで、オープンできないデータが存在しています。

＊全国合成レーダは以下に解説があります。
[１kmメッシュ全国合成レーダーGPV](http://www.jmbsc.or.jp/jp/online/file/f-online30100.html)

このデータを用いたDeep LearningプログラムをKerasで作りたいので、pythonで扱えるようにしたいところです。
そこで気象庁フォーマットをサポートしているwgrib2というツールを使用して、pythonで扱えるnetCDF形式に変換しました。


## pygribでオープンできない全国合成レーダ
全国合成レーダは[京都大学生存圏研究所](http://database.rish.kyoto-u.ac.jp/arch/jmadata/)からダウンロードしました。
このデータファイルをpygribでオープンすると以下のエラーが表示されます。

```python
>>> import pygrib
>>> grb=pygrib.open("Z__C_RJTD_20240331230000_RDR_JMAGPV_Ggis1km_Prr05lv_ANAL_grib2.bin")
ECCODES ERROR   :  Unable to find template productDefinition from grib2/template.4.50008.def 
ECCODES ERROR   :  Unable to find template productDefinition from grib2/template.4.50008.def 
```

これは気象庁が使用している定義がpygribの中に存在していないためと思われます。


## wgrib2を用いた変換

### wgrib2

wgrib2については、@ysomeiさんが下記に纏めておられます。
アメリカの「大気海洋庁（NOAA）」が提供しているツールです。
[今更ですが、気象データを解析してみた](https://qiita.com/ysomei/items/12d6622bed030f6ac793)

色々なオプションがあって様々なことができるのですが、例えば以下のようにGRIBファイルの内容をCSVに出力するとことができます。
~~できたCSVファイルをのぞくと、データを指定するラベルに`var0_1_201`と`surface`が入っています。~~
できたCSVファイルをのぞくと、データを指定するラベルに`var0_1_203`と`surface`が入っています。
最後の３要素が、経度、緯度、データ（「データ代表値」）です。

```bash:
> wgrib2 -V GRIB_formatted_file -csv tmp.csv
> cat tmp.csv
略
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.419,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.431,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.444,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.456,46.0125,0
"2024-03-31 23:00:00","2024-03-31 23:00:00","var0_1_203","surface",141.469,46.0125,0
後略
```

### 合成レーダのデータ形式（データ代表値とレベル値）
合成レーダの「データ代表値」は降雨強度であるmm/hが単位です。
元のGRIB2形式のファイルには、直接降雨強度の数値がセットされてはおらず、「レベル値」という0から251の値がセットされています。
例えば「レベル値」が「37」というのは、降雨強度が「7.0mm/h以上、7.5mm/h未満」ということを意味しています。
その「データ代表値」は「7.75mm/h」であるということです。

つまり、「レベル値」から「データ代表値」への変換テーブルが存在していることになりますが、wgrib2はそこもやってくれるようです（wgrib2を解読した訳ではありませんが）。

### NetCDF形式への変換
wgrib2を用いて、GRIB2からnetCDF形式へ変換するには下記のオプションを使います。

```bash:wgrib.sh
wgrib2 GRIB_formatted_file -netcdf GRIB_formatted_file.nc
```
これでnetCDF形式ファイルが出来上がります。
ただし、圧縮が解除されるためか、猛烈にファイルサイズが大きくなりますので注意が必要です。

```bash:ls
103424   Z__C_RJTD_20240331230000_RDR_JMAGPV_Ggis1km_Prr05lv_ANAL_grib2.bin
34454836 Z__C_RJTD_20240331230000_RDR_JMAGPV_Ggis1km_Prr05lv_ANAL_grib2.bin.nc
```

## pythonでnetCDF4を用いたデータの読み込み

このファイルはpythonモジュールnetCDF4を用いると下記のように取り扱うことができます。

```python:netcdf4_open.py
from netCDF4 import Dataset
(中略）

nc = Dataset( fn , 'r' )
xx = nc.variables['longitude'][:] #
yy = nc.variables['latitude' ][:] #
data = nc.variables['var0_1_203_surface'][:]
lons, lats = np.meshgrid( xx , yy )
        
_save_file = savefilenamelist[ptr_sfile]

drawmap( data[0,:,:] , levels_rai , precip_colormap , _save_file ,_flag_show , "RAIN" , norm)
```

ファイルに含まれているデータを特定するために、`var0_1_203_surface`を指定しています。
wgrib2でGRIBファイルの中身を参照した際には、`var0_1_203`と`surface`に分かれて見えていましたが、NetCDFファイルに変換した際に`var0_1_203_surface`というラベルに変換されるようです。
下記のように、netcdfファイルを情報出力させると、`variables`の最後にこのラベルが出てきています。

```Python:DumpNetCDF
>>> from netCDF4 import Dataset
>>> nc=Dataset("www.nc")
>>> nc
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    Conventions: COARDS
    History: created by wgrib2
    GRIB2_grid_template: 0
    dimensions(sizes): latitude(3360), longitude(2560), time(1)
    variables(dimensions): float64 latitude(latitude), float64 longitude(longitude), float64 time(time), float32 var0_1_203_surface(time,latitude,longitude)
    groups:
```

### 可視化

可視化するところのpythonソースです。
matplotlibを用いています。

手前味噌ではありますが、可視化のあれこれは、[気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（2) 入力となる気象データ可視化の話](https://qiita.com/m-taque/items/988b08185097dca5f5b5 "気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（2) 入力となる気象データ可視化の話")に記載しております。
可視化するカラーマップに白色のレンジがあるので、下記の可視化では地図そのものに少し色をつけています。
色のついた地図とコンターマップの重ね合わせの順番をzorderで指定しています。

```python:draw_routine
def drawmap( values , levels , cmap , _save_filename , flag_show , name ,norm) :
        fig,ax = plt.subplots(figsize=(5,5))
        plt.subplots_adjust(left=0.02, right=0.98, top=0.98, bottom=0.02)
     
        m = Basemap(projection='stere', llcrnrlat=25, urcrnrlat=47, llcrnrlon=125, urcrnrlon=150, lat_0=60,  lon_0=140, resolution='i' )
        
        m.drawparallels(np.arange(-80.,81.,10.))
        m.drawmeridians(np.arange(-180.,181.,10.))
        m.drawcoastlines()
        m.drawmapboundary(fill_color='#eeeeee' , zorder=-1)
        m.fillcontinents(color='#ceceee' , zorder=-1)
        
        x , y = m(lons, lats)

        m.contourf( x , y , values , levels , cmap=cmap , zorder=0 , norm=norm )

        fig.savefig(_save_filename)
        if flag_show == 'yes' :
                plt.show()
        plt.close()

```

下図のような画像が作成できます(2019年9月8日15時UTCのデータ）。

![gsrd.v5.201909081500.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/af0ff51f-5504-368f-e021-e8301b91eb64.png)

##まとめ

GRIB2形式だけどpygribで扱えないファイルをnetCDF形式に変換してpythonから扱えるようにした方法を記載しました。