Rounding error in geopotential_to_height with large dataset

[sgdecker]

My code:

import pickle

import numpy as np

import metpy
import metpy.calc as mpcalc
from metpy.units import units

print(metpy.__version__)


# Simple case works as expected
phi = [2218.999755859375, 202024.84375, 201592.84375] * units('J/kg')
z = mpcalc.geopotential_to_height(phi)

a = phi[-1].m
b = z[-1].m

print(f'Phi: {a:.3f}  z: {b:.3f}')


# Much larger dataset: Rounding error
with open('phi.pickle', 'rb') as f:
    phi = pickle.load(f)

phi = np.asarray(phi) * units('J/kg')
z = mpcalc.geopotential_to_height(phi)

a = phi[-1,-1,-1].m
b = z[-1,-1,-1].m

print(f'Phi: {a:.3f}  z: {b:.3f}')

The output:

0.10.0+45.g19deb8b7
Phi: 201592.844  z: 20597.248
Phi: 201592.844  z: 20597.000

• Expected output

I expect the value of z to be unchanged in the two scenarios, rather than rounded off to the nearest integer.

The pickle file is available here:
https://drive.google.com/file/d/1o7vxZfUH2HnUp7ghAB3HdPBMV90ZADse/view?usp=sharing

[jthielen]

At least when I tested this out (on the current master branch), this appears to be a dtype issue. In the first case, the numpy array had float64, whereas in the second, it had float32. Adding the dtype argument to np.asarray gives the expected output:

Updated Code Sample

import pickle

import numpy as np

import metpy
import metpy.calc as mpcalc
from metpy.units import units

print(metpy.__version__)


# Simple case
phi = [2218.999755859375, 202024.84375, 201592.84375] * units('J/kg')
z = mpcalc.geopotential_to_height(phi)

a = phi[-1].m
b = z[-1].m

print(f'Phi: {a:.3f}  z: {b:.3f}')


# Much larger dataset
with open('phi.pickle', 'rb') as f:
    phi = pickle.load(f)

phi = np.asarray(phi, dtype=np.float64) * units('J/kg')
z = mpcalc.geopotential_to_height(phi)

a = phi[-1,-1,-1].m
b = z[-1,-1,-1].m

print(f'Phi: {a:.3f}  z: {b:.3f}')

Output

0.8.0+271.gd3c29401
Phi: 201592.844  z: 20597.248
Phi: 201592.844  z: 20597.248

[sgdecker]

OK, so the "large dataset" is a red herring. It still seems there is excessive loss of precision, though, as I'd expect float32 to be good for six digits of accuracy (or at least for the sixth digit to be off by no more than one).

Revised program:

import numpy as np

import metpy
import metpy.calc as mpcalc
from metpy.units import units

print(metpy.__version__)


for p in range(201590, 201600):
    phi32 = np.array([p], dtype=np.float32) * units('J/kg')
    z32 = mpcalc.geopotential_to_height(phi32)
    phi64 = np.array([p], dtype=np.float64) * units('J/kg')
    z64 = mpcalc.geopotential_to_height(phi64)
    
    a = phi32[0].m
    b = z32[0].m
    c = z64[0].m
    d = abs(b - c)
    e = np.asarray(c, dtype=np.float32)
    
    print(f'Phi: {a:.3f}  z32: {b:.3f}  z64: {c:.6f}  d: {d:.3f}  e: {e:.6f}')

Output:

0.10.0+45.g19deb8b7
Phi: 201590.000  z32: 20597.000  z64: 20596.956894  d: 0.043  e: 20596.957031
Phi: 201591.000  z32: 20597.000  z64: 20597.059397  d: 0.059  e: 20597.058594
Phi: 201592.000  z32: 20597.000  z64: 20597.161900  d: 0.162  e: 20597.162109
Phi: 201593.000  z32: 20597.000  z64: 20597.264403  d: 0.264  e: 20597.263672
Phi: 201594.000  z32: 20597.500  z64: 20597.366906  d: 0.133  e: 20597.367188
Phi: 201595.000  z32: 20597.500  z64: 20597.469408  d: 0.031  e: 20597.468750
Phi: 201596.000  z32: 20597.500  z64: 20597.571911  d: 0.072  e: 20597.572266
Phi: 201597.000  z32: 20597.500  z64: 20597.674414  d: 0.174  e: 20597.673828
Phi: 201598.000  z32: 20597.500  z64: 20597.776917  d: 0.277  e: 20597.777344
Phi: 201599.000  z32: 20597.500  z64: 20597.879420  d: 0.379  e: 20597.878906

Expected output:
All values of d less than 0.1, or, in other words, for values of z32 to be much closer to the values of e than they are.

[jthielen]

Given the operations involved in calculating height from geopotential, I'm not sure if this loss of precision with float32 is avoidable: (EDIT: as pointed out by @dopplershift, algebraic reformulation helps greatly)

@dopplershift Any thoughts?

[dopplershift]

The naive formulation (i.e. direct from Hobbs) has some pretty terrible floating point aspects to it, including two almost catastrophic subtraction cancellations plus reciprocal of fractions. Just bad. I've reworked and gotten much better results. PR incoming...

[sgdecker]

Thanks for taking care of this!