ecCodes学习笔记：GRIB2文件插值与裁剪

August 04, 2018 (最后修改: March 08, 2021)

使用 ecCodes 和 scipy 实现对 GRIB2 文件的裁剪与插值

ecCodes 提供 Python API，可以方便地解码 GRIB2 数据。
GRIB2 文件由多个消息构成，循环使用 codes_grib_new_from_file() 方法，逐个读取 GRIB2 消息。

with open(file_path, 'rb') as f:
    while 1:
        grib_message = eccodes.codes_grib_new_from_file(f)
        if grib_message is None:
            break

ecCodes 内部使用 NumPy 保存数据，支持 SciPy 技术栈。 scipy 的 interpolation 包提供多种插值方法，支持从一维到n维的各类数据。 GRIB2 数据可以看做二维数据，因此可以使用 scipy 实现裁剪与插值。

获取源数据信息

使用 ecCodes 获取 GRIB 数据的网格信息

left_lon = eccodes.codes_get(grib_message, 'longitudeOfFirstGridPointInDegrees')
right_lon = eccodes.codes_get(grib_message, 'longitudeOfLastGridPointInDegrees')
lon_step = eccodes.codes_get(grib_message, 'iDirectionIncrementInDegrees')
nx = eccodes.codes_get(grib_message, 'Ni')

top_lat = eccodes.codes_get(grib_message, 'latitudeOfFirstGridPointInDegrees')
bottom_lat = eccodes.codes_get(grib_message, 'latitudeOfLastGridPointInDegrees')
lat_step = eccodes.codes_get(grib_message, 'jDirectionIncrementInDegrees')
ny = eccodes.codes_get(grib_message, 'Nj')

orig_values = eccodes.codes_get_values(grib_message)

GRIB2 数据点保存顺序由 scanning mode 值决定。 NWPC 生成的 GRIB2 数据 scanning mode 都是 0，也就是经度从小到大，纬度从大到小，同纬度的点排在一起（按行排列）。下面只考虑该类型的数据。

gribdata 实现

scipy.interpolate.griddata(
    points, 
    values, 
    xi, 
    method='linear', 
    fill_value=nan,
    rescale=False
)

gribdata() 方法支持对任意纬度数据的插值，用法比较简单。
首先要准备插值需要的源数据点坐标和目标数据点坐标。

points

源数据每个点的坐标数据，按纬度分开，每个维度一个数组。

# TODO: meshgrid()...
lons = []
lats = []
for y in np.arange(top_lat, bottom_lat, -lat_step):
    for x in np.arange(left_lon, right_lon, lon_step):
        lons.append(x)
        lats.append(y)

目标数据每个点的坐标，每个坐标对应插值后的一个点。注意坐标范围要在源数据的范围内，超出范围的数据会使用 fill_value 填充。

points = []
for y in np.arange(top_lat, bottom_lat, -lat_step):
    for x in np.arange(left_lon, right_lon, lon_step):
       points.append([x, y])

下面调用 gribdata()，就可以返回插值和裁剪后的数据数组。

target_values = griddata(
  (orig_lons, orig_lats),
  orig_values, 
  target_points, 
  method='linear'
)

修改数据范围需要修改上一步提取的与数据范围相关的 key，并将数据更新为插值后的数据。

eccodes.codes_set(grib_message, 'longitudeOfFirstGridPointInDegrees', left_lon)
eccodes.codes_set(grib_message, 'longitudeOfLastGridPointInDegrees', right_lon)
eccodes.codes_set(grib_message, 'iDirectionIncrementInDegrees', lon_step)
eccodes.codes_set(grib_message, 'Ni', nx)

eccodes.codes_set(grib_message, 'latitudeOfFirstGridPointInDegrees', top_lat)
eccodes.codes_set(grib_message, 'latitudeOfLastGridPointInDegrees', bottom_lat)
eccodes.codes_set(grib_message, 'jDirectionIncrementInDegrees', lat_step)
eccodes.codes_set(grib_message, 'Nj', ny)

eccodes.codes_set_values(grib_message, target_values)
eccodes.codes_write(grib_message, output_file)

对 regular_ll 网格类型的优化

gribdata 针对任意的n维数据进行插值，所以效率不高。 GRIB2 数据通常都是二维网格数据，scipy 提供针对这类数据的插值类 RegularGridInterpolator。

class scipy.interpolate.RegularGridInterpolator(
    points, 
    values, 
    method='linear', 
    bounds_error=True, 
    fill_value=nan
)

这里的 points 和上面略有不同

points

每个纬度的坐标值：

y_array = np.arange(bottom_lat, top_lat, lat_step)
x_array = np.arange(left_lon, right_lon, lon_step)

values

原始 values 为一维数组，需要 reshape 为二维矩阵。
下面构建插值类，调用生成的对象得到插值后的结果。

target_function = RegularGridInterpolator(
    (orig_y_array, orig_x_array), 
    orig_values.reshape(ny, nx)
)

points = []
for y in np.arange(bottom_lat, top_lat, lat_step):
    for x in np.arange(left_lon, right_lon, lon_step):
        points.append([y, x])

target_values = target_function(target_xy_points)