How to apply custom calculation between two IRIS cubes (GRIB files)? Considering also using xarray

Question

I am trying to do some calculation between two iris cubes (GRIB files), here it is what I'm trying to achieve:

First cube: ERA5-Land dataset, downloaded from official site via cdsapi API routine, cropped to custom Lat and Lon, in this example, I have only 2m air temperature, in celsius, hourly, for 3 days:

print(air_temperature)

air_temperature / (celsius)         (time: 72; latitude: 18; longitude: 27)
    Dimension coordinates:
        time                             x             -              -
        latitude                         -             x              -
        longitude                        -             -              x
    Auxiliary coordinates:
        forecast_period                  x             -              -
    Scalar coordinates:
        height                      2 m
        originating_centre          European Centre for Medium Range Weather Forecasts

Then, I have a series of sampling points at gives coordinates:

## Sample points coordinates
ws_latitudes = np.array([40.64, 41.19, 41.11, 41.19, 40.86, 40.93, 40.83, 40.25, 40.79, 40.56, 41.42, 41.42, 41.02, 41.24, 40.64, 40.13, 41.33, 40.61])
ws_longitudes = np.array([14.54, 15.13, 14.82, 13.83, 15.28, 14.02, 15.03, 15.66, 14.16, 15.23, 13.88, 15.04, 14.34, 14.47, 14.83, 15.45, 14.33, 14.97])

ws_samplepoints = [("latitude", ws_latitudes), ("longitude", ws_longitudes)]

The other cube (GRIB file) is a 2D cube ("timeless") of elevation:

I've downloaded ERA-Land geopontential GRIB2 file from here: https://confluence.ecmwf.int/display/CKB/ERA5-Land%3A+data+documentation#ERA5Land:datadocumentation-parameterlistingParameterlistings

geopotential = "geo_1279l4_0.1x0.1.grib2"
geopot_cube = iris.load_cube(geopotential)
print(geopot_cube)
geopotential / (m2 s-2)             (latitude: 1801; longitude: 3600)
    Dimension coordinates:
        latitude                             x                -
        longitude                            -                x
    Scalar coordinates:
        forecast_period             0 hours
        forecast_reference_time     2013-08-09 12:00:00
        time                        2013-08-09 12:00:00
    Attributes:
        GRIB_PARAM                  GRIB2:d000c003n004
        centre                      'European Centre for Medium Range Weather Forecasts'
z, Geopotential, m**2 s**-2

Then, to convert the geopotential to elevation, I've divided by 9.80665 m/s^2

elev_cube = geopot_cube / 9.80665
elev_cube.rename("Elevation")
elev_cube.units = "m"
print(elev_cube)
Elevation / (m)                     (latitude: 1801; longitude: 3600)
    Dimension coordinates:
        latitude                             x                -
        longitude                            -                x
    Scalar coordinates:
        forecast_period             0 hours
        forecast_reference_time     2013-08-09 12:00:00
        time                        2013-08-09 12:00:00
    Attributes:
        GRIB_PARAM                  GRIB2:d000c003n004
        centre                      'European Centre for Medium Range Weather Forecasts'

The resulting cube has been cropped to the same lat and lon as air temperature above (probably not necessary):

area_slicer = iris.Constraint(longitude=lambda v: 13.45 <= v <= 16.14, latitude=lambda v: 39.84 <= v <= 41.6)
elevcube_slice = elev_cube.extract(area_slicer)
print(elevcube_slice)
Elevation / (m)                     (latitude: 18; longitude: 27)
    Dimension coordinates:
        latitude                             x              -
        longitude                            -              x
    Scalar coordinates:
        forecast_period             0 hours
        forecast_reference_time     2013-08-09 12:00:00
        time                        2013-08-09 12:00:00
    Attributes:
        GRIB_PARAM                  GRIB2:d000c003n004
        centre                      'European Centre for Medium Range Weather Forecasts'

Now here is the point: having these two cubes, I have to calculate a new temperature value at every sample points given the linear equation:

where:

= temperature to calculate at given coordinates sample points;

= temperature read from the first GRIB file (2m air temperature) at sample points coordinates

= sample point elevation

= elevation from second GRIB file at sample points coordinates

as temperature/meter

How could I achieve this? Even when I try to do very simple math between the two cubes, for example a simple multiplication: print(air_temperature * elevcube_slice)

I have this error: ValueError: Coordinate 'latitude' has different points for the LHS cube 'air_temperature' and RHS cube 'Elevation'.

To double check, both cubes have same CS:

cselev = elevcube_slice.coord_system()
cstemperature = air_temperature.coord_system()
print(cselev, cstemperature)
GeogCS(6371229.0) GeogCS(6371229.0)

I've also considered to switch to xarray if it is possible and suggested, probably working with xarray dataset is easier?

Answer 1

Iris is strict about metadata and will fail loudly when they don't match in operations you try to do.

The error you get tells you what's going on: ValueError: Coordinate 'latitude' has different points for the LHS cube 'air_temperature' and RHS cube 'Elevation'.

So you can investigate and compare your left and right hand sides with cube.coord('latitude').points

Xarray on the other hand is not strict about metadata and assumes you know what you are doing, i.e. will also do operations that will give wrong results.

Both have their merits. I'm siding with xarray for reading and analysing files. And iris when writing files.