You are viewing an old version of this page. View the current version.

Compare with Current View Page History

Version 1 Next »

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(m,n) coefficients in the discrete representation of the data on a grid by a continuous function expressed as a truncated series of spherical harmonics:

\[ \begin{eqnarray*} A(\lambda,\mu,\eta,t) & = & \sum_{m=0}^{\mbox{T}} \sum_{|n|\leq m}^{\mbox{T}} X(m,n) (\eta,t) P_{m}^{n}(\mu) e^{im\lambda} \end{eqnarray*} \]

where μ = sinθ with λ  the longitude and θ the latitude of the grid point, T is the triangular spectral truncation number, Pm n are the associated Legendre polynomials of the first kind,  and  eimλ are the Fourier functions. 

In the GRIB binary data section, the complex X(m,n) coefficients are stored for m ≥ 0 as pairs of real numbers Re(X(m,n)) and Im(X(m,n)) 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(0,1)), Im(X(0,1), Re(X(0,2)), Im(X(0,2)), … , Re(X(0,T)), Im(X(0,T)), Re(X(1,1)), Im(X(1,1)), Re(X(1,2)), Im(X(1.2)), … , Re(X(1,T)), Im(X(1,T)), ... ,Re(X(T-1,T)), Im(X(T-1,T)), Re(X(T,T)), Im(X(T,T)).

Only the coefficients for n ≥ 0 are stored.  The coefficients for n < 0 are the complex conjugates of the corresponding coefficient with n > 0 :  X(m,-n)= X* (m,n)

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

  1. Get the J key (also known as pentagonalResolutionParameterJ) which for ECMWF data where triangular truncation is used, provides the spectral truncation number, T.
  2. 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.
 
from __future__ import print_function
import traceback
import sys

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 m=%d\t n=%d %.10f\t%.10f" % (i, m, n, 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())


See also:



  • No labels