Cube creation for grib2 dump

[Marston]

I have a case study to create a grib2 from a xarray dataArray. But I'm having trouble creating a proper cube with the information needed to write to grib2. I have not yet tried to dump the cube to grib2 because I'm stuck in the following areas:

import xarray as xr
import numpy as np
import pandas as pd
import iris
from iris.coords import DimCoord
from iris.time import PartialDateTime
from iris.cube import Cube

def build_cube(da, varbname, units):
   """ 
   Build an Iris cube from a dataArray. Doing this manually due to lack of proper metadata in the dataArray.
   """

    dt =da.coords['time']
    pdt = PartialDateTime(year=dt.year, month=dt.month, day=dt.day, hour=dt.hour, minute=dt.minute)
    time = DimCoord(pdt,
                    standard_name='time')
    lev = DimCoord(da.coords['lev'].values, standard_name='lev', units='Height MSL')
    latitude = DimCoord(da.coords['lat'].values,
                        standard_name='latitude',
                        units='degrees')
    longitude = DimCoord(da.coords['lon'].values,
                         standard_name='longitude',
                         units='degrees')
    cube = Cube(da.values, 
                standard_name=varbname.values,
                long_name=varbname.values,
                var_name=varbname.values,
                units='K',  
                attributes=None, # Will add variable attributes here. 
                dim_coords_and_dims=[ (time, 0), (lev, 1), (latitude, 2), (longitude, 3) ] 
    )    
    return out

1.) I'm having trouble understanding Cube dimension time and how to construct it.
2.) Where to add the grib2 information in the Cube.
3.) What Cube needs to write the grib2 file?

[pp-mo]

1.) I'm having trouble understanding Cube dimension time and how to construct it.

I've not looked at your code in any great detail, but you seem to be doing everything for yourself.

Certainly, building a coordinate out of a PartialDateTime is not going to work.
PartialDateTime is only provided to assist comparison operations. Iris requires its coordinates to be encoded as pure time offsets (floating-point numbers), with the calendar and base-time/date encoded in the 'units' parameter.
E.G. something like DimCoord(day_numbers, standard_name='time', units=cf_units.Unit('days since 1970-01-01', calendar='standard'))

Anyway ... you should certainly try out the xarray DataArray.to_iris method, as mentioned here in xarray docs on 'other libraries'
Although that is a bit limited, and may not deliver exactly what you want.
There is another initiative in the offing to make this translation a bit easier and more comprehensive, but it is still a little way off.

---

2.) Where to add the grib2 information in the Cube.
3.) What Cube needs to write the grib2 file?

In principle, once you have a suitable Iris cube it will just save to a grib file on command.
But obviously it relies on being encoded in the right way.
As your data appears to be 2-D with no vertical coordinate, there shouldn't be a lot of complexity in that.
The "phenomenon type" is generally translated from the cube 'standard_name', and possibly its cell-methods, into grib parameter code.
However, if there is no translation in the table, you can add a GRIB_PARAM attribute to control this on output.
Unfortunately there is no user documentation on this at present, but you can see in the code how this is implemented. More advice available if you need it !

[Marston]

@pp-mo

Thank you for the response and the links. I really like the idea of Iris to grib2 (as much as it is not my favorite file format), and I'm trying to learn to encode the dataArray so that I can obtain a proper Cube or do the necessary content enrichment before I dump to grib2. I was trying to build the Cube from scratch to learn but without good docs, that is a much harder way.
I have all the info to achieve this, but I cannot find good examples, as you pointed out. I truly appreciate your advice. Below is a screenshot of my current Cube. As you can see it is 3D having 31 levels. I know grib will store these as individual records, but I need to encode each record properly. I can use your help with this. Some pseudo code perhaps on how to encode the Cude properly or modify it by adding the missing information? For example, I see from the iris-grib code that the GRIB_PARAM string can be added to the Cube's global properties. But how do I label, the vertical levels as having height MSL and not isobaric, e.g?

[Marston]

@pp-mo

So I've read and explored and came up a little short but I think with a little help I can get a working example across the finish line. Here's the idea:

• I read in a netcdf4 file into Xarray. I tried iris.fileformats.netcdf.loader.load_cubes but this resulted in an object that makes no sense to me and I could not go any further.
• I select a 3D array and loop over the levels.
• Add/append each level to a grib2 file where each file will contain all the records for the above 3D variable.
  Here is the code I'm developing:

data = xr.open_dataset(ifile, engine='netcdf4', group='galwem')
# see image below for: print(data) 
RE_PARSE_FOURNUMS = re.compile(4 * r'[^\d]*(\d*)')
gribEdition = 2
for i in range(1): # just take one variable, aka feature for now
    name = data.names3d.isel(features3d=i)
    varb = data.variables3d.isel(features3d=i)
    varb = varb.rename({'lat':'Y', 'lon':'X'})
    # using a local LUT for grib info on each variable
    df = gribInfo.where(gribInfo.Parameter == name).dropna()
    key = df.SearchKey.item()
    print(key)  # 0-2-10-L100-1
    print(RE_PARSE_FOURNUMS.match(key).groups()) # ('0', '2', '10', '100')
    key = f'{gribEdition}-{key}'
    for lev in range(1): 
        qslice = varb.isel(time=0, msllevels=lev)
        qslice.attrs={'long_name':str(name.values), 'GRIB_PARAM': key}
        print(name.values, qslice) 
        # (array('Absolute_vorticity', dtype=object),
        # <xarray.DataArray 'variables3d' (Y: 721, X: 1440)>
        # [1038240 values with dtype=float64]
        # Coordinates:
        # time   datetime64[ns] 2021-09-01T18:00:00
        # * Y   (Y) float32 90.0 89.75 89.5 89.25 ... -89.25 -89.5 -89.75 -90.0
        # * X   (X) float32 -180.0 -179.8 -179.5 -179.2 ... 179.2 179.5 179.8
        # msllevels  int32 110000
        # Attributes:
        # long_name:   Absolute_vorticity
        # GRIB_PARAM:  2-0-2-10-L100-1)
        
        # Convert to cube
        cube = qslice.to_iris()
        # trying to edit the cube
        #cube.coord(axis='X').coord_system=iris.coord_systems.GeogCS(654321)
        #cube.coord(axis='Y').coord_system=iris.coord_systems.GeogCS(654321)   
        #cube.add_aux_coord(iris.coords.DimCoord(0, standard_name='forecast_period', units='hours'))
        #cube.add_aux_coord(iris.coords.DimCoord(0, "height", units="ft"))
        #cube.data.max() 
        #cube.data=ma.masked_invalid(cube.data)
        
        print(cube)
        
        # Absolute_vorticity / (unknown)      (latitude: 721; longitude: 1440)
        # Dimension coordinates:
        # latitude                             x               -
        # longitude                            -               x
        # Scalar coordinates:
        # msllevels                   110000
        # time                        2021-09-01 18:00:00
        # Attributes:
        # GRIB_PARAM                  '2-0-2-10-L100-1'
        
        # test writing to grib2
        iris.save(cube, 'test.grib2')

        # File ~/mambaforge-pypy3/envs/contrail/lib/python3.11/site-packages/iris_grib/__init__.py:801, in save_pairs_from_cube(cube)
        # 799 y_coords = cube.coords(axis='y', dim_coords=True)
        # 800 if len(x_coords) != 1 or len(y_coords) != 1:
        # --> 801     raise TranslationError("Did not find one (and only one) x or y coord")
        # 803 # Save each latlon slice2D in the cube
        # 804 for slice2D in cube.slices([y_coords[0], x_coords[0]]):
        # TranslationError: Did not find one (and only one) x or y coord

So here is where I am. I have changed the dimensions but the conversion to cube converts the lat and lon dims from x, y to latitude and longitude. Not sure how to fix this. Lots of scattered code out there. Could use some help/advice.