Tests broken with updates to Pint constants and unit defaults

[jthielen]

So, I pointed my local Pint repo to the master branch instead of a tagged version or feature branch as I had been doing, and discovered that numerous tests are failing (see this gist)! Most are due calculations falling outside tolerance (which looks to be due to hgrecco/pint#811), but there is this particularly strange one

pint.errors.DimensionalityError: Cannot convert from 'inch_Hg ** 0.190267' ([mass] ** 0.190267 / [length] ** 0.190267 / [time] ** 0.380533) to 'hectopascal ** 0.190267' ([mass] ** 0.190267 / [length] ** 0.190267 / [time] ** 0.380533)

(link)

this one from a change of the default name for degrees Celsius

AssertionError: assert 'degree_Celsius' == 'degC'

(link),

and one that was is skipped for 0.9 but isn't fixed yet in Pint (see hgrecco/pint#751)

TypeError: Neither Quantity object nor its magnitude (0)supports indexing

(link).

I'm thinking most of these should be pretty straightforward to fix...but seeing all this after #997 and #1111 makes me think: is it possible to set up an allowed-failure CI build that tests against the master branch of pint to ensure these things get caught?

[dopplershift]

Pint is simple enough to build that we could add installing pint from git to our existing set of dev builds.

[jthielen]

In trying to resolve the broken tests that are due to the changed pint constants, most can be fixed by adjusting the tolerances to less stringent, but still reasonable values.

However, the geopotential <-> height conversion functions seem to be very sensitive to this change:

_________________________ test_height_to_geopotential __________________________

    def test_height_to_geopotential():
        """Test conversion from height to geopotential."""
        height = units.Quantity([0, 1000, 2000, 3000], units.m)
        geopot = height_to_geopotential(height)
        assert_array_almost_equal(geopot, units.Quantity([0., 9817, 19632,
>                                 29443], units('m**2 / second**2')), 0)

metpy/calc/tests/test_basic.py:293: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

actual = array([    0.        ,  9818.14474002, 19633.20840782, 29445.19245349])
desired = array([    0.,  9817., 19632., 29443.]), decimal = 0

    def assert_array_almost_equal(actual, desired, decimal=7):
        """Check that arrays are almost equal, including units.

        Wrapper around :func:`numpy.testing.assert_array_almost_equal`
        """
        actual, desired = check_and_drop_units(actual, desired)
>       numpy.testing.assert_array_almost_equal(actual, desired, decimal)
E       AssertionError: 
E       Arrays are not almost equal to 0 decimals
E       
E       Mismatch: 25%
E       Max absolute difference: 2.19245349
E       Max relative difference: nan
E        x: array([    0.,  9818., 19633., 29445.])
E        y: array([    0.,  9817., 19632., 29443.])

metpy/testing.py:160: AssertionError
______________________ test_height_to_geopotential_32bit _______________________

    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)

metpy/calc/tests/test_basic.py:303: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

actual = array([201604.31, 201605.28, 201606.27, 201607.23, 201608.22, 201609.19,
       201610.17, 201611.14, 201612.12, 201613.1 , 201614.06],
      dtype=float32)
desired = array([201590.42, 201591.39, 201592.38, 201593.34, 201594.31, 201595.3 ,
       201596.27, 201597.25, 201598.22, 201599.2 , 201600.17],
      dtype=float32)
decimal = 2

    def assert_almost_equal(actual, desired, decimal=7):
        """Check that values are almost equal, including units.

        Wrapper around :func:`numpy.testing.assert_almost_equal`
        """
        actual, desired = check_and_drop_units(actual, desired)
>       numpy.testing.assert_almost_equal(actual, desired, decimal)
E       AssertionError: 
E       Arrays are not almost equal to 2 decimals
E       
E       Mismatch: 100%
E       Max absolute difference: 13.90625
E       Max relative difference: 6.898136e-05
E        x: array([201604.31, 201605.28, 201606.27, 201607.23, 201608.22, 201609.19,
E              201610.17, 201611.14, 201612.12, 201613.1 , 201614.06],
E             dtype=float32)
E        y: array([201590.42, 201591.39, 201592.38, 201593.34, 201594.31, 201595.3 ,
E              201596.27, 201597.25, 201598.22, 201599.2 , 201600.17],
E             dtype=float32)

metpy/testing.py:151: AssertionError
______________________ test_geopotential_to_height_32bit _______________________

    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)

metpy/calc/tests/test_basic.py:321: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

actual = array([20595.531, 20595.635, 20595.736, 20595.84 , 20595.941, 20596.045,
       20596.146, 20596.25 , 20596.354, 20596.455], dtype=float32)
desired = array([20596.957, 20597.059, 20597.162, 20597.266, 20597.365, 20597.469,
       20597.57 , 20597.674, 20597.777, 20597.88 ], dtype=float32)
decimal = 2

    def assert_almost_equal(actual, desired, decimal=7):
        """Check that values are almost equal, including units.

        Wrapper around :func:`numpy.testing.assert_almost_equal`
        """
        actual, desired = check_and_drop_units(actual, desired)
>       numpy.testing.assert_almost_equal(actual, desired, decimal)
E       AssertionError: 
E       Arrays are not almost equal to 2 decimals
E       
E       Mismatch: 100%
E       Max absolute difference: 1.4257812
E       Max relative difference: 6.922291e-05
E        x: array([20595.53, 20595.63, 20595.74, 20595.84, 20595.94, 20596.04,
E              20596.15, 20596.25, 20596.35, 20596.46], dtype=float32)
E        y: array([20596.96, 20597.06, 20597.16, 20597.27, 20597.37, 20597.47,
E              20597.57, 20597.67, 20597.78, 20597.88], dtype=float32)

metpy/testing.py:151: AssertionError

These mismatches seem about as bad as, if not worse than, the 32-bit destructive cancellation errors seen in #1075, for which specific strict-tolerance tests were implemented in #1082. What should be done about these?

[dopplershift]

@jthielen Sounds like some digging into those calculations is needed. Are you just doing that by installing pint master, at least to reproduce those particular failures?

[jthielen]

Installing pint master is the easiest way to reproduce this set of failures, and that's how I produced the gist in the initial comment. (When I've been going through this, however, I've been using my __array_function__ implementation branch at andrewgsavage/pint#6 to try fixing these issues at the same time as #1111.)

A quick look back at the calculation shows that the only thing that appears to have changed is G. In pint v0.9, its value was:

6.67384e-11 m^3 kg^-1 s^-2

whereas in current pint master, it is

6.67430e-11 m^3/(kg s^2)

If I alter the definition of G in metpy/constants.py to

G = units.Quantity(6.67384e-11, 'm^3 / kg / s^2')

the geopotential <-> height conversion functions go back to passing.

[dopplershift]

I am deeply disturbed by those tests failing, producing IMO large differences, due to a change in G in the 4th decimal place.

[jthielen]

This issue has reappeared with the continued updates to pint and MetPy's test suite. Right now (other than the infinite recursion with powers, which an upstream issue: hgrecco/pint#941, and the changes needed below), the following tests are failing:

• Not sure what is wrong:
  • test_isentropic_pressure_p_increase_rh_out (returning all nans)
• Precision too high, should be ~3 decimals
  • test_multiple_lfs_wide
  • test_multiple_el_wide
  • test_muliple_el_most_cape
  • test_muliple_lfc_most_cape
  • test_el_lfc_most_cape_bottom
  • test_lcl_grid_surface_LCLs

There were two places where it was assumed that units were stripped, and now they are not. The following diff shows the basics of the changes needed:

diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index b7028aaa..fe0b706c 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -359,7 +359,7 @@ def lcl(pressure, temperature, dewpt, max_iters=50, eps=1e-5):
 
     # np.isclose needed if surface is LCL due to precision error with np.log in dewpoint.
     # Causes issues with parcel_profile_with_lcl if removed. Issue #1187
-    lcl_p = np.where(np.isclose(lcl_p, pressure), pressure, lcl_p) * pressure.units
+    lcl_p = np.where(np.isclose(lcl_p, pressure.m), pressure.m, lcl_p) * pressure.units
 
     return lcl_p, dewpoint(vapor_pressure(lcl_p, w)).to(temperature.units)
 
diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py
index 7910c658..6507b0fa 100644
--- a/src/metpy/calc/tools.py
+++ b/src/metpy/calc/tools.py
@@ -777,6 +777,10 @@ def lat_lon_grid_deltas(longitude, latitude, **kwargs):
     if latitude.ndim < 2:
         longitude, latitude = np.meshgrid(longitude, latitude)
 
+    # pyproj requires ndarrays, not Quantities
+    longitude = np.asarray(longitude)
+    latitude = np.asarray(latitude)
+
     geod_args = {'ellps': 'sphere'}
     if kwargs:
         geod_args = kwarg

Although, on second thought, the lat_lon_grid_deltas could probably use more robust unit-handling along the way.

When would it be best to get a PR for these in?

[dopplershift]

@jthielen I think it'd be good to get one in for 0.12. Also, if we need any more discussion, probably best to do it in a follow-up issue since this one was already marked as closed in 0.11.