Support saving to GDT 40

[lbdreyer]

Currently iris-grib supports loading GDT 40 (see here) but not saving it.

We would like to write out data on an octahedral reduced gaussian grid.

I can provide some example files of what "good" looks like.

[trexfeathers]

Thanks @lbdreyer, we would love to implement this (previously discussed in email 2026-03-11). Do you have any deadlines we should be aware of?

[lbdreyer]

Unfortunately, the best I can give you is that it is needed for PS48. Hopefully, that is enough to go on?

I have already provided @pp-mo with the patch to iris-grib that I have been using so far in case that helps with development

[pp-mo]

Update 2026-06-24

Now working on this!

For the sake of being definite, it seems the key purpose (usecase) here is actually round-tripping data on a reduced-gaussian grid.
Some things I think we know:

• there is nothing 'special' about the octahedral grid, as far as the grib message is concerned :
  it just means a particular choice of initial latitude offset and number-of-points per parallel. Detail here
• we need to save on GDT40, as we already load
• with the test data, we already observed another problem, in that it loads without a vertical coordinate
  • this is either because the source data is missing surface-presure fields, or possibly because iris-grib can't handle a 1D reference field (TBD)
• in loading, we convert to an iris cube with 1 dimension, lat+lon aux-coords mapping it, and a specific pattern of points within those coords
  • E.G. LIKE lats = [3*-30, 5*-10, 5*+10, 3*+30] and lons = [0, 120, 240, 0, 72, 144, 216, 288, 0, 72, 144, 216, 288, 0, 120, 240]
• in configuring a GDT40 message, the actual latitude points are not specified, but entirely determined by the choice of 'N' (octets 68-71)
  • this means that in saving from Iris the exact values could be "lost", and ideally we would check they exactly match the expected values
  • .. this is also true of longitudes, but there the calculation is very simple.
  • specifically in the case of latitudes, the calculation is not simple : you have to solve a D.E. to get them
    • like this code in eccodes
    • for load, we have hitherto chosen not to do that in iris-grib, but instead use the ecccodes "distinctLatitudes" computed key to calculate values in eccodes
• there is also a possible problem with saving to regular gaussian grids :
  • this presently 'works', but saves as an non-regular lat/lon grid (GDT3.4)
  • .. so if we need it to save as it was loaded, more control may be needed

Currently I still have some outstanding unresolved problems with fixing this:

• cannot manage to create a sample message from the "GR1B2" template, but apparently need to start from the "reduced_gg_sfc_grib2"
  • as seen here (for those that can view that)
• even when using the "reduced_gg_sfc_grib2" template as a basis, I still cannot use codes_get_array(.., "pl") on an in-memory message
  • .. though it does works if I save the same message, unchanged, to a file and then load it back in.

[sawom666]

Re: Starting from the sample GRIB2.tmpl and creating a reduced Gaussian grid
You need to do the following (shown here using grib_filter)

% cat > f <<EOF
  set gridType = "reduced_gg";
  set pl = {4 ,5 ,6}; # just some dummy values
  # Of course you also have to set other related keys e.g., N, numberOfDataPoints etc
  write;
EOF
% grib_filter -o rgg.grib2 f samples/GRIB2.tmpl
% grib_ls rgg.grib2

Now it should work. You can do the above using C, Fortran or Python too
:)

[sawom666]

Note: Once you have encoded your GRIB message, it is always a good idea to check it i.e., run a set of basic sanity checks.
This can be done via the key isMessageValid. For example

% grib_get -p isMessageValid my.grib

If this returns 1, then the basic checks have passed. If it returns 0, you will see some error messages posted to stderr

[sawom666]

Re: Modifying and accessing the PL array from the sample file reduced_gg_sfc_grib2.tmpl
Let's try that in Python. We will use the in-memory version, (not saving to disk and re-reading)

import eccodes
gid = eccodes.codes_grib_new_from_samples("reduced_gg_sfc_grib2")
print( "BEFORE ", eccodes.codes_get_array(gid, "pl") )
eccodes.codes_set_array(gid, "pl", [21,22,33])
print( "NOW ", eccodes.codes_get_array(gid, "pl") )
eccodes.codes_release(gid)

Works for me ;)

[pp-mo]

@sawom666 thanks for your very fast responses !!

Re: Modifying and accessing the PL array from the sample file reduced_gg_sfc_grib2.tmpl
Let's try that in Python. ...

Yes, that does work!

But our existing code creates all saved messages from a simple "GRIB2" template, and that appears not to work...
E.G. ...

import eccodes

# NOTE: see below
# grib = eccodes.codes_grib_new_from_samples("GRIB2")
grib = eccodes.codes_grib_new_from_samples("reduced_gg_sfc_grib2")

eccodes.codes_set(grib, "gridDefinitionTemplateNumber", 40)

eccodes.codes_set(grib, "shapeOfTheEarth", 6)
eccodes.codes_set(grib, "latitudeOfFirstGridPoint", 89284228)  # in micro degrees
eccodes.codes_set(grib, "longitudeOfFirstGridPoint", 0)

eccodes.codes_set(grib, "N", 96)
eccodes.codes_set(grib, "Nj", 192)
eccodes.codes_set_missing(grib, "Ni")

eccodes.codes_set(grib, "scanningMode", 0)

eccodes.codes_set_long(grib, 'basicAngleOfTheInitialProductionDomain', 0)
eccodes.codes_set_long(grib, 'resolutionAndComponentFlags', 0)
eccodes.codes_set_missing(grib, 'iDirectionIncrement')

eccodes.codes_set(grib, "numberOfOctectsForNumberOfPoints", 2)
eccodes.codes_set(grib, "interpretationOfNumberOfPoints", 1)

# Create desired latitude sequence for O96
#  --> [20, 24, ... , 396, 400, 400, 396, ... 24, 20]
pl = list(range(20, 401, 4))
pl = pl + pl[::-1]
print(len(pl), pl)
eccodes.codes_set_array(grib, "pl", pl)

That works fine as written -->

192 [20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68, 72, 76, 80, 84, 88, 92, 96, 100, 104, 108, 112, 116, 120, 124, 128, 132, 136, 140, 144, 148, 152, 156, 160, 164, 168, 172, 176, 180, 184, 188, 192, 196, 200, 204, 208, 212, 216, 220, 224, 228, 232, 236, 240, 244, 248, 252, 256, 260, 264, 268, 272, 276, 280, 284, 288, 292, 296, 300, 304, 308, 312, 316, 320, 324, 328, 332, 336, 340, 344, 348, 352, 356, 360, 364, 368, 372, 376, 380, 384, 388, 392, 396, 400, 400, 396, 392, 388, 384, 380, 376, 372, 368, 364, 360, 356, 352, 348, 344, 340, 336, 332, 328, 324, 320, 316, 312, 308, 304, 300, 296, 292, 288, 284, 280, 276, 272, 268, 264, 260, 256, 252, 248, 244, 240, 236, 232, 228, 224, 220, 216, 212, 208, 204, 200, 196, 192, 188, 184, 180, 176, 172, 168, 164, 160, 156, 152, 148, 144, 140, 136, 132, 128, 124, 120, 116, 112, 108, 104, 100, 96, 92, 88, 84, 80, 76, 72, 68, 64, 60, 56, 52, 48, 44, 40, 36, 32, 28, 24, 20]

Process finished with exit code 0

but if I replace the "reduced_gg_sfc_grib2" with a plain "GRIB2", then -->

/data/apps/sss/environments/default-2026_02_12/bin/python /home/users/patrick.peglar/worklog/period_20260525_iris3x16/dev/grib_gdt40_save/template_gdt40_minimal.py 
Traceback (most recent call last):
  File "/home/users/patrick.peglar/worklog/period_20260525_iris3x16/dev/grib_gdt40_save/template_gdt40_minimal.py", line 30, in <module>
    eccodes.codes_set_array(grib, "pl", pl)
  File "/data/apps/sss/environments/default-2026_02_12/lib/python3.12/site-packages/gribapi/gribapi.py", line 2174, in grib_set_array
    grib_set_long_array(msgid, key, value)
  File "/data/apps/sss/environments/default-2026_02_12/lib/python3.12/site-packages/gribapi/gribapi.py", line 1293, in grib_set_long_array
    GRIB_CHECK(lib.grib_set_long_array(h, key.encode(ENC), inarray, len(inarray)))
  File "/data/apps/sss/environments/default-2026_02_12/lib/python3.12/site-packages/gribapi/gribapi.py", line 226, in GRIB_CHECK
    errors.raise_grib_error(errid)
  File "/data/apps/sss/environments/default-2026_02_12/lib/python3.12/site-packages/gribapi/errors.py", line 381, in raise_grib_error
    raise ERROR_MAP[errid](errid)
gribapi.errors.KeyValueNotFoundError: Key/value not found
192 [20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68, 72, 76, 80, 84, 88, 92, 96, 100, 104, 108, 112, 116, 120, 124, 128, 132, 136, 140, 144, 148, 152, 156, 160, 164, 168, 172, 176, 180, 184, 188, 192, 196, 200, 204, 208, 212, 216, 220, 224, 228, 232, 236, 240, 244, 248, 252, 256, 260, 264, 268, 272, 276, 280, 284, 288, 292, 296, 300, 304, 308, 312, 316, 320, 324, 328, 332, 336, 340, 344, 348, 352, 356, 360, 364, 368, 372, 376, 380, 384, 388, 392, 396, 400, 400, 396, 392, 388, 384, 380, 376, 372, 368, 364, 360, 356, 352, 348, 344, 340, 336, 332, 328, 324, 320, 316, 312, 308, 304, 300, 296, 292, 288, 284, 280, 276, 272, 268, 264, 260, 256, 252, 248, 244, 240, 236, 232, 228, 224, 220, 216, 212, 208, 204, 200, 196, 192, 188, 184, 180, 176, 172, 168, 164, 160, 156, 152, 148, 144, 140, 136, 132, 128, 124, 120, 116, 112, 108, 104, 100, 96, 92, 88, 84, 80, 76, 72, 68, 64, 60, 56, 52, 48, 44, 40, 36, 32, 28, 24, 20]

Process finished with exit code 1

So, any ideas why, and do I just need to set something else ?

It's important because all of our current save code does start with grib_new_from_samples("GRIB2")
and changing that makes the save code rather tangled, as it breaks the neat separation of concerns in the existing code structure.

[sawom666]

I will take a look

[sawom666]

OK. I have rewritten your script and now it works. This time we start from GRIB2.tmpl and modify it accordingly

import eccodes

grib = eccodes.codes_grib_new_from_samples("GRIB2")

eccodes.codes_set(grib, "gridType", "reduced_gg")

eccodes.codes_set(grib, "shapeOfTheEarth", 6)
eccodes.codes_set(grib, "latitudeOfFirstGridPoint", 89284228)  # in micro degrees
eccodes.codes_set(grib, "longitudeOfFirstGridPoint", 0)
eccodes.codes_set(grib, "N", 96)
eccodes.codes_set(grib, "Nj", 192)
eccodes.codes_set_missing(grib, "Ni")
eccodes.codes_set(grib, "scanningMode", 0)
eccodes.codes_set_long(grib, 'basicAngleOfTheInitialProductionDomain', 0)
eccodes.codes_set_long(grib, 'resolutionAndComponentFlags', 0)
eccodes.codes_set_missing(grib, 'iDirectionIncrement')
eccodes.codes_set(grib,` "interpretationOfNumberOfPoints", 1)

# Create desired latitude sequence for O96
#  --> [20, 24, ... , 396, 400, 400, 396, ... 24, 20]
pl = list(range(20, 401, 4))
pl = pl + pl[::-1]
print(len(pl), pl)
eccodes.codes_set_array(grib, "pl", pl)
### save
OUTPUT = "out.grib2"
fout = open(OUTPUT, "wb")
eccodes.codes_write(grib, fout)
fout.close()

Some points:
1/ You can directly set the key gridType rather than the cryptic gridDefinitionTemplateNumber
2/ No need to set numberOfOctectsForNumberOfPoints

Hope this helps
:)

[pp-mo]

OK. I have rewritten your script and now it works. This time we start from GRIB2.tmpl and modify it accordingly

OK thanks for this -- it does work!

But I did find that it doesn't work if again I replace
eccodes.codes_set(grib, "gridType", "reduced_gg")
with
eccodes.codes_set(grib, "gridDefinitionTemplateNumber", 40)
Then I get
gribapi.errors.KeyValueNotFoundError: Key/value not found

The reason I want to do that is because in iris-grib we try AFAP to avoid the use of computed keys and only access grib data in the terms belonging to the definitive WMO spec.
So, I would say that although "gridDefinitionTemplateNumber" is cryptic from a Python programmer POV, it is the correct term which relates directly to what you see in a file (e.g. from "grib_dump -O").

In this case, though, it clearly isn't just a name alias : the 'gridType' setting operation is doing something else inside eccodes which then allows it to set the pl values : in fact, it seems to indicate that 'pl' is not a valid key, as if the message has not been fully "switched" to the desired template definition ?

My 'other' problem

While this code (and the very similar code here in my draft PR) do allow the message to save,
I also find that I still can't perform a eccodes.codes_get_array(grib, "distinctLatitudes") call.
I get gribapi.errors.WrongGridError: Grid description is wrong or inconsistent.

Possibly this can be fixed by setting addtional keys, but I have found that if I write the message to a file, and then read it back, I can perform that operation. So, possibly, this just can't be done with an in-memory message created from a template ??

The reason I would like to do that, is to check that it is 'correct' to write a given Iris cube as a GDT40 message.
At present, I am checking that the data dimensions and coordinate structure "imply" that it can be saved to GDT40.
But since Iris cubes have full coordinate information, it is then still possible that if we write to GRIB GDT40 and read back, the latitude values might then not be the same.
It's tricky because computing the values is highly non-trivial, so it would seem best to rely on eccodes for that (as detailed here).

So, It would be great if you can cast some light on these problems!