Interpolating from pressure to height levels

This notebook demonstrates how to interpolate pressure level data to height levels.

import earthkit.meteo.vertical.array as vertical

Getting the data

We get some sample pressure level data for two points.

from earthkit.meteo.utils.sample import get_sample

DATA = get_sample("vertical_pl_data")
p = DATA.p  # pressure [Pa]
t = DATA.t  # temperature [K]
z = DATA.z  # geopotential [m2/s2]
sp = DATA.p_surf  # surface pressure [Pa]
zs = DATA.z_surf  # surface geopotential [m2/s2]

p, t.shape

(array([100000.,  85000.,  70000.,  50000.,  40000.,  30000.,  20000.,
         15000.,  10000.]),
 (9, 2))

Using interpolate_pressure_to_height_levels

We can use interpolate_pressure_to_height_levels() to carry out the interpolation. The example below interpolates temperature to geometric heights above the ground.

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

t_res = vertical.interpolate_pressure_to_height_levels(
    t,  # data to interpolate
    target_h,
    z,
    zs,
    h_type="geometric",
    h_reference="ground",
    interpolation="linear",
)

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

1000.0 [294.20429573 299.22387254]
5000.0 [271.02124509 272.90306903]

Heights above sea level

The reference height can also be the sea level.

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

# zs is not used in interpolate_pressure_to_height_levels()
# when h_reference="sea", so we can pass None for it
t_res = vertical.interpolate_pressure_to_height_levels(
    t,  # data to interpolate
    target_h,
    z,
    None,
    h_type="geometric",
    h_reference="sea",
    interpolation="linear",
)

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

1000.0 [298.45168008 298.99485246]
5000.0 [274.03261154 272.7133842 ]

Geopotential heights

The height type can be geopotential height.

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

t_res = vertical.interpolate_pressure_to_height_levels(
    t,  # data to interpolate
    target_h,
    z,
    zs,
    h_type="geopotential",
    h_reference="ground",
    interpolation="linear",
)

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

1000.0 [294.20228758 299.22210157]
5000.0 [270.99379632 272.87589766]

Using aux levels

In the input data the lowest pressure level is 1000 hPa. If we look at the surface pressure values we can see that the first point on 1000 hPa is lying below the surface (100000 > 94984), while the second point is lying above it (100000 < 101686).

sp

array([ 94984.34375, 101686.34375])

This means that if we interpolate to a height just above the surface the interpolation will result in missing value in the second point.

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

t_res = vertical.interpolate_pressure_to_height_levels(
    t,  # data to interpolate
    target_h,
    z,
    zs,
    h_type="geometric",
    h_reference="ground",
    interpolation="linear",
)

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

2.0 [302.09183704          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 = [2.0]  # geometric height (m), above the ground

t_res = vertical.interpolate_pressure_to_height_levels(
    t,  # data to interpolate
    target_h,
    z,
    zs,
    h_type="geometric",
    h_reference="ground",
    aux_bottom_h=0.0,
    aux_bottom_data=[304.0, 306.0],
    interpolation="linear",
)

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

2.0 [302.09183704 305.9996225 ]

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 e.g. the geometric height on the pressure levels geometric_height_from_geopotential(). This height data then can be used to interpolate multiple input data arrays.

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

# geometric height is a non-linear function of geopotential. So we need to
# do the subtraction as follows:
h = vertical.geometric_height_from_geopotential(z) - vertical.geometric_height_from_geopotential(zs)

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,  # data to interpolate
    h,
    target_h,
    interpolation="linear",
)

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

1000.0 [294.20429573 299.22387254]
5000.0 [271.02124509 272.90306903]

Heights above sea level

The reference height can also be the sea level.

# compute geometric height above sea level
h_sea = vertical.geometric_height_from_geopotential(z)

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

t_res = vertical.interpolate_monotonic(
    t,  # data to interpolate
    h_sea,
    target_h,
    interpolation="linear",
)

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

1000.0 [298.45168008 298.99485246]
5000.0 [274.03261154 272.7133842 ]

Using aux levels

For the second input data point, as described above, the lowest pressure level is lying above the surface so interpolation to a height just above the surface will result in a missing value.

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

t_res = vertical.interpolate_monotonic(
    t,  # data to interpolate
    h,
    target_h,
    interpolation="linear",
)

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

2.0 [302.09183704          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 = [2.0]  # geometric height (m), above the ground

t_res = vertical.interpolate_monotonic(
    t,  # data to interpolate
    h,
    target_h,
    aux_min_level_coord=0.0,
    aux_min_level_data=[304.0, 306.0],
    interpolation="linear",
)

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

2.0 [302.09183704 305.9996225 ]