Inaccurate method for calculating GRIB1 coordinates

[tjtg]

When opening GRIB1 data defined on a 1/24 degree (ie. 0.041666... degree) resolution regular latitude-longitude (Plate Carree) grid, iris-grib calculates longitude coordinates incorrectly.
This seems to be due to a combination of low numerical precision in GRIB1 metadata and an inaccurate calculation method in iris-grib.

In GRIB1, the metadata for these latitude-longitude grid coordinates is expressed using 3-octet (24 bit) integers in units of 10^-3 degrees. (for reference, see GRIB 1 guide, table C sundry grid definitions, latitude/longitude grids).
For example, the 0.0416... degree grid spacing will be rounded off to 42 millidegrees, an error of about 1%.
GRIB2 and GRIB3 have a more accurate representation using 4-octet integers and units of 10^-6 degrees.

The iris-grib method _regular_longitude_common uses Ni (the number of i-direction points), iDirectionIncrementInDegrees (the spacing between grid points in the i-direction) and longitudeOfFirstGridPointInDegrees to calculate the longitude coordinates in degrees. The "InDegrees" coordinates are provided by eccodes/gribapi as floating-point numbers so that they're independent of the GRIB version in use.

This calculation method based on multiplying the increment re-uses the rounded-off grid spacing, resulting in a significantly different last grid point coordinate to that specified in the GRIB message.
A more accurate calculation method would be to use longitudeOfFirstGridPointInDegrees and longitudeOfLastGridPointInDegrees with numpy linspace to generate the series of coordinates, then perhaps check for consistency with iDirectionIncrementInDegrees.

There is some similar code for processing latitude earlier in the file, but I first encountered this issue with longitudes.

[pp-mo]

Hi @tjtg thanks for raising this !
We are looking into the implications. Will update in due course

[ESadek-MO]

@SciTools/peloton This has been quiet for a while, is this still relevant @tjtg?