lfc() fails with parcel_temperature_profile argument

[sgdecker]

Here is some code triggering the error:

import numpy as np
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import metpy

print(metpy.__version__)

P = np.array([ 1024.95703125,  1016.61474609,  1005.33056641,   991.08544922,   973.4163208,
               951.3381958,    924.82836914,   898.25482178,   873.46124268,   848.69830322,
               823.92553711,   788.49304199,   743.44580078,   700.50970459,   659.62017822,
               620.70861816,   583.69421387,   548.49719238,   515.03826904,   483.24401855,
               453.0418396,    424.36477661,   397.1505127,    371.33441162,   346.85922241,
               323.66995239,   301.70935059,   280.92651367,   261.27053833,   242.69168091,
               225.14237976,   208.57781982,   192.95333862,   178.22599792,   164.39630127,
               151.54336548,   139.68635559,   128.74923706,   118.6588974,    109.35111237,
               100.76405334,    92.84288025,    85.53556824,    78.79430389,    72.57549286,
               66.83885193,    61.54678726,    56.66480637,    52.16108322]) * units('hPa')

T = np.array([  6.00750732,   5.14892578,   4.177948,     3.00268555,   1.55535889,  -0.25527954,
                -1.93988037,  -3.57766724,  -4.40600586,  -4.19238281,  -3.71185303,  -4.47943115,
                -6.81280518,  -8.08685303,  -8.41287231, -10.79302979, -14.13262939, -16.85784912,
                -19.51675415, -22.28689575, -24.99938965, -27.79664612, -30.90414429, -34.49435425,
                -38.438797, -42.27981567, -45.99230957, -49.75340271, -53.58230591, -57.30686951,
                -60.76026917, -63.92070007, -66.72470093, -68.97846985, -70.4264679, -71.16407776,
                -71.53797913, -71.64375305, -71.52735901, -71.53523254, -71.61097717, -71.92687988,
                -72.68682861, -74.129776,   -76.02471924, -76.88977051, -76.26008606, -75.90351868,
                -76.15809631]) * units('degC')

Td = np.array([4.50012302,   3.42483997,   2.78102994,   2.24474645,   1.593485,   0.94408154,
               -3.8044982,   -3.55629468,  -9.7376976,  -10.2950449,  -9.67498302, -10.30486488,
               -8.70559597,  -8.71669006, -12.66509628, -18.6697197,  -23.00351334, -29.46240425,
               -36.82178497, -41.68824768, -44.50320816, -48.54426575, -52.50753403, -51.09564209,
               -48.92690659, -49.97380829, -51.57516098, -52.62096405, -54.24332809, -57.09109879,
               -60.5596199,  -63.93486404, -67.07530212, -70.01263428, -72.9258728, -76.12271881,
               -79.49847412, -82.2350769,  -83.91127014, -84.95665741, -85.61238861, -86.16391754,
               -86.7653656,  -87.34436035, -87.87495422, -88.34281921, -88.74453735, -89.04680634,
               -89.26436615]) * units('degC')

mlP, mlT, mlTd = mpcalc.mixed_parcel(P, T, Td)
print(mlP, mlT, mlTd)

ml_profile = mpcalc.parcel_profile(P, mlT, mlTd)
print(ml_profile)
print(len(P), len(T), len(Td), len(ml_profile))

mllfc_p, mllfc_T = mpcalc.lfc(P, T, Td, parcel_temperature_profile=ml_profile)
print(mllfc_p, mllfc_T)

I expect information about the MLLFC to appear, but instead, this is my output:

0.9.1+6.gb7f97991
1024.95703125 hectopascal 5.743482661750363 degC 2.08279974748133 degC
[278.89348266 278.24322133 277.35755297 276.22928809 274.81361678 273.52555096 272.08325952 270.57324142 269.10199475 267.56801846 265.96421247 263.5399307  260.21335861 256.7576317  253.17725629 249.48027885 245.67713199 241.78102546 237.80684667 233.77057953 229.68754004 225.57206704 221.43686749 217.29167908 213.14477576 209.0024604  204.86851485 200.74607368 196.63681225 192.54139059 188.45988434 184.39212913 180.33722545 176.29398272 172.27346048 168.31396424 164.44241275 160.65710232 156.95530012 153.33548532 149.79516513 146.33278761 142.94622388 139.63343896 136.39268027 133.22221317 130.12009053 127.08448424 124.11358965] kelvin
49 49 49 49
Traceback (most recent call last):
  File "lfc_bug.py", line 47, in <module>
    mllfc_p, mllfc_T = mpcalc.lfc(P, T, Td, parcel_temperature_profile=ml_profile)
  File "/home/decker/src/git_repos/metpy/metpy/xarray.py", line 381, in wrapper
    return func(*args, **kwargs)
  File "/home/decker/src/git_repos/metpy/metpy/units.py", line 305, in wrapper
    return func(*args, **kwargs)
  File "/home/decker/src/git_repos/metpy/metpy/calc/thermo.py", line 383, in lfc
    return x[0], y[0]
  File "/home/decker/local/anaconda3/envs/devel/lib/python3.6/site-packages/pint/quantity.py", line 1281, in __getitem__
    value = self._magnitude[key]
IndexError: index 0 is out of bounds for axis 0 with size 0

I also tried the following variation, but it failed in the same way:

import numpy as np
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import metpy

print(metpy.__version__)

P = np.array([ 1024.95703125,  1016.61474609,  1005.33056641,   991.08544922,   973.4163208,
               951.3381958,    924.82836914,   898.25482178,   873.46124268,   848.69830322,
               823.92553711,   788.49304199,   743.44580078,   700.50970459,   659.62017822,
               620.70861816,   583.69421387,   548.49719238,   515.03826904,   483.24401855,
               453.0418396,    424.36477661,   397.1505127,    371.33441162,   346.85922241,
               323.66995239,   301.70935059,   280.92651367,   261.27053833,   242.69168091,
               225.14237976,   208.57781982,   192.95333862,   178.22599792,   164.39630127,
               151.54336548,   139.68635559,   128.74923706,   118.6588974,    109.35111237,
               100.76405334,    92.84288025,    85.53556824,    78.79430389,    72.57549286,
               66.83885193,    61.54678726,    56.66480637,    52.16108322]) * units('hPa')

T = np.array([  6.00750732,   5.14892578,   4.177948,     3.00268555,   1.55535889,  -0.25527954,
                -1.93988037,  -3.57766724,  -4.40600586,  -4.19238281,  -3.71185303,  -4.47943115,
                -6.81280518,  -8.08685303,  -8.41287231, -10.79302979, -14.13262939, -16.85784912,
                -19.51675415, -22.28689575, -24.99938965, -27.79664612, -30.90414429, -34.49435425,
                -38.438797, -42.27981567, -45.99230957, -49.75340271, -53.58230591, -57.30686951,
                -60.76026917, -63.92070007, -66.72470093, -68.97846985, -70.4264679, -71.16407776,
                -71.53797913, -71.64375305, -71.52735901, -71.53523254, -71.61097717, -71.92687988,
                -72.68682861, -74.129776,   -76.02471924, -76.88977051, -76.26008606, -75.90351868,
                -76.15809631]) * units('degC')

Td = np.array([4.50012302,   3.42483997,   2.78102994,   2.24474645,   1.593485,   0.94408154,
               -3.8044982,   -3.55629468,  -9.7376976,  -10.2950449,  -9.67498302, -10.30486488,
               -8.70559597,  -8.71669006, -12.66509628, -18.6697197,  -23.00351334, -29.46240425,
               -36.82178497, -41.68824768, -44.50320816, -48.54426575, -52.50753403, -51.09564209,
               -48.92690659, -49.97380829, -51.57516098, -52.62096405, -54.24332809, -57.09109879,
               -60.5596199,  -63.93486404, -67.07530212, -70.01263428, -72.9258728, -76.12271881,
               -79.49847412, -82.2350769,  -83.91127014, -84.95665741, -85.61238861, -86.16391754,
               -86.7653656,  -87.34436035, -87.87495422, -88.34281921, -88.74453735, -89.04680634,
               -89.26436615]) * units('degC')

mlP, mlT, mlTd = mpcalc.mixed_parcel(P, T, Td)
print(mlP, mlT, mlTd)

P[0], T[0], Td[0] = mlP, mlT, mlTd

P2, T2, Td2, ml_profile = mpcalc.parcel_profile_with_lcl(P, T, Td)
print(len(P2), len(T2), len(Td2), len(ml_profile))

mllfc_p, mllfc_T = mpcalc.lfc(P2, T2, Td2, parcel_temperature_profile=ml_profile)
print(mllfc_p, mllfc_T)

[sgdecker]

A clue:
find_intersections finds one intersection for this data, at a pressure of 1009.2 hPa, but the computed LCL is at 1001.5 hPa. The inequality here thus returns False, triggering the bug.

My hypothesis is that the calls to lcl in this function should use parcel_temperature_profile[0] rather than temperature[0].

Will make that change and see what happens!

[sgdecker]

Nope, the LCL is now 1005 hPa. Looks like there needs to be code to set the LFC to the LCL when pLFC > pLCL.

[sgdecker]

I have code now that doesn't crash for this example, but breaks a test. I'm not sure the test is correct, though. I'll submit a PR for feedback.