Reshaping a cube

[Peter Watson]

Hello,

I have a cube of seasonal forecast data with the dimensions described below - 'time' is an auxiliary coordinate with 'Forecast_start_time' and 'forecast_month' (i.e. the lead time of the forecast in months) being the temporal dimension coordinates. However, other scripts that I have need 'time' to be a dimension coordinate. Is there any way I can reshape the cube to have 'time' as a dim coord and 'Forecast_start_time' and 'forecast_month' as aux coords? All of the values of 'time' are unique (i.e. none of the months in the forecasts are repeated), so I don't see a reason in principle why it could not be used as a dimension.

low_type_cloud_area_fraction / (1)  (forecast_month: 6; Forecast_start_time: 2; latitude: 72; longitude: 144)

     Dimension coordinates:

          forecast_month                                     x                       -                                 -                      -

          Forecast_start_time                              -                       x                                 -                      -

          latitude                                                  -                       -                                 x                      -

          longitude                                               -                       -                                 -                      x

     Auxiliary coordinates:

          forecast_period                                    x                       x                                -                       -

          time                                                      x                       x                                -                       -

Thanks very much for any help you can give.

[Peter Watson]

I take the lack of responses to mean there is not a simple way to do this. So I wrote a function to do it myself, basically by making a new cube from a data array and coordinates made by reshaping those from the original cube. Then I thought it wouldn't be too much effort to generalise it, and I've shared what I have below, in case anyone else ever wants to do something similar. I've tested it with various transpositions of the cube I printed above and with the coordinate being promoted having between 1 and 3 dimensions, and it seems to work fine. It would be useful to know if anyone finds this helpful, in case it could be worth adding it into Iris.

def promote_coord(cube,coord_to_promote):
    """
    Promotes an AuxCoord on the cube to a DimCoord. The AuxCoord may be
    associated with multiple cube dimensions. Other AuxCoords must be
    associated either with the same set of cube dimensions or with a
    single dimension. The DimCoords associated with those of the
    AuxCoord get demoted to an AuxCoord. (Note this function needs to be
    modified to work with cubes that have dimensions without a DimCoord.
    Also, the promoted coordinate is made the leading coordinate.)
   
    Args:


    * cube
        An instance of :class:`iris.cube.Cube`


    * coord_to_promote:
        The :attr:`standard_name`, :attr:`long_name`, or
        :attr:`var_name` of an instance of an instance of
        :class:`iris.coords.AuxCoord`.
   
    For example::


        >>> print cube

        low_type_cloud_area_fraction / (1)  (forecast_month: 6; Forecast_start_time: 2; latitude: 72; longitude: 144)
        Dimension coordinates:
          forecast_month                           x                       -            -              -
          Forecast_start_time                      -                       x            -              -
          latitude                                 -                       -            x              -
          longitude                                -                       -            -              x
        Auxiliary coordinates:
          forecast_period                          x                       x            -              -
          time                                     x                       x            -              -

        Scalar coordinates:
          Ens member: 0      
        >>> cube2=promote_coord(cube, 'time')
        >>> print cube2
        low_type_cloud_area_fraction / (1)  (time: 12; latitude: 72; longitude: 144)
        Dimension coordinates:
          time                           x             -              -
          latitude                       -             x              -
          longitude                      -             -              x
        Auxiliary coordinates:
          Forecast_start_time            x             -              -
          forecast_month                 x             -              -
          forecast_period                x             -              -
        Scalar coordinates:
          Ens member: 0
    """
   
    dim_coord_names=[coord.name() for coord in cube.dim_coords]
    coords_to_demote=cube.coord_dims(coord_to_promote)
   
    #Make the coords to demote the leading cube dimensions - this makes it simple to reshape the cube data as needed below.
    other_dims=tuple([i for i in range(len(dim_coord_names)) if i not in coords_to_demote])
    transpose_indices=coords_to_demote+other_dims
    cube.transpose(transpose_indices)
   
    coords_to_demote_new_inds=cube.coord_dims(coord_to_promote)
    other_dims_new_inds=tuple([i for i in range(len(dim_coord_names)) if i not in coords_to_demote_new_inds])


    #The shape the final cube will have    
    new_shape=(np.product(cube.shape[:len(coords_to_demote)]),)+cube.shape[len(coords_to_demote):]


    #Making the new leading coordinate
    new_points=cube.coord(coord_to_promote).points.transpose(cube.coord_dims(coord_to_promote)).reshape(new_shape[0])  #transpose to make sure the points dimensions correspond to the cube dimensions before reshaping.
   
    #Get indices for sorting the coordinate values into ascending or descending order, if they are not already ordered. The new DimCoord coordinate values must be monotonic and non-repeating.
    assert not np.any([new_points[i] in [new_points[j] for j in range(len(new_points)) if j!=i] for i in range(len(new_points))]), coord_to_promote+' contains repeating coordinate values and cannot be made a DimCoord.'
   
    if not(np.all(new_points[1:]>new_points[:-1]) or np.all(new_points[1:]<new_points[:-1])):  #checking if coordinate values are not already all ascending/descending
        sort_inds=np.argsort(new_points)
        new_points=new_points[sort_inds]
    else:
        sort_inds=range(len(new_points))
       
    if cube.coord(coord_to_promote).bounds is not None:
        bounds_transpose_dims=cube.coord_dims(coord_to_promote)+(len(cube.coord_dims(coord_to_promote)),)
        new_bounds=cube.coord(coord_to_promote).bounds.transpose(bounds_transpose_dims).reshape(new_shape[0],2)
        new_bounds=new_bounds[sort_inds,:]
    else:
        new_bounds=None
       
    promoted_coord=iris.coords.DimCoord(
        new_points,
        standard_name=cube.coord(coord_to_promote).standard_name,
        long_name=cube.coord(coord_to_promote).long_name,
        var_name=cube.coord(coord_to_promote).var_name,
        units=cube.coord(coord_to_promote).units,
        bounds=new_bounds,
        attributes=cube.coord(coord_to_promote).attributes,
        coord_system=cube.coord(coord_to_promote).coord_system,
        #Cannot copy the circular attribute because, as an AuxCoord, coord_to_promote does not have one
        )
   
    #New indices of DimCoords that will remain as DimCoords (going from 1 to number of other DimCoords)
    new_other_dims=range(1,len(other_dims_new_inds)+1)
   
    #Making full list of tuples of new DimCoords and the index of their dimension on the final cube
    new_dim_coords=[(promoted_coord,0)]+[(cube.dim_coords[i],j) for i,j in zip(other_dims_new_inds,new_other_dims)]
   
    #Making new AuxCoords
    new_aux_coords=[]
   
    #First, create AuxCoords from DimCoords that are being demoted with element-to-element correspondence with the coord being promoted
    for i in range(len(coords_to_demote)):
   
        coord=cube.dim_coords[i]
       
        #First tile the coordinate values for each value of the previous dimensions to be demoted
        new_points=np.tile(coord.points, np.product(cube.shape[:i]))
       
        #Then repeat the values for all subsequent dimensions to be demoted
        new_points=np.repeat(new_points, np.product(cube.shape[i+1:len(coords_to_demote)]))
       
        #Sort in the same way the promoted dimension was sorted, so the values match up.
        new_points=new_points[sort_inds]
       
        #Get new bounds in a similar way
        if coord.bounds is not None:
            new_bounds=np.tile(coord.bounds, (np.product(cube.shape[:i]),1))
            new_bounds=np.repeat(new_bounds, np.product(cube.shape[i+1:len(coords_to_demote)]), axis=0)
            new_bounds=new_bounds[sort_inds,:]
        else:
            new_bounds=None
           
        new_aux_coord=iris.coords.AuxCoord(
            new_points,
            standard_name=coord.standard_name,
            long_name=coord.long_name,
            var_name=coord.var_name,
            units=coord.units,
            bounds=new_bounds,
            attributes=coord.attributes,
            coord_system=coord.coord_system,
            )  
       
        new_aux_coords+=[(new_aux_coord,0)]  #this coord is associated with the first dimension of the new cube
   
    #Now create new AuxCoords from existing AuxCoords - if they correspond to any of the demoted dimensions, then they need to be reshaped. Else, copy them from the existing cube.
    for coord in cube.aux_coords:
        if coord.name()!=coord_to_promote:
            assert cube.coord_dims(coord.name()) in [coords_to_demote_new_inds]+[(i,) for i in coords_to_demote_new_inds]+[()], 'Cannot proceed as '+coord.name()+' has more than one dimension and is associated with different dimension coordinates to those associated with '+coord_to_promote
            coord=cube.coord(coord.name())
           
            if cube.coord_dims(coord.name())==coords_to_demote_new_inds:
                coord_dims=cube.coord_dims(coord)
                new_points=coord.points.transpose(coord_dims).reshape(new_shape[0])
                new_points=new_points[sort_inds]
   
                if coord.bounds is not None:
                    bounds_transpose_dims=coord_dims+(len(coord_dims),)
                    new_bounds=coord.bounds.transpose(bounds_transpose_dims).reshape(new_shape[0],2)
                    new_bounds=new_bounds[sort_inds,:]
                else:
                    new_bounds=None
   
                new_dim=0
           
            elif cube.coord_dims(coord.name()) in coords_to_demote_new_inds:
                coords_to_demote_ind=list(coords_to_demote_new_inds).index(cube.coord_dims(coord.name()))
               
                #Get new points array in similar way as for DimCoords that are being demoted
                new_points=np.tile(coord.points, np.product(cube.shape[:coords_to_demote_ind]))
                new_points=np.repeat(new_points, np.product(cube.shape[coords_to_demote_ind+1:len(coords_to_demote_new_inds)]))
                new_points=new_points[sort_inds]
               
                if coord.bounds is not None:
                    new_bounds=np.tile(coord.bounds, (np.product(cube.shape[:coords_to_demote_ind]),1))
                    new_bounds=np.repeat(new_bounds, np.product(cube.shape[coords_to_demote_ind+1:len(coords_to_demote_new_inds)]), axis=0)
                    new_bounds=new_bounds[sort_inds,:]
                else:
                    new_bounds=None
           
                new_dim=0
           
            else:
                new_points=coord.points
                new_bounds=coord.bounds
                new_dim=()
                   
            new_aux_coord=iris.coords.AuxCoord(
                new_points,
                standard_name=coord.standard_name,
                long_name=coord.long_name,
                var_name=coord.var_name,
                units=coord.units,
                bounds=new_bounds,
                attributes=coord.attributes,
                coord_system=coord.coord_system,
                )  
           
            new_aux_coords+=[(new_aux_coord, new_dim)]
   
    #New cube data
    new_data=cube.data.reshape(new_shape)
    new_data=new_data[sort_inds,...]
   
    cube=iris.cube.Cube(  
        new_data,
        standard_name=cube.standard_name,
        long_name=cube.long_name,
        var_name=cube.var_name,
        units=cube.units,
        attributes=cube.attributes,
        cell_methods=cube.cell_methods,
        dim_coords_and_dims=new_dim_coords,
        aux_coords_and_dims=new_aux_coords,
        aux_factories=cube.aux_factories  #not sure if this line will always work properly, as I'm not familiar with cube.aux_factories
        )


    return cube

[marqh]

Hi Peter

This is an interesting problem space, thank you for sharing it with us.  I have looked into a generic cube.reshape operation and found too many issues with it to be confident in implementation.

This is a more constrained case, which could well be the limit which allows some operations in some circumstances.

I'm interested in pursuing this implementation and seeing if it could have a life.  I wonder what the test cases would be.  Which cases should work and which cases should fail.

If we could come up with a few of those, we might have a chance of proposing an enhancement

all the best
mark

[Peter Watson]

Thanks for your input Mark. I suppose I stuck to thinking about it in terms of promoting coordinates rather than reshaping cubes because that was sufficient for solving my problem. Adding functionality to promote multi-dimensional AuxCoords allows the aspect of reshaping to occur where cube dimensions are combined, so that the number of cube dimensions decreases - you can add an AuxCoord associated with the DimCoords to be combined with values that give the desired ordering of points, and then promote it. Increasing the number of cube dimensions is harder I think.

I don't see why promoting coords couldn't be made to work quite generally - in my mind, it's basically a matter of reshaping the cube's data array and reshaping the points and bounds arrays of the cube's coordinates to maintain the correspondence with the coord values of the promoted coord, and copying all of the other attributes. I think my function above could be extended to work for cubes with AuxCoords that have multiple dimensions that are not the same as those of the coordinate being promoted, but I just didn't need that myself and couldn't find time to add that in. I'm very far from being an expert on Iris, so probably there are lots of potential problems that I am unaware of - however, my function seems to work fine on cubes I have produced from my data and so I would have thought it would work for many use cases, as I don't think I'm doing anything special. Are there special cases that you think would fail?

I'm not quite sure what test cases would be useful (I'm a scientist rather than a developer really, so I'm not sure what is the best way to go about these things). For promoting coordinates, obviously the function needs to work for cubes and coordinates spanning the possible combinations of: cubes with 1 to N dimensions (for a not too big choice of N for testing); cubes with AuxCoords with 1 to M dimensions with M<=N; and cubes with AuxCoords associated with different combinations of the DimCoords on each cube. It should fail, as my function does, if the coordinate being promoted does not have values that can be rearranged into a monotonic sequence - this may restrict the allowable datatype of the coord values.

Cheers,

Peter