注意（2023/08/04）

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

grib2とxarrayの相性問題

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

cfgribエンジンという革命

しかし，いつの間にかxarrayでcfgribエンジン(?)を選択することにより，gribを扱えるようになっていた．最近のMetpyMondaysでも詳細に説明されている．

import xarray as xr

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

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

なお，cfgribは独立したpythonパッケージであり，xarrayやpygrib同様，事前にインストールが必要．pipやconda forgeなど．

surfaceデータの落とし穴？

早速，JRA-55を読み込んでみる．まずは気圧面データを読み込んでみる．

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権限のあるディレクトリで実行すればエラーは出てこない）

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...

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

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の読み込みがスキップされている．

ds


<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が見当たらない．．

調べると，heightAboveGroundの要素が複数あり（気温，湿度などは2m，風は10m）そのことが原因で読み込みに失敗する模様．

結局の解決はcfgribの公式ドキュメントへ辿り着き，キーワードでフィルタリングをすることで10m風の要素を取り出すことができた．'topLevel'を10することがポイント．

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...

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