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

@@ -410,22 +410,21 @@ def height_to_geopotential(height):
     >>> from metpy.constants import g, G, me, Re
     >>> import metpy.calc
     >>> from metpy.units import units
-    >>> height = np.linspace(0,10000, num = 11) * units.m
+    >>> height = np.linspace(0, 10000, num=11) * units.m
     >>> geopot = metpy.calc.height_to_geopotential(height)
     >>> geopot
-    <Quantity([     0.           9817.46806283  19631.85526579  29443.16305888
-    39251.39289118  49056.54621087  58858.62446525  68657.62910064
-    78453.56156253  88246.42329545  98036.21574306], 'meter ** 2 / second ** 2')>
+    <Quantity([     0.           9817.46806283  19631.85526579  29443.16305887
+    39251.39289118  49056.54621087  58858.62446524  68657.62910064
+    78453.56156252  88246.42329544  98036.21574305], 'meter ** 2 / second ** 2')>
 
     Notes
     -----
     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
@@ -448,12 +447,12 @@ def geopotential_to_height(geopot):
     >>> from metpy.constants import g, G, me, Re
     >>> import metpy.calc
     >>> from metpy.units import units
-    >>> height = np.linspace(0,10000, num = 11) * units.m
+    >>> height = np.linspace(0, 10000, num=11) * units.m
     >>> geopot = metpy.calc.height_to_geopotential(height)
     >>> geopot
-    <Quantity([     0.           9817.46806283  19631.85526579  29443.16305888
-    39251.39289118  49056.54621087  58858.62446525  68657.62910064
-    78453.56156253  88246.42329545  98036.21574306], 'meter ** 2 / second ** 2')>
+    <Quantity([     0.           9817.46806283  19631.85526579  29443.16305887
+    39251.39289118  49056.54621087  58858.62446524  68657.62910064
+    78453.56156252  88246.42329544  98036.21574305], 'meter ** 2 / second ** 2')>
     >>> height = metpy.calc.geopotential_to_height(geopot)
     >>> height
     <Quantity([     0.   1000.   2000.   3000.   4000.   5000.   6000.   7000.   8000.
@@ -464,10 +463,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

@@ -283,13 +283,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.'