注意(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...
一応これで使える.もっと上手な取り回し方はありそう.