Incorrect values parsed for Complex Packing due to missing 2nd Scale Factor

Simple Description on the Bug

When parsing data at a constant Isobaric level that is encoded with Complex Packing, there are some messages that have a "Missing" 2nd Scale Factor that cannot be parsed resulting in incorrect values being produced.

Steps to Reproduce

Attempt to parse a GRIB file with a supposedly missing 2nd Scale Factor. I have posted an example GRIB file that displays this here: https://drive.google.com/drive/folders/1Zxs21uamCSF9LVUWHYO-dE-A9rzb91S_?usp=sharing. The error_grib_file contains the troublesome data and the correct_test_file contains GRIB data that can be properly parsed.

Note that I am only parsing messages that are at Isobaric layers (Code Table 4-5 100), not those that are altitude layers. The .idx index file shows which messages are of which type. The specific message that I am testing on is number 295 in the index file.

The following code explains how I was comparing the data from the test file to a different file that should have fairly similar data. Note that this was copy-pasted then adapted, so it will likely require some changes to function properly.

    #[test]
    fn test_grib_format() {
        // read through every message in a grib file and parse the first few values from each
        // message
        let p = "error_grib_file";
        let f = std::fs::File::open(&p).unwrap();
        let rd = std::io::BufReader::new(f);
        let f1 = grib::from_reader(rd).unwrap();

        let p = "correct_test_file";
        let f = std::fs::File::open(&p).unwrap();
        let rd = std::io::BufReader::new(f);
        let f2 = grib::from_reader(rd).unwrap();

        // set up decoder
        let f1 = f1
            .iter()
            .map(|lyr| {
                let ((i, j), fm) = lyr;
                let strn = format!(
                    "\n\nIncorrect\n\n{}, {}\n{}\n{:?}\n{:?}\n{:?}\n{:?}\n{:?}",
                    i,
                    j,
                    fm.describe(),
                    fm.indicator(),
                    fm.grid_def(),
                    fm.prod_def(),
                    fm.prod_def().forecast_time().unwrap().to_string(),
                    fm.repr_def(),
                );
                let dc = Grib2SubmessageDecoder::from(fm).unwrap();

                let dv = dc.dispatch().unwrap();

                let val = dv.take(10).collect::<Vec<_>>();

                (strn, val)
            })
            .collect::<Vec<_>>();

        // set up decoder
        let f2 = f2
            .iter()
            .map(|lyr| {
                let ((i, j), fm) = lyr;
                let strn = format!(
                    "\n\nCorrect\n\n{}, {}\n{}\n{:?}\n{:?}\n{:?}\n{:?}\n{:?}",
                    i,
                    j,
                    fm.describe(),
                    fm.indicator(),
                    fm.grid_def(),
                    fm.prod_def(),
                    fm.prod_def().forecast_time().unwrap().to_string(),
                    fm.repr_def(),
                );
                let dc = Grib2SubmessageDecoder::from(fm).unwrap();

                let dv = dc.dispatch().unwrap();

                let val = dv.take(10).collect::<Vec<_>>();

                (strn, val)
            })
            .collect::<Vec<_>>();

/*
Code to pair up messages together for printing purposes
*/

        f2.into_iter().zip(f1.into_iter()).for_each(|(f, l)| {
            println!("{}", f.0);
            println!("{:?}", f.1);
            println!("{}", l.0);
            println!("{:?}", l.1);
        });

        // ensure the test fails
        panic!()
    }

Expected Behavior

Parse the values properly.

Actual Behavior

Produce values that are wildly inaccurate and nowhere near the values produced by the file that is parsed correctly. Here is example output from the aforementioned message that I tested with:

Correct

Grid:                                   Latitude/longitude
  Number of points:                     1038240
Product:                                Analysis or forecast at a horizontal level or in a horizontal layer at a point in time
  Parameter Category:                   Momentum
  Parameter:                            u-component of wind
  Generating Proceess:                  Forecast
  Forecast Time:                        0
  Forecast Time Unit:                   Hour
  1st Fixed Surface Type:               Isobaric surface
  1st Scale Factor:                     0
  1st Scaled Value:                     20000
  2nd Fixed Surface Type:               code '255' is not implemented
  2nd Scale Factor:                     0
  2nd Scaled Value:                     0
Data Representation:                    Grid point data - complex packing and spatial differencing
  Number of represented values:         1038240

Indicator { discipline: 0, total_length: 581327 }
GridDefinition { payload: [0, 0, 15, 215, 160, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 160, 0, 0, 2, 209, 0, 0, 0, 0, 255, 255, 255, 255, 5, 93, 74, 128, 0, 0, 0, 0, 48, 133, 93, 74, 128, 21, 113, 89, 112, 0, 3, 208, 144, 0, 3, 208, 144, 0] }
ProdDefinition { payload: [0, 0, 0, 0, 2, 2, 2, 0, 81, 0, 0, 0, 1, 0, 0, 0, 0, 100, 0, 0, 0, 78, 32, 255, 0, 0, 0, 0, 0] }
"0 [h]"
ReprDefinition { payload: [0, 15, 215, 160, 0, 3, 195, 167, 184, 23, 0, 0, 0, 1, 9, 0, 1, 0, 98, 88, 209, 154, 255, 255, 255, 255, 0, 0, 105, 56, 0, 4, 0, 0, 0, 1, 1, 0, 0, 0, 71, 7, 2, 2] }
[-24.74382, -24.64382, -24.64382, -24.543821, -24.44382, -24.44382, -24.34382, -24.24382, -24.24382, -24.14382]

Incorrect

Grid:                                   Latitude/longitude
  Number of points:                     1038240
Product:                                Analysis or forecast at a horizontal level or in a horizontal layer at a point in time
  Parameter Category:                   Temperature
  Parameter:                            Temperature
  Generating Proceess:                  Forecast
  Forecast Time:                        0
  Forecast Time Unit:                   Hour
  1st Fixed Surface Type:               Isobaric surface
  1st Scale Factor:                     0
  1st Scaled Value:                     30
  2nd Fixed Surface Type:               code '255' is not implemented
  2nd Scale Factor:                     Missing
  2nd Scaled Value:                     255
Data Representation:                    Grid point data - complex packing and spatial differencing
  Number of represented values:         1038240

Indicator { discipline: 0, total_length: 675339 }
GridDefinition { payload: [0, 0, 15, 215, 160, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 160, 0, 0, 2, 209, 0, 0, 0, 0, 255, 255, 255, 255, 133, 93, 74, 128, 0, 0, 0, 0, 48, 5, 93, 74, 128, 21, 113, 89, 112, 0, 3, 208, 144, 0, 3, 208, 144, 64] }
ProdDefinition { payload: [0, 0, 0, 0, 0, 0, 2, 255, 16, 0, 0, 0, 1, 0, 0, 0, 0, 100, 0, 0, 0, 0, 30, 255, 255, 0, 0, 0, 255] }
"0 [h]"
ReprDefinition { payload: [0, 15, 215, 160, 0, 3, 70, 174, 52, 0, 0, 0, 0, 2, 8, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 89, 0, 4, 0, 0, 0, 1, 1, 0, 0, 0, 64, 6, 2, 2] }
[252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28]

Additional Context

Additionally, the Parameter Category and Parameter are seemingly incorrectly parsed as they should both be Momentum, u-component of wind and the incorrect file produces Temperature, Temperature.

@LafeWessel Thank you for reporting. Since I am a little confused, let me ask a question: what does "incorrect values" mean? Does "Missing" seem to be incorrect? Or, decoded values are incorrect?

Note

I am sorry if you are familiar with fixed surfaces in GRIB2, but I would like to explain a little about it here.

Fixed surfaces

Each "layer" (or submessage) in GRIB2 has one or two fixed surfaces. Since most layers express values at just one altitude level, they have only one fixed surface (e.g. humidity at 700hPa, temperature at 30cm depth, etc.). For such layers, the second fixed surface is missing.

However, there are cases that it would be nice for layers to have the second fixed surface. Fox example, some values are defined for layers with a certain thickness. Actually, I have not seen so much data that has a second fixed surface, so I asked ChatGPT to explain its usage to me; according to ChatGPT, the second fixed surface has the following uses (I am not sure how much of a lie this explanation contains).

Me:

Each message in GRIB2 files contains 1 or 2 fixed surfaces. In what cases is the second fixed surface used?

ChatGPT:

In GRIB2 files, each message can indeed contain one or two fixed surfaces. The second fixed surface is typically used to represent a secondary level or layer of the atmosphere or Earth's surface that is relevant to the meteorological data being stored.

The second fixed surface is often employed in the following scenarios:

Vertical interpolation: In some cases, meteorological variables may need to be interpolated or derived between two fixed surfaces to provide a more complete representation of the atmosphere. For example, temperature data may be provided at both 1000 hPa (first fixed surface) and 850 hPa (second fixed surface) to allow for vertical interpolation and calculation of temperature values at different pressure levels.

Layer variables: Certain meteorological variables are defined for specific layers or thicknesses in the atmosphere. In such cases, the first fixed surface represents the lower boundary of the layer, while the second fixed surface represents the upper boundary. This allows for the representation of variables like the thickness of a layer, cloud cover within a specific layer, or precipitation accumulation over a particular layer.

Earth's surface and pressure levels: The second fixed surface can also be used in combination with the Earth's surface or pressure levels. For example, a message might include data at the Earth's surface (e.g., surface temperature) as the first fixed surface, and then provide data at a specific pressure level (e.g., 500 hPa) as the second fixed surface.

The inclusion of a second fixed surface in a GRIB2 message enhances the flexibility and completeness of the data by accommodating various vertical representations and calculations. It allows for more comprehensive analysis and modeling of the atmosphere or Earth's surface based on different levels, layers, or transitions between fixed surfaces.

Expressions of fixed surfaces

When a layer does not have the second fixed surface, that fixed surface is expressed as "missing". That status is expressed by setting code 255 to the type of that fixed surface. So, the following means 2nd fixed surface is missing. (I think it would be nice to express the type as "Missing" instead of "code '255' is not implemented"...)

  2nd Fixed Surface Type:               code '255' is not implemented

When a fixed surface is missing, the scale factor and scaled value are usually set to 0xff and 0xffffffff, but there are cases where the scale factor and scaled value are both set to 0. Anyway, following two layers both have a missing second surface:

  2nd Fixed Surface Type:               code '255' is not implemented
  2nd Scale Factor:                     0
  2nd Scaled Value:                     0

  2nd Fixed Surface Type:               code '255' is not implemented
  2nd Scale Factor:                     Missing
  2nd Scaled Value:                     255

Scale factor and scaled value of fixed surface

With a fixed surface type, scale factor and scaled value are used like this:

975 hPa (975 * 10 ** (-(-2)) Pa (The unit of isobaric surface is defined to be Pa))

  1st Fixed Surface Type:               Isobaric surface
  1st Scale Factor:                     -2
  1st Scaled Value:                     975

Mean sea level:

  1st Fixed Surface Type:               Mean sea level
  1st Scale Factor:                     Missing
  1st Scaled Value:                     Missing

1.5 m high above the ground (15 * 10 ** -1 m)

  1st Fixed Surface Type:               Specified height level above ground
  1st Scale Factor:                     1
  1st Scaled Value:                     15

Thanks for the clarification! I confess I do not know that much about the Fixed Surfaces - though I do now.

My report is specifically about incorrect values being parsed from the messages, not a missing Fixed Surface. Though I don't know how the parsing works internally, my thought was that the 255 from the 2nd Scaled Value was "polluting" the values that were parsed from the message. Here's my reasoning:

From the correct message, we see:

  1st Fixed Surface Type:               Isobaric surface
  1st Scale Factor:                     0
  1st Scaled Value:                     20000
  2nd Fixed Surface Type:               code '255' is not implemented
  2nd Scale Factor:                     0
  2nd Scaled Value:                     0

which produces the following (first ten) values:

[-24.74382, -24.64382, -24.64382, -24.543821, -24.44382, -24.44382, -24.34382, -24.24382, -24.24382, -24.14382]

For the incorrect message, we see:

  1st Fixed Surface Type:               Isobaric surface
  1st Scale Factor:                     0
  1st Scaled Value:                     30
  2nd Fixed Surface Type:               code '255' is not implemented
  2nd Scale Factor:                     Missing
  2nd Scaled Value:                     255

which produces the following (first ten) values:

[252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28]

Both of these messages should produce roughly the same data, but they are clearly very far off from one another. My thought was that the 255 from the 2nd Scaled Value from the incorrect message caused the values parsed to be incorrect as they are rather close to 255. Whereas the correct message has 0 for a 2nd Scaled Value and produces the correct values.

I would guess that the 2nd Scale Factor is missing, but my hunch was/is that the placeholder of 255 for the missing 2nd Scaled Factor causes different values to be parsed.

Does that make any more sense, or have I misunderstood what your comment says/how it works internally?

Thank you. I understand the intent of your report. However, so far it doesn't sound like there is a problem with the library. Here are two reasons why I think so.

Section 4 is not used in decoding

First, the information contained in Section 4 (Product Definition Section), such as fixed surfaces, should not affect decoding. A GRIB2 message consists of following sections:

Section 0 (Indicator Section) b"GRIB", GRIB version, message size, etc.
Section 1 (Identification Section) Information about the data source, such as issuing authority, reference time, etc.
Section 2 (Local Use Section) Information defined by each local authorities.
Section 3 (Grid Definition Section) Definition of grid system. The type of grid system (e.g., equirectangular, Lambert, etc.) and the parameters that specifically define each grid system (e.g., longitude and latitude at the four corners, number of grid points in the latitudinal and longitudinal directions, etc.).
Section 4 (Product Definition Section) Definition of the meaning of the included values. For example, meteorological elements, fixed surface (altitude level), forecast time, etc.
Section 5 (Data Representation Section) Definition of what representation is used for including grid point values. Simply put, it defines the compression method and the parameter values used to compress the data with that method.
Section 6 (Bit-map Section) A bitmap representing a series of missing values in the data (array). The information about whether each grid point value is a missing value or not is separated from the array of grid point values and stored in this section, and only the grid point values other than missing values are compressed using the method defined in Section 5 and stored in Section 7.
Section 7 (Data Section) The main body of data compressed by the method defined in Section 5. Some parameters are not included in Section 5, but are included here.
Section 8 (End Section) b "7777".

Of these 9 sections, only Sections 3, 5, 6, and 7 are required to decode grid point values. Information on meteorological elements, fixed surfaces, and forecast times are included in Section 4, function as data attributes and are necessary to understand the meaning of the data, but are not used in the decoding process.

Values

Second, please note the following points in your data.

  Parameter:                            u-component of wind
(snip)
[-24.74382, -24.64382, -24.64382, -24.543821, -24.44382, -24.44382, -24.34382, -24.24382, -24.24382, -24.14382]

  Parameter:                            Temperature
(snip)
[252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28, 252.28]

The u component of the wind is positive for the east and negative if the wind blows to the west. Also, it is not surprising that wind speeds of -24 m/s can be observed regularly at high altitudes.

The temperature cannot be lower than absolute zero, and it is reasonable to assume that the meteorological data is within 273±100 K. Therefore, the value of 252 K is not particularly strange.

This is very helpful, thank you. After further poking around in the code that I have, I realized a couple things:

The messages were not matching up how I anticipated, so for the example messages that I expected to be "u-component", only one of them actually was. None were being read incorrectly.
The lat-lon grid I was using is flipped on the y-axis for the second set of data I have, so the logic I used to parse it out was showing correct data when I tested at 0 latitude but was incorrect anywhere else.

Thus, the bug is on my end in how I consume the GRIB data, there is nothing wrong with your library! Thanks for all the help, sorry for wasting your time.

noritada / grib-rs