Convert profiles vertical coordinate from Hybrid-Pressure to pressure

%matplotlib inline
# for figures in notebook

# import & initialize epygram
import epygram
epygram.init_env()

import os
INPUTS_DIR = os.path.join('..', 'inputs')

r = epygram.resources.meta_resource(os.path.join(INPUTS_DIR, 'ICMSHAROM+0022'), 'r', 'CL')  # CL = CombineLevels
# in a "CL" meta-resource, fields are to be grabbed through a GRIB2 nomenclature
t_as_grib2 = dict(discipline=0, parameterCategory=0, parameterNumber=0,  # this is temperature
                  typeOfFirstFixedSurface=119)  # hybrid-pressure levels

# read as 3D
t3d = r.readfield(t_as_grib2)
t3d.sp2gp()

# extract profiles
profiles = []
for position in [(6.1,46.2), (6.8,45.8), (6,46.5)]:
    profgeo = t3d.geometry.make_profile_geometry(*position)  # make a V1D geometry at that position
    profiles.append(t3d.extract_subdomain(profgeo))  # extract subdomain corresponding to the geometry of the profile

# plot with level as vertical coord
fig = ax = None
for p in profiles:
    fig, ax = p.plotfield(labels='Temperature @ ({}, {})'.format(p.geometry.grid['longitudes'][0],
                                                                 p.geometry.grid['latitudes'][0]),
                          over=(fig, ax))  # superpose profiles

and now on a Pressure z-axis

# for the conversion, we need Surface pressure
ps = epygram.formats.resource(r.resource.container.abspath, 'r').readfield('SURFPRESSION')
# Ps is spectral and ln(Ps)
ps.sp2gp()
ps.operation('exp')

from epygram.geometries.VGeometry import hybridP2pressure

for p in profiles:
    ps_local = ps.getvalue_ll(p.geometry.grid['longitudes'][0],
                              p.geometry.grid['latitudes'][0])
    # create the P vertical coordinate
    p.geometry.vcoordinate = hybridP2pressure(p.geometry.vcoordinate,
                                              ps_local,
                                              vertical_mean='geometric')

fig = ax = None
for p in profiles:
    fig, ax = p.plotfield(labels='Temperature @ ({}, {})'.format(p.geometry.grid['longitudes'][0],
                                                                 p.geometry.grid['latitudes'][0]),
                          over=(fig, ax))  # superpose profiles