Updates to pint and xarray usage for future version compatability

There were *many* places in the code base where passing pint Quantities
to numpy functions was assumed to return a plain ndarray. This behavior
will no longer be the case after pint implements __array_function__ in
its next version, and so, this modifies usage to explicitly request
magnitude where it is assumed that quantities are stripped.

Also, makes some miscelaneous changes with regards to the split between
scalar and sequence quantity classes.

Updates test tolerances to account for changes to base
constants from pint v0.9 defaults to v0.10 defaults (except for
geopotential <-> height conversions which are problematic and resolved
elsewhere).

Finally, resolves issues with xarray v0.13.0 (no implicit casting to
ndarray or inplace argument anymore).

examples/cross_section.py

@@ -42,8 +42,7 @@
 ##############################
 # Get the cross section, and convert lat/lon to supplementary coordinates:
 
-cross = cross_section(data, start, end)
-cross.set_coords(('lat', 'lon'), True)
+cross = cross_section(data, start, end).set_coords(('lat', 'lon'))
 print(cross)
 
 ##############################

metpy/calc/basic.py

@@ -931,8 +931,9 @@ def altimeter_to_station_pressure(altimeter_value, height):
     # N-Value
     n = (mpconsts.Rd * gamma / mpconsts.g).to_base_units()
 
-    return ((altimeter_value ** n - ((p0 ** n * gamma * height) / t0)) ** (1 / n) + (
-            0.3 * units.hPa))
+    return ((altimeter_value ** n
+             - ((p0.to(altimeter_value.units) ** n * gamma * height) / t0)) ** (1 / n)
+            + 0.3 * units.hPa)
 
 
 @exporter.export

metpy/calc/indices.py

@@ -49,10 +49,10 @@ def precipitable_water(dewpt, pressure, bottom=None, top=None):
     dewpt = dewpt[sort_inds]
 
     if top is None:
-        top = np.nanmin(pressure) * pressure.units
+        top = np.nanmin(pressure.magnitude) * pressure.units
 
     if bottom is None:
-        bottom = np.nanmax(pressure) * pressure.units
+        bottom = np.nanmax(pressure.magnitude) * pressure.units
 
     pres_layer, dewpt_layer = get_layer(pressure, dewpt, bottom=bottom, depth=bottom - top)
 
@@ -109,7 +109,8 @@ def mean_pressure_weighted(pressure, *args, **kwargs):
     # function. Said integral works out to this function:
     pres_int = 0.5 * (layer_p[-1].magnitude**2 - layer_p[0].magnitude**2)
     for i, datavar in enumerate(args):
-        arg_mean = np.trapz(layer_arg[i] * layer_p, x=layer_p) / pres_int
+        arg_mean = np.trapz((layer_arg[i] * layer_p).magnitude,
+                            x=layer_p.magnitude) / pres_int
         ret.append(arg_mean * datavar.units)
 
     return ret
@@ -163,7 +164,8 @@ def bunkers_storm_motion(pressure, u, v, heights):
     # Take the cross product of the wind shear and k, and divide by the vector magnitude and
     # multiply by the deviaton empirically calculated in Bunkers (2000) (7.5 m/s)
     shear_cross = concatenate([shear[1], -shear[0]])
-    rdev = shear_cross * (7.5 * units('m/s').to(u.units) / np.hypot(*shear))
+    shear_mag = np.hypot(*(arg.magnitude for arg in shear)) * shear.units
+    rdev = shear_cross * (7.5 * units('m/s').to(u.units) / shear_mag)
 
     # Add the deviations to the layer average wind to get the RM motion
     right_mover = wind_mean + rdev

metpy/calc/kinematics.py

@@ -21,7 +21,7 @@
 
 
 def _stack(arrs):
-    return concatenate([a[np.newaxis] for a in arrs], axis=0)
+    return concatenate([a[np.newaxis] if iterable(a) else a for a in arrs], axis=0)
 
 
 def _is_x_first_dim(dim_order):

metpy/calc/tests/test_basic.py

@@ -396,7 +396,7 @@ def test_add_height_to_pressure():
 def test_add_pressure_to_height():
     """Test the height at pressure above height calculation."""
     height = add_pressure_to_height(110.8286757 * units.m, 100 * units.hPa)
-    assert_almost_equal(height, 988.0028867 * units.meter, 5)
+    assert_almost_equal(height, 988.0028867 * units.meter, 3)
 
 
 def test_sigma_to_pressure():

metpy/calc/tests/test_calc_tools.py

@@ -275,8 +275,8 @@ def get_bounds_data():
 def test_get_bound_pressure_height(pressure, bound, hgts, interp, expected):
     """Test getting bounds in layers with various parameter combinations."""
     bounds = _get_bound_pressure_height(pressure, bound, heights=hgts, interpolate=interp)
-    assert_array_almost_equal(bounds[0], expected[0], 5)
-    assert_array_almost_equal(bounds[1], expected[1], 5)
+    assert_array_almost_equal(bounds[0], expected[0], 4)
+    assert_array_almost_equal(bounds[1], expected[1], 4)
 
 
 def test_get_bound_invalid_bound_units():
@@ -363,8 +363,8 @@ def test_get_layer(pressure, variable, heights, bottom, depth, interp, expected)
     """Test get_layer functionality."""
     p_layer, y_layer = get_layer(pressure, variable, heights=heights, bottom=bottom,
                                  depth=depth, interpolate=interp)
-    assert_array_almost_equal(p_layer, expected[0], 5)
-    assert_array_almost_equal(y_layer, expected[1], 5)
+    assert_array_almost_equal(p_layer, expected[0], 4)
+    assert_array_almost_equal(y_layer, expected[1], 4)
 
 
 def test_greater_or_close():

metpy/calc/tests/test_kinematics.py

@@ -1086,8 +1086,8 @@ def test_q_vector_without_static_stability(q_vector_data):
                           [8.6038280e-14, 4.6968342e-13, -4.6968342e-13, -8.6038280e-14]])
                 * units('m^2 kg^-1 s^-1'))
 
-    assert_almost_equal(q1, q1_truth, 20)
-    assert_almost_equal(q2, q2_truth, 20)
+    assert_almost_equal(q1, q1_truth, 18)
+    assert_almost_equal(q2, q2_truth, 18)
 
 
 def test_q_vector_with_static_stability(q_vector_data):
@@ -1110,8 +1110,8 @@ def test_q_vector_with_static_stability(q_vector_data):
                           [4.4714390e-08, 2.3490489e-07, -2.2635441e-07, -4.0003646e-08]])
                 * units('kg m^-2 s^-3'))
 
-    assert_almost_equal(q1, q1_truth, 14)
-    assert_almost_equal(q2, q2_truth, 14)
+    assert_almost_equal(q1, q1_truth, 12)
+    assert_almost_equal(q2, q2_truth, 12)
 
 
 @pytest.fixture
@@ -1726,5 +1726,5 @@ def test_q_vector_4d(data_4d):
                            [-1.06658890e-13, -2.19817426e-13, -8.35968065e-14, 1.88190788e-13],
                            [-2.27182863e-13, -2.74607819e-13, -1.10587309e-13,
                             -3.88915866e-13]]]]) * units('m^2 kg^-1 s^-1')
-    assert_array_almost_equal(q1, q1_truth, 20)
-    assert_array_almost_equal(q2, q2_truth, 20)
+    assert_array_almost_equal(q1, q1_truth, 18)
+    assert_array_almost_equal(q2, q2_truth, 18)

metpy/calc/tests/test_thermo.py

@@ -710,8 +710,8 @@ def test_cape_cin():
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0])
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 75.7340825 * units('joule / kilogram'), 6)
-    assert_almost_equal(cin, -89.8179205 * units('joule / kilogram'), 6)
+    assert_almost_equal(cape, 75.7340825 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -89.8179205 * units('joule / kilogram'), 2)
 
 
 def test_cape_cin_no_el():
@@ -721,8 +721,8 @@ def test_cape_cin_no_el():
     dewpoint = np.array([19., -11.2, -10.8, -10.4]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC')
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 0.08610409 * units('joule / kilogram'), 6)
-    assert_almost_equal(cin, -89.8179205 * units('joule / kilogram'), 6)
+    assert_almost_equal(cape, 0.08610409 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -89.8179205 * units('joule / kilogram'), 2)
 
 
 def test_cape_cin_no_lfc():
@@ -732,8 +732,8 @@ def test_cape_cin_no_lfc():
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]).to('degC')
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 0.0 * units('joule / kilogram'), 6)
-    assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 6)
+    assert_almost_equal(cape, 0.0 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 2)
 
 
 def test_find_append_zero_crossings():
@@ -980,8 +980,8 @@ def test_surface_based_cape_cin():
     temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     cape, cin = surface_based_cape_cin(p, temperature, dewpoint)
-    assert_almost_equal(cape, 75.7340825 * units('joule / kilogram'), 6)
-    assert_almost_equal(cin, -136.607809 * units('joule / kilogram'), 6)
+    assert_almost_equal(cape, 75.7340825 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, -136.607809 * units('joule / kilogram'), 2)
 
 
 def test_most_unstable_cape_cin_surface():
@@ -990,8 +990,8 @@ def test_most_unstable_cape_cin_surface():
     temperature = np.array([22.2, 14.6, 12., 9.4, 7., -38.]) * units.celsius
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     mucape, mucin = most_unstable_cape_cin(pressure, temperature, dewpoint)
-    assert_almost_equal(mucape, 75.7340825 * units('joule / kilogram'), 6)
-    assert_almost_equal(mucin, -136.607809 * units('joule / kilogram'), 6)
+    assert_almost_equal(mucape, 75.7340825 * units('joule / kilogram'), 2)
+    assert_almost_equal(mucin, -136.607809 * units('joule / kilogram'), 2)
 
 
 def test_most_unstable_cape_cin():
@@ -1151,7 +1151,7 @@ def test_wet_bulb_temperature():
     """Test wet bulb calculation with scalars."""
     val = wet_bulb_temperature(1000 * units.hPa, 25 * units.degC, 15 * units.degC)
     truth = 18.34345936 * units.degC  # 18.59 from NWS calculator
-    assert_almost_equal(val, truth)
+    assert_almost_equal(val, truth, 5)
 
 
 def test_wet_bulb_temperature_1d():
@@ -1162,7 +1162,7 @@ def test_wet_bulb_temperature_1d():
     val = wet_bulb_temperature(pressures, temperatures, dewpoints)
     truth = [21.4449794, 16.7368576, 12.0656909] * units.degC
     # 21.58, 16.86, 12.18 from NWS Calculator
-    assert_array_almost_equal(val, truth)
+    assert_array_almost_equal(val, truth, 5)
 
 
 def test_wet_bulb_temperature_2d():
@@ -1178,7 +1178,7 @@ def test_wet_bulb_temperature_2d():
              [20.5021631, 15.801218, 11.1361878]] * units.degC
     # 21.58, 16.86, 12.18
     # 20.6, 15.9, 11.2 from NWS Calculator
-    assert_array_almost_equal(val, truth)
+    assert_array_almost_equal(val, truth, 5)
 
 
 def test_static_stability_adiabatic():
@@ -1229,8 +1229,8 @@ def test_lfc_not_below_lcl():
                           8.6, 8.1, 7.6, 7., 6.5, 6., 5.4]) * units.degC
     lfc_pressure, lfc_temp = lfc(levels, temperatures, dewpoints)
     # Before patch, LFC pressure would show 1000.5912165339967 hPa
-    assert_almost_equal(lfc_pressure, 811.8263397 * units.mbar, 6)
-    assert_almost_equal(lfc_temp, 6.4992871 * units.celsius, 6)
+    assert_almost_equal(lfc_pressure, 811.8263397 * units.mbar, 3)
+    assert_almost_equal(lfc_temp, 6.4992871 * units.celsius, 3)
 
 
 def test_multiple_lfcs():
@@ -1258,10 +1258,10 @@ def test_multiple_lfcs():
     lfc_pressure_bottom, lfc_temp_bottom = lfc(levels, temperatures, dewpoints,
                                                which='bottom')
     lfc_pressure_all, _ = lfc(levels, temperatures, dewpoints, which='all')
-    assert_almost_equal(lfc_pressure_top, 705.4346277 * units.mbar, 6)
-    assert_almost_equal(lfc_temp_top, 4.8922235 * units.degC, 6)
-    assert_almost_equal(lfc_pressure_bottom, 884.1954356 * units.mbar, 6)
-    assert_almost_equal(lfc_temp_bottom, 13.9597571 * units.degC, 6)
+    assert_almost_equal(lfc_pressure_top, 705.4346277 * units.mbar, 3)
+    assert_almost_equal(lfc_temp_top, 4.8922235 * units.degC, 3)
+    assert_almost_equal(lfc_pressure_bottom, 884.1954356 * units.mbar, 3)
+    assert_almost_equal(lfc_temp_bottom, 13.9597571 * units.degC, 3)
     assert_almost_equal(len(lfc_pressure_all), 2, 0)
 
 
@@ -1289,10 +1289,10 @@ def test_multiple_els():
     el_pressure_top, el_temp_top = el(levels, temperatures, dewpoints)
     el_pressure_bottom, el_temp_bottom = el(levels, temperatures, dewpoints, which='bottom')
     el_pressure_all, _ = el(levels, temperatures, dewpoints, which='all')
-    assert_almost_equal(el_pressure_top, 228.0575059 * units.mbar, 6)
-    assert_almost_equal(el_temp_top, -56.8123126 * units.degC, 6)
-    assert_almost_equal(el_pressure_bottom, 849.7942185 * units.mbar, 6)
-    assert_almost_equal(el_temp_bottom, 12.4233265 * units.degC, 6)
+    assert_almost_equal(el_pressure_top, 228.0575059 * units.mbar, 3)
+    assert_almost_equal(el_temp_top, -56.8123126 * units.degC, 3)
+    assert_almost_equal(el_pressure_bottom, 849.7942185 * units.mbar, 3)
+    assert_almost_equal(el_temp_bottom, 12.4233265 * units.degC, 3)
     assert_almost_equal(len(el_pressure_all), 2, 0)
 
 
@@ -1303,8 +1303,8 @@ def test_cape_cin_custom_profile():
     dewpoint = np.array([19., -11.2, -10.8, -10.4, -10., -53.2]) * units.celsius
     parcel_prof = parcel_profile(p, temperature[0], dewpoint[0]) + 5 * units.delta_degC
     cape, cin = cape_cin(p, temperature, dewpoint, parcel_prof)
-    assert_almost_equal(cape, 1443.505086499895 * units('joule / kilogram'), 6)
-    assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 6)
+    assert_almost_equal(cape, 1443.505086499895 * units('joule / kilogram'), 2)
+    assert_almost_equal(cin, 0.0 * units('joule / kilogram'), 2)
 
 
 def test_parcel_profile_below_lcl():

metpy/calc/thermo.py

@@ -264,7 +264,7 @@ def dt(t, p):
         frac = ((mpconsts.Rd * t + mpconsts.Lv * rs)
                 / (mpconsts.Cp_d + (mpconsts.Lv * mpconsts.Lv * rs * mpconsts.epsilon
                                     / (mpconsts.Rd * t * t)))).to('kelvin')
-        return frac / p
+        return (frac / p).magnitude
 
     if ref_pressure is None:
         ref_pressure = pressure[0]
@@ -287,14 +287,14 @@ def dt(t, p):
 
     if ref_pressure > pressure.min():
         # Integrate downward in pressure
-        pres_down = np.append(ref_pressure, pressure[(ref_pres_idx - 1)::-1])
-        trace_down = si.odeint(dt, temperature.squeeze(), pres_down.squeeze())
+        pres_down = np.append(ref_pressure.m, pressure[(ref_pres_idx - 1)::-1].m)
+        trace_down = si.odeint(dt, temperature.m.squeeze(), pres_down.squeeze())
         ret_temperatures = np.concatenate((ret_temperatures, trace_down[:0:-1]))
 
     if ref_pressure < pressure.max():
         # Integrate upward in pressure
-        pres_up = np.append(ref_pressure, pressure[ref_pres_idx:])
-        trace_up = si.odeint(dt, temperature.squeeze(), pres_up.squeeze())
+        pres_up = np.append(ref_pressure.m, pressure[ref_pres_idx:].m)
+        trace_up = si.odeint(dt, temperature.m.squeeze(), pres_up.squeeze())
         ret_temperatures = np.concatenate((ret_temperatures, trace_up[1:]))
 
     if pres_decreasing:
@@ -628,7 +628,7 @@ def _parcel_profile_helper(pressure, temperature, dewpt):
     temp_lower = dry_lapse(press_lower, temperature)
 
     # If the pressure profile doesn't make it to the lcl, we can stop here
-    if _greater_or_close(np.nanmin(pressure), press_lcl.m):
+    if _greater_or_close(np.nanmin(pressure.m), press_lcl.m):
         return (press_lower[:-1], press_lcl, np.array([]) * press_lower.units,
                 temp_lower[:-1], temp_lcl, np.array([]) * temp_lower.units)
 
@@ -1449,17 +1449,17 @@ def cape_cin(pressure, temperature, dewpt, parcel_profile):
 
     # CAPE
     # Only use data between the LFC and EL for calculation
-    p_mask = _less_or_close(x, lfc_pressure) & _greater_or_close(x, el_pressure)
-    x_clipped = x[p_mask]
-    y_clipped = y[p_mask]
+    p_mask = _less_or_close(x.m, lfc_pressure) & _greater_or_close(x.m, el_pressure)
+    x_clipped = x[p_mask].magnitude
+    y_clipped = y[p_mask].magnitude
     cape = (mpconsts.Rd
             * (np.trapz(y_clipped, np.log(x_clipped)) * units.degK)).to(units('J/kg'))
 
     # CIN
     # Only use data between the surface and LFC for calculation
-    p_mask = _greater_or_close(x, lfc_pressure)
-    x_clipped = x[p_mask]
-    y_clipped = y[p_mask]
+    p_mask = _greater_or_close(x.m, lfc_pressure)
+    x_clipped = x[p_mask].magnitude
+    y_clipped = y[p_mask].magnitude
     cin = (mpconsts.Rd
            * (np.trapz(y_clipped, np.log(x_clipped)) * units.degK)).to(units('J/kg'))
 
@@ -1498,7 +1498,7 @@ def _find_append_zero_crossings(x, y):
     y = y[sort_idx]
 
     # Remove duplicate data points if there are any
-    keep_idx = np.ediff1d(x, to_end=[1]) > 1e-6
+    keep_idx = np.ediff1d(x.magnitude, to_end=[1]) > 1e-6
     x = x[keep_idx]
     y = y[keep_idx]
     return x, y
@@ -1638,7 +1638,7 @@ def _isen_iter(iter_log_p, isentlevs_nd, ka, a, b, pok):
     slices = [np.newaxis] * ndim
     slices[axis] = slice(None)
     slices = tuple(slices)
-    pres = np.broadcast_to(pres[slices], temperature.shape) * pres.units
+    pres = np.broadcast_to(pres[slices].magnitude, temperature.shape) * pres.units
 
     # Sort input data
     sort_pres = np.argsort(pres.m, axis=axis)
@@ -1904,7 +1904,7 @@ def mixed_layer(p, *args, **kwargs):
     ret = []
     for datavar_layer in datavars_layer:
         actual_depth = abs(p_layer[0] - p_layer[-1])
-        ret.append((-1. / actual_depth.m) * np.trapz(datavar_layer, p_layer)
+        ret.append((-1. / actual_depth.m) * np.trapz(datavar_layer.m, p_layer.m)
                    * datavar_layer.units)
     return ret
 
@@ -2047,7 +2047,7 @@ def thickness_hydrostatic(pressure, temperature, **kwargs):
 
     # Take the integral (with unit handling) and return the result in meters
     return (- mpconsts.Rd / mpconsts.g * np.trapz(
-        layer_virttemp.to('K'), x=np.log(layer_p / units.hPa)) * units.K).to('m')
+        layer_virttemp.m_as('K'), x=np.log(layer_p.m_as('hPa'))) * units.K).to('m')
 
 
 @exporter.export
@@ -2306,7 +2306,7 @@ def static_stability(pressure, temperature, axis=0):
     """
     theta = potential_temperature(pressure, temperature)
 
-    return - mpconsts.Rd * temperature / pressure * first_derivative(np.log(theta / units.K),
+    return - mpconsts.Rd * temperature / pressure * first_derivative(np.log(theta.m_as('K')),
                                                                      x=pressure, axis=axis)