Cube operations with different coordinates

[Nicholas Tyrrell]

I haven't got a specific problem as such, but rather an ongoing issue that I'd like to get the community's advice on.

I find the use of metadata by Iris to be very useful and powerful, but this is an example something I do all the time:

diff_cube = ctrl_cube.copy()
diff_cube = ctrl_cube.data - exp1_cube.data

Where ctrl_cube is a model control run, and exp1_cube is a perturbation experiment and I want to see the difference. Each cube has the same lat, lon, pressure, and, say, 12 months of data (e.g. a model climatology) but the time coordinate is different so I need to do a hack to make it work. The cubes usually have matching auxiliary time coordinates, e.g. months and seasons.

Another common example is calculating the bias in a model run compared to ERA interim. In this case the shape of the cubes is the same, but all the coordinates differ, including attributes.

bias_cube = ctrl_cube.copy()
bias_cube = ctrl_cube.data - erai_cube.data

So my question is; is there a better way to do this, and how do other people deal with this? I assume these type of operations are very common...

An ideal solution for me would be if there was some way to choose to an Auxiliary coordinate in each cube that matches, and use that as the common coord for the operation. I realise that breaks with the philosophy of Iris regarding the strict use of metadata, so in the meantime I just have to break the rules myself! 

[Simon C. Peatman]

Nicholas,

I don't understand exactly what's happening in your code snippets.  In the first line you create a copy of a Cube, but then you overwrite this with an array.  Unless the second line was supposed to say diff_cube.data = ... ?

Another option is to copy a coordinate from one Cube to another.  For example:

exp1_cube.remove_coord(exp1_cube.coord(axis='T'))
exp1_cube.add_dim_coord(ctrl_cube.coord(axis='T').copy(), time_position)
diff_cube = ctrl_cube - exp1_cube

(I don't think this is particularly better or worse than your way, just different.)

Also, if the only problem is clashing attributes then you could just strip out the offending attributes, if you're happy to lose some metadata. iris.experimental.equalise_cubes.equalise_attributes does this for Cube attributes but not for coordinates.  Below is an extended version which handles coordinates as well.  (I haven't rigorously tested this but I've used it quite a lot and it seems to work well.)

def my_equalise_attributes(cubes):
    """
    Delete cube attributes that are not identical over all cubes in a group.

    This function simply deletes any attributes which are not the same for
    all the given cubes.  The cubes will then have identical attributes.  The
    given cubes are modified in-place.

    Arguments:

    * cubes (iterable of :class:`iris.cube.Cube`):
        A collection of cubes to compare and adjust.

    """

    # Work out which attributes are identical across all the cubes.
    common_keys = cubes[0].attributes.keys()
    dtype = cubes[0].dtype
    for cube in cubes[1:]:
        cube_keys = cube.attributes.keys()
        common_keys = [
            key for key in common_keys
            if key in cube_keys
            and np.all(cube.attributes[key] == cubes[0].attributes[key])]
        # Match data type and fill value to those of first cube.
        if cube.dtype != dtype:
            cube._my_data = cube._my_data.astype(dtype)

    # Remove all the other attributes.
    for cube in cubes:
        for key in cube.attributes.keys():
            if key not in common_keys:
                del cube.attributes[key]

    # Iterate for each co-ordinate in first cube.
    for coord0 in cubes[0].coords():
        coord_name = coord0.standard_name
        if coord_name is not None:
            dtype_points = coord0.points.dtype
            dtype_bounds = getattr(coord0.bounds, 'dtype', None)
            long_name = coord0.long_name
            var_name = coord0.var_name
            circular = getattr(coord0, 'circular', None)
            # Work out which attributes are identical (for this co-ordinate)
            # across all the cubes.
            common_coord_keys = coord0.attributes.keys()
            for cube in cubes[1:]:
                coord = cube.coord(coord_name)
                cube_coord_keys = coord.attributes.keys()
                common_coord_keys = [key for key in common_coord_keys if key in
                                     cube_coord_keys and
                                     np.all(coord.attributes[key] ==
                                            coord0.attributes[key])]
                # Match names and data types to those of coord in first cube.
                if coord.points.dtype != dtype_points:
                    coord.points = coord.points.astype(dtype_points,
                                                       copy=False)
                d = getattr(coord.bounds, 'dtype', None)
                if d != dtype_bounds and d is not None:
                    coord.bounds = coord.bounds.astype(dtype_bounds,
                                                       copy=False)
                if coord.long_name != long_name:
                    coord.long_name = long_name
                if coord.var_name != var_name:
                    coord.var_name = var_name
                if circular is not None and \
                        getattr(coord, 'circular', None) != circular:
                    coord.circular = circular

            # Remove all the other attributes (for this co-ordinate).
            for cube in cubes:
                coord = cube.coord(coord_name)
                for key in coord.attributes.keys():
                    if key not in common_coord_keys:
                        del coord.attributes[key]

[Nicholas Tyrrell]

Hi Simon

Ah, apologies, I forgot a ".data" in my snippet of code. So it should read

diff_cube = ctrl_cube.copy()
diff_cube.data = ctrl_cube.data - exp1_cube.data

So I replace the data and use the rest of the coordinates and attributes from one of the original cubes.

Thank you for posting that code, I think that looks like it could be very useful. I have tried removing attributes and replacing coordinates before, but in the end I found just copying a cube was a simpler, if less "Iris-y" way to do things. 

[RuthC]

There is an open issue about these coordinate checks here:

https://github.com/SciTools/iris/issues/2599

I believe that the intention is to make them less strict, but the core developers simply haven't had time to work on this.  In the meantime, I think all you can do is work around it by either modifying the offending metadata or just working on the cube.data as you say.

I wrote a function that will print all the differences between two coordinates, which helps me to figure out what a workaround should look like for a given case.

https://gist.github.com/rcomer/443dc06dbef2d804ebf703b44aba5bc5#file-coord_comparison-py

Ruth