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