How to create an extra dimension and then merge

[Laura Burgin]

I've got several files containing monthly mean temperature on a lat/lon grid for different pressure levels that I'd like to merge into a single cube so I can take cross sections through it.

I've tried copying some code from here:
https://github.com/SciTools/iris/issues/484
Which I think does what I want and works fine:

>>> cube.add_dim_coord(DimCoord(np.arange(3), standard_name='latitude', units='degrees'), 0)
>>> cube.add_dim_coord(DimCoord(np.arange(4), standard_name='longitude', units='degrees'), 1)
>>> other = cube.copy()
>>> cube.add_aux_coord(AuxCoord(np.array([10], dtype=np.int), standard_name='height', units='m'))
>>> other.add_aux_coord(AuxCoord([20.], bounds=[[15., 25.]], standard_name='height', units='m'))
>>> merged = iris.cube.CubeList([cube, other]).merge()[0]
>>> print merged
air_temperature / (K)               (height: 2; latitude: 3; longitude: 4)
     Dimension coordinates:
          height                           x            -             -
          latitude                         -            x             -
          longitude                        -            -             x
>>> merged.coord('height')
DimCoord(array([ 10.,  20.]), standard_name='height', units=Unit('m'))

However when I try it with my data, I don't end up with a new pressure level coordinate:

>>> ERAInt_mm_1000 = iris.load_cube('/data/local/aplb/MARIUS/ERAInterim/1980-2010.t1000.monmean.nc')

>>> ERAInt_mm_1000.add_aux_coord(AuxCoord(np.array([1000], dtype=np.int), standard_name='air_pressure', units='hPa'))

>>> del ERAInt_mm_1000.attributes['history']

>>> print ERAInt_mm_1000

>>> ERAInt_mm_850 = iris.load_cube('/data/local/aplb/MARIUS/ERAInterim/1980-2010.t850.monmean.nc')

>>> ERAInt_mm_850.add_aux_coord(AuxCoord(np.array([850], dtype=np.int), standard_name='air_pressure', units='hPa'))

>>> del ERAInt_mm_850.attributes['history']

 

>>> ERAInt_mm = iris.cube.CubeList([ERAInt_mm_1000, ERAInt_mm_850]).merge()[0]

>>> print ERAInt_mm_1000

>>> print ERAInt_mm_850

>>> print ERAInt_mm

air_temperature / (K)               (time: 12; latitude: 241; longitude: 480)
     Dimension coordinates:
          time                           x             -               -
          latitude                       -             x               -
          longitude                      -             -               x
     Scalar coordinates:
          air_pressure: 1000 hPa
     Attributes:
          CDI: Climate Data Interface version 1.7.0 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.7.0 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.0
air_temperature / (K)               (time: 12; latitude: 241; longitude: 480)
     Dimension coordinates:
          time                           x             -               -
          latitude                       -             x               -
          longitude                      -             -               x
     Scalar coordinates:
          air_pressure: 850 hPa
     Attributes:
          CDI: Climate Data Interface version 1.7.0 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.7.0 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.0
air_temperature / (K)               (time: 12; latitude: 241; longitude: 480)
     Dimension coordinates:
          time                           x             -               -
          latitude                       -             x               -
          longitude                      -             -               x
     Scalar coordinates:
          air_pressure: 1000 hPa
     Attributes:
          CDI: Climate Data Interface version 1.7.0 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.7.0 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.0

Any pointers? Thanks!

[Andrew Dawson]

My guess would be that the two cubes have some metadata inconsistency that cause them not to be merged. Can you try calling merge_cube() instead? This will fail, but I think it should give you an error message describing the problem (perhaps not in enough detail to immediately resolve, but should give us a direction to work in):

iris.cube.CubeList([ERAInt_mm_1000, ERAInt_mm_850]).merge_cube()

[Laura Burgin]

Thanks for the quick response Andrew!

Yep, it reported an error do with with time:
iris.exceptions.MergeError: failed to merge into a single cube.
  Coordinates in cube.dim_coords differ: time.
 
I tried two of the other pressure level files and they merged no problem. I'll go back and check if there's an issue with the original data.

Thanks

[Andrew Dawson]

Iris is extremely sensitive to coordinate metadata (perhaps too sensitive for many applications). I'd suggest you print out the time coordinate of each cube and check for any differences in metadata, it is often attributes of the coordinate that are the problem, and you can easily sort this out within iris:

print(ERAInt_mm_1000.coord('time'))
print(ERAInt_mm_850.coord('time'))