Cubes won't merge but will concatenate after promoting scalar coordinate to singleton dimension

[david-bentley]

🐛 Bug Report

Loading some NetCDF files fail to merge into a single cube, but promoting the relevant scalar coordinates and concatenating does collapse into a single cube.

How To Reproduce

I have some NetCDF files which I would expect to merge into a single cube when loaded, but the cubes fail to merge. For example

>>> from pathlib import Path

>>> import iris

>>> p = Path("/path/to/files/")
>>> files = list(p.glob("*nc"))

>>> def callback(cube, field, filename):
...     del cube.attributes["history"]

>>> name_con = iris.Constraint(cube_func=lambda cube: cube.name().startswith("air_temperature"))
>>> coord_con = iris.Constraint(cube_func=lambda cube: any(coord.name() == "pressure"
...                                                        for coord in cube.dim_coords))

>>> cubes = iris.load(files, name_con & coord_con, callback=callback)
>>> print(cubes)
0: air_temperature / (K)               (time: 2; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
1: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)

If I load the cubes individually and then attempt to manually merge them I get a MergeError as you might expect:

>>> # now load raw and attempt to coerce them into a single cube
>>> cubes = iris.load_raw(files, name_con & coord_con, callback=callback)
>>> print(cubes)
0: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
1: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
2: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)

>>> # attempt to merge into a single cube
>>> print(cubes.merge_cube())
Traceback (most recent call last):
  ...
  File ".../lib/python3.12/site-packages/iris/cube.py", line 348, in merge_cube
    proto_cube.register(c, error_on_mismatch=True)
  File ".../lib/python3.12/site-packages/iris/_merge.py", line 1301, in register
    match = coord_payload.match_signature(
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".../lib/python3.12/site-packages/iris/_merge.py", line 271, in match_signature
    raise iris.exceptions.MergeError(msgs)
iris.exceptions.MergeError: failed to merge into a single cube.
  Coordinates in cube.aux_coords (non-scalar) differ: air_temperature status_flag.

Metadata wise the coordinates look compatible

>>> # examine the status_flag coordinate
>>> for cube in cubes:
...     print(cube.coord("air_temperature status_flag"))
AuxCoord :  air_temperature status_flag / (no_unit)
    points: <lazy>
    shape: (33, 970, 1042)
    dtype: int8
    standard_name: 'air_temperature status_flag'
    var_name: 'flag'
    attributes:
        flag_meanings  'above_surface_pressure below_surface_pressure'
        flag_values    [0 1]
AuxCoord :  air_temperature status_flag / (no_unit)
    points: <lazy>
    shape: (33, 970, 1042)
    dtype: int8
    standard_name: 'air_temperature status_flag'
    var_name: 'flag'
    attributes:
        flag_meanings  'above_surface_pressure below_surface_pressure'
        flag_values    [0 1]
AuxCoord :  air_temperature status_flag / (no_unit)
    points: <lazy>
    shape: (33, 970, 1042)
    dtype: int8
    standard_name: 'air_temperature status_flag'
    var_name: 'flag'
    attributes:
        flag_meanings  'above_surface_pressure below_surface_pressure'
        flag_values    [0 1]

If I promote the relevant scalar coordinate to a singleton dimension (along with other coordinates spanning that dimension) I can then successfully concatenate into a single cube

>>> # attempt to concatenate the cubes by promoting the scalar dimension to a
>>> # singleton dimension
>>> cubes = iris.cube.CubeList(iris.util.new_axis(cube, "time",
...                                               expand_extras=(cube.coord("forecast_period"),
...                                                              cube.coord("air_temperature status_flag")))
...                            for cube in cubes)
>>> print(cubes.concatenate_cube())
air_temperature / (K)               (time: 3; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
    Dimension coordinates:
        time                             x            -                            -                             -
        pressure                         -            x                            -                             -
        projection_y_coordinate          -            -                            x                             -
        projection_x_coordinate          -            -                            -                             x
    Auxiliary coordinates:
        forecast_period                  x            -                            -                             -
        air_temperature status_flag      x            x                            x                             x
    Scalar coordinates:
        forecast_reference_time     2025-10-01 00:00:00
    Attributes:
        Conventions                 'CF-1.7, UKMO-1.0'
        institution                 'Met Office'
        least_significant_digit     np.int64(1)
        mosg__forecast_run_duration 'PT54H'
        mosg__grid_domain           'uk_extended'
        mosg__grid_type             'standard'
        mosg__grid_version          '1.7.0'
        mosg__model_configuration   'uk_det'
        source                      'Met Office Unified Model'
        title                       'UKV Model Forecast on UK 2 km Standard Grid'
        um_version                  '13.1'

which is what I would expect to happen. It is unclear to me why the status_flag coordinate should prevent a merge but not prevent concatenation?

Expected behaviour

On load there would be a single cube with a time dimension of length 3.

Environment

• OS & Version: RHEL7
• Iris Version: 3.13.1

Additional context

I get the same result with all of the COMBINE_POLICY options:

# try the different load policies
>>> for policy in ("legacy", "recommended", "default", "comprehensive"):
...     iris.COMBINE_POLICY.set(policy)
...     cubes = iris.load(files, name_con & coord_con, callback=callback)
...     print(policy)
...     print(cubes)
legacy
0: air_temperature / (K)               (time: 2; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
1: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
recommended
0: air_temperature / (K)               (time: 2; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
1: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
default
0: air_temperature / (K)               (time: 2; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
1: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
comprehensive
0: air_temperature / (K)               (time: 2; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
1: air_temperature / (K)               (pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)

[stephenworsley]

from @SciTools/peloton :
We believe this is expected behaviour. The points of air_temperature_status_flag look like they differ on their points, if this is the case, merge is designed to fail, it will only expand the dimensions of scalar coordinates, not multidimensional AuxCoords. The promoting of the scalar coordinates here is absolutely the right way to solve this by explicitly signaling that the dimensions of air_temperature_status_flag are supposed to be expanded. If this was to be implicitly inferred from whether this coordinate varies over your cubes, then you will end up with inconsistent behaviour in cases where coordinates are frequently similar but sporadically different (as might be expected with status flags).

[david-bentley]

from @SciTools/peloton : We believe this is expected behaviour. The points of air_temperature_status_flag look like they differ on their points, if this is the case, merge is designed to fail, it will only expand the dimensions of scalar coordinates, not multidimensional AuxCoords. The promoting of the scalar coordinates here is absolutely the right way to solve this by explicitly signaling that the dimensions of air_temperature_status_flag are supposed to be expanded. If this was to be implicitly inferred from whether this coordinate varies over your cubes, then you will end up with inconsistent behaviour in cases where coordinates are frequently similar but sporadically different (as might be expected with status flags).

I think in this case the status_flag is inherently time dependent though, because it depends on the pressure at the given timestep. It may be that the first two timesteps happen to have the same values, but this is not true in general.

I'd certainly argue that (from a weather-timescale perspective) the status_flag is more likely to be time-dependent than orography, which can now span the time dimension when loaded.

[david-bentley]

Or looking at this from another perspective, here is (a portion of) the ncdump -h of one of the files:

netcdf \20251001T0000Z-PT0000H00M-temperature_on_pressure_levels {
dimensions:
	pressure = 33 ;
	projection_y_coordinate = 970 ;
	projection_x_coordinate = 1042 ;
	bnds = 2 ;
variables:
	float air_temperature(pressure, projection_y_coordinate, projection_x_coordinate) ;
		air_temperature:least_significant_digit = 1LL ;
		air_temperature:standard_name = "air_temperature" ;
		air_temperature:units = "K" ;
		air_temperature:ancillary_variables = "flag" ;
		air_temperature:grid_mapping = "lambert_azimuthal_equal_area" ;
		air_temperature:coordinates = "flag forecast_period forecast_reference_time time" ;
...
	byte flag(pressure, projection_y_coordinate, projection_x_coordinate) ;
		flag:standard_name = "air_temperature status_flag" ;
		flag:flag_meanings = "above_surface_pressure below_surface_pressure" ;
		flag:flag_values = 0b, 1b ;
	int forecast_period ;
		forecast_period:units = "seconds" ;
		forecast_period:standard_name = "forecast_period" ;
	int64 forecast_reference_time ;
		forecast_reference_time:units = "seconds since 1970-01-01 00:00:00" ;
		forecast_reference_time:standard_name = "forecast_reference_time" ;
		forecast_reference_time:calendar = "gregorian" ;
	int64 time ;
		time:units = "seconds since 1970-01-01 00:00:00" ;
		time:standard_name = "time" ;
		time:calendar = "gregorian" ;

So air_temperature is 3D, with dimensions pressure, x and y, and there are some associated coordinates (flag, forecast_period, forecast_reference_time, time) which do not span the dimensions of the data. The two cubes which do get merged give

air_temperature / (K)               (time: 2; pressure: 33; projection_y_coordinate: 970; projection_x_coordinate: 1042)
    Dimension coordinates:
        time                             x            -                            -                             -
        pressure                         -            x                            -                             -
        projection_y_coordinate          -            -                            x                             -
        projection_x_coordinate          -            -                            -                             x
    Auxiliary coordinates:
        forecast_period                  x            -                            -                             -
        air_temperature status_flag      -            x                            x                             x
    Scalar coordinates:
        forecast_reference_time     2025-10-01 00:00:00
    Attributes:
        Conventions                 'CF-1.7, UKMO-1.0'
        institution                 'Met Office'
        least_significant_digit     np.int64(1)
        mosg__forecast_run_duration 'PT54H'
        mosg__grid_domain           'uk_extended'
        mosg__grid_type             'standard'
        mosg__grid_version          '1.7.0'
        mosg__model_configuration   'uk_det'
        source                      'Met Office Unified Model'
        title                       'UKV Model Forecast on UK 2 km Standard Grid'
        um_version                  '13.1'

Of those coordinates associated with air_temperature but not with any of its dimensions, time and forecast_period get promoted to span a new dimension, whereas forecast_reference_time and status_flag do not get promoted to span this new dimension.

How do time and forecast_period get promoted and forecast_reference_time and status_flag do not, if not being implicitly inferred from whether or not the coordinate value varies over the cubes?

[trexfeathers]

How do time and forecast_period get promoted and forecast_reference_time and status_flag do not, if not being implicitly inferred from whether or not the coordinate value varies over the cubes?

time and forecast_period started as scalar coordinates, which are the core of merging logic - Iris looks for scalar coordinates which are different in every Cube being merged; these can be easily combined into a 1-dimensional coordinate in the resulting Cube. And of course one of these coordinates (time in this case) is used as the basis for the new dimension.

forecast_reference_time is also scalar, but is identical across all Cubes, so remains scalar - there is only one value.

status_flag is more complicated, as it begins as multi-dimensional and also varies across the Cubes. This sits at the very edge of what Iris can handle, and you can see in this case it is not handled (I believe if it were different on every Cube we might have been fine). To state the obvious: Iris is not perfect and never will be, which is why it is important that there are ways for users to help it handle complex cases like this. We may be able to fix this in future, if the cost-benefit matches available resource, but the workarounds are there if we don't manage it.

[david-bentley]

For reference, here is the behaviour for a scalar vs non-scalar flag coordinate, based on the code at the bottom of the comment.

For a scalar status flag coordinate we get

flag value always the same
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
merged
0: air_temperature / (K)               (time: 2; height: 5; latitude: 6; longitude: 7)


flag value always different
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
merged
0: air_temperature / (K)               (time: 2; height: 5; latitude: 6; longitude: 7)


True
False
True
False
True
flag value a mixture of values
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
2: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
3: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
4: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
merged
0: air_temperature / (K)               (time: 5; height: 5; latitude: 6; longitude: 7)

and for a non-scalar status flag coordinate

flag value always the same
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
merged
0: air_temperature / (K)               (time: 2; height: 5; latitude: 6; longitude: 7)


flag value always different
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
merged
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)


False
True
False
False
True
flag value a mixture of values
0: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
2: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
3: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
4: air_temperature / (K)               (height: 5; latitude: 6; longitude: 7)
merged
0: air_temperature / (K)               (time: 3; height: 5; latitude: 6; longitude: 7)
1: air_temperature / (K)               (time: 2; height: 5; latitude: 6; longitude: 7)

The scalar case looks like it handles everything the same, or everything different, or a mixture. In the non-scalar case, the values which are the same are merged together, but even when every value is different the dimension is not promoted.

from datetime import datetime, timedelta
import random

import cf_units
import iris.coords
import iris.cube
import numpy as np

SCALAR = True

def create_cube(reference_time: datetime, forecast_period: int, flag: bool,
                n_z: int = 5, n_y: int = 6, n_x: int = 7) -> iris.cube.Cube:
    """
    Create dummy air temperature cube.

    Arguments:
        reference_time: the forecast time
        forecast_period: the number of hours from the forecast reference time
        flag: a flag indicating if something is True or False

    Keyword Arguments:
        n_z: the size of the vertical dimension
        n_y: the size of the latitude dimension
        n_x: the size of the longitude dimension
    """
    # dim coords
    height_coord = iris.coords.DimCoord(1 + np.arange(n_z), standard_name="height",
                                        units="m")
    lat_coord = iris.coords.DimCoord(np.linspace(-90, 90, n_y),
                                     standard_name="latitude", units="degrees")
    lon_coord = iris.coords.DimCoord(np.linspace(-180, 180, n_x),
                                     standard_name="longitude", units="degrees")
    dim_coords = [(height_coord, 0), (lat_coord, 1), (lon_coord, 2)]

    # aux coords
    time_units = cf_units.Unit("hours since 1970-01-01", calendar="gregorian")
    frt_coord = iris.coords.AuxCoord(time_units.date2num(reference_time),
                                     standard_name="forecast_reference_time",
                                     units=time_units)
    fp_coord = iris.coords.AuxCoord(forecast_period,
                                    standard_name="forecast_period",
                                    units="hours")
    t_coord = iris.coords.AuxCoord(time_units.date2num(reference_time + timedelta(hours=forecast_period)),
                                   standard_name="time", units=time_units)
    if SCALAR:
        flag_coord = iris.coords.AuxCoord(int(flag), standard_name="status_flag", units=1)
        aux_coords = [(frt_coord, None), (fp_coord, None), (t_coord, None), (flag_coord, None)]
    else:
        flag_coord = iris.coords.AuxCoord(int(flag) * np.ones((n_z, n_y, n_x)),
                                          standard_name="status_flag", units=1)
        aux_coords = [(frt_coord, None), (fp_coord, None), (t_coord, None), (flag_coord, (0, 1, 2))]
    

    # cube
    cube = iris.cube.Cube(np.random.random((5, 6, 7)),
                          standard_name="air_temperature", units="K",
                          dim_coords_and_dims=dim_coords,
                          aux_coords_and_dims=aux_coords)

    return cube


if __name__ == "__main__":
    dt = datetime(2025, 11, 14, 12)
    cubes = iris.cube.CubeList()
    # flag always the same
    for i in range(2):
        cubes.append(create_cube(dt, i, True))
    print("flag value always the same")
    print(cubes)
    print("merged")
    print(cubes.merge())
    print("\n")

    # flag always different
    cubes = iris.cube.CubeList()
    for i, f in zip(range(2), (True, False)):
        cubes.append(create_cube(dt, i, f))
    print("flag value always different")
    print(cubes)
    print("merged")
    print(cubes.merge())
    print("\n")

    # flag sometimes same, sometimes different
    cubes = iris.cube.CubeList()
    for i in range(5):
        flag = random.choice((True, False))
        print(flag)
        cubes.append(create_cube(dt, i, flag))
    print("flag value a mixture of values")
    print(cubes)
    print("merged")
    print(cubes.merge())