GRIB学习笔记：样例文件分析 - 常量要素场
目录
本文通过分析一个GRIB2消息介绍所有数值都相等的常量要素场。
常量要素场（Constant Field）
如果某个要素场的所有数值都是同一个值，那么仅需要保存一个数值就能代表整个要素场。 虽然GRIB2使用的JPEG 2000压缩算法会考虑重复出现的数值，但如果我们只保存一个数值，就能避免对数据进行压缩和解压缩，也能节省一部分空间。
所以，ecCodes在处理常量要素场时，仅使用referenceValue保存数据。
示例
样例文件来自中国气象局数值预报中心GRAPES全球模式的预报结果，选用2019年8月17日00时次000时效的净长波辐射（Net long wave radiation flux (W m-2)）。
grib_dump的输出如下所示
***** FILE: constant_field.grib2
#============== MESSAGE 1 ( length=205 ) ==============
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 = 205
====================== SECTION_1 ( length=21, padding=0 ) ======================
1-4 section1Length = 21
5 numberOfSection = 1
6-7 centre = 38 [Beijing (RSMC) (common/c-11.table) ]
8-9 subCentre = 0
10 tablesVersion = 4 [Version implemented on 7 November 2007 (grib2/tables/1.0.table) ]
11 localTablesVersion = 1 [Unknown code table entry (grib2/tables/4/1.1.table) ]
12 significanceOfReferenceTime = 0 [Analysis (grib2/tables/4/1.2.table) ]
13-14 year = 2019
15 month = 8
16 day = 7
17 hour = 0
18 minute = 0
19 second = 0
20 productionStatusOfProcessedData = 0 [Operational products (grib2/tables/4/1.3.table) ]
21 typeOfProcessedData = 0 [Analysis products (grib2/tables/4/1.4.table) ]
====================== SECTION_3 ( length=72, padding=0 ) ======================
1-4 section3Length = 72
5 numberOfSection = 3
6 sourceOfGridDefinition = 0 [Specified in Code table 3.1 (grib2/tables/4/3.0.table) ]
7-10 numberOfDataPoints = 1036800
11 numberOfOctectsForNumberOfPoints = 0
12 interpretationOfNumberOfPoints = 0 [There is no appended list (grib2/tables/4/3.11.table) ]
13-14 gridDefinitionTemplateNumber = 0 [Latitude/longitude (Also called equidistant cylindrical, or Plate Carree) (grib2/tables/4/3.1.table) ]
15 shapeOfTheEarth = 6 [Earth assumed spherical with radius of 6,371,229.0 m (grib2/tables/4/3.2.table) ]
16 scaleFactorOfRadiusOfSphericalEarth = 0
17-20 scaledValueOfRadiusOfSphericalEarth = 0
21 scaleFactorOfEarthMajorAxis = 0
22-25 scaledValueOfEarthMajorAxis = 0
26 scaleFactorOfEarthMinorAxis = 0
27-30 scaledValueOfEarthMinorAxis = 0
31-34 Ni = 1440
35-38 Nj = 720
39-42 basicAngleOfTheInitialProductionDomain = 0
43-46 subdivisionsOfBasicAngle = 0
47-50 latitudeOfFirstGridPoint = 89875000
51-54 longitudeOfFirstGridPoint = 0
55 resolutionAndComponentFlags = 48 [00110000]
56-59 latitudeOfLastGridPoint = -89875000
60-63 longitudeOfLastGridPoint = 359750000
64-67 iDirectionIncrement = 250000
68-71 jDirectionIncrement = 250000
72 scanningMode = 0 [00000000]
====================== SECTION_4 ( length=58, padding=0 ) ======================
1-4 section4Length = 58
5 numberOfSection = 4
6-7 NV = 0
8-9 productDefinitionTemplateNumber = 8 [Average, accumulation, extreme values or other statistically processed values at a horizontal level or in a horizontal layer in a continuous or non-continuous time interval (grib2/tables/4/4.0.table) ]
10 parameterCategory = 5 [Long-wave Radiation (grib2/tables/4/4.1.0.table) ]
11 parameterNumber = 5 [Net long wave radiation flux (W m-2) (grib2/tables/4/4.2.0.5.table) ]
12 typeOfGeneratingProcess = 0 [Analysis (grib2/tables/4/4.3.table) ]
13 backgroundProcess = 0
14 generatingProcessIdentifier = 15
15-16 hoursAfterDataCutoff = 0
17 minutesAfterDataCutoff = 0
18 indicatorOfUnitOfTimeRange = 1 [Hour (grib2/tables/4/4.4.table) ]
19-22 forecastTime = 0
23 typeOfFirstFixedSurface = 1 [Ground or water surface (grib2/tables/4/4.5.table) ]
24 scaleFactorOfFirstFixedSurface = MISSING
25-28 scaledValueOfFirstFixedSurface = MISSING
29 typeOfSecondFixedSurface = 255 [Missing (grib2/tables/4/4.5.table) ]
30 scaleFactorOfSecondFixedSurface = MISSING
31-34 scaledValueOfSecondFixedSurface = MISSING
35-36 yearOfEndOfOverallTimeInterval = 2019
37 monthOfEndOfOverallTimeInterval = 8
38 dayOfEndOfOverallTimeInterval = 7
39 hourOfEndOfOverallTimeInterval = 0
40 minuteOfEndOfOverallTimeInterval = 0
41 secondOfEndOfOverallTimeInterval = 0
42 numberOfTimeRange = 1
43-46 numberOfMissingInStatisticalProcess = 0
47 typeOfStatisticalProcessing = 1 [Accumulation (grib2/tables/4/4.10.table) ]
48 typeOfTimeIncrement = 2 [Successive times processed have same start time of forecast, forecast time is incremented (grib2/tables/4/4.11.table) ]
49 indicatorOfUnitForTimeRange = 1 [Hour (grib2/tables/4/4.4.table) ]
50-53 lengthOfTimeRange = 0
54 indicatorOfUnitForTimeIncrement = 255 [Missing (grib2/tables/4/4.4.table) ]
55-58 timeIncrement = 0
====================== SECTION_5 ( length=23, padding=0 ) ======================
1-4 section5Length = 23
5 numberOfSection = 5
6-9 numberOfValues = 1036800
10-11 dataRepresentationTemplateNumber = 40 [JPEG2000 Packing (grib2/tables/4/5.0.table) ]
12-15 referenceValue = 0
16-17 binaryScaleFactor = -10
18-19 decimalScaleFactor = 0
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
====================== SECTION_6 ( length=6, padding=0 ) ======================
1-4 section6Length = 6
5 numberOfSection = 6
6 bitMapIndicator = 255 [A bit map does not apply to this product (grib2/tables/4/6.0.table) ]
====================== SECTION_7 ( length=5, padding=0 ) ======================
1-4 section7Length = 5
5 numberOfSection = 7
====================== SECTION_8 ( length=4, padding=0 ) ======================
1-4 7777 = 7777通过上面的dump文本可以看到，常量要素场的Section 7只有5个字节，没有数据内容。
我们知道，GRIB 2的实际数据与referenceValue，binaryScaleFactor和decimalScaleFactor相关。那么常量要素场的值是否与这三个都有关系呢？下面做一个简单的测试。
测试
使用grib_set修改上面文件中的数据值。
grib_set -d 1 constant_field.grib2 constant_field2.grib2使用grib_dump查看新生成的grib2文件。
***** FILE: constant_field2.grib2
#============== MESSAGE 1 ( length=205 ) ==============
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 = 205
====================== SECTION_1 ( length=21, padding=0 ) ======================
1-4 section1Length = 21
5 numberOfSection = 1
6-7 centre = 38 [Beijing (RSMC) (common/c-11.table) ]
8-9 subCentre = 0
10 tablesVersion = 4 [Version implemented on 7 November 2007 (grib2/tables/1.0.table) ]
11 localTablesVersion = 1 [Unknown code table entry (grib2/tables/4/1.1.table) ]
12 significanceOfReferenceTime = 0 [Analysis (grib2/tables/4/1.2.table) ]
13-14 year = 2019
15 month = 8
16 day = 7
17 hour = 0
18 minute = 0
19 second = 0
20 productionStatusOfProcessedData = 0 [Operational products (grib2/tables/4/1.3.table) ]
21 typeOfProcessedData = 0 [Analysis products (grib2/tables/4/1.4.table) ]
====================== SECTION_3 ( length=72, padding=0 ) ======================
1-4 section3Length = 72
5 numberOfSection = 3
6 sourceOfGridDefinition = 0 [Specified in Code table 3.1 (grib2/tables/4/3.0.table) ]
7-10 numberOfDataPoints = 1036800
11 numberOfOctectsForNumberOfPoints = 0
12 interpretationOfNumberOfPoints = 0 [There is no appended list (grib2/tables/4/3.11.table) ]
13-14 gridDefinitionTemplateNumber = 0 [Latitude/longitude (Also called equidistant cylindrical, or Plate Carree) (grib2/tables/4/3.1.table) ]
15 shapeOfTheEarth = 6 [Earth assumed spherical with radius of 6,371,229.0 m (grib2/tables/4/3.2.table) ]
16 scaleFactorOfRadiusOfSphericalEarth = 0
17-20 scaledValueOfRadiusOfSphericalEarth = 0
21 scaleFactorOfEarthMajorAxis = 0
22-25 scaledValueOfEarthMajorAxis = 0
26 scaleFactorOfEarthMinorAxis = 0
27-30 scaledValueOfEarthMinorAxis = 0
31-34 Ni = 1440
35-38 Nj = 720
39-42 basicAngleOfTheInitialProductionDomain = 0
43-46 subdivisionsOfBasicAngle = 0
47-50 latitudeOfFirstGridPoint = 89875000
51-54 longitudeOfFirstGridPoint = 0
55 resolutionAndComponentFlags = 48 [00110000]
56-59 latitudeOfLastGridPoint = -89875000
60-63 longitudeOfLastGridPoint = 359750000
64-67 iDirectionIncrement = 250000
68-71 jDirectionIncrement = 250000
72 scanningMode = 0 [00000000]
====================== SECTION_4 ( length=58, padding=0 ) ======================
1-4 section4Length = 58
5 numberOfSection = 4
6-7 NV = 0
8-9 productDefinitionTemplateNumber = 8 [Average, accumulation, extreme values or other statistically processed values at a horizontal level or in a horizontal layer in a continuous or non-continuous time interval (grib2/tables/4/4.0.table) ]
10 parameterCategory = 5 [Long-wave Radiation (grib2/tables/4/4.1.0.table) ]
11 parameterNumber = 5 [Net long wave radiation flux (W m-2) (grib2/tables/4/4.2.0.5.table) ]
12 typeOfGeneratingProcess = 0 [Analysis (grib2/tables/4/4.3.table) ]
13 backgroundProcess = 0
14 generatingProcessIdentifier = 15
15-16 hoursAfterDataCutoff = 0
17 minutesAfterDataCutoff = 0
18 indicatorOfUnitOfTimeRange = 1 [Hour (grib2/tables/4/4.4.table) ]
19-22 forecastTime = 0
23 typeOfFirstFixedSurface = 1 [Ground or water surface (grib2/tables/4/4.5.table) ]
24 scaleFactorOfFirstFixedSurface = MISSING
25-28 scaledValueOfFirstFixedSurface = MISSING
29 typeOfSecondFixedSurface = 255 [Missing (grib2/tables/4/4.5.table) ]
30 scaleFactorOfSecondFixedSurface = MISSING
31-34 scaledValueOfSecondFixedSurface = MISSING
35-36 yearOfEndOfOverallTimeInterval = 2019
37 monthOfEndOfOverallTimeInterval = 8
38 dayOfEndOfOverallTimeInterval = 7
39 hourOfEndOfOverallTimeInterval = 0
40 minuteOfEndOfOverallTimeInterval = 0
41 secondOfEndOfOverallTimeInterval = 0
42 numberOfTimeRange = 1
43-46 numberOfMissingInStatisticalProcess = 0
47 typeOfStatisticalProcessing = 1 [Accumulation (grib2/tables/4/4.10.table) ]
48 typeOfTimeIncrement = 2 [Successive times processed have same start time of forecast, forecast time is incremented (grib2/tables/4/4.11.table) ]
49 indicatorOfUnitForTimeRange = 1 [Hour (grib2/tables/4/4.4.table) ]
50-53 lengthOfTimeRange = 0
54 indicatorOfUnitForTimeIncrement = 255 [Missing (grib2/tables/4/4.4.table) ]
55-58 timeIncrement = 0
====================== SECTION_5 ( length=23, padding=0 ) ======================
1-4 section5Length = 23
5 numberOfSection = 5
6-9 numberOfValues = 1036800
10-11 dataRepresentationTemplateNumber = 40 [JPEG2000 Packing (grib2/tables/4/5.0.table) ]
12-15 referenceValue = 1
16-17 binaryScaleFactor = -10
18-19 decimalScaleFactor = 0
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
====================== SECTION_6 ( length=6, padding=0 ) ======================
1-4 section6Length = 6
5 numberOfSection = 6
6 bitMapIndicator = 255 [A bit map does not apply to this product (grib2/tables/4/6.0.table) ]
====================== SECTION_7 ( length=5, padding=0 ) ======================
1-4 section7Length = 5
5 numberOfSection = 7
====================== SECTION_8 ( length=4, padding=0 ) ======================
1-4 7777 = 7777可以看到，referenceValue 从 0 变成 1，即数据值只与 referenceValue 有关，与scale factor无关。
解码
首先判断Section 7中是否有数据值，如果没有数据，则认为该消息是一个常量要素场。
如果section 7只有5个字节，则判断该消息是常量要素场。
bool GribSection7::parseFile(std::FILE* file, bool header_only) {
const auto buffer_length = section_length_ - 5;
if (buffer_length == 0) {
return true;
}
// ...skip...
return true;
}解码数据值时，如果读取的数据为空，则按照常量要素场解码。
bool DataValuesProperty::decodeValues(GribMessageHandler* container) {
// constant field has no data values
if (raw_value_bytes_.empty()) {
return decodeConstantFields(container);
} else {
return decodeNormalFields(container);
}
}解码常量要素场时，将referenceValue作为数据值。
bool DataValuesProperty::decodeConstantFields(GribMessageHandler* container) {
data_count_ = container->getLong("numberOfValues");
const auto reference_value = static_cast<float>(container->getDouble("referenceValue"));
code_values_.resize(data_count_);
std::fill(std::begin(code_values_), std::end(code_values_), reference_value);
return true;
}编码
编码前先根据实际的数据得到编码需要的各个数值。如果需要保存的data_count_为0，则按照常量要素场编码数据。
bool DataValuesProperty::encodeValues(GribMessageHandler* container) {
calculate(container);
raw_value_bytes_.clear();
if (data_count_ == 0) {
return encodeConstantFields(container);
} else {
return encodeNormalFields(container);
}
}根据数据的最大值和最小值判断是否为常量要素场。
// algorithm is from NCEP wgrib2 (grib2/g2clib-1.4.0/jpcpack.c)
void DataValuesProperty::calculate(GribMessageHandler* container) {
// ...skip ...
const auto [min_value_iter, max_value_iter] = std::minmax_element(std::begin(code_values_), std::end(code_values_));
const auto min_value = (*min_value_iter) * decimal_scale * bscale;
const auto max_value = (*max_value_iter) * decimal_scale * bscale;
const auto max_number_step = static_cast<long>(max_value-min_value);
// ...skip...
auto reference_value = static_cast<float>(min_value);
// if max_value and min_value are not equal, use data values.
// or use empty data values for constant field.
if(max_value!=min_value && max_number_step != 0) {
const auto temp = std::log((static_cast<double>(max_number_step) + 1)) / log2;
bits_per_value = std::ceil(temp);
} else {
data_count_ = 0;
bits_per_value = 0;
reference_value = *min_value_iter;
}
container->setDouble("referenceValue", reference_value);
container->setLong("bitsPerValue", bits_per_value);
}对于常量要素场，将referenceValue设为数据值，将bitsPerValue设为0。
bool DataValuesProperty::encodeConstantFields(GribMessageHandler* container) {
const auto reference_value = code_values_[0];
const auto bits_per_value = 0;
container->setDouble("referenceValue", reference_value);
container->setLong("bitsPerValue", bits_per_value);
raw_value_bytes_.clear();
data_count_ = 0;
return true;
}参考
上述代码摘取自perillaroc/nwpc-codes-cpp。
