Add support for rotated lat/lon (Harmonie)

[bartdevriendt]

Added support for rotated lat lon grid as used in the harmonie weather model

[nmangue]

Thank you for your contribution ! Is it possible for you to share a grib file to add a test ?

[bartdevriendt]

Hi,

I added a test, and fixed some values in the wrong scale.

[nmangue]

Your changes add support to decode the grid definition. But it doesn't make it possible to read the grib file you added. To do so, we need to add :

• the computation of the point coordinates from the définition
• the decoding of the point values since they are encoded with a "complex" packing (only the simple one is currently implemented)

Would you like to work on these features ? If so I can explain how to implement them.

[bartdevriendt]

Hi,

I can give it a try.

[nmangue]

The first point is that this GRIB file uses a complex/second-order packing instead of a "simple" one. This currently raises an exception in Grib1BinaryDataSection.cs@89. The code checks for grid-point data an simple packing flags (see Code table 11 in the GRIB 1 specifications). Instead of that we need to branch the execution and run a distinct procedure to decode the second-order packing data. The decoding still needs to ensure that the file contains grid-point data. But based on the packing flag, it should call the right method. You can use the NetCDF code as a reference for the complex packing decoding : readExtendedComplexPacking. The algorithm is intimidating but it should be adaptable to NGrib without to many changes. To check that the code is valid you can call reader.ReadRecordRawData(record) and it should return without error and the excepted values (that you can retrieve using a reference tool such as Panoply)

The second point is the grid definition in itself. The decoding is implemented in Grib1GridDefinitionSection.EnumerateGridPoints. The current implementation doesn't support a Scan mode other than North to South and West to East (ScanMode != 0). The dataset has a scan mode of 64 meaning that the points are in the +j direction (South to North) according to Code table 8. You have to change the sign of yStep based on this flag (see yStep definition) :

var yStep = ((ScanMode & 64) == 0 ? -1 : 1) * Dy;

It should be enough to get the rotated coordinates. NB : The returned coordinates as expressed in the projection defined in the GRIB file.

Please let me know if you need more guidance.

Reference : GRIB edition 1 WMO Specifications

[bartdevriendt]

Hi,

I updated the code:

• fixed the scanmode direction for the y direction
• fixed the grid reading skipping wrong number of bytes
• added the complex reading code, but put in comment as it is not needed. The bug in the grid reading caused the byte position to be incorrect causing the simple packing flag to be incorrect (harmonie thus is using simple packing)
• added a test to validate the data reading

[nmangue]

Hello, everything looks fine. Thank you again for your work. I will publish an new version with your changes.