GRIB笔记:使用cfgrib加载GRIB文件

目录

cfgrib是 ECMWF 开发的 GRIB Python接口,支持Unidata’s Common Data Model v4,符合CF Conventions。 高层API接口为 xarray 提供GRIB解码引擎。 底层访问和解码由 ECMWF 的ecCodes库实现。

功能

cfgrib正在开发中,处于Beta版本的功能有:

  • 支持xarray使用engine="cfgrib"读取GRIB文件。
  • 使用cfgrib.open_datasets能读取大部分GRIB 1和2文件,包括包含不同层次类型的文件。
  • 支持Linux、MacOS和Windows,唯一的依赖是ecCodes的C库。
  • 所有支持的平台都可以使用conda-forge包安装。
  • 延迟和高效读取数据,节省内存占用和磁盘访问。
  • 允许使用dask进行大于内存的分布式处理。
  • 支持将坐标转换为不同的数据模型和命名约定。
  • 支持将GRIB文件的索引写入磁盘,以在打开时保存全文件扫描。

安装

推荐使用conda安装。

conda install -c conda-forge eccodes
conda install -c conda-forge cfgrib

如果使用中国气象局的CMA-PI高性能计算机,则需要手动安装。

推荐使用apps/anaconda3/2019.10环境。从PyPi网站中下载cfgrib,attrs和cffi三个包的wheel文件, 将这三个包安装到本地用户目录。 例如使用下面的命令安装cfgrib的预编译包。

pip install cfgrib-0.9.8.1-py2.py3-none-any.whl --user

准备

从 ECMWF 官网下载 ERA5 的示例文件:

wget http://download.ecmwf.int/test-data/cfgrib/era5-levels-members.grib

首先载入需要用到的一些库。

import xarray as xr
import cfgrib

读取示例文件

读取示例文件,返回xarray.Dataset

data_set = xr.open_dataset(
    "era5-levels-members.grib",
    engine="cfgrib",
)
data_set
<xarray.Dataset>
Dimensions:        (isobaricInhPa: 2, latitude: 61, longitude: 120, number: 10, time: 4)
Coordinates:
  * number         (number) int64 0 1 2 3 4 5 6 7 8 9
  * time           (time) datetime64[ns] 2017-01-01 ... 2017-01-02T12:00:00
    step           timedelta64[ns] ...
  * isobaricInhPa  (isobaricInhPa) int64 850 500
  * latitude       (latitude) float64 90.0 87.0 84.0 81.0 ... -84.0 -87.0 -90.0
  * longitude      (longitude) float64 0.0 3.0 6.0 9.0 ... 351.0 354.0 357.0
    valid_time     (time) datetime64[ns] ...
Data variables:
    z              (number, time, isobaricInhPa, latitude, longitude) float32 ...
    t              (number, time, isobaricInhPa, latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            1
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2020-03-16T13:44:28 GRIB to CDM+CF via cfgrib-0....

可以通过设置backend_kwargs中的read_keys,增加需要读取的 GRIB key。

data_set = xr.open_dataset(
    "era5-levels-members.grib", 
    engine="cfgrib",
    backend_kwargs={
        "read_keys": [
            "experimentVersionNumber",
        ],
    })
data_set.t.attrs['GRIB_experimentVersionNumber']
'0001'

过滤GRIB文件

xr.open_dataset仅能读取有同样shortName且能放到同一个超立方体(hypercube)中,即同样的要素、相同的层次类型。

GRAPES GFS模式每个时效生成一个单一的GRIB2文件,包含所有要素场。

下面代码展示使用xr.open_dataset直接读取GRIB 2文件会抛出异常。 cfgrib 无法有效生成 GRAPES 模式的索引,所以将backend_kwargs中的indexpath设为空字符串,关掉 GRIB 文件索引功能。

data_set = xr.open_dataset(
    "gmf.gra.2020031500000.grb2",
    engine="cfgrib",
    backend_kwargs={
        "indexpath": "",
    },
)
data_set
...
DatasetBuildError: multiple values for unique key, try re-open the file with one of:
    filter_by_keys={'typeOfLevel': 'surface'}
    filter_by_keys={'typeOfLevel': 'nominalTop'}
    filter_by_keys={'typeOfLevel': 'heightAboveGround'}
    filter_by_keys={'typeOfLevel': 'cloudBase'}
    filter_by_keys={'typeOfLevel': 'atmosphere'}
    filter_by_keys={'typeOfLevel': 'meanSea'}
    filter_by_keys={'typeOfLevel': 'isobaricInhPa'}
    filter_by_keys={'typeOfLevel': 'isobaricInPa'}
    filter_by_keys={'typeOfLevel': 'depthBelowLandLayer'}
    filter_by_keys={'typeOfLevel': 'heightAboveGroundLayer'}

backend_kwargs中的filter_by_keys中使用 GRIB Key 设置过滤条件。

下面的代码筛选文件中等压面层的温度场。

data_set = xr.open_dataset(
    "gmf.gra.2020031500000.grb2",
    engine="cfgrib",
    backend_kwargs={
        "indexpath": "",
        "filter_by_keys": {
            "shortName": "t",
            "typeOfLevel": "isobaricInhPa"
        },
    },
)
data_set
<xarray.Dataset>
Dimensions:        (isobaricInhPa: 36, latitude: 720, longitude: 1440)
Coordinates:
    time           datetime64[ns] ...
    step           timedelta64[ns] ...
  * isobaricInhPa  (isobaricInhPa) int64 1000 975 950 925 900 850 ... 5 4 3 2 1
  * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    valid_time     datetime64[ns] ...
Data variables:
    t              (isobaricInhPa, latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             babj
    GRIB_centreDescription:  Beijing 
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             Beijing 
    history:                 2020-03-16T06:26:07 GRIB to CDM+CF via cfgrib-0....

数据集中t变量就是包含36个层次的温度场。 选择850hPa的温度场。

data_set.t.sel(isobaricInhPa=850)
<xarray.DataArray 't' (latitude: 720, longitude: 1440)>
[1036800 values with dtype=float32]
Coordinates:
    time           datetime64[ns] ...
    step           timedelta64[ns] ...
    isobaricInhPa  int64 850
  * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    valid_time     datetime64[ns] ...
Attributes:
    GRIB_paramId:                             130
    GRIB_shortName:                           t
    GRIB_units:                               K
    GRIB_name:                                Temperature
    GRIB_cfName:                              air_temperature
    GRIB_cfVarName:                           t
    GRIB_dataType:                            an
    GRIB_missingValue:                        9999
    GRIB_numberOfPoints:                      1036800
    GRIB_typeOfLevel:                         isobaricInhPa
    GRIB_NV:                                  0
    GRIB_stepUnits:                           1
    GRIB_stepType:                            instant
    GRIB_gridType:                            regular_ll
    GRIB_gridDefinitionDescription:           Latitude/longitude 
    GRIB_Nx:                                  1440
    GRIB_iDirectionIncrementInDegrees:        0.25
    GRIB_iScansNegatively:                    0
    GRIB_longitudeOfFirstGridPointInDegrees:  0.0
    GRIB_longitudeOfLastGridPointInDegrees:   359.75
    GRIB_Ny:                                  720
    GRIB_jDirectionIncrementInDegrees:        0.25
    GRIB_jPointsAreConsecutive:               0
    GRIB_jScansPositively:                    0
    GRIB_latitudeOfFirstGridPointInDegrees:   89.875
    GRIB_latitudeOfLastGridPointInDegrees:    -89.875
    long_name:                                Temperature
    units:                                    K
    standard_name:                            air_temperature

使用xarray自带的绘图功能。

data_set.t.sel(isobaricInhPa=850).plot()

自动过滤

cfgrib提供open_datasets函数用于自动选择合适的filter_by_keys并返回所有有效的xarray.Dateset列表。

下面代码载入整个时效的文件,运行速度很慢。

data_set = cfgrib.open_datasets(
    "gmf.gra.2020031500000.grb2",
    engine="cfgrib",
    backend_kwargs={
        "indexpath": "",
    },
)
data_set

截取一段返回结果,可以看到返回结果中将层次类型为isobaricInhPa且有30个层次的变量放到同一个xarray.Dataset中。

[
  ...
 <xarray.Dataset>
 Dimensions:        (isobaricInhPa: 30, latitude: 720, longitude: 1440)
 Coordinates:
     time           datetime64[ns] 2020-03-15
     step           timedelta64[ns] 00:00:00
   * isobaricInhPa  (isobaricInhPa) int64 1000 975 950 925 900 ... 70 50 30 20 10
   * latitude       (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
   * longitude      (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
     valid_time     datetime64[ns] 2020-03-15
 Data variables:
     paramId_0      (isobaricInhPa, latitude, longitude) float32 ...
     q              (isobaricInhPa, latitude, longitude) float32 ...
     vo             (isobaricInhPa, latitude, longitude) float32 ...
     d              (isobaricInhPa, latitude, longitude) float32 ...
     r              (isobaricInhPa, latitude, longitude) float32 ...
     papt           (isobaricInhPa, latitude, longitude) float32 ...
     dpt            (isobaricInhPa, latitude, longitude) float32 ...
     clwmr          (isobaricInhPa, latitude, longitude) float32 ...
     icmr           (isobaricInhPa, latitude, longitude) float32 ...
     rwmr           (isobaricInhPa, latitude, longitude) float32 ...
     snmr           (isobaricInhPa, latitude, longitude) float32 ...
     grle           (isobaricInhPa, latitude, longitude) float32 ...
     wz             (isobaricInhPa, latitude, longitude) float32 ...
     ccl            (isobaricInhPa, latitude, longitude) float32 ...
 Attributes:
     GRIB_edition:            2
     GRIB_centre:             babj
     GRIB_centreDescription:  Beijing 
     GRIB_subCentre:          0
     Conventions:             CF-1.7
     institution:             Beijing ,
  ...
]

高级特性

cfgrib的engine支持xarray的所有只读特性,例如:

  • 使用xarray.open_mddataset将多个GRIB文件合并到一个单一的dataset。
  • 使用dask处理大于内存的数据集。
  • 使用dask.distributed进行分布式处理。

后续会研究如何使用这些特性。

另外cfgrib还支持写入GRIB文件等特性。

参考

ecmwf/cfgrib

Using cfgrib to load GRIB data into an xarray dataset