ecCodes学习笔记：如何关掉ecCodes的常量场特性

June 05, 2023 (最后修改: June 05, 2023)

简述

设置环境变量 ECCODES_GRIB_LARGE_CONSTANT_FIELDS 可以关掉 ecCodes 的常量场 (constant field) 特性。

export ECCODES_GRIB_LARGE_CONSTANT_FIELDS=1

背景

CMA-GFS 提供 GRIB2 全球预报产品 (grib2-orig)，为了降低数据传输量，从全球预报产品中抽取四分之一区域生成 GRIB2 东北半球产品 (grib2-ne)。在 CMA-GFS V4.0 升级前 grib2-ne 使用 CEMC 内部开发的 grib_post.exe 制作，但 V4.0 分辨率从 25 KM 升到 12.5 KM，串行程序运行速度显著变慢。为了提高产品生成速度，实现了 python + dask + eccodes 的并行抽取版本，将产品制作时间从十多分钟降低为两分钟，详情见参考文章4。但使用 wgrib2 工具读取新算法生成的 grib2-ne 数据会出现警告信息，笔者检查数据发现有告警信息的要素场都是常量场。经过上网搜索发现问题出在 wgrib2 不支持常量编码方式，通过查看 ecCodes 源码最终发现通过设置环境变量 ECCODES_GRIB_LARGE_CONSTANT_FIELDS 可以关掉 ecCodes 的常量场特性。

本文首先介绍什么是常量场，并给出 CMA-GFS 常量场示例，然后介绍 wgrib2 和 g2ctl.pl 告警信息，并说明告警实际上不影响对数据的使用，最后通过分析源码介绍如何关掉 ecCodes 的常量场特性。

什么是常量场

常量场 (constant field) 是所有数据点值都相同的要素场。

WMO 技术手册在 GRIB2 规则章节中有如下注释（参考文章1）：

For grid-point values, complex packing features are intended to reduce the whole size of the GRIB message (data compression without loss of information with respect to simple packing).

说明复杂打包特性的目标是减少 GRIB 消息的整体大小。

手册在介绍 GRIB2 的数据段 (Section 7) 模板 Data template 7.2 – Grid point data – complex packing 时有如下注释：

For groups with a constant value, associated field width is 0, and no incremental data are physically present.

如果数据值都相同，该模板长度为 0，也就是说对于常量场，数据段中不需要定义模板。

ecCodes 在对常量场在进行编码时，通过仅存储一个数据值来显著减少消息大小（参考文章2）。

数据值存储在 referenceValue 中并编码为浮点数（32 位 IEEE 浮点数）
bitsPerValue 设置为 0，数据段设为空
解码时使用数据值填充数组，数据大小由 numberOfDataPoints 确定

下面以 CMA-GFS V4.0 的一个实际例子说明常量场。

常量场数据样例

以 2023 年 5 月 29 日 00 时次预报 024 小时数据为例说明，文件名 gmf.gra.2023052900024.grb2。

使用 grib_copy 抽取其中第 742 个要素场 (DEPR 20hPa)。

grib_copy -w count=742 gmf.gra.2023052900024.grb2 742.grib2

在 JupyterLab 中使用 cemc-oper/reki 库加载要素场并绘图

import matplotlib.pyplot as pyplot
from reki.format.grib.eccodes import load_field_from_file

orig_path = "some/path/to/742.grib2"
orig_field = load_field_from_file(orig_path, count=1)
orig_field.plot()

图形如下所示：

CMA-GFSV4.0 2023052900 时次 024 时效 20hPa DEPR 全球

可以看到，全球场不是常量场，但仅裁剪东北半球四分之一数据就是一个常量场。

同样的代码显示裁剪后的要素场：

CMA-GFSV4.0 2023052900 时次 024 时效 20hPa DEPR 东北半球

常量场数据对比

对比不同编码方式的 grib2-ne，其中

742.v1.grib2：使用常量场特性
742.v2.grib2：使用常规编码方式

使用 grib_dump -O 打印两个文件的信息，省略部分相同的输出信息，主要展示 SECTION 5 和 SECTION 7 的信息。

742.v1.grib2 常量场特性

***** FILE: 742.v1.grib2
#==============   MESSAGE 1 ( length=181 )                 ==============
1-4       identifier = GRIB
5-6       reserved = MISSING
7         discipline = 0 [Meteorological products (grib2/tables/4/0.0.table) ]
8         editionNumber = 2
9-16      totalLength = 181
...skip...
======================   SECTION_5 ( length=23, padding=0 )    ======================
1-4       section5Length = 23
5         numberOfSection = 5
6-9       numberOfValues = 259560
10-11     dataRepresentationTemplateNumber = 40 [JPEG2000 Packing (grib2/tables/4/5.0.table) ]
12-15     referenceValue = 30.0021
16-17     binaryScaleFactor = 0
18-19     decimalScaleFactor = 3
20        bitsPerValue = 0
21        typeOfOriginalFieldValues = 0 [Floating point (grib2/tables/4/5.1.table) ]
22        typeOfCompressionUsed = 0 [Lossless (grib2/tables/4/5.40.table) ]
23        targetCompressionRatio = 255
...skip...
======================   SECTION_7 ( length=5, padding=0 )     ======================
1-4       section7Length = 5
5         numberOfSection = 7
======================   SECTION_8 ( length=4, padding=0 )     ======================
1-4       7777 = 7777

742.v2.grib2 常规编码

***** FILE: 742.v2.grib2
#==============   MESSAGE 1 ( length=329 )                 ==============
1-4       identifier = GRIB
5-6       reserved = MISSING
7         discipline = 0 [Meteorological products (grib2/tables/4/0.0.table) ]
8         editionNumber = 2
9-16      totalLength = 329
...skip...
======================   SECTION_5 ( length=23, padding=0 )    ======================
1-4       section5Length = 23
5         numberOfSection = 5
6-9       numberOfValues = 259560
10-11     dataRepresentationTemplateNumber = 40 [JPEG2000 Packing (grib2/tables/4/5.0.table) ]
12-15     referenceValue = 30.0021
16-17     binaryScaleFactor = 0
18-19     decimalScaleFactor = 0
20        bitsPerValue = 15
21        typeOfOriginalFieldValues = 0 [Floating point (grib2/tables/4/5.1.table) ]
22        typeOfCompressionUsed = 0 [Lossless (grib2/tables/4/5.40.table) ]
23        targetCompressionRatio = 255
...skip...
======================   SECTION_7 ( length=153, padding=0 )   ======================
1-4       section7Length = 153
5         numberOfSection = 7
6-153     codedValues = (259560,148) {
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01,
3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01, 3.0002111435e+01
... 259460 more values
} # data_jpeg2000_packing codedValues
======================   SECTION_8 ( length=4, padding=0 )     ======================
1-4       7777 = 7777

使用 VSCode 对比上面两个输出结果，可以更明显看出两种编码方式的区别

启用常量场特性功能的 v1 数据（左）与常规编码的 v2 数据（右）grib_dump 命令输出结果对比

正如 ecCodes 文档描述，启用常量场特性的 v1 数据的数据段 (Section 7) 仅包含头信息，没有具体数据，数据描述段 (SECTION 5) 的 bitsPerValue 值为 0。而常规编码的 v2 数据的数据段有 JPEG 2000 压缩后的数据，bitsPerValue 值为 15。

值得注意的是，v2 数据虽然有 259560 个点，但压码后的数据段仅占 148 字节。这是因为数据段都是相同的数据，JPEG 压缩对这类数据的压缩效果较好。所以从 v1 到 v2 数据大小的增量较小（totalLength 值：182 vs 329）。

从概念上看，采用常量场特性的 v1 数据直接省略数据段，符合 WMO GRIB2 编码的标准，但业务上广泛应用的 wgrib2 工具对使用常量场特性的要素场会报警告。

问题

在展示问题之前需要首先加载 CMA-PI HPC 预安装的 WGRIB2 和 GrADS 软件包。

module load compiler/intel/composer_xe_2017.2.174
module load mathlib/wgrib2/2.0.6/intel
module mathlib/grads/2.2.0/gnu

wgrib2 告警

wgrib2 可以正常读取常规编码的 v2 文件：

wgrib2 -V 742.v2.grib2

1:0:vt=2023053000:20 mb:24 hour fcst:DEPR Dew Point Depression (or Deficit) [K]:
    ndata=259560:undef=0:mean=30.0021:min=30.0021:max=30.0021
    grid_template=0:winds(N/S):
        lat-lon grid:(721 x 360) units 1e-06 input WE:NS output WE:SN res 48
        lat 89.937500 to 0.187500 by 0.250000
        lon 0.000000 to 180.000000 by 0.250000 #points=259560

但 wgrib2 读取使用常量场特性的 v1 文件则会给出 Warning 报警。

wgrib2 -V 742.v1.grib2

Warning: g2lib/g2clib jpeg encode/deocde may differ from WMO standard, use -g2clib 0 for WMO standard
1:0:vt=2023053000:20 mb:24 hour fcst:DEPR Dew Point Depression (or Deficit) [K]:
    ndata=259560:undef=0:mean=30.0021:min=30.0021:max=30.0021
    grid_template=0:winds(N/S):
        lat-lon grid:(721 x 360) units 1e-06 input WE:NS output WE:SN res 48
        lat 89.937500 to 0.187500 by 0.250000
        lon 0.000000 to 180.000000 by 0.250000 #points=259560

告警原因原因详见参考文章3。

虽然给出了一个告警信息，但从 mean，min，max 三个信息看 wgrib2 还是能识别出 referenceValue 值作为数据值。

g2ctl.pl 告警

g2ctl.pl 是 GrADS 提供的命令行工具，可以为 grib2 文件生成 CTL 描述，方便在 GrADS 中加载 GRIB2 格式数据。

为 742.v1.grib2 生成 CTL 文件，将标准输出打印到 742.v1.ctl 中：

g2ctl.pl -verf 742.v1.grib2 > 742.v1.ctl

Warning: g2lib/g2clib jpeg encode/deocde may differ from WMO standard, use -g2clib 0 for WMO standard
Warning: g2lib/g2clib jpeg encode/deocde may differ from WMO standard, use -g2clib 0 for WMO standard

虽然命令行输出了两行 Warning 警告，但生成的 742.v1.ctl 中保存了数据的正确描述：

dset ^742.v1.grib2
index ^742.v1.grib2.idx
undef 9.999E+20
title 742.v1.grib2
* produced by g2ctl v0.0.8.8
* command line options: -verf 742.v1.grib2
* griddef=1:0:(721 x 360):grid_template=0:winds(N/S): lat-lon grid:(721 x 360) units 1e-06 input WE:NS output WE:SN res 48 lat 89.937500 to 0.187500 by 0.250000 lon 0.000000 to 180.000000 by 0.250000 #points=259560:winds(N/S)

dtype grib2
ydef 360 linear 0.187500 0.25
xdef 721 linear 0.000000 0.250000
tdef 1 linear 00Z30may2023 1mo
zdef 1 linear 1 1
vars 1
DEPR20mb   0,100,2000   0,0,7 ** 20 mb Dew Point Depression (or Deficit) [K]
ENDVARS

使用 gribmap 生成 GRIB2 索引文件 742.v1.grib2.idx：

gribmap -i 742.v1.ctl

使用 GrADS 画图：

ga-> open 742.v1.ctl
Scanning description file:  742.v1.ctl
Data file 742.v1.grib2 is open as file 1
LON set to 0 180
LAT set to 0.1875 89.9375
LEV set to 1 1
Time values set: 2023:5:30:0 2023:5:30:0
E set to 1 1
ga-> d DEPR20mb
Constant field.  Value = 30.0021

GrADS 绘制使用常量场特性的 v1 数据

可以看到 GrADS 可以加载 v1 版本的常量要素场。

虽然 wgrib2 和 g2ctl.pl 看起来都支持常量场特性，但有 warning 告警总是不太合适的，为了保持编码方式的一致性，还是要寻找关闭 ecCodes 常量场特性的方法。

解决方案

因为互联网上 ecCodes 相关资料较少，所以笔者首先通过研究源代码找到可行的解决方案，然后尝试不同的解决方案是否有效。

源码分析

笔者查看了 ecCodes 2.30.0 的源代码。(ecCodes的源代码非常难理解，用 C 语言模拟 C++ 对象继承！)

在 grib_accessor_class_data_simple_packing.cc 文件的 pack_double 函数中有如下代码与常量场特性有关：

large_constant_fields = grib_producing_large_constant_fields(gh, self->edition);
if (large_constant_fields) {
    if ((err = grib_set_long_internal(gh, self->binary_scale_factor, 0)) != GRIB_SUCCESS)
        return err;
    if ((err = grib_set_long_internal(gh, self->decimal_scale_factor, 0)) != GRIB_SUCCESS)
        return err;
    if (bits_per_value == 0) {
        if ((err = grib_set_long_internal(gh, self->bits_per_value, 16)) != GRIB_SUCCESS)
            return err;
    }
    return GRIB_SUCCESS;
}
else {
    bits_per_value = 0;
    if ((err = grib_set_long_internal(gh, self->bits_per_value, bits_per_value)) != GRIB_SUCCESS)
        return err;
    return GRIB_CONSTANT_FIELD;
}

上述代码首先调用 grib_producing_large_constant_fields() 函数返回是否需要编码常量场：

如果返回是，则正常编码，并返回 GRIB_SUCCESS
如果返回否，则启用常量场特性，将 bitsPerValue 设为 0，不编码数据段，并返回 GRIB_CONSTANT_FIELD。

grib_producing_large_constant_fields 函数代码在 grib_util.cc 文件中定义：

int grib_producing_large_constant_fields(grib_handle* h, int edition)
{
    // First check if the transient key is set
    grib_context* c                 = h->context;
    long produceLargeConstantFields = 0;
    if (grib_get_long(h, "produceLargeConstantFields", &produceLargeConstantFields) == GRIB_SUCCESS &&
        produceLargeConstantFields != 0) {
        return 1;
    }

    if (c->gribex_mode_on == 1 && edition == 1) {
        return 1;
    }

    // Finally check the environment variable via the context
    return c->large_constant_fields;
}

通过三种方式来确定编码常量场的数组：

produceLargeConstantFields 不为 0
gribex_model_on 启用且 GRIB 版本为 1
large_constant_fields 值为 1

而 large_constant_fields 的值在 grib_context.cc 文件中的 grib_context_get_default 函数中定义：

grib_context* grib_context_get_default()
{
    // ...skip...

        write_on_fail                       = codes_getenv("ECCODES_GRIB_WRITE_ON_FAIL");
        bufrdc_mode                         = getenv("ECCODES_BUFRDC_MODE_ON");
        bufr_set_to_missing_if_out_of_range = getenv("ECCODES_BUFR_SET_TO_MISSING_IF_OUT_OF_RANGE");
        bufr_multi_element_constant_arrays  = getenv("ECCODES_BUFR_MULTI_ELEMENT_CONSTANT_ARRAYS");
        grib_data_quality_checks            = getenv("ECCODES_GRIB_DATA_QUALITY_CHECKS");
        large_constant_fields               = codes_getenv("ECCODES_GRIB_LARGE_CONSTANT_FIELDS");
        no_abort                            = codes_getenv("ECCODES_NO_ABORT");
        debug                               = codes_getenv("ECCODES_DEBUG");
        gribex                              = codes_getenv("ECCODES_GRIBEX_MODE_ON");
        ieee_packing                        = codes_getenv("ECCODES_GRIB_IEEE_PACKING");
        io_buffer_size                      = codes_getenv("ECCODES_IO_BUFFER_SIZE");
        log_stream                          = codes_getenv("ECCODES_LOG_STREAM");
        no_big_group_split                  = codes_getenv("ECCODES_GRIB_NO_BIG_GROUP_SPLIT");
        no_spd                              = codes_getenv("ECCODES_GRIB_NO_SPD");
        keep_matrix                         = codes_getenv("ECCODES_GRIB_KEEP_MATRIX");
        file_pool_max_opened_files          = getenv("ECCODES_FILE_POOL_MAX_OPENED_FILES");

    // ...skip...
    return &default_grib_context;
}

large_constant_fields 的值来自环境变量 ECCODES_GRIB_LARGE_CONSTANT_FIELDS。

根据上述分析，可以通过两种方式关闭 ecCodes 的常量场特性：

为环境变量 ECCODES_GRIB_LARGE_CONSTANT_FIELDS 设置非 0 值
为 GRIB Key produceLargeConstantFields 设置非 0 值

实现

最简单的实现方案是在启动程序前设置环境变量 ECCODES_GRIB_LARGE_CONSTANT_FIELDS

export ECCODES_GRIB_LARGE_CONSTANT_FIELDS=1
srun --mpi=pmi2 \
    ${python_interpreter_path} -m reki_data_tool.postprocess.grid.gfs.ne dask-v1 \
    --input-file-path=${input_file_path} \
    --longitude="0:180" \
    --latitude="89.9375:0" \
    --output-file-path=${output_file_path} \
    --engine=mpi

第二种方案需要在编码前设置 produceLargeConstantFields 值，下面是对 reki.format.grib.eccodes.operator.extract_region() 函数代码的修改：

eccodes.codes_set_long(message, "produceLargeConstantFields", 1)
eccodes.codes_set_double_array(message, "values", values)

笔者倾向于使用第一种方案，已在 CMA-GFS V4.0 业务系统 2023060100 时次开始将环境变量添加到 grib2-ne 任务脚本中，修复 2023052200 时次 V4.0 更新引入的小 BUG。

总结

从本次 CMA-GFS V4.0 升级工作中得到的一个经验教训是对业务系统产品算法变更要进行完整的测试，不只要测试数据本身，还要测试数据处理工具链的兼容性。随着 CMA 系列数值模式产品的用户越来越多，任何一个非兼容性改动都可能给用户对接带来额外的工作量。

更改数据产品算法需要进行详细的测试，包括：

完整性：要素场不能多也不能少。新算法生成数据的元信息描述清单必须与旧算法生成数据保持一致，比如 wgrib2 的输出保持一致
正确性：最佳方式是新数据与原数据保持二进制一致，但 CMA 模式 GRIB2 数据使用 JPEG 2000 有损压缩，无法保持二进制一致，所以需要寻找合理的对比方法来检查数据是否一致以及数据是否正确
兼容性：确保新数据能兼容针对旧算法数据开发的程序，可以无缝对接。如果引入了不兼容部分，需要谨慎评估由此带来的额外工作量

在 CMA-GFS V4.0 升级修改数据任务算法过程中，笔者只关注了前两项，还没有找到合适的方法对比两个要素场是否一致，完全忽略了兼容性，导致新数据与旧数据在常量场上的编码方式不同，即便不影响使用，也没有达到保持一致的目的。

随着模式分辨率的提高，模式输出的数据量将显著增长，目前大部分串行后处理程序都面临运行时间增长的问题，需要进行并行化改造。在更新数据产品制作算法的同时也要保证前后数据的一致性。笔者后续将进一步寻找比较两个 GRIB2 要素场是否一致的合适方法。

参考

参考文章

WMO 技术手册：Manual on Codes - International Codes, Volume I.2, Annex II to the WMO Technical Regulations: Part B – Binary Codes, Part C – Common Features to Binary and Alphanumeric Codes
ECMWF 官网文档：What is a constant field - ecCodes GRIB FAQ
WGRIB2 官网文档：wgrib2: -g2clib
《Dask应用：批量转换GRIB2文件中的要素场》

工具软件

数据准备工具库：cemc-oper/reki

数据产品算法库：cemc-oper/reki-data-tool