View Source

Some ECMWF data produced by the IFS is stored in GRIB with gridType=sh to indicate that the values are stored as spherical harmonic coefficients. These are the X(n,m) coefficients in the discrete representation of the data on a grid at time, t, and vertical coordinate, η, by a continuous function expressed as a truncated series of spherical harmonics:

\begin{eqnarray*}
A(\lambda,\mu,\eta,t) & = &  \sum_{m=-\mbox{T}}^{\mbox{T}} \, \sum_{n=|m|}^{\mbox{T}}  X(n,m) (\eta,t) \overline{P}_{n}^{m}(\mu) e^{im\lambda}
\end{eqnarray*}

where μ = sinθ with λ the longitude and θ the latitude of the grid point, T the triangular spectral truncation number and e^imλ are the Fourier functions. The normalised associated Legendre polynomials of the first kind of degree n and order m are denoted by

\overline{P}_n^m

with the normalisation defined by:

\begin{eqnarray*}
\frac{1}{2} \int_{-1}^{1}\, \{\overline{P}_{n}^{m}(\mu)\}^2 \,d\mu = 1\, .
\end{eqnarray*}

In the GRIB binary data section, the complex X(n,m) coefficients are stored for m ≥ 0 as pairs of real numbers Re(X(n,m)) and Im(X(n,m)) ordered with n increasing from m to T, first for m = 0 and then for m = 1, 2, . . . T.

Hence, the ecCodes values array contains the coefficients in the order Re(X(0,0)), Im(X(0,0)), Re(X(1,0)), Im(X(1,0), Re(X(2,0)), Im(X(2,0)), … , Re(X(T,0)), Im(X(T,0)), Re(X(1,1)), Im(X(1,1)), Re(X(2,1)), Im(X(2,1)), … , Re(X(T,1)), Im(X(T,1)), ... ,Re(X(T,T-1)), Im(X(T,T-1)), Re(X(T,T)), Im(X(T,T)).

Only the coefficients for m ≥ 0 are stored. The coefficients for m < 0 are obtained from the complex conjugates of the corresponding coefficient with m > 0 as X(n,-m)= X^*(n,m) / (-1)^m

The spherical harmonic (n, m) wave number space is shown in the figure.

User Documentation > How to access the data values of a spherical harmonic field in GRIB - ecCodes GRIB FAQ > nm_coeff_pm.png

The figure depicts the spherical harmonic (n, m) wave number space for a triangular truncation, T.

Only the coefficients for m and n represented by the filled red dots are stored in the GRIB Binary Data Section.

When the complex coefficients are extracted as the ecCodes values array, they are returned as real and imaginary pairs, starting from m=0, n=0 (i.e., at the apex of the inverted triangle) and moving upwards in direction of increasing n until n=T. It then moves to the next value of m and continues in this way, reading from the n = m line, upwards in the direction of increasing n until the final pair with n = m = T is reached.

To extract the real and imaginary pairs of the coefficients from the values array returned by ecCodes needs the following two steps:

Get the J key (also known as pentagonalResolutionParameterJ) which for ECMWF data where triangular truncation is used, provides the spectral truncation number, T.
Read the complex coefficients as pairs of real numbers from the values array starting at m=0

The following Python script illustrates the algorithm:

# Copyright 2020 ECMWF.
#
# This software is licensed under the terms of the Apache Licence Version 2.0
# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
#
# In applying this licence, ECMWF does not waive the privileges and immunities
# granted to it by virtue of its status as an intergovernmental organisation
# nor does it submit to any jurisdiction.

import sys
import traceback

from eccodes import *
 
INPUT = "../../data/spherical_model_level.grib2"
VERBOSE = 1  # verbose error reporting
 
def example():
    f = open(INPUT, "rb")

    while 1:
        gid = codes_grib_new_from_file(f)
        if gid is None:
            break

        # Get the values. This will be a numpy array
        values = codes_get_array(gid,"values")
 
        # Store the real and imaginary parts of the coefficients in arrays a and b, respectively
        # The real parts 'a' stored in every second element starting at element 0
        a = values[0::2]
        # The imaginary parts 'b' stored in every second element starting at element 1
        b = values[1::2]
  
        codes_release(gid)

        # Loop through the values and print the m and n indices together with the corresponding
        # real and imaginary parts of the coefficients  
        m = 0
        n = 0
        for i in range(len(a)):
            print("%d\t n=%d\t m=%d %.10f\t%.10f" % (i, n, m, a[i], b[i]))
            n += 1
            if n > T:
                m += 1
                n = m 
 
    f.close()

def main():
    try:
        example()
    except CodesInternalError as err:
        if VERBOSE:
            traceback.print_exc(file=sys.stderr)
        else:
            sys.stderr.write(err.msg + "\n")
 
        return 1


if __name__ == "__main__":
    sys.exit(main())

Related articles