Python处理气象信息grib，grib2文件

Posted 2023-01-09 程序媛一枚~

写这篇博客源于博友的提问，

1. 效果图如下：

2. 安装

pip install pygrib
pip install gributils

GRIB 是世界气象组织 （World Meterological Organization WMO） 用于分发网格数据的标准。该模块包含一个python程序员接口，用于使用ECWMF GRIB API C库读取和写入GRIB网格（版本1和2）。

# 导入包
import pygrib

# 打开grib文件，并创建grib消息迭代器
grbs = pygrib.open('sampledata/flux.grb')

# 有seek, tell, read, readline and close方法，不同于文件的区别是，文件以字节码跳跃，而这里以grib消息跳跃
grbs.seek(2)

grbs.tell()

# 返回一个list有接下来的 N条消息，这里是N=1条
grb = grbs.read(1)[0]  # read returns a list with the next N (N=1 in this case) messages.

# 打印grib对象消息的概括信息
print(grb)

grbs.tell()

# 打印文件清单
grbs.seek(0)

for grb in grbs:
    print(grb)
# 1:Precipitation rate:kg m**-2 s**-1 (avg):regular_gg:surface:level 0:fcst time 108-120 hrs (avg):from 200402291200
# 2:Surface pressure:Pa (instant):regular_gg:surface:level 0:fcst time 120 hrs:from 200402291200
# 3:Maximum temperature:K (instant):regular_gg:heightAboveGround:level 2 m:fcst time 108-120 hrs:from 200402291200
# 4:Minimum temperature:K (instant):regular_gg:heightAboveGround:level 2 m:fcst time 108-120 hrs:from 200402291200


# 查找具有匹配名称的第一条 grib 消息：
grb = grbs.select(name='Maximum temperature')[0]

# 使用提取数据值values,（grb.keys（） 将返回可用键的列表）：
maxt = grb.values  # same as grb['values']

# The data is returned as a numpy array, or if missing values or a bitmap
# are present, a numpy masked array.  Reduced lat/lon or gaussian grid
# data is automatically expanded to a regular grid. Details of the internal
# representation of the grib data (such as the scanning mode) are handled
# automatically.
# 数据以numpy array格式返回，如果缺少值或位图，则返回一个 numpy masked array。
# 缺少lat纬度/lon经度或高斯网格数据自动扩展为常规网格。grib数据的内部细节的表示（如扫描模式）会自动处理。
print(maxt.shape, maxt.min(), maxt.max())
# (94, 192) 223.7 319.9

# 获取网格的经纬度
lats, lons = grb.latlons()

print(lats.shape, lats.min(), lats.max(), lons.shape, lons.min(), lons.max())
# (94, 192) -88.5419501373 88.5419501373  0.0 358.125

# 获取第2条grib信息
grb = grbs.message(2)  # same as grbs.seek(1); grb=grbs.readline()

print(grb)
# 2:Surface pressure:Pa (instant):regular_gg:surface:level 0:fcst time 120 hrs:from 200402291200

# 提取经纬度子集，如北美范围的经纬度
data, lats, lons = grb.data(lat1=20, lat2=70, lon1=220, lon2=320)

print(data.shape, lats.min(), lats.max(), lons.min(), lons.max())
# (26, 53) 21.904439458 69.5216630593 221.25 318.75

# 根据key重赋值
grb['forecastTime'] = 240
grb.dataDate = 20100101

# 获取与编码消息关联的二进制字符串：
msg = grb.tostring()

# 关闭grib文件
grbs.close()

# 写入新的grib文件
grbout = open('test.grb', 'wb')
grbout.write(msg)
grbout.close()

print(pygrib.open('test.grb').readline())
# 1:Surface pressure:Pa (instant):regular_gg:surface:level 0:fcst time 240 hrs:from 201001011200

参考

• https://github.com/seminar2012/pygrib
• https://jswhit.github.io/pygrib/api.html#example-usage
• https://www.heywhale.com/mw/project/61d8592d9d17520017a8d815
• https://cloud.tencent.com/developer/article/1618335
• https://www.bilibili.com/read/cv17061670