Loading Grib data - Iris won't use ensemble member number as dimension

[Pagw]

Hello,

I'm trying to load some ECMWF ensemble forecast data of a particular variable in Grib format, with different messages corresponding to forecasts starting from different dates, different lead times of forecasts and different ensemble members. If I do iris.load_cube(file) with file being the file name, I get the following error:

ConstraintMismatchError: failed to merge into a single cube.
  Duplicate 'UNKNOWN LOCAL PARAM 139.128' cube, with scalar coordinates time=Cell(point=-86232.0, bound=None), forecast_period=Cell(point=696, bound=None), originating_centre=Cell(point='European Centre for Medium Range Weather Forecasts', bound=None)

If I create files with subsets of the Grib messages from the original file, this error only occurs for files with more than one ensemble member. The ensemble members are labelled with Grib key perturbationNumber, and Iris does not seem to be detecting this as information that differentiates the Grib messages. I've attached a Grib file with two messages containing data for the same forecast start time and verification time, but different ensemble members, and output from the grib_api grib_dump function on this file, so you can see what the data look like.

Is there a way to fix this and have perturbationNumber as a cube dimension?

[Pagw]

I came across the following page that explains how to create a dimension for the ensemble member number based on the filename, where the data for different ensemble members is separated into different files. Is there a way to do a similar thing based on the Grib key perturbationNumber within the file messages?

http://scitools.org.uk/iris/docs/latest/examples/graphics/lagged_ensemble.html

[Andrew Dawson]

Hi Peter

There is a mention of the perturbationParameter key in the GRIB load rules, which should result in a DimCoord called 'ensemble_member'. It looks like this isn't happening in your case. I notice this is only activated for GRIB2 files, are you working with GRIB1 or GRIB2? This may need to be fixed in Iris, as it happens the GRIB loading is being completely overhauled over the next few weeks, so now is a good time to raise this issue.

In any case, you might be able to work-around this using a load-time callback as in the example you linked to. Notice a callback takes 3 arguments: cube, field, and filename. You would need to access the perturbationNumber attribute of the field argument. The field argument should be the internal representation of a GRIB message with access to the keys via GIRB API.

I haven't tested any of this, so post back if you have problems with this.

[Pagw]

Thanks for the reply Andrew. Yes I'm working with GRIB1 files - quite a lot of ECMWF datasets use GRIB1 so having compatibility with this would be helpful.

I don't know myself how to structure a callback to do what is required. In the callback in the link I gave, it looks to me like the cube needs to have been created before the callback can work, since it needs to evaluate cube.coords('realization'), but in my case it seems like Iris won't make the cube. Do you or anyone else have any idea about how to code up such a callback? If not, I guess I will have to use another method to read in my data.

[Andrew Dawson]

The callback is applied to each sub-cube before a merge is attempted, so it should work.

The iris load process will construct a cube for each individual GRIB message, and then attempt to combine the cubes as best as it can. Usually this does not result in failure, but in your case it fails to distinguish ensemble members, so they appear to be duplicates. The failure to create your cube is actually the failure to make sense of multiple sub-cubes that have already been created. By adding an ensemble dimension you will allow iris to distinguish these cubes and the merge should work. It might not be exactly what you want the first time round, but if you take it step by step you should end up with something useful.

[Pagw]

Thanks again Andrew. I don't see how to get Iris to read the perturbationNumber variable so that it can be used to create a dimension for the ensemble members, though - I don't see it amongst any of the cube coordinates or metadata even when I just load a grib file with a single message. I've tried  the following on a Grib file containing a single message, based on what is at the webpage I linked to previously:

def ens_callback(cube, field, filename):
        import iris.coords
        ens_member_coord = iris.coords.AuxCoord(np.int32(perturbationNumber), 'ens_member')
        cube.add_aux_coord(ens_member_coord)

cube=iris.load_cube(fname,iris.Constraint(ens_member=lambda value:True),callback=ens_callback)

However, this of course fails because Iris does not know what perturbationNumber is. I have no idea myself how to get around this.

[Andrew Dawson]

As I said before, you need to use the perturbationNumber attribute of the field argument:

def ens_callback(cube, field, filename):

    ens_member_coord = iris.coords.AuxCoord(field.perturbationNumber, 'ens_member')

    ...

[Pagw]

Thanks a lot, it seems to be working now - though, for the information of anyone else reading this, I had to modify the argument to iris.coords.AuxCoord slightly so that that part reads iris.coords.AuxCoord(np.int32(field.perturbationNumber), long_name='ens_member') .

[Pagw]

Also the iris.Constraint part of the argument to iris.load_cube() does not seem to be necessary.