...
| Mathdisplay |
|---|
\begin{eqnarray*}
A(\lambda,\mu,\eta,t) & = & \sum_{mn=0}^{\mbox{T}} \sum_{|n|\leq mm = -n}^{\mbox{T}n} X(n,m,n) (\eta,t) P_{mn}^{nm}(\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, Pnm n are the associated Legendre polynomials of the first kind, and eimλ are the Fourier functions. The order of the summations can be interchanged to provide an equivalent expression:
| Mathdisplay |
|---|
\begin{eqnarray*}
A(\lambda,\mu,\eta,t) & = & \sum_{m=-T}^{\mbox{T}} \sum_{n = |m|}^{T} X(n,m) (\eta,t) P_{n}^{m}(\mu) e^{im\lambda}
\end{eqnarray*} |
In the GRIB binary data section, the complex X(n,m,n) coefficients are stored for m ≥ 0 as pairs of real numbers Re(X(n,m,n)) and Im(X(n,m,n)) ordered with n increasing from m to T, first for m = 0 and then for m = 1, 2, . . . T.
...
| Note | ||
|---|---|---|
| ||
Only the coefficients for n m ≥ 0 are stored. The coefficients for n m < 0 are the complex conjugates of the corresponding coefficient with n m > 0 : X(mn,-nm)= X* (n,m,n) |
To extract the real and imaginary pairs of the coefficients from the values array returned by ecCodes needs the following two steps:
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
# 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 mn=%d\t nm=%d %.10f\t%.10f" % (i, mn, nm, 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()) |
...