python读取函数 python读取grib

前面已经写过如何在windows下安装wgrib，见这个博文 实际应用中，发现grib2文件的读法还是有挺多需要细写的，这里以国家气象信息中心的三维云实况产品（Z_NAFP_C_BABJ_20210830203447_P_3DCloudA_RT_CHN_0P05_HOR-CCP3-2021083020.GRB2）为例。

方法一：用wgrib将grib2格式转换为csv格式读取

查看文件信息，在命令提示符下或者anaconda下，输入命令

wgrib2 -v0 Z_NAFP_C_BABJ_20210830203447_P_3DCloudA_RT_CHN_0P05_HOR-CCP3-2021083020.GRB2

得到结果

从结果看，共有43层，我们需要的数据名称是CDCC。在python下输入命令

import os
file_grb2 = 'Z_NAFP_C_BABJ_20210830203447_P_3DCloudA_RT_CHN_0P05_HOR-CCP3-2021083020.GRB2'
file_csv=file_grb2+'.csv'
data_match = 'CDCC'
cmd1 = 'wgrib2 %s -match \'%s\' -csv %s' % (file_grb2,data_match,file_csv)
success = os.system(cmd1)

得到的csv文件很大，格式是这样的

可以用读csv文件的方法来读取

import pandas as pd
data = pd.read_csv(file_csv, sep=',',names = ['time0','time1','data_type','height','lon','lat','value'],lineterminator="\n" )

读取到的数据是dataframe格式的。

dataframe格式筛选和取值就很方便。

value = data[(data['height']=='50 mb') & (data['lon']==70) & (data['lat']==60)]
print(value)

但是因为grib2文件本身就大，转换为csv文件后就更大了，是原grib2文件的几十倍大。所以，文件转换和数据读取的耗时都很长，还容易死机。消耗的计算时间和存储空间都很大，不建议使用。

方法二：用wgrib将grib2格式转换为nc格式读取

file_nc = file_grb2+'.nc'
cmd1 = 'wgrib2 %s -netcdf %s' % (file_grb2,file_nc)
success = os.system(cmd1)

得到的nc文件比grib2文件大一倍多。

然后开始读取nc文件

import netCDF4
from netCDF4 import Dataset
nc_obj=Dataset(file_nc)
print(nc_obj)

nc_obj的输出结果是这样的

可以看到nc_obj有很多属性，可以输出查看

print(nc_obj.dimensions)
print(nc_obj.variables)

nc_obj.属性是dict字典形式的。

如果要取用云量数据，可以variables属性里用字典取用的方式来取。

print(nc_obj.variables['CDCC_50mb'])

这样就取到了50mb高度上的云量

结果是（time,latitude,longitude）三个维度的，其中时间只有一维。将经纬度值转换为格点序号，这里取用纬度序号为1，经度序号为2。

print(nc_obj.variables['CDCC_50mb'][0][1][2])

这样就取到了某个高度、某个纬度、某个经度的云量值。

方法三：用xarray直接读取grib2格式

import xarray as xr
import cfgrib
ds1 = xr.open_dataset(file_grb2,engine='cfgrib')
print(ds1)

结果如下

那么， 提出变量ccl

print(ds1.variables['ccl'])

结果是

发现它是一个三维array，维度是（高度，纬度，经度），计算出高度、纬度、经度的序号，并提取数值。

print(ds1.variables['ccl'][0][1][2].values)