How to apply custom calculation between two cubes (GRIB files)?

[mcrimaldi]

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:

$T_{at sample points}= T_{airtempGRIB}+(z_{sp}-z_{gb})*0.005$

where:
$T_{at sample points}$ = temperature to calculate at given coordinates sample points
$T_{airtempGRIB}$ = temperature read from the first GRIB file (2m air temperature) at sample points coordinates
$z_{sp}$ = sample point elevation
$z_{gb}$ = elevation from second GRIB file at sample points coordinates
$0.005$ 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)

Thank you so much for helping

[stephenworsley]

What are the points of the latitude for each of these cubes? i.e. what do you get when you call:

print(air_temperature.coord("latitude").points)
print(elevcube_slice.coord("latitude").points)

I suspect that if your cubes are not defined on the same grids, it may be appropriate to interpolate onto sample points before doing your calculations so that your cubes are defined on the same coordinates. See this interpolation function https://scitools-iris.readthedocs.io/en/latest/generated/api/iris/analysis/trajectory.html#iris.analysis.trajectory.interpolate .

I'm also curious where your two different elevation values $z_{sp}$ and $z_{gb}$ derive from since elev_cube seems to be the only source for elevation.

[mcrimaldi]

I double checked all the lat and lon coordinates, rounded all to 1st decimal and changed the crop to be both of same dimension:

print(air_temperature.coord("latitude").points)
print(elevcube_slice.coord("latitude").points)

print(air_temperature.coord("longitude").points)
print(elevcube_slice.coord("longitude").points)

[41.6 41.5 41.4 41.3 41.2 41.1 41.  40.9 40.8 40.7 40.6 40.5 40.4 40.3
 40.2 40.1 40.  39.9 39.8]
[41.6 41.5 41.4 41.3 41.2 41.1 41.  40.9 40.8 40.7 40.6 40.5 40.4 40.3
 40.2 40.1 40.  39.9 39.8]
[13.5 13.6 13.7 13.8 13.9 14.  14.1 14.2 14.3 14.4 14.5 14.6 14.7 14.8
 14.9 15.  15.1 15.2 15.3 15.4 15.5 15.6 15.7 15.8 15.9 16.  16.1]
[13.5 13.6 13.7 13.8 13.9 14.  14.1 14.2 14.3 14.4 14.5 14.6 14.7 14.8
 14.9 15.  15.1 15.2 15.3 15.4 15.5 15.6 15.7 15.8 15.9 16.  16.1]

I can now run the trajectory interpolation you suggested but I don't understand the result and how I can perform the equation $T_{at sample points}= T_{airtempGRIB}+(z_{sp}-z_{gb})*0.005$ as in my first post:

zgb = trajectory.interpolate(elev_cube, ws_samplepoints)

print(zgb)
print(zgb.coord("longitude").points)
print(zgb.coord("latitude").points)

Elevation / (m)                     (-- : 18)
    Auxiliary 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'
[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]
[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]

[pp-mo]

@mcrimaldi I double checked all the lat and lon coordinates, rounded all to 1st decimal ... I can now run the trajectory interpolation you suggested but I don't understand the result

If you have two cubes with all-very-close coordinate points, then I think there is very little harm in simply changing the coordinates of one to match the other : Interpolation or regridding seems unnecessary and unhelpful in such a case.
I think I would try to do it a reasonably generic + safe way, something like this (NB untested!!) ...

for axis in ('x', 'y'):
    dimcoord_1 = cube_1.coord(axis=axis, dim_coords=True)
    dimcoord_2 = cube_2.coord(axis=axis, dim_coords=True)
    # Ensure they are basically the same, before proceeding.
    assert np.allclose(dimcoord_1.points, dimcoord_2.points)
    assert dimcoord_1.has_bounds() == dimcoord_2.has_bounds()
    if dimcoord_1.has_bounds()
        assert np.allclose(dimcoord_1.bounds, dimcoord_2.bounds)
    # change cube#2 to use coords from cube#1 
    cube_2.replace_coord(dimcoord_1)  

Or, if you're sure of your facts, then just cube_2.replace_coord(cube_1.coord(axis=axis)) will do it.

[mcrimaldi]

@mcrimaldi I double checked all the lat and lon coordinates, rounded all to 1st decimal ... I can now run the trajectory interpolation you suggested but I don't understand the result

If you have two cubes with all-very-close coordinate points, then I think there is very little harm in simply changing the coordinates of one to match the other : Interpolation or regridding seems unnecessary and unhelpful in such a case. I think I would try to do it a reasonably generic + safe way, something like this (NB untested!!) ...

@pp-mo Both cubes (air_temperature and elevcube_slice) are cropped to same sampling area (from a vector file) used to download GRIB files via ERA5 API:

def takearea():
    """ Take max min of coordinates as tuple"""
    global area

    supported_drivers['KML'] = 'rw' # Enable fiona to read KML files
    area_shp = gpd.read_file('./shps/area.kml', driver='KML') # read KML file fo sampling area
    area_tup = area_shp.total_bounds # take total bounds as tuple
    area = list(map(lambda item: round((item), 1), area_tup)) # convert coordinates tuple to list rounding to 1 decimal
    order = [3, 0, 1, 2] # rearrange the coordinates order to be ERA5 compliant
    area = [area[i] for i in order]

    return area

So, as the coordinates prints above show, I am pretty sure coordinates are the same.
By the way, my problem here is to apply the equation $T_{at sample points}= T_{airtempGRIB}+(z_{sp}-z_{gb})*0.005$ to calculate new temperature from sample points elevation and elevation from 'elevcube_slice'.

I was thinking about create a third cube made by sampling points coordinates and elevations, then doing the maths between these three cubes, but as usual I don't know how and if it is possible.

Maybe a cube with same lat and lon dimension of the other twos with all masked data except for the sampling points?

Thank you all for helping me!

[pp-mo]

Ah, I think I didn't read the whole spec, but was only talking about the the problem with nearly-equal lat+lon coordinates.
So I missed the need to extract at specific (non-grid) station locations anyway... Interpolation makes perfect sense for that, so you may as well apply the same to both the inputs, and not worry about the slight difference.

So I think you now have a suitable Zgb cube, with 18 points corresponding to the station locations.

In order to do arithmetic, you need this and the other two cubes to all share a "station dimension" (length 18).
That means also interpolating the other gridded data inputs (air-temperature and surface-elevation) from lats and lons onto the fixed station points -- to get a "air_temperature_at_stations(time: 72; station: 18)" and a "sample_elevation_at_stations(station: 18)".

Thus, for the main air-temp data, you will want something like
T_airtemp_grib = trajectory.interpolate(air_temperature, ws_samplepoints)

I'm not totally sure what 'Zsp' needs to be, but I think the 'air_temperature' data has a constant 'height' == 2m, so these are all surface (or near-surface) values. In which case, you need to get Zsp by interpolating from a standard orography (i.e. surface elevation) field, which with luck will be available in the same file.
So possibly, all that is needed is orography = iris.load(grib_file, 'surface_altitude') + 2.0 ??
You would then do the same interpolation, i.e. Zsp = trajectory.interpolate(orography, ws_samplepoints).

Finally, in order to do the cube arithmetic, I think you will need all 3 cubes to have a matching "station" dimension so that Iris can identify and line up the 'station' dimension common to all of them. So it is probably also necessary to 'label' the station dimension by adding a matching dimension coordinate to all 3 of them.
From your above table, the "ID" column is monotonic, so that can be the value of the 'station' dim-coord.
Assuming that is always simply 1..18, it can be create as
station_coord = iris.coords.DimCoord(np.arange(1,19), long_name='station'),
and added to all 3 of the cubes like this :
cube.add_dim_coord(station_coord, axis=cube.ndim-1)
(assuming the stations is the last dim in each cube)

So I think we then have one 2-dimensional cube T_airtemp_grib(time: 72; station:18), and two 1-D cubes Zsp(station: 18) and Zgb(station: 18). Arithmetic between these should then succeed !
It can also fail due to additional metadata which differs between the cubes, typically either auxiliary-coordinates or attributes : Offending information can usually be just removed, until they agree.

[pp-mo]

(re-reading, slight mistake :
"cube.add_dim_coord(station_coord, axis=cube.ndim-1)"
should be
"cube.add_dim_coord(station_coord, data_dim=cube.ndim-1)"
I got the kwarg name wrong.
)

[mcrimaldi]

@pp-mo thanks, I think I'm almost there...

I found another little "problem". I thought there was the elevation data (Elevation column in the table) in ws_samplepoints but no, there are only the lat and lon coordinates:

# Latitudes from file
lat_file = "lat_file.csv"
ws_latitudes = np.genfromtxt(lat_file, delimiter=",")

# Longitudes from file
lon_file = "lon_file.csv"
ws_longitudes = np.genfromtxt(lon_file, delimiter=",")
# #ws_longitudes = np.loadtxt(lon_file, delimiter=",")

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

[('latitude',
  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])),
 ('longitude',
  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]))]

the $z_{sp}$ is exactly the elevation of each weather station (Elevation column in the table).
So, I could create a cube with Latitudes, Longitudes and Elevations of each record of the table?
Because if i create a numpy.array with all elevations, I will have a 1D vector without the coordinates.
I'm feeling we almost solved my issue, but probably I am missing something or there is something else I don't understand.

the T_airtemp_grib = trajectory.interpolate(air_temperature, ws_samplepoints) is what I need because I suppose is the cube with the temperatures of downloaded ERA-5 GRIB file (air_temperature) at each weather stations coordinate (ws_samplepoints);

$z_{gb}$ is the elevation at weather station coordinates taken from elevcube_slice that is another GRIB file, so I think Zgb = trajectory.interpolate(elevcube_slice, ws_samplepoints) is what I need.

The final step, as i wrote before, is the creation of $z_{sp}$ with all elevations at weather station coordinates. Then, I hope, I can do maths as you explained.