Incorrect missing value handling when converting geopotential height (gh) to geopotential (z) in GRIB files

[zhenkunl]

What happened?

Description:
I encountered an issue when converting 850hPa geopotential height (gh) to geopotential (z) using the formula z = gh * 9.80665 in NCEP GFS forecast GRIB files. After writing the converted data back to a GRIB file and reading it again, some values appear as 3.4028234663852886e+38 (IEEE 32-bit float max), while they were originally close to 9999 (likely missing values) before the conversion.

The issue persists across:

• earthkit-data (which depends on eccodes-python)
• eccodes-python directly
• eccodes C library

This suggests a problem in the underlying eccodes C library when handling missing values during data transformation and write-back.

Attachments:
I've attached a sample GRIB file (gh850.grib) containing the 850hPa geopotential height field from NCEP GFS that reproduces the issue.

Additional notes:
The issue appears specifically when:

• The original GRIB file uses a missing value marker (e.g., 9999)
• After numerical transformation, the missing value handling seems to break
• The output file contains IEEE max float values instead of correct values which close to 9999

I suspect this might be related to how eccodes handles missing value indicators during data modification. After further investigation, I tested earlier versions of the eccodes C library and found that version 2.26 processes the conversion correctly without generating the anomalous 3.4028234663852886e+38 values. This suggests that a regression was introduced in later versions of eccodes. Any insight or fix would be appreciated.

What are the steps to reproduce the bug?

Steps to reproduce:

1. Read a GRIB file containing 850hPa geopotential height (gh) from NCEP GFS.
2. Convert gh to z by multiplying by 9.80665.
3. Write the modified field to a new GRIB file.
4. Read the new file and observe that values originally near 9999 become 3.4028234663852886e+38.

Minimal reproducible example:
Using earthkit-data

import earthkit.data as ekd
# Read input GRIB file
gh = ekd.from_source("file", "gh850.grib")
f = gh[0]
# Convert gh to z
values = f.to_numpy() * 9.80665
f_new = f.clone(values=values, metadata=f.metadata().override(shortName="z"))
# Write back to new GRIB file
f_new.to_target("file", "z850.grib")

# Read back and check values
z = ekd.from_source("file", "z850.grib")
f = z[0]
md = f.metadata()
print(md["missingValue"], md["level"], f.values.max(), f.values.min())
# 9999 850 3.4028234663852886e+38 9277.9560546875

Version

v2.44.0

Platform (OS and architecture)

ubuntu 22.04/Mac OS 15.7.3

Relevant log output

In [1]: import earthkit.data as ekd

In [2]: import numpy as np

In [3]: z = ekd.from_source("file", "z850.grib")

In [4]: f = z[0]

In [5]: values = f.to_numpy()

In [6]: indices = np.where(values > 10000000.0)

In [7]: indices
Out[7]:
(array([140, 602, 604, 617, 618, 619, 638, 647, 668]),
 array([1245,  544,  484,  491, 1399,   12,   22, 1398,  802]))

In [8]: gh = ekd.from_source("file", "gh850.grib")

In [9]: f = gh[0]

In [10]: values_orig = f.to_numpy()

In [11]: values_new = f.to_numpy() * 9.80665

In [12]: values_orig[indices]
Out[12]:
array([1019.62425, 1019.62425, 1019.62425, 1019.62425, 1019.62425,
       1019.62425, 1019.60825, 1019.60825, 1019.62425])

In [13]: values_new[indices]
Out[13]:
array([9999.09815126, 9999.09815126, 9999.09815126, 9999.09815126,
       9999.09815126, 9999.09815126, 9998.94124486, 9998.94124486,
       9999.09815126])

Accompanying data

The sample GRIB file is gh850.tar.gz

[shahramn]

I tried this with the ecCodes tools ( not using earthkit, python etc)
First of all the input file gh850.grib has no missing values. We can check this via wgrib2 as well as ecCodes:

% wgrib2 -stats gh850.grib 
1:0:ndata=1038240:undef=0:mean=1400.73:min=946.088:max=1649.45:cos_wt_mean=1451.27

Note how the "undef" field is zero

% grib_ls -jp numberOfDataPoints,numberOfCodedValues,numberOfMissing,missingValuesPresent gh850.grib 
    "numberOfDataPoints": 1038240,
    "numberOfCodedValues": 1038240,
    "numberOfMissing": 0,
    "missingValuesPresent": 0

so no missing values.

[shahramn]

Now let's use the ecCodes tool grib_set to multiply all values by 9.80665

grib_set -s scaleValuesBy=9.80665 gh850.grib modified.grib

Now check the statistics before and after

% grib_ls -jn statistics -P file gh850.grib modified.grib 
  {
    "file": "gh850.grib",
    "max": 1649.45,
    "min": 946.088,
    "avg": 1400.73,
    "sd": 145.291,
    "skew": -0.780557,
    "kurt": -0.496769,
  },
  {
    "file": "modified.grib",
    "max": 16175.5,
    "min": 9277.96,
    "avg": 13736.5,
    "sd": 1424.78,
    "skew": -0.780535,
    "kurt": -0.496832,

So no problem there

[sandorkertesz]

Now let's use the ecCodes tool grib_set to multiply all values by 9.80665

grib_set -s scaleValuesBy=9.80665 gh850.grib modified.grib

Now check the statistics before and after

% grib_ls -jn statistics -P file gh850.grib modified.grib 
  {
    "file": "gh850.grib",
    "max": 1649.45,
    "min": 946.088,
    "avg": 1400.73,
    "sd": 145.291,
    "skew": -0.780557,
    "kurt": -0.496769,
  },
  {
    "file": "modified.grib",
    "max": 16175.5,
    "min": 9277.96,
    "avg": 13736.5,
    "sd": 1424.78,
    "skew": -0.780535,
    "kurt": -0.496832,

So no problem there

Please check the actual (numpy) values and you will see the problem! If in the earthkit code above I also print:

print(md["max"])

I get

16175.4560546875

though the values max is float max.

So it indicates that there is a discrepancy in ecCodes between the values and the statistics namespace.

[shahramn]

I see it now. Yes there is indeed a problem

% grib_ls -jp file,numberOfDataPoints,numberOfCodedValues,numberOfMissing,missingValuesPresent \
    gh850.grib modified.grib
 {
    "file": "gh850.grib",
    "numberOfDataPoints": 1038240,
    "numberOfCodedValues": 1038240,
    "numberOfMissing": 0,
    "missingValuesPresent": 0
 },
 {
    "file": "modified.grib",
    "numberOfDataPoints": 1038240,
    "numberOfCodedValues": 1038240,
    "numberOfMissing": 9,
    "missingValuesPresent": 1
}

The scaling has introduced 9 dummy missing values! This is a bug.
Thank you for reporting this.
We will work on a fix

[zhenkunl]

There is no missing values in original gh850.grib since it is a forecast field. I am just wondering why values those are close to 9999 will be converted to float max. I attempted to investigate and fix the issue myself, but after searching through the entire eccodes repository codebase, I couldn't find any occurrence of the value 3.4028234663852886e+38 (likely FLT_MAX in grib_optimize_decimal_factor()?).
For comparison, I uploaded another 850hPa geopotential height GRIB file (gh850-2.grib.tar.gz) and performed the exact same operations using the same code, and it processed correctly without any issues.

[zhenkunl]

I also observed an intriguing phenomenon: when I assigned the entire array to a value close to 9999 (e.g., 9999.09815126), write the data to a GRIB file, then read it back, I noticed some inexplicable changes in the data values. Moreover, each time I perform the same operation, the resulting data values are different. Below is the code I used:

In [1]: import earthkit.data as ekd

In [2]: import numpy as np

In [3]: gh = ekd.from_source("file", "gh850.grib")

In [4]: f = gh[0]

In [5]: values = f.to_numpy()

In [6]: new_values = np.full_like(values, 9999.09815126)

In [7]: f_new = f.clone(values=new_values, metadata=f.metadata().override(shortName="z"))

In [8]: f_new.to_target("file", "z850.grib")

In [9]: z = ekd.from_source("file", "z850.grib")

In [10]: f = z[0]

In [11]: md = f.metadata()

In [12]: print(md["missingValue"], md["level"], f.values.max(), f.values.min())
9999 850 3.602879701896397e+16 0.0

In [13]: new_values.max()
Out[13]: np.float64(9999.09815126)

In [14]: new_values.min()
Out[14]: np.float64(9999.09815126)

In [15]: f.values
Out[15]:
array([7.77091700e+15, 2.57113733e+16, 9.00719925e+15, ...,
       2.76573308e+16, 2.77089347e+16, 2.74275605e+16], shape=(1038240,))

[zhenkunl]

When performing this test, all the values read back from the GRIB file differ from the original 9999.09815126 or 3.4028234663852886e+38. This non-reproducible, non-deterministic behavior when handling values near the missing value threshold suggests there might be issues with how certain floating-point values are encoded, stored, or decoded in the GRIB format, particularly around the 9999 range which is commonly used as a missing value indicator in meteorological data.

[ZhangquanZhangQ]

I encountered similar issues to the author's, and I tried using libraries like pygrib, earthkit, isir, and cfgrib, all of which exhibited comparable problems. Of course, it might be because these libraries all rely on eccodes at their core.

[shahramn]

In the meantime, you could (as a temporary measure) convert these GRIB files to use a different packing scheme. For example CCSDS
This is how you can do this:

grib_set -w isGridded=1 -r -s packingType=grid_ccsds gh850.grib gh850.ccsds.grib

Now apply the scaling and it will work ( no missing values issue )

grib_set -s scaleValuesBy=9.80665  gh850.ccsds.grib  gh850.ccsds.z.grib

Also see https://confluence.ecmwf.int/pages/viewpage.action?pageId=146658047

[zhenkunl]

In the meantime, you could (as a temporary measure) convert these GRIB files to use a different packing scheme. For example CCSDS This is how you can do this:

grib_set -w isGridded=1 -r -s packingType=grid_ccsds gh850.grib gh850.ccsds.grib

Now apply the scaling and it will work ( no missing values issue )

grib_set -s scaleValuesBy=9.80665  gh850.ccsds.grib  gh850.ccsds.z.grib

Also see https://confluence.ecmwf.int/pages/viewpage.action?pageId=146658047

Yes, it works for me. Thanks for your prompt reply. Looking forward to your fix for this issue.