Fieldlist: Computing values on hybrid levels¶
This notebook demonstrates how to compute various quantities on hybrid levels using GRIB fieldlist data. Fieldlist support requires the usage of earthkit-data.
[1]:
import earthkit.data as ekd
import numpy as np
import earthkit.meteo.vertical.fieldlist as vertical
Getting the data¶
The input is GRIB data on 137 IFS model levels for a single step. We fetch the file and inspect its contents.
[2]:
fl = ekd.from_source("sample", "tq_ml137_europe.grib2").to_fieldlist()
fl.head()
[2]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | lnsp | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 1 | hybrid | None | regular_ll |
| 1 | z | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 1 | hybrid | None | regular_ll |
| 2 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 1 | hybrid | None | regular_ll |
| 3 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 2 | hybrid | None | regular_ll |
| 4 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 3 | hybrid | None | regular_ll |
The file contains:
t – temperature (K) on 137 hybrid full-levels
q – specific humidity (kg/kg) on 137 hybrid full-levels
lnsp – logarithm of surface pressure on model level 1 (used to recover surface pressure)
z – surface geopotential (m²/s²) on model level 1
Please note that although lnsp and z are valid for the surface in our data they are available on hybrid level 1. This is how IFS forecast and analysis data is stored in ECMWF’s MARS archive (from where our the data is originating).
We extract each quantity below. The log surface pressure is converted to surface pressure in Pa.
[3]:
# LNSP → surface pressure [Pa]
lnsp = fl.sel({"parameter.variable": "lnsp", "vertical.level": 1})[0]
sp = lnsp.set(values=np.exp(lnsp.values))
t = fl.sel({"parameter.variable": "t"}) # temperature, 137 levels
q = fl.sel({"parameter.variable": "q"}) # specific humidity, 137 levels
zs = fl.sel({"parameter.variable": "z", "vertical.level": 1})[0] # surface geopotential
t.tail()
[3]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 133 | hybrid | None | regular_ll |
| 1 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 134 | hybrid | None | regular_ll |
| 2 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 3 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 4 | t | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
The A and B coefficients that define the hybrid levels are embedded in the GRIB metadata and are inferred automatically by all fieldlist vertical functions — there is no need to provide them explicitly.
If you do need to inspect these coefficient you can get them from any GRIB field on a hybrid level as follows:
[4]:
# get pv array from raw GRIB metadata, this is not yet available as high-level field metadata
pv = lnsp.get("metadata.pv")
# get A and B coefficients from pv
coeff_num = int(len(pv) / 2)
A, B = pv[:coeff_num], pv[coeff_num:]
len(A), len(B), A[:5], B[:5]
[4]:
(138,
138,
array([0. , 2.00036502, 3.10224104, 4.66608381, 6.82797718]),
array([0., 0., 0., 0., 0.]))
Computing pressure on hybrid levels¶
:py:meth:~earthkit.meteo.vertical.fieldlist.pressure_on_hybrid_levels computes the pressure on all full hybrid levels from the surface pressure.
[5]:
p_full = vertical.pressure_on_hybrid_levels(sp, output="full")
p_full.tail()
[5]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 133 | hybrid | None | regular_ll |
| 1 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 134 | hybrid | None | regular_ll |
| 2 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 3 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 4 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
The output is a FieldList sorted by model level number in ascending order (level 1 at the top of the atmosphere, level 137 near the surface).
[6]:
# the first value from the lowest 3 levels(=135,136,137) above the surface
p_full[-3:].values[:, 0]
[6]:
array([101538.53712552, 101815.60759216, 102069.07594384])
Using the output keyword argument, additional pressure-related quantities can be computed in a single call. See pressure_on_hybrid_levels() for details.
[7]:
p_full, p_half, alpha, delta = vertical.pressure_on_hybrid_levels(
sp, alpha_top="ifs", output=["full", "half", "alpha", "delta"]
)
len(p_full), len(p_half), len(alpha), len(delta)
[7]:
(137, 138, 137, 137)
It is possible to request only a contiguous subset of levels via the levels kwarg. The example below computes the pressure on the 3 lowest model levels.
[8]:
p_full_sub = vertical.pressure_on_hybrid_levels(sp, levels=[135, 136, 137], output="full")
p_full_sub.ls()
[8]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 1 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 2 | pres | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
Computing geopotential thickness¶
relative_geopotential_thickness_on_hybrid_levels() computes the geopotential thickness (m²/s²) between the surface and each full hybrid level.
[9]:
z_thickness = vertical.relative_geopotential_thickness_on_hybrid_levels(t, q, sp)
z_thickness.tail()
[9]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 133 | hybrid | None | regular_ll |
| 1 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 134 | hybrid | None | regular_ll |
| 2 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 3 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 4 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
[10]:
# the first value from the lowest 3 levels(=135,136,137) above the surface
z_thickness[-3:].values[:, 0]
[10]:
array([503.24894095, 288.97561477, 93.33643824])
A subset of levels can be used by pre-filtering the input FieldLists. The level range must be contiguous and must include the bottom-most level.
[11]:
# use only the 3 lowest model levels
t_sub = t.sel({"vertical.level": [135, 136, 137]})
q_sub = q.sel({"vertical.level": [135, 136, 137]})
z_thickness_sub = vertical.relative_geopotential_thickness_on_hybrid_levels(t_sub, q_sub, sp)
z_thickness_sub.ls()
[11]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 1 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 2 | rel_geopot_thick | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
Computing geopotential on hybrid levels¶
geopotential_on_hybrid_levels() computes the geopotential (m²/s²) on all full hybrid levels. It requires the surface geopotential zs.
[12]:
z = vertical.geopotential_on_hybrid_levels(t, q, zs, sp)
z.tail()
[12]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | z | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 133 | hybrid | None | regular_ll |
| 1 | z | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 134 | hybrid | None | regular_ll |
| 2 | z | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 3 | z | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 4 | z | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
[13]:
# the first value from the lowest 3 levels(=135,136,137) above the surface
z[-3:].values[:, 0]
[13]:
array([433.73099661, 219.45767043, 23.81849391])
Computing height on hybrid levels¶
:py:meth:`~earthkit.meteo.vertical.fieldlist.height_on_hybrid_levels` computes the height (m) on all full hybrid levels. The height type and reference level are controlled by ``h_type`` and ``h_reference``: - **h_type**: ``"geometric"`` (default) or ``"geopotential"`` - **h_reference**: ``"ground"`` (default) or ``"sea"``Geometric height above the ground¶
[14]:
h = vertical.height_on_hybrid_levels(t, q, zs, sp, h_type="geometric", h_reference="ground")
h.tail()
[14]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 133 | hybrid | None | regular_ll |
| 1 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 134 | hybrid | None | regular_ll |
| 2 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 3 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 4 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
[15]:
# the first value from the lowest 3 levels(=135,136,137) above the surface
h[-3:].values[:, 0]
[15]:
array([51.31740957, 29.46738267, 9.51766097])
Geometric height above sea level¶
[16]:
h_sea = vertical.height_on_hybrid_levels(t, q, zs, sp, h_type="geometric", h_reference="sea")
h_sea.tail()
[16]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 133 | hybrid | None | regular_ll |
| 1 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 134 | hybrid | None | regular_ll |
| 2 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 3 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 4 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
[17]:
# the first value from the lowest 3 levels(=135,136,137) above the surface
h_sea[-3:].values[:, 0]
[17]:
array([44.22855996, 22.37853306, 2.42881137])
Using a subset of hybrid levels¶
A subset of levels can be used by pre-filtering the input FieldLists. The example below computes the geometric height above the ground for the 3 lowest model levels.
[18]:
h_sub = vertical.height_on_hybrid_levels(t_sub, q_sub, zs, sp, h_type="geometric", h_reference="ground")
h_sub.ls()
[18]:
| parameter.variable | time.valid_datetime | time.base_datetime | time.step | vertical.level | vertical.level_type | ensemble.member | geography.grid_type | |
|---|---|---|---|---|---|---|---|---|
| 0 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 135 | hybrid | None | regular_ll |
| 1 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 136 | hybrid | None | regular_ll |
| 2 | h | 2019-06-02 12:00:00 | 2019-06-02 12:00:00 | 0 days | 137 | hybrid | None | regular_ll |
[ ]: