Loading a cube giving different AuxCoord behaviour?

[A-Gainford]

I’m still quite new to using iris so this might not even be the right question to ask, apologies. I am working with data which spans six hours, with three sets of data for each hour. I have some working code which reads in a given number of raw .pp files as:

rain_constraint = iris.AttributeConstraint(STASH='m01s04i201')
snow_constraint = iris.AttributeConstraint(STASH='m01s04i202')
ens_rain_cube = iris.load_cube(filename_list, rain_constraint, callback=ensemble_callback_oper)
ens_snow_cube = iris.load_cube(filename_list, snow_constraint, callback=ensemble_callback_oper)
ens_full_cube = ens_rain_cube + ens_snow_cube

Where the filename_list spans data from the entire six hour period. The callback simply adds a realization coord to the cube if it is not present in the input data. After adding the cubes for rain and snow together, it is able to construct a cube which looks something like this:

unknown / (kg m-2)                  (-- : 15; realization: 18; grid_latitude: 546; grid_longitude: 421)
    Dimension coordinates:
        realization                     -                x                  -                    -
        grid_latitude                   -                -                  x                    -
        grid_longitude                  -                -                  -                    x
    Auxiliary coordinates:
        time                            x                x                  -                    -
    Scalar coordinates:
        forecast_reference_time     2021-07-20 00:00:00
    Attributes:
        source                      Data from Met Office Unified Model
        um_version                  11.5

The time AuxCoord spans both the realization and the unnamed dimension, as desired.

I am trying to do the same for a sample of a much larger dataset, where keeping the full .pp files around would be impractical. So instead, after retrieving the data for the individual times, I have extracted the relevant fields using similar code as above and saved the resulting cube as smaller .nc files.

Now, however, when trying to load data for multiple times into a single cube using iris.load_cube(fnms), iris gives an error with:

ConstraintMismatchError: failed to merge into a single cube. Coordinates in cube.aux_coords (non-scalar) differ: time.

After trying to load just the three sets of data for a single time, I find that the cube looks something like this:

unknown / (kg m-2)                  (realization: 3; -- : 126; grid_latitude: 752; grid_longitude: 740)
    Dimension coordinates:
        realization                             x       -                   -                    -
        grid_latitude                           -       -                   x                    -
        grid_longitude                          -       -                   -                    x
    Auxiliary coordinates:
        forecast_period                         -       x                   -                    -
        time                                    -       x                   -                    -
    Scalar coordinates:
        forecast_reference_time     2019-06-28 19:00:00
    Attributes:
        Conventions                 CF-1.7
        source                      Data from Met Office Unified Model
        um_version                  11.5

where, now, it appears that the time coord is not able to be expanded into the dimension which has 126 values.

I have also tried to load individual cubes at different times into a single cube using CubeList.merge_cube() but run into the same issue. Is there any way that I can recreate the same behaviour during both loads?

Any help or suggestions would be very much appreciated. Many thanks.

[pp-mo]

You should be able to avoid that message "failed to merge into a single cube", if you used "iris.load" instead of "iris.load_cube".

I think what is going on is that the individual netcdf files contain multiple timepoints, so that they have a time dimension -- like your example immediately above

unknown / (kg m-2)     (realization: 3; -- : 126; grid_latitude: 752; grid_longitude: 740)

So, the key problem is that to combine several such cubes into a single one, you need to "concatenate" them, but iris loading can only "merge" cubes.
This is because "merge" works to create a new dimension from scalar coordinates, whereas "concatenate" combines cubes on an existing dimension -- see the distinction explained here.

( By contrast, PP files don't themselves contain a time dimension : each 2D field loads as a separate "raw cube", and the time dimensio is contructed by merging these, during loading )

To be frank, this is an annoying and ultimately unnecessary distinction, but it is where we are at present with our legacy code.
I know that work is currently under way to fix this, but it will still be a while coming, we have a very longstanding issue about it !

In the meantime, you may be able to use something like iris.load(files, constraint, callback=x).concatenate_cube().

Does that help at all ?

[A-Gainford]

Thanks for your suggestions and explanations.

I added a callback to remove any unnecessary coords which prevent concatenating ("forecast_period" and "forecast_reference_time") but it still returned an error saying that it couldn't return a single cube.

It seems like the problem was the dimensions of the coord('time').points.data and coord('time').bounds.data. The previous cube which was loaded directly from .pp files had shapes (18,15) for the coord('time').points/bounds.data, meaning that there was an array of 15 time values associated with each of the 18 realizations. This seems to be the key difference between the two ways of loading - concatenate_cube() was not able to construct a new AuxCoord for 'time' which was expanded over the 'realization' dimension.

So what I needed was for the shape of the coord('time').points.data for each cube in the CubeList returned by iris.load() to have the shape (x,3), where x is the number of values in a single coord('time').points or bounds data array, and 3 is the number of realizations loaded into each cube (which all have the same data for coord('time')).

I appear to have fixed this by forcing the points and bounds in 'time' to be the correct shape for my data (using 3 copies of the input points and bounds for each of the 6 cubes in the CubeList) using the following function:

def expand_time_coord_dimension(cubelist):
    for cube in cubelist:
        cube_time = cube.coord('time')
        cube_new_time = iris.coords.AuxCoord(
            points = [cube_time.points, cube_time.points, cube_time.points],
            bounds = [cube_time.bounds, cube_time.bounds, cube_time.bounds],
            standard_name = cube_time.standard_name,
            units = cube_time.units,
            var_name = cube_time.var_name,
        )
        cube.remove_coord('time')
        cube.add_aux_coord(cube_new_time, data_dims=[0,1])

So now I'm able to construct the cube that I want using:

ens_cube = iris.load(fnms, callback=remove_forecast_period_and_ref_time_callback)
expand_time_coord_dimension(ens_cube)
concat_cube = ens_cube.concatenate_cube()

which gives the following cube:

unknown / (kg m-2)                  (realization: 6; -- : 126; grid_latitude: 752; grid_longitude: 740)
    Dimension coordinates:
        realization                             x       -                   -                    -
        grid_latitude                           -       -                   x                    -
        grid_longitude                          -       -                   -                    x
    Auxiliary coordinates:
        time                                    x       x                   -                    -
    Attributes:
        Conventions                 CF-1.7
        source                      Data from Met Office Unified Model
        um_version                  11.5

and I can add a common 'forecast_reference_time' back to the cube to make it exactly the same as the .pp cube (ignoring the different dimensions assoicated with 'realization' and 'unknown', which shouldn't present an issue).

[A-Gainford]

Actually I have jumped the gun a bit with the answer here: this worked for the few files which I tested but it does not work for a full set of 18 realizations. (Not sure whether this should be a seperate question or not)

I'm not completely sure whether this is the issue or not, but I think the realization data is now a problem. In the forecast outputs, the '0' member is sometimes run out of order with the rest of the realizations: e.g, the data at T=22:00 will have realizations [15, 16, 17] and the data at T=23:00 will have realizations [0, 18, 19].

When loading data directly from .pp files with load_cube(), the output cube is able to rearrange this to be monotonic in the realization. concatenate() does not seem to be able to do this, and can only return two cubes: one with realizations [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 20, 21, 22], and the other with realizations [0, 18, 19].

Is there any way to tell concatenate() to rearrange the data so that it can return a single cube?

[trexfeathers]

Is there any way to tell concatenate() to rearrange the data so that it can return a single cube?

concatenate() does indeed need help in this case. We are working on an experimental alternative, so this may not be a problem in future.

[rcomer]

One difference between your two examples is that only the second has a “forecast_period” coordinate. This seems odd to me if they come from the same data set. I think Iris gets the “forecast_period” from the lbft header on the pp file.

[pp-mo]

Archiving "answered" Q+As