MetPy broken by pint 0.9

[markmuetz]

I was successfully calculating CAPE and CIN with an older version of metpy. The profiles are taken from a mean state of an idealized cloud resolving model. On upgrade to a newer version (v0.10.0), I started seeing a problem (see Actual output below). I don't think there's a problem with my profiles as it was working previously. I've worked out where the bug was introduced - see the expected and actual output commits. The diff between the commits looked fairly extensive though so I couldn't work out what the problem was!

# coding: utf-8
import numpy as np
import metpy.calc as mpcalc
from metpy.units import units

p_profile = np.array([ 99731.40672841,  99012.93988483,  98287.88833531,  97546.8524674 ,
    96777.36391641,  95964.53852765,  95092.95684841,  94149.52201132,
    93125.46324786,  92017.03931377,  90823.55098039,  89543.68226431,
    88172.29988771,  86699.1136601 ,  85111.1941357 ,  83398.47584964,
    81559.42072823,  79604.8563692 ,  77557.7790759 ,  75449.20161866,
    73312.07613362,  71175.14623916,  69059.73471423,  66980.00666125,
    64944.56262355,  62957.06524521,  61017.06243663,  59125.39164667,
    57281.97692881,  55486.35911778,  53737.6394924 ,  52034.68456921,
    50376.26195021,  48761.2365222 ,  47188.53315282,  45657.07499016,
    44165.77122667,  42713.50509724,  41299.23639793,  39922.00321203,
    38580.86339437,  37274.89576525,  36003.23062364,  34765.01566383,
    33559.42640456,  32385.67271967,  31243.00589608,  30130.71333169,
    29048.10915457,  27994.54529059,  26969.41383295,  25972.13470436,
    25002.15734616,  24058.95931997,  23142.03833883,  22250.90317508,
    21385.06413239,  20544.03870877,  19727.34792228,  18934.56937955,
    18165.222981  ,  17419.72499893,  16697.63397071,  15998.39190248,
    15321.11449205,  14664.17223577,  14023.91945431,  13393.32474551,
    12761.08598215,  12110.50906806,  11420.61415495,  10671.79230543,
    9851.26937387,   8965.08588726,   8038.62814417,   7106.86959838,
    6205.94601236,   5366.5115902 ,   4609.25652035,   3945.73599008,
    3373.40907093,   2883.69777091,   2466.05682486,   2110.23178103,
    1806.63466476,   1547.5486207 ,   1325.09426688,   1137.88414121,
    979.78621433,    845.76894618,    730.8180016 ,    634.74196791,
    554.20262164,    485.85214384,    427.233084  ,    376.19262967,
    331.4337417 ,    291.91256049,    257.01114311])

T_profile = np.array([ 297.88216024,  297.26764764,  296.65196701,  296.03670196,
    295.39354615,  294.72073495,  294.00409131,  293.24006615,
    292.42457574,  291.56661565,  290.67420312,  289.74522486,
    288.77161911,  287.74363449,  286.6552643 ,  285.50326008,
    284.28942922,  283.0232586 ,  281.7408227 ,  280.49193362,
    279.29388391,  278.06855099,  276.72343495,  275.31697755,
    273.95124832,  272.53171796,  270.84519572,  269.28773695,
    267.81665268,  266.40268414,  264.99541812,  263.55980528,
    262.07484043,  260.5509311 ,  258.99223251,  257.39246525,
    255.74414459,  254.03528329,  252.27274778,  250.46484962,
    248.60996504,  246.70679961,  244.75953215,  242.76608658,
    240.72712625,  238.64492011,  236.52352912,  234.36643591,
    232.17562005,  229.95590655,  227.71308036,  225.45059003,
    223.17285164,  220.88446062,  218.58858736,  216.28651421,
    213.97720399,  211.66094137,  209.33702652,  207.01898461,
    204.69056742,  202.60783781,  200.56208789,  198.56303812,
    196.65661834,  194.92303101,  193.37052771,  192.02912928,
    191.07334475,  190.49447264,  190.10126068,  189.99390591,
    189.86897907,  190.65551486,  192.71116526,  195.2892728 ,
    197.79736664,  200.10150056,  202.35481869,  205.49557229,
    208.39005467,  210.76135549,  212.59158953,  214.04803363,
    214.97290906,  215.81884392,  215.28176629,  219.44640798,
    223.53130664,  227.41802446,  229.02162464,  237.54145616,
    246.92800635,  254.70113502,  260.85998109,  263.66823025,
    264.84128327,  264.19585255,  263.4474366 ])

Td_profile = np.array([  1.84993136e+01,   1.83990128e+01,   1.80350756e+01,
    1.78648581e+01,   1.76550614e+01,   1.74167816e+01,
    1.70949046e+01,   1.66780604e+01,   1.61550825e+01,
    1.55049912e+01,   1.47250689e+01,   1.38661200e+01,
    1.29599474e+01,   1.20161241e+01,   1.10249275e+01,
    9.97909842e+00,   8.85500100e+00,   7.62307849e+00,
    6.23219929e+00,   4.63567267e+00,   2.91344167e+00,
    1.29529983e+00,  -9.16774322e-02,  -1.41783117e+00,
    -2.81993735e+00,  -4.08960092e+00,  -5.12625444e+00,
    -6.39948701e+00,  -7.85507497e+00,  -9.45793521e+00,
    -1.11417929e+01,  -1.28568797e+01,  -1.45949085e+01,
    -1.63995801e+01,  -1.82750345e+01,  -2.01964984e+01,
    -2.21988241e+01,  -2.42426426e+01,  -2.63251545e+01,
    -2.84502144e+01,  -3.06267387e+01,  -3.28925518e+01,
    -3.52345816e+01,  -3.75635832e+01,  -3.98531663e+01,
    -4.21356269e+01,  -4.43935769e+01,  -4.65902982e+01,
    -4.87199139e+01,  -5.08274384e+01,  -5.29529931e+01,
    -5.50926325e+01,  -5.72436566e+01,  -5.93384711e+01,
    -6.13493619e+01,  -6.33398863e+01,  -6.53582740e+01,
    -6.73937865e+01,  -6.93380560e+01,  -7.12644617e+01,
    -7.32610179e+01,  -7.51962068e+01,  -7.70251001e+01,
    -7.87856463e+01,  -8.04500545e+01,  -8.19548546e+01,
    -8.32098918e+01,  -8.42502644e+01,  -8.50160692e+01,
    -8.55792584e+01,  -8.60612449e+01,  -8.65586201e+01,
    -8.71538847e+01,  -8.78787417e+01,  -8.87433130e+01,
    -8.97310365e+01,  -9.08195544e+01,  -9.19865565e+01,
    -9.32082321e+01,  -9.44569536e+01,  -9.57126923e+01,
    -9.69620366e+01,  -9.81969865e+01,  -9.94166767e+01,
    -1.00619768e+02,  -1.01804824e+02,  -1.02977761e+02,
    -1.04132423e+02,  -1.05265633e+02,  -1.06381712e+02,
    -1.07479662e+02,  -1.08658343e+02,  -1.10052896e+02,
    -1.11889563e+02,  -1.14756163e+02,           np.nan,
    np.nan,              np.nan,              np.nan])

print(mpcalc.surface_based_cape_cin(p_profile * units('Pa'), T_profile * units('K'), Td_profile * units('degC')))

• Problem description: I believe you should be able to calculate CAPE/CIN for the given profiles.

• Expected output (eb60959):

/home/markmuetz/anaconda3/envs/metpy_v10.0_test_minimal/lib/python3.6/site-packages/pint/quantity.py:1377: UnitStrippedWarning: The unit of the quantity is stripped.
  warnings.warn("The unit of the quantity is stripped.", UnitStrippedWarning)
(<Quantity(1547.307203463148, 'joule / kilogram')>, <Quantity(-6.785224296243305, 'joule / kilogram')>)

• Actual output (abcb7e9 and later)

/home/markmuetz/anaconda3/envs/metpy_v10.0_test_minimal/lib/python3.6/site-packages/pint/quantity.py:1377: UnitStrippedWarning: The unit of the quantity is stripped.
  warnings.warn("The unit of the quantity is stripped.", UnitStrippedWarning)
Traceback (most recent call last):
  File "cape_calc.py", line 92, in <module>
    print(mpcalc.surface_based_cape_cin(p_profile * units('Pa'), T_profile * units('K'), Td_profile * units('degC')))
  File "/home/markmuetz/projects/MetPy/metpy/xarray.py", line 381, in wrapper
    return func(*args, **kwargs)
  File "/home/markmuetz/projects/MetPy/metpy/units.py", line 305, in wrapper
    return func(*args, **kwargs)
  File "/home/markmuetz/projects/MetPy/metpy/calc/thermo.py", line 1652, in surface_based_cape_cin
    p, t, td, profile = parcel_profile_with_lcl(pressure, temperature, dewpoint)
  File "/home/markmuetz/projects/MetPy/metpy/xarray.py", line 381, in wrapper
    return func(*args, **kwargs)
  File "/home/markmuetz/projects/MetPy/metpy/units.py", line 305, in wrapper
    return func(*args, **kwargs)
  File "/home/markmuetz/projects/MetPy/metpy/calc/thermo.py", line 514, in parcel_profile_with_lcl
    new_temp = _insert_lcl_level(pressure, temperature, p_lcl)
  File "/home/markmuetz/projects/MetPy/metpy/calc/thermo.py", line 552, in _insert_lcl_level
    interp_temp = interpolate_1d(lcl_pressure, pressure, temperature)
  File "/home/markmuetz/projects/MetPy/metpy/xarray.py", line 381, in wrapper
    return func(*args, **kwargs)
  File "/home/markmuetz/anaconda3/envs/metpy_v10.0_test_minimal/lib/python3.6/site-packages/pint/registry_helpers.py", line 194, in wrapper
    result = func(*new_values, **kw)
  File "/home/markmuetz/projects/MetPy/metpy/interpolate/one_dimension.py", line 119, in interpolate_1d
    variables = [arr[sorter] for arr in args]
  File "/home/markmuetz/projects/MetPy/metpy/interpolate/one_dimension.py", line 119, in <listcomp>
    variables = [arr[sorter] for arr in args]
TypeError: tuple indices must be integers or slices, not tuple

• Versions:
  • Python 3.6.6
  • `python -c 'import metpy; print(metpy.version)' : 0.8.0+101.gabcb7e9 (first commit with bug)

[am-thyst]

Thanks @markmuetz - thought I was going crazy. Could've sworn it was working a couple of weeks ago. Time to downgrade metpy until the bug is fixed?

edit: downgrading metpy just brought up another error. Keeping metpy as 0.10.0 and downgrading pint to 0.8.1 as mentioned in #998 seems to fix it with a quick test:

from metpy.calc.thermo import surface_based_cape_cin
from metpy.units import units
t = [288, 284, 278] * units.kelvin
dp = [285, 278, 269] * units.kelvin
p = [1013, 950, 900] * units.hPa

surface_based_cape_cin(p, t, dp)
(<Quantity(21.53649906623308, 'joule / kilogram')>, <Quantity(-8.118571674744551, 'joule / kilogram')>)

Now to see if it works on my entire dataset

[markmuetz]

Thanks @am-thyst - downgrading pint to 0.8.1 worked. For me this was: conda install pint=0.8.1.

[dopplershift]

Leaving this open to track our problem with pint 0.9.

[jthielen]

For whatever reason, I got side-tracked into trying to figure out what was wrong with MetPy and pint 0.9. From what I could tell, I got it traced down to two separate issues:

• New usage of inspect.signature in pint messes with *args when used with units.wraps decorator (this breaks the interpolation_1d functions)
• Single value quantities are currently evaluated as iterable (see iter for single value.... hgrecco/pint#751, this breaks metpy.calc.kinematics._check_and_flip, and in turn a good part of the kinematics tests for dim_order='xy')

I have a branch up with workarounds for both these issues (https://github.com/jthielen/MetPy/tree/pint-0.9). Almost all the tests are passing for me locally, just two matplotlib tests are failing on the second issue.

With that in mind, should we push for hgrecco/pint#751 to get fixed upstream? Also, is the workaround for the units.wraps issue okay, or should a better fix be implemented here or upstream?

[dopplershift]

@jthielen So I have fixes for the pint units.wraps issue, I was just needing to get to tests so I could push upstream. Then my plan was to make MetPy monkeypatch pint to fix. I want to fix the issue in pint, but I'm undecided if the best path forward for MetPy is my monkeypatch or your change to no longer use the problematic wrapper.

I didn't see, what's your solution for the other issue?

[jthielen]

Cool! I simply wrapped np.iterable (in cbook.py) to have a special case for Quantities:

def iterable(value):
    """Determine if value can be iterated over."""
    # Special case for pint Quantities
    if hasattr(value, 'magnitude'):
        value = value.magnitude
    return np.iterable(value)

[dopplershift]

@jthielen Can you put in a PR with your fixes for pint? That PR could also remove the restriction against pint 0.9 then.

I'm thinking we need a 0.10.1 right now with fixes for everything upstream that's broken us, esp. xarray.

[jthielen]

Sure! I should be able to get it up soon.