Problem parsing negative longitude values in Lambert Azimuthal Equal Area grid specs

[JohnHalleyGotway]

I'm writing to inform you of a potential bug in the grib2c library revealed via testing with the wgrib2 utility.

Here's the GRIB2 file to demonstrate the issue: ukv_agl_temperature_1.5_12.grib2.gz

Running wgrib2 -v prints the specification of the Lambert Azimuthal Equal Area grid on which this data resides.

wgrib2 -V ukv_agl_temperature_1.5_12.grib2 
    Lambert Azimuthal Equal Area grid: (1042 x 970) input WE:SN output WE:SN res 48
    Lat1 44.517153 Lon1 2164.600777 Cen Lon -2.500000 Std Par 54.900000
    Dx 2000.000000 m Dy 2000.000000 m mode 48

The problem is in this Lon1 2164.600777 value. That should actually be -17.11713 instead!

This can be seen by running the following wgrib2 command:

wgrib2 ukv_agl_temperature_1.5_12.grib2 -get_byte 3 43 4 -get_int 3 43 1

1:0:3-43=129,5,47,201:3-43=-17117129

The GRIB2 Template 3.140 (for LAEA grids), indicates that "LO1" is stored in bytes 43 to 46. The -get_int option for bytes 43 to 46 correctly returns a value of -17117129 in micro-degrees, or -17.117129 in degrees.

Note that I also confirmed this by running the grib_dump utility from ECMWF:

grib_dump ukv_agl_temperature_1.5_12.grib2

longitudeOfFirstGridPointInDegrees = -17.1171;

The issue seems to be that the grib2c library isn't correctly parsing the leading negative sign when populating the igdstmpl array. This issue may also exist for the central longitude value of -2.5. While wgrib2 does correctly report its negative value, the MET code on which I work did not. You can see a comment about this issue in this line of the MET code.

FYI, my current workaround in the MET code is calling the following parse_int4() utility to reprocess the bytes checking for the leading negative sign.

int parse_int4(g2int i) {
   unsigned char c[4];
   unsigned long n = i;

   c[0] = (n >> 24) & 0xFF;
   c[1] = (n >> 16) & 0xFF;
   c[2] = (n >> 8) & 0xFF;
   c[3] = n & 0xFF;

   //  convert unsigned char to signed integer
   int i_val;
   if(c[0] & 0x80) {
      i_val = -(((c[0] & 0x7f) << 24) + (c[1] << 16) + (c[2] << 8) + c[3]);
   }
   else {
      i_val = (c[0] << 24) + (c[1] << 16) + (c[2] << 8) + c[3];
   }

   return(i_val);
}

[edwardhartnett]

@webisu is this the problem you were referring to between g2c and wgrib2?

[webisu]

The grib2 standard sets the convention that longitudes are between 0 and 360 degrees. So Sec 3, Octets 43-46 and Octets 51-54 should be interpreted as an unsigned integer. Strictly speaking, the original grib file is encoded incorrectly. However, there is no drawback if one interprets the values to be a 4-byte signed integer, So wgrib2 v3.1.2 interprets "lon1" and "Cen lon" to be a signed integers. How other libraries handle poorly encoded grib files is up to their maintainers.

[JohnHalleyGotway]

@webisu thanks for the followup. So perhaps this issue belongs to wgrib2 more so than the grib2c library?

Although, I did discuss this with @GwenChen-NOAA at the METplus NOAA User Telecon this morning. She pointed me to a more recent version of wgrib2 (v3.1.4b2) on Hera with which I could test. And indeed, the problem is resolved in that version.

# On Hera
wget https://github.com/NOAA-EMC/NCEPLIBS-g2c/files/12135232/ukv_agl_temperature_1.5_12.grib2.gz
gunzip  ukv_agl_temperature_1.5_12.grib2.gz
source /home/Wesley.Ebisuzaki/wgrib2_mods.sh
/home/Wesley.Ebisuzaki/bin/wgrib2.v3.1.4b2 -V ukv_agl_temperature_1.5_12.grib2

Produces this output:

    Lambert Azimuthal Equal Area grid: (1042 x 970) input WE:SN output WE:SN res 48
    Lat1 44.517153 Lon1 -17.117129 Cen Lon -2.500000 Std Par 54.900000
    Dx 2000.000000 m Dy 2000.000000 m mode 48

Note that Lon1 -17.117129 replaces Lon1 2164.600777.

So I'd say that this issue can probably just be closed as being no longer relevant... but will leave it to the repo managers to do so.