Merge pull request #1144 from jthielen/pint-patch-10

Updates to pint/xarray master branch (and xarray v0.13) compatability

.travis.yml

@@ -61,9 +61,25 @@ matrix:
       env: TASK="docs" PRE="--pre"
     - python: nightly
       env: PRE="--pre"
+    - python: 3.7
+      env:
+        - TASK="coverage"
+        - VERSIONS="git+git://github.com/hgrecco/pint@master#egg=pint git+git://github.com/pydata/xarray@master#egg=xarray"
+    - python: 3.7
+      env:
+        - TASK="docs"
+        - VERSIONS="git+git://github.com/hgrecco/pint@master#egg=pint git+git://github.com/pydata/xarray@master#egg=xarray"
   allow_failures:
     - python: "3.7-dev"
     - python: nightly
+    - python: 3.7
+      env:
+        - TASK="coverage"
+        - VERSIONS="git+git://github.com/hgrecco/pint@master#egg=pint git+git://github.com/pydata/xarray@master#egg=xarray"
+    - python: 3.7
+      env:
+        - TASK="docs"
+        - VERSIONS="git+git://github.com/hgrecco/pint@master#egg=pint git+git://github.com/pydata/xarray@master#egg=xarray"
 
 before_install:
   # We hard-code the sphinx_rtd_theme to lock in our build with patch for

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)