代码收藏家技术教程 2024-08-05

Python读取grib数据获取变量推荐姿势

前情提要

最近使用的EC和GFS预报数据给的都是grib2格式的，之前用惯nc格式的，python读取grib2数据的时候还走了些弯路，看到很多博客上给的教程其实不能满足我的需求，现在搞明白了分享一下

pygrib安装

第一个问题就是我电脑上pygrib安装都折腾了一阵子
我电脑上直接pip是装不上去的

pip3 install pygrib

使用conda的时候最好使用的是conda-forge源的，用这个装成功了

conda install -c conda-forge pygrib

pygrib使用

遍历一下查看数据

以一个GFS的预报数据为例，我们使用pygrib读取之后，得到的是一个可迭代对象，可以拿来循环

import pygrib
grbs = pygrib.open('gfs.t00z.pgrb2.0p25.f024')
for grb in grbs:
    print(grb)

我们来看一下遍历出的grbmessage对象打印出来是什么

1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 24 hrs:from 202407210000
2:Cloud mixing ratio:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 24 hrs:from 202407210000
3:Ice water mixing ratio:kg kg**-1 (instant):regular_ll:hybrid:level 1:fcst time 24 hrs:from 202407210000

可以发现，每一个grbMessage对象可以是二维格点数据，同时单位，气压层等维度信息也是有的
以“510:Temperature:K (instant):regular_ll:isobaricInhPa:level 95000 Pa:fcst time 24 hrs:from 202407210000”为例，表示的意思如下

510: 这是记录的索引或序号。它标识这是文件中的第 510 个记录。
Temperature: 这是变量的名称。在这里，它表示“温度”。
K (instant): 这是变量的单位和时间类型。
K 表示开尔文，是温度的单位。
instant 表示这是瞬时值。
regular_ll: 这是网格类型，regular_ll 表示规则的纬度/经度网格。
isobaricInhPa: 这是变量所属的层。在这里，它表示等压面。
level 95000 Pa: 这是变量的层次信息。这里表示95000帕斯卡 (Pa) 等压面（950 hPa 等压面）。
fcst time 24 hrs: 这是预报时间。这里表示从初始时间起 24 小时后的预报。
from 202407210000: 这是初始时间，表示 2024 年 7 月 21 日 00:00。

获取变量

我下载的明明有41个高度层的信息，但是这样便利出来每一条数据都只是一个二维变量，我最终想获得的温度数据应该是三维数据，这意味着我需要将遍历的数据中为Temperature的调出来将他们拼成三维数组，这样才获取了一个变量完整的数据

可以看到，遍历出的数据以Temperature开头的就有56个，这是将非等压面的一些奇怪的Temperature也算上的总数

因此我们需要从中获取grbs中真正是气压层的Temperature

import pygrib
import numpy as np

def select_variable(filepath, variable, typeOfLevel):
    # 打开GRIB2文件
    grbs = pygrib.open(filepath)
    for grb in grbs:
        print(grb)
    temperature_records = grbs.select(name='Temperature', typeOfLevel='isobaricInhPa')
    # 获取记录数量
    num_records = len(temperature_records)
    
    # 假设所有记录的维度相同，获取第一个记录的维度
    # 这里使用第一个记录来确定数据的维度
    first_record = temperature_records[0].values
    data_shape = first_record.shape
    
    # 创建一个三维数组来存储温度数据
    select_data = np.zeros((num_records, *data_shape))
    
    # 提取每个记录的数据并存储到三维数组中
    for i, record in enumerate(temperature_records):
        select_data[i, :, :] = record.values
    return select_data


filepath = '../download_res/ncep/20240721/gfs.t00z.pgrb2.0p25.f024'
variable = 'Temperature'
typeOfLevel='isobaricInhPa'
T = select_variable(filepath, variable, typeOfLevel)
print(np.shape(T))

这个函数和案例就展示了如何从grib数据中获取需要的变量
这里需要解释的是typeOfLevel=‘isobaricInhPa’，这是因为这个数据中只要是有不同高度层的数据，他的typeOfLevel就是isobaricInhPa，因为除了高度层的Temperature可以看到这种不是描述等压面温度的Temperature，如果不加上typeOfLevel='isobaricInhPa就会把他们也筛选出来