grib2フォーマットの気象データをxarrayで読み込む方法と注意

### 注意（2023/08/04）

以下の方法は経度がズレるバグがある，と報告を受けました．確かにdsの出力ではlongitudeは-1.25から開始していますが，データ自体は0から開始しています．修正されたら更新します．


## grib2とxarrayの相性問題

気象データはgrib1やgrib2というフォーマットで提供されることが多い．しかし，MetpyMondaysで[過去に紹介された](https://www.youtube.com/watch?v=yLoudFv3hAY)ように，gribフォーマットとxarrayの相性が良くなく，pygribを使って複雑なループを回すか，wgribをインストールし利用するか（macとの相性が悪い），という苦行を強いられていた．

## cfgribエンジンという革命
しかし，いつの間にかxarrayでcfgribエンジン(?)を選択することにより，gribを扱えるようになっていた．最近のMetpyMondaysでも[詳細に説明されている](https://www.youtube.com/watch?v=LYWfD_sFpK0)．

```python3
import xarray as xr

ds = xr.open_dataset('file.grb', engine='cfgrib')
```



(個人的な話で恐縮だが，これにより現在扱うすべてのバイナリデータをxararyで扱えるようになり，精神が安定した．)

なお，cfgribは独立したpythonパッケージであり，xarrayやpygrib同様，事前にインストールが必要．[pip](https://pypi.org/project/cfgrib/0.8.4.5/)や[conda forge](https://anaconda.org/conda-forge/cfgrib)など．


## surfaceデータの落とし穴？

早速，[JRA-55](https://jra.kishou.go.jp/JRA-55/index_ja.html)を読み込んでみる．まずは気圧面データを読み込んでみる．

```python3
ds = xr.open_dataset('/DATA/JRA-55/Daily/anl_p125/2010/201001/anl_p125_hgt.2010010100', engine='cfgrib')
```

```
Can't create file '/DATA/JRA-55/Daily/anl_p125/2010/201001/anl_p125_hgt.2010010100.923a8.idx'
Traceback (most recent call last):
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/messages.py", line 534, in from_indexpath_or_filestream
    with compat_create_exclusive(indexpath) as new_index_file:
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/contextlib.py", line 119, in __enter__
    return next(self.gen)
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/messages.py", line 500, in compat_create_exclusive
    fd = os.open(path, os.O_WRONLY | os.O_CREAT | os.O_EXCL)
PermissionError: [Errno 13] Permission denied: '/DATA/JRA-55/Daily/anl_p125/2010/201001/anl_p125_hgt.2010010100.923a8.idx'
Can't read index file '/DATA/JRA-55/Daily/anl_p125/2010/201001/anl_p125_hgt.2010010100.923a8.idx'
Traceback (most recent call last):
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/messages.py", line 544, in from_indexpath_or_filestream
    index_mtime = os.path.getmtime(indexpath)
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/genericpath.py", line 55, in getmtime
    return os.stat(filename).st_mtime
FileNotFoundError: [Errno 2] No such file or directory: '/DATA/JRA-55/Daily/anl_p125/2010/201001/anl_p125_hgt.2010010100.923a8.idx'
```

indexファイルを作れないというエラーが返ってくるが，読み込めている．（write権限のあるディレクトリで実行すればエラーは出てこない）

```python3
ds
```

```
<xarray.Dataset>
Dimensions:        (isobaricInhPa: 37, latitude: 145, longitude: 288)
Coordinates:
    number         int64 ...
    time           datetime64[ns] ...
    step           timedelta64[ns] ...
  * isobaricInhPa  (isobaricInhPa) float64 1e+03 975.0 950.0 ... 3.0 2.0 1.0
  * latitude       (latitude) float64 90.0 88.75 87.5 ... -87.5 -88.75 -90.0
  * longitude      (longitude) float64 -1.25 0.0 1.25 2.5 ... 355.0 356.2 357.5
    valid_time     datetime64[ns] ...
Data variables:
    gh             (isobaricInhPa, latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            1
    GRIB_centre:             rjtd
    GRIB_centreDescription:  Japanese Meteorological Agency - Tokyo
    GRIB_subCentre:          241
    Conventions:             CF-1.7
    institution:             Japanese Meteorological Agency - Tokyo
    history:                 2023-02-25T23:22 GRIB to CDM+CF via cfgrib-0.9.1...
```

ところが，地表面データを読み込む際は奇妙なことが起こる．

```python3
ds = xr.open_dataset('/DATA/JRA-55/Daily/anl_surf125/2010/201001/anl_surf125.2010010100', engine='cfgrib')
```

```
Can't create file '/DATA/JRA-55/Daily/anl_surf125/2010/201001/anl_surf125.2010010100.923a8.idx'
Traceback (most recent call last):
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/messages.py", line 534, in from_indexpath_or_filestream
    with compat_create_exclusive(indexpath) as new_index_file:
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/contextlib.py", line 119, in __enter__
    return next(self.gen)
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/messages.py", line 500, in compat_create_exclusive
    fd = os.open(path, os.O_WRONLY | os.O_CREAT | os.O_EXCL)
PermissionError: [Errno 13] Permission denied: '/DATA/JRA-55/Daily/anl_surf125/2010/201001/anl_surf125.2010010100.923a8.idx'
Can't read index file '/DATA/JRA-55/Daily/anl_surf125/2010/201001/anl_surf125.2010010100.923a8.idx'
Traceback (most recent call last):
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/messages.py", line 544, in from_indexpath_or_filestream
    index_mtime = os.path.getmtime(indexpath)
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/genericpath.py", line 55, in getmtime
    return os.stat(filename).st_mtime
FileNotFoundError: [Errno 2] No such file or directory: '/DATA/JRA-55/Daily/anl_surf125/2010/201001/anl_surf125.2010010100.923a8.idx'
skipping variable: paramId==165 shortName='u10'
Traceback (most recent call last):
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/dataset.py", line 676, in build_dataset_components
    dict_merge(variables, coord_vars)
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/dataset.py", line 607, in dict_merge
    raise DatasetBuildError(
cfgrib.dataset.DatasetBuildError: key present and new value is different: key='heightAboveGround' value=Variable(dimensions=(), data=2.0) new_value=Variable(dimensions=(), data=10.0)
skipping variable: paramId==166 shortName='v10'
Traceback (most recent call last):
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/dataset.py", line 676, in build_dataset_components
    dict_merge(variables, coord_vars)
  File "/mnt/data1/kasuga/miniconda/lib/python3.9/site-packages/cfgrib/dataset.py", line 607, in dict_merge
    raise DatasetBuildError(
cfgrib.dataset.DatasetBuildError: key present and new value is different: key='heightAboveGround' value=Variable(dimensions=(), data=2.0) new_value=Variable(dimensions=(), data=10.0)
```

先ほどと同じエラーも出ているが，よく読むと`u10`と`v10`の読み込みがスキップされている．

```python3
ds
```

```python3

<xarray.Dataset>
Dimensions:            (latitude: 145, longitude: 288)
Coordinates:
    number             int64 ...
    time               datetime64[ns] ...
    step               timedelta64[ns] ...
    surface            float64 ...
  * latitude           (latitude) float64 90.0 88.75 87.5 ... -87.5 -88.75 -90.0
  * longitude          (longitude) float64 -1.25 0.0 1.25 ... 355.0 356.2 357.5
    valid_time         datetime64[ns] ...
    meanSea            float64 ...
    heightAboveGround  float64 ...
Data variables:
    sp                 (latitude, longitude) float32 ...
    msl                (latitude, longitude) float32 ...
    t2m                (latitude, longitude) float32 ...
    pt                 (latitude, longitude) float32 ...
    depr               (latitude, longitude) float32 ...
    sh2                (latitude, longitude) float32 ...
    r2                 (latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            1
    GRIB_centre:             rjtd
    GRIB_centreDescription:  Japanese Meteorological Agency - Tokyo
    GRIB_subCentre:          241
    Conventions:             CF-1.7
    institution:             Japanese Meteorological Agency - Tokyo
    history:                 2023-02-25T23:26 GRIB to CDM+CF via cfgrib-0.9.1...
```

そして，確かに，`Data variables`に地表10m西風`u10`と北風`v10`が見当たらない．．

[調べる](https://github.com/ecmwf/cfgrib/issues/245)と，`heightAboveGround`の要素が複数あり（気温，湿度などは2m，風は10m）そのことが原因で読み込みに失敗する模様．

結局の解決は[cfgribの公式ドキュメント](https://pypi.org/project/cfgrib/0.8.4.5/)へ辿り着き，キーワードでフィルタリングをすることで10m風の要素を取り出すことができた．`'topLevel'`を`10`することがポイント．

```python3
ds = xr.open_dataset(f, engine='cfgrib',
                     backend_kwargs={'filter_by_keys': {'typeOfLevel': 'heightAboveGround', 'topLevel': 10}})
```

```
<xarray.Dataset>
Dimensions:            (latitude: 145, longitude: 288)
Coordinates:
    number             int64 ...
    time               datetime64[ns] ...
    step               timedelta64[ns] ...
    heightAboveGround  float64 ...
  * latitude           (latitude) float64 90.0 88.75 87.5 ... -87.5 -88.75 -90.0
  * longitude          (longitude) float64 0.0 1.25 2.5 ... 356.2 357.5 358.8
    valid_time         datetime64[ns] ...
Data variables:
    u10                (latitude, longitude) float32 ...
    v10                (latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            1
    GRIB_centre:             rjtd
    GRIB_centreDescription:  Japanese Meteorological Agency - Tokyo
    GRIB_subCentre:          241
    Conventions:             CF-1.7
    institution:             Japanese Meteorological Agency - Tokyo
    history:                 2023-02-25T23:39 GRIB to CDM+CF via cfgrib-0.9.1...
```

一応これで使える．もっと上手な取り回し方はありそう．