BUG: Fixup geopotential<->height calculation (Fixes Unidata#1075)

Previous implementations was a naive implementation that directly
translated formulas from source material. It included a lot of division
and subtraction, resulting in catastrophic cancellation and severe loss
of floating point precision. New formulation is just an algebraic
reformulation that avoids these problems.

metpy/calc/basic.py

@@ -420,10 +420,9 @@ def height_to_geopotential(height):
     Derived from definition of geopotential in [Hobbs2006]_ pg.14 Eq.1.8.
 
     """
-    # Calculate geopotential
-    geopot = mpconsts.G * mpconsts.me * ((1 / mpconsts.Re) - (1 / (mpconsts.Re + height)))
-
-    return geopot
+    # Direct implementation of formula from Hobbs yields poor numerical results (see
+    # gh-1075), so was replaced with algebraic equivalent.
+    return (mpconsts.G * mpconsts.me / mpconsts.Re) * (height / (mpconsts.Re + height))
 
 
 @exporter.export
@@ -462,10 +461,10 @@ def geopotential_to_height(geopot):
     Derived from definition of geopotential in [Hobbs2006]_ pg.14 Eq.1.8.
 
     """
-    # Calculate geopotential
-    height = (((1 / mpconsts.Re) - (geopot / (mpconsts.G * mpconsts.me))) ** -1) - mpconsts.Re
-
-    return height
+    # Direct implementation of formula from Hobbs yields poor numerical results (see
+    # gh-1075), so was replaced with algebraic equivalent.
+    scaled = geopot * mpconsts.Re
+    return scaled * mpconsts.Re / (mpconsts.G * mpconsts.me - scaled)
 
 
 @exporter.export

metpy/calc/tests/test_basic.py

@@ -281,13 +281,34 @@ def test_height_to_geopotential():
                               29443], units('m**2 / second**2')), 0)
 
 
+# See #1075 regarding previous destructive cancellation in floating point
+def test_height_to_geopotential_32bit():
+    """Test conversion to geopotential with 32-bit values."""
+    heights = np.linspace(20597, 20598, 11, dtype=np.float32) * units.m
+    truth = np.array([201590.422, 201591.391, 201592.375, 201593.344,
+                      201594.312, 201595.297, 201596.266, 201597.250,
+                      201598.219, 201599.203, 201600.172], dtype=np.float32) * units('J/kg')
+    assert_almost_equal(height_to_geopotential(heights), truth, 2)
+
+
 def test_geopotential_to_height():
     """Test conversion from geopotential to height."""
     geopotential = units.Quantity([0, 9817.70342881, 19632.32592389,
                                   29443.86893527], units('m**2 / second**2'))
     height = geopotential_to_height(geopotential)
     assert_array_almost_equal(height, units.Quantity([0, 1000, 2000, 3000], units.m), 0)
 
+
+# See #1075 regarding previous destructive cancellation in floating point
+def test_geopotential_to_height_32bit():
+    """Test conversion from geopotential to height with 32-bit values."""
+    geopot = np.arange(201590, 201600, dtype=np.float32) * units('J/kg')
+    truth = np.array([20596.957, 20597.059, 20597.162, 20597.266,
+                      20597.365, 20597.469, 20597.570, 20597.674,
+                      20597.777, 20597.881], dtype=np.float32) * units.m
+    assert_almost_equal(geopotential_to_height(geopot), truth, 2)
+
+
 # class TestIrrad(object):
 #    def test_basic(self):
 #        'Test the basic solar irradiance calculation.'