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

python读取、使用grib文件，对grib数据进行插值

简介：

由于实际工作需要，需要对grib文件中的气象数据进行插值处理。了解到在python中有pygrib库可以提取文件，且根据EARA提供的气象数据在其官网上可以直接利用python代码提取。（cds气象数据提取地为Copernicus Climate Data Store | Copernicus Climate Data Store）于是决定使用python来进行插值处理。（pygrib文件调用api说明网址为：API — pygrib documentation）

1.pygrib提取grib气象数据

从cds网站上爬取数据需要在电脑上进行cdsapi处理，这部分可以在部分博主的相关博客，一下为提取对流降雨数据的示例代码，运行后可以得到一个grib文件。

import cdsapi

c = cdsapi.Client()

c.retrieve(
    'reanalysis-era5-single-levels',
    {
        'product_type': 'reanalysis',
        'format': 'grib',
        'variable': [
            'convective_precipitation',
        ],
        'year': '2019',
        'month': [
            '01', '02', '03',
        ],
        'day': [
            '01', '02', '03',
            '04', '05', '06',
            '07', '08', '09',
            '10', '11', '12',
            '13', '14', '15',
            '16', '17', '18',
            '19', '20', '21',
            '22', '23', '24',
            '25', '26', '27',
            '28', '29', '30',
            '31',
        ],
        'time': [
            '00:00', '01:00', '02:00',
            '03:00', '04:00', '05:00',
            '06:00', '07:00', '08:00',
            '09:00', '10:00', '11:00',
            '12:00', '13:00', '14:00',
            '15:00', '16:00', '17:00',
            '18:00', '19:00', '20:00',
            '21:00', '22:00', '23:00',
        ],
    },
    'download.grib')

2.pygrib库读取grib文件

考虑到grib文件含有经纬度信息、时间信息及变量数据是一个多维数据，同时在python中NPZ文件是一种用于存储和加载多个NumPy数组的文件格式。NPZ文件是NumPy库提供的一种压缩文件格式，可以将多个数组保存在一个单独的文件中。NPZ文件是通过使用NumPy的savez函数创建的，该函数接受一个文件名参数和要保存的多个数组作为输入。它将这些数组保存在一个压缩的NPZ文件中。NPZ文件可以通过NumPy的load函数加载。加载NPZ文件时，会返回一个类似于字典的对象，其中包含已保存的数组，并且可以通过键来访问每个数组。因此决定将由pygrib提取出的数据保存在npz格式文件中。

import pygrib
import numpy as np
from tqdm import tqdm

# 定义要读取的GRIB文件路径
file_path = 'ERA5_data.grib'

# 使用pygrib打开GRIB文件
grbs = pygrib.open(file_path)

# 创建一个空的字典来存储数据
data_dict = {}
timestamps = []

# 创建一个空的三维数组，用于存储数据
data_array = np.empty((205, 253, 5989))

# 获取消息记录总数
num_records = grbs.messages

# 读取每个消息记录
for i, grb in enumerate(tqdm(grbs, total=num_records, desc='Updating Data')):
    name = grb.name

    # 提取数据和经纬度
    data, lats, lons = grb.data(lat1=3, lat2=54, lon1=73, lon2=136)

    # 获取有效时间
    validity_date = grb.validityDate
    validity_time = grb.validityTime
    timestamp = f"{validity_date} {validity_time}"
    timestamps.append(timestamp)

    # 将数据存储为NumPy数组
    data_array[:, :, i] = data

# 保存data_array为npz文件
data_dict = {
    'data_array': data_array,
    'timestamps': timestamps
}
np.savez('data_array.npz', **data_dict)

在这里提取时间时可以根据自身对时间格式的需要将时间修改后存入npz文件中。由于所需要的数据生成的grib文件数据量较大，使用tqdm来显示循环的进度即显示代码运行进度。

3.插值处理

3.1数据准备

目前需要在两段时间中进行气象数据的线性插值（cds数据为每一小时预测），需要考虑一种格式将数据很好的保存下来方便提取进行线性插值，由于气象数据为在经纬度网格上平面展开的二维模式，当加入时间后便可以理解为一个多维数组，层数为时间于是可以有：（这里将时间单独保存了，将时间理解为了指示标记作用）

import numpy as np
from tqdm import tqdm

#加载保存的字典文件
data = np.load('data_array.npz', allow_pickle=True)

# 从 .npz 文件中提取数据数组
data_array = data['data_array']
timestamps = data['timestamps']
connective_precipitation = np.empty((205, 253, 1498))
# 该多维数组的形状大小根据grib文件数据格式决定

# 创建字典存储时间戳数据
t1 = []
for j in range(1498):
    connective_precipitation[:,:,j] = data_array[:,:,j]
    t1.append(timestamps[j])

# 保存变量为npz文件
data_dict = {
    'connective_precipitation': connective_precipitation,
    't1': t1
}
np.savez('connective_precipitation.npz', **data_dict)

3.2建立插值函数

对每段时间采用15分钟插值一次，如上所说时间作为标记作用因此没有使用strftime来单独进行时间格式转换及时间计算，于是在插值时将时间作为了百进制。同时因为是简单的线性插值选择了构建简单易于理解的一次函数。代码如下：

import numpy as np
from tqdm import tqdm

data = np.load('connective_precipitation.npz', allow_pickle=True)

# 从字典中提取数据数组和时间戳数组
connective_precipitation = data['connective_precipitation']
t1 = data['t1']

#建立线性插值函数
def interpolate_arrays(arr_list):
    # 确保输入数组的形状相同
    assert len(set(arr.shape for arr in arr_list)) == 1, "输入数组的形状不相同"

    # 获取数组的形状和数组数量
    rows, cols = arr_list[0].shape

    # 创建一个空数组，用于存储插值结果
    interpolated_segments = np.empty((rows, cols, 3))

    # 对每个元素进行分段线性插值
    for k in range(3):
        interpolated_segment = (arr_list[1] - arr_list[0]) * 0.25 * (k + 1) + arr_list[0]
        interpolated_segments[:, :, k] = interpolated_segment

    # 返回插值结果数组
    return interpolated_segments

# 数组拼接
cp_intered = np.empty((   ,     ,    ))  # 数组大小根据需要计算
for i in tqdm(range(1497), desc='Processing'):
    arr_list = [connective_precipitation[:, :, i], connective_precipitation[:, :, i + 1]]
    interpolated_arr = interpolate_arrays(arr_list)

    for j in range(0,   , 4):
        cp_intered[:, :, j] = connective_precipitation[:, :, i]
        cp_intered[:, :, j + 1:j + 4] = interpolated_arr

    cp_intered[:, :,     ] = connective_precipitation[:, :, 1497]

print(cp_intered.shape)
# 保存插值后的数据
np.savez('cp_intered.npz', cp_intered=cp_intered)

如对时间显示确有需要，可以再次前提上多加入对时间格式的处理，正如前所述我这里的时间仅作为了标记作用使用，对上述代码中的时间格式处理可以使用：

        # 正则表达式模式，匹配数字和字符串部分
        pattern = r'(\d+)\s(.*)'

        # 提取数字和字符串部分
        matches = [re.match(pattern, ti)]

        # 分别获取数字和字符串部分
        numbers = [int(matches[0].group(2))]
        strings = [matches[0].group(1)]

        # 加法操作
        numbers_plus_15 = np.array(numbers) + 15 * (k + 1)

如对显示的格式不满意可以进一步选择datatime库的strptime或strftime来对时间格式进行修改。