Correct way to change cube and coordinates data points without changing attributes

[Denis Sergeev]

Hi all,

I would really appreciate if anyone could explain me the correct method to create a cube using coordinates and metadata of an existing cube.

I have a Met Office Unified Model output file in .pp format, which I open in Iris as a 4D cube with dimensions of "time", "pressure", "grid_latitude" and "grid_longitude".

Now, my idea is to 'unstagger' the cube data, i.e. to interpolate the values into cell centres simply by averaging neighbouring points:

new_data = 0.5*(original_data[:-1, ...] + original_data[1:, ...])

and obviously average the corresponding coordinate points (e.g. along longitude) in the same manner.

In my script I do this by extracting the cube's data and coordinate points, and after averaging I create a new cube by calling iris.cube.Cube with the same attributes as the old cube, except for the data array and 
dim_coords_and_dims (I also do not pass the "aux_coords_and_dims", as I do not need them), which I also reproduce by using iris.coords.DimCoord instance with reduced number of points along the averaging axis (to match the shape of the data array). In the iris.coords.DimCoord function I also pass all the attributes of the coordinate of the original cube. The only thing that I change is the name of horizontal coordinates: "grid_latitude" to "latitude" and "grid_longitude" to "longitude". 

Having done all this, I write the new cube to a NetCDF file.

However, my problem is that when I call new_cube.coords() function an error is returned:

CoordinateNotFoundError: 'Expected to find exactly 1  coordinate, but found none.'

while calling the same function in the original cube works fine. 

Using iris.analysis.cartography.rotate_winds() function with the recreated cubes returns an error too:

/path/to/iris/aux_factory.pyc in updated(self, new_coord_mapping)

    222         for key, coord in self.dependencies.iteritems():

    223             if coord:

--> 224                 coord = new_coord_mapping[id(coord)]

    225             new_dependencies[key] = coord

    226         return type(self)(**new_dependencies)

KeyError: 57483088

And finally, calling new_cube.coord_system() also does not work, it just returns None, while it should return RotatedGeogCS(10.5, 154.8, ellipsoid=GeogCS(6371229.0)).

Can anyone tell me what the problem is? Perhaps it is a wrong method of creating a cube all along?

Thanks,

Denis

[RuthC]

Hi Denis,

I don't know why your cube creation isn't working, but would the following function do what you want?
iris.analysis.interpolate.linear()

Ruth

[Denis Sergeev]

Thanks, Ruth!

Apparently, I was trying to reinvent the wheel, but you've opened my eyes...

Cheers,

Denis