Grib to Netcdf conversion

Hi Marston,

The command
grib_set -v -s editionNumber=1 ICMGGgruz+201701 ICMGGgruz+201701.grb ICMGGgruz+201701

is a 'trick' to make any decoding software think it's dealing with GRIB-1 instead of GRIB-2. This normally works even though the rest of the GRIB messages are in GRIB-2. It's sometimes needed for older software that doesn't know how to handle GRIB-2.

However, one limitation with GRIB-1 is only 8 bits are available to encode the number of levels. This is why ECMWF use GRIB-2 for the upper air data, because 137 model levels is too many for GRIB-1. The 276 comes from 2*(137+1) which is the number of A & B coefficients on the model half-levels.

You might find you don't need to do this step if the software you use for the netcdf conversion is using ECMWF grib_api or eccodes or understands GRIB-2

  Glenn

This is good thing to know but I there is more to this than revealed here.
If I skip this grib_set part then:

cdo -v -O -R copy ICMGGgruz+201701 ICMGGgruz+201701_R.grb

Error (decodeIS) : GRIB version 2 unsupported!

OpenMP:  num_procs = 32  max_threads = 1

This thread should be read together with the "how to article" on grib conversion. I've posted some pertinent  issues I've experienced when using CDO and grib_api to convert grib to netcdf AND getting the data on a regular lat lon.

If I use grib_api 1.23.0 to convert the grib to netcdf:

grib_to_netcdf -k 2 -D NC_FLOAT -o ICMGGgruz+201701.nc ICMGGgruz+201701

the file conversion just freezes.
I've tried this grib_api only conversion to netcdf but the problem becomes that I only get a RGG netcdf, which I cannot use with NCL, Python, or with CDO.

Just an update:

grib_to_netcdf -k 2 -D NC_FLOAT -o ICMGGgruz+000000.nc ICMGGgruz+000000
grib_to_netcdf: Version 1.23.0
grib_to_netcdf: Processing input file 'ICMGGgruz+000000'.
grib_to_netcdf: Found 145 GRIB fields in 1 file.
grib_to_netcdf: Ignoring key(s): method, type, stream, refdate, hdate
grib_to_netcdf: Creating netCDF file 'ICMGGgruz+000000.nc'
grib_to_netcdf: NetCDF library version: 4.3.2 of Feb  3 2015 10:33:48 $
grib_to_netcdf: Creating large (64 bit) file format.
GRIB_API ERROR   :  First GRIB is not on a regular lat/lon grid or on a regular Gaussian grid. Exiting.

So this kind of describes the circle I've been going around when it comes to OPenIFS grib to netcdf.
It seems CDO is the only regular software out there that can convert from RGG to RLL. It is ironic that even grib_api seems to require that the grib file is RLL and yet DOES NOT provide this function at all. CDO have been using cgribex, which calls 2 emos functions, to carry out this conversion, BUT cgribex cannot handle the grib types in OpenIFS grib files. So grib_api is not a full replacement for cgribex in CDO.

Hi Marston,

Issues

To summarise first:

1. grib_to_netcdf does not do any interpolation from reduced Gaussian (varying number of longitude points per latitude) to regular grid and is not meant to. As the man page for grib_to_netcdf states, it only works with regularly gridded data (same number of longitude points on each latitude).

The ECMWF software to do the interpolation would be EMOSLIB. We'll only consider CDO here but I can add more about EMOSLIB if you'd like.

2. CDO's -R option, to convert from reduced Gaussian grid to regular grid, only works with GRIB 1 format.  This seems to be a limitation of CDO, because if as you say they use EMOSLIB, the EMOSLIB functions work with GRIB 1 and 2.

3. The trick of setting 'edition=1' in GRIB-2 messages in order to use 'cdo -R' doesn't work with 137 levels as the total number of half levels can't be encoded into the 8 bit value provided by GRIB-1.

Solution

Assume a gridpoint GRIB file with 137 level data in GRIB-2 messages: ICMGGgruzINIUA.

When trying to set the GRIB edition number to 1, it's the 'pv' GRIB parameter that causes the problem (nothing to do with potential vorticity!). This is an array holding the half level values of the model's A & B coefficients that define the location of the levels.

The workaround is to delete this array while setting the edition number to 1, like this:

grib_set -s deletePV=1,edition=1 ICMGGgruzINIUA ICMGGgruzINIUA_noPV.grib1

The 'deletePV' option here is known as a 'concept' rather than being an entry in the grib file. It ensures the PV array is deleted correctly.

Now cdo will be able to do the conversion to a regular lat-lon grid:

cdo -R copy ICMGGgruzINIUA_noPV.grib1 ICMGGgruzINIUA_noPV_ll.grib1

However, since these are really GRIB-2 messages, it is best to change the edition number back, otherwise the parameter names may not be recognised when converting to netcdf. This avoids the use of -t ecmwf option to cdo which is not recommended for GRIB_2 messages:

grib_set -s edition=2 ICMGGgruzINIUA_noPV_ll.grib1 ICMGGgruzINIUA_noPV_ll

Now there's a choice in converting to netCDF.

If cdo has been compiled with the ECMWF eccodes or grib_api libraries then this will work:

cdo -f nc copy  ICMGGgruzINIUA_noPV_ll ICMGGgruzINIUA_noPV_ll.nc

or you can use grib_to_netcdf from the ECMWF eccodes/grib_api distribution if cdo is not using eccodes/grib_api:

grib_to_netcdf -D NC_FLOAT -o ICMGGgruzINIUA_noPV_ll.nc ICMGGgruzINIUA_noPV_ll

where the -D NC_FLOAT option ensures the parameters are encoded as real numbers instead of scaled integers with an offset.

Both these options work for me with no errors. I'm using cdo 1.7.2 compiled with grib_api 1.15.0, cgribex 1.7.5 and netCDF 4.4.1.  For grib_to_netcdf I'm using eccodes/2.1.0

Limitation
Because the 'pv' array was deleted from the GRIB messages, it means the resulting netCDF file will not have any record of the half level model coordinates.

This may not be important but if it is, it's better to use EMOSLIB to do the reduced grid to regular grid conversion. There is a way of injecting the half levels back into the regularly gridded GRIB file using cdo but it's not particularly straightforward.

Please let me know if this works for you. I can then update the information on the converting grib to netcdf page.

Thanks for highlighting this problem.

And many thanks to Paul Dando for spending time looking into this (and other options).

   Glenn

Hi Guys,

Thanks for this work around. I tested it this morning on small sample and it worked. I will need the A and B levels later on, but I can get these offline, or can I?

A few things:

1.) With regards to CDO, they do not automatically update the grib table in cgribex, therefore as of the latest CDO version, conversion of OpenIFS grib to netcdf using CDO results in the names being mislabeled. Having spoken to CDO about this, there recommendation is that one should turn of cgribex and use grib_api as the grib decoder engine. This is works well for the variable name translations but results in other problems. See CDO issues: #7813 and #7797.

2.) I think a solution might be to convert the EMOS routines: qu2reg3 and rowina3 to python routines, read the grib files using ECMWF python interface to grib_api and proceed from there on. If you could give me a brief explanation of how these 2 EMOSLIB routines work, I would appreciate it.

Hi Marston,

Glad it worked.

To extract the A & B coefficients, you can either use cdo:

cdo vct in.grib > pv.txt

or with grib_api, you need to use a grib_filter command as it's not possible to do it with 'grib_get':

cat<<EOF>print_pv
print '[pv!1%.6f]' ;
EOF
grib_filter print_pv in.grib

The A's are in the first L+1 elements and the B's in elements L+2 to 2L+2. Or use a Fortran or C code to read the grib directly with the grib_api library.

On your points 1 & 2 above:

   Cheers,  Glenn (with thanks again to Paul Dando)

Hi Glenn,

Thanks again for the help. It's kind of pain to have so many workarounds for each grib type. But I'm glad to be able to move on from this. With regards to the need to control the transformation, I'm trying to track down some artifacts in my COSP output from FullPOS. The issue is that I do not do my post-processing in fortran. I use NCL/Python and that is why I'm so dependent upon tools like CDO.
I do not have the time left in my project to move my processing chain to fortran so I'll shelf this idea
Good to know that I can still access the A/B values.

Hi Glenn,

This is hopefully a final update on this issue. The CDO guys have put forward a solution that is elegant and a template for the future:

The CDO option -R to convert RGG to RLL is implemented in CDO cgribex library and works only for GRIB1 data. There is a CDO operator with the same feature but not the same performance:

cdo setgridtype,regular infile outfile

Hi Marston,

You were asking about using ECMWF EMOS to do the interpolation to a regular grid.

Assuming a recent version of the EMOS library, there are two tools that can be used.

In the 'bin' directory there is a command 'emos_tool' which can be used as:

emos_tool --regular=256 [--area="40/-10/31/20"]  in.grib out.grib

The grid specified for the --regular option follows the grid naming convention for EMOSLIB here: Reduced Gaussian Grids.

Another command can be found in the 'tools' directory. This allows the interpolation function to be specified:

tools/int --INTOUT:gridname=F256 --INTOUT:area="40/-10/31/20" --input=1.grib --output=2.grib --intf2

Here 'F256' means 'full-grid' not reduced, and --intf2 is the EMOSLIB interpolation function to be used.

Hope that helps,

                    Glenn

Hi Glenn,

Thanks for this update

/M