Cube creation

[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?

[lbdreyer]

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

Iris stores the time points as numbers since an epoch, which is the standard NetCDF (and CF) representation for datetimes. For example, you could create a coordinate with points=(6, 12, 18) and units='hours since 1970-01-01 00:00:00' which would represent time points of 1970-01-01 06:00:00, 1970-01-01 12:00:00, 1970-01-01 18:00:00.

Xarray typically converts the numbers into numpy.datetimes64 (and discards the units), although you can switch this off by including decode_times=False when you open the file (e.g. with xr.open_dataset).

However, if you need to convert already loaded numpy.datetimes64 objects back into numbers, you can do this as follows:

1. Set up an appropriate unit

from cf_units import Unit, CALENDAR_STANDARD
time_unit = Unit('hours since 1970-01-01 00:00:00', CALENDAR_STANDARD) # you can use the same unit

2. Get the values of your time variable, as datetime objects:

datetimes = da['time'].dt.date .values

3. Convert the datetime objects to numbers

time_points = time_unit.date2num(datetimes)

4. And then you could create your time coordinate:

time = DimCoord(time_points, standard_name='time', units=time_unit )

2.) Where to add the grib2 information in the Cube.

Iris follows CF as its data model, so you will need to represent your metadata following these conventions. Then at save time, Iris will then convert this metadata into GRIB constructs.

As an example, below we create a cube:

from cf_units import Unit, CALENDAR_STANDARD
import iris
from iris.coord_systems import GeogCS
from iris.coords import DimCoord
from iris.cube import Cube
import numpy as np

time = DimCoord([12], standard_name='time', units=Unit('hours since 1970-01-01 00:00:00', CALENDAR_STANDARD) )
height = DimCoord([10], standard_name='height', units='m')
# Set up a coordinate system for our latitude and longitude coords
cs = GeogCS(6371229)
latitude = DimCoord([0, 30, 60], standard_name='latitude', units='degrees', coord_system=cs)
longitude = DimCoord([-60, 0, 60],standard_name = 'longitude', units='degrees', coord_system=cs)
cube = Cube(np.ones((1,1,3,3)), standard_name='air_temperature', units='K', dim_coords_and_dims=[ (time, 0), (height, 1), (latitude, 2), (longitude, 3)])

When we print out the cube it looks like:

air_temperature / (K)               (time: 1; height: 1; latitude: 3; longitude: 3)
    Dimension coordinates:
        time                             x          -            -             -
        height                           -          x            -             -
        latitude                         -          -            x             -
        longitude                        -          -            -             x

We then save out the cube:

iris.save(cube, 'temp_data.grib2')

If you have a look at a grib dump of the file:
$ grib_dump -O temp_data.grib2

You will see that:

• As the cube had a coordinate system of GeogCS, a gridDefinitionTemplateNumber = 0 (Latitude/longitude) has been used
• As the cube represents a simple forecast (no realisation coordinate, no special statistical processing) productDefinitionTemplateNumber = 0 (Analysis or forecast at a horizontal level or in a horizontal layer at a point in time) is used
• As the cube had a standard_name='air_temperature', this is used to set: discipline=0 (Meteorological products), parameterCategory = 0 (Temperature), and parameterNumber = 0 (Temperature (K))
• As the z-axis coordinate is height, Iris sets typeOfFirstFixedSurface = 103 (Specified height level above ground (m))

3.) What Cube needs to write the grib2 file?

As demonstrated above, you should be able to save out your Cube as follows
iris.save(cube, 'temp_data.grib2')
Note the extension .grib2 used in the filename. This is used by iris to know that you want it to be saved out to GRIB, as opposed to a different file format such as netCDF. More information on saving can be found here: https://scitools-iris.readthedocs.io/en/latest/userguide/saving_iris_cubes.html

It is worth noting that Iris doesn't support all GRIB product/grib definition templates (there are quite a few!) so whether it can save out your data or not will depend on what is contained.