Interpolating from hybrid to height levels

import earthkit.meteo.vertical.array as vertical

Hybrid levels

Hybrid levels are terrain following at the bottom and isobaric at the top of the atmosphere with a transition in between. This is the vertical coordinate system in the IFS model.

For details about the hybrid levels see IFS Documentation CY47R3 - Part IV Physical processes, Chapter 2, Section 2.2.1.. Also see the /how-tos/hybrid_levels.ipynb notebook for the related computations in earthkit-meteo.

Getting the data

First, we get some sample data containing 137 IFS model levels (hybrid full-levels) for two points. This data is in ascending order with regards to the model level number. So the first level is model level 1 (i.e. the top), while the last one is model level 137 (i.e. the bottom).

from earthkit.meteo.utils.sample import get_sample

DATA = get_sample("vertical_hybrid_data")
sp = DATA.p_surf  # surface pressure [Pa]
zs = DATA.z_surf  # surface geopotential [m2/s2]
t = DATA.t  # temperature [K]
q = DATA.q  # specific humidity [kg/kg]

sp.shape, zs.shape, t.shape, q.shape

((2,), (2,), (137, 2), (137, 2))

t[::40]

array([[197.06175232, 198.54808044],
       [230.29292297, 219.00386047],
       [234.56352234, 229.85160828],
       [264.58778381, 292.04481506]])

Next, get the hybrid level definitions using hybrid_level_parameters().

A, B = vertical.hybrid_level_parameters(137, model="ifs")
A.shape, B.shape

((138,), (138,))

Using interpolate_hybrid_to_height_levels

We can use interpolate_hybrid_to_height_levels() for hybrid to height level interpolation. The example below interpolates temperature to geometric heights above the ground.

target_h = [1000.0, 5000.0]  # # geometric height, above the ground

t_res = vertical.interpolate_hybrid_to_height_levels(
    t,  # data to interpolate
    target_h,
    t,
    q,
    zs,
    sp,
    A,
    B,
    h_type="geometric",
    h_reference="ground",
    interpolation="linear",
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

1000.0 [262.36938948 291.44520344]
5000.0 [236.77461002 265.49528592]

Heights above sea level

The reference height can also be the sea level.

target_h = [1000.0, 5000.0]  # # geometric height, above the sea

t_res = vertical.interpolate_hybrid_to_height_levels(
    t,  # data to interpolate
    target_h,
    t,
    q,
    zs,
    sp,
    A,
    B,
    h_type="geometric",
    h_reference="sea",
    interpolation="linear",
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

1000.0 [265.83447529 291.04194846]
5000.0 [239.80992741 264.96290891]

Geopotential heights

The height type can be geopotential height.

target_h = [1000.0, 5000.0]  # # geopotential height, above the ground

t_res = vertical.interpolate_hybrid_to_height_levels(
    t,  # data to interpolate
    target_h,
    t,
    q,
    zs,
    sp,
    A,
    B,
    h_type="geopotential",
    h_reference="ground",
    interpolation="linear",
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

1000.0 [262.36579436 291.44471712]
5000.0 [236.7517288  265.47139844]

Using a subset of levels

It is possible to use only a subset of the model levels in the input data. This can significantly speed up the computations and reduce memory usage.

In this case the model level range in the input must be contiguous and include the bottom-most level. The following example only uses the lowest 50 model levels above the surface.

# using the 50 lowest model levels as input
# Please note we still need to use the full A and B arrays
t_res = vertical.interpolate_hybrid_to_height_levels(
    t[-50:],  # data to interpolate
    target_h,
    t[-50:],
    q[-50:],
    zs,
    sp,
    A,
    B,
    h_type="geometric",
    h_reference="ground",
    interpolation="linear",
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

1000.0 [262.36938948 291.44520344]
5000.0 [236.77461002 265.49528592]

Using aux levels

Since, the lowest (full) model level is always above the surface the interpolation fails if the target height is between the surface and the lowest model level. We can compute the height of the lowest model level in our data:

vertical.height_on_hybrid_levels(t[-1:], q[-1:], zs, sp, A, B, h_type="geometric", h_reference="ground")

array([[ 9.34500108, 10.21640974]])

As a consequence, interpolation to e.g. 5 m above the ground does not work:

target_h = [5.0]  # # geometric height, above the ground

t_res = vertical.interpolate_hybrid_to_height_levels(
    t,  # data to interpolate
    target_h,
    t,
    q,
    zs,
    sp,
    A,
    B,
    h_type="geometric",
    h_reference="ground",
    interpolation="linear",
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

5.0 [nan nan]

To overcome this problem we can define “aux” levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the aux_bottom_* kwargs.

target_h = [5.0]  # # geometric height, above the ground

t_res = vertical.interpolate_hybrid_to_height_levels(
    t,  # data to interpolate
    target_h,
    t,
    q,
    zs,
    sp,
    A,
    B,
    h_type="geometric",
    h_reference="ground",
    aux_bottom_h=0.0,
    aux_bottom_data=[280.0, 300.0],
    interpolation="linear",
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

5.0 [274.04815857 296.50007348]

Using interpolate_monotonic

We can also use interpolate_monotonic() to carry out the interpolation. This is a generic method and we need an extra step to compute the height on (full) model levels using height_on_hybrid_levels(). This height data then can be used to interpolate multiple input data arrays.

First, compute the geometric height above the ground on all the model levels in the input data.

h = vertical.height_on_hybrid_levels(t, q, zs, sp, A, B, h_type="geometric", h_reference="ground")
h.shape, h[::40]

((137, 2),
 array([[81782.73207519, 80656.03851153],
        [23758.58686088, 24154.77147989],
        [ 9064.8524558 , 10238.94032239],
        [  654.70913487,   736.54960026]]))

Next, interpolate the temperature to the target height levels.

target_h = [1000.0, 5000.0]  # geometric height (m), above the ground

t_res = vertical.interpolate_monotonic(t, h, target_h, interpolation="linear")

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

1000.0 [262.36938948 291.44520344]
5000.0 [236.77461002 265.49528592]

Repeat the same computation for specific humidity.

q_res = vertical.interpolate_monotonic(q, h, target_h, interpolation="linear")

for hv, rv in zip(target_h, q_res):
    print(hv, rv)

1000.0 [0.00112569 0.00748045]
5000.0 [7.55903873e-05 1.37677882e-03]

Heights above sea level

The reference height can also be the sea level. This time we need to pass the proper surface geopotential to height_on_hybrid_levels().

h_sea = vertical.height_on_hybrid_levels(t, q, zs, sp, A, B, h_type="geometric", h_reference="sea")

target_h_sea = [1000.0, 5000.0]  # geometric height, above the sea

t_res = vertical.interpolate_monotonic(t, h_sea, target_h_sea, interpolation="linear")

for hv, rv in zip(target_h_sea, t_res):
    print(hv, rv)

1000.0 [265.83447529 291.04194846]
5000.0 [239.80992741 264.96290891]

Using a subset of levels

We can also use a subset of levels. This can significantly speed up the computations and reduce memory usage.

In this case the model level range in the input data must be contiguous and include the bottom-most level. The example below only uses the lowest 50 model levels above the surface.

# compute height on the 50 lowest model levels
# Please note we still need to use the full A and B arrays.
h_sub = vertical.height_on_hybrid_levels(t[-50:], q[-50:], zs, sp, A, B, h_type="geometric", h_reference="ground")

# interpolate to height levels
t_res = vertical.interpolate_monotonic(t[-50:], h_sub, target_h, interpolation="linear")


for hv, rv in zip(target_h, t_res):
    print(hv, rv)

1000.0 [262.36938948 291.44520344]
5000.0 [236.77461002 265.49528592]

Using aux levels

As described above, the lowest (full) model level is always above the surface so the interpolation fails if the target height is between the surface and the lowest model level.

target_h = [5.0]  # geometric height, above the ground

t_res = vertical.interpolate_monotonic(t, h, target_h, interpolation="linear")

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

5.0 [nan nan]

To overcome this problem we can define “aux” levels with a prescribed input data. As a demonstration, we set the temperature values on the surface in all the input points with the aux_min_level_* kwargs.

target_h = [5.0]  # geometric height, above the ground

t_res = vertical.interpolate_monotonic(
    t, h, target_h, aux_min_level_coord=0, aux_min_level_data=[280.0, 300.0], interpolation="linear"
)

for hv, rv in zip(target_h, t_res):
    print(hv, rv)

5.0 [274.04815857 296.50007348]