Bug in shading CAPE/CIN or operator error?

[jsillin]

Hey MetPy friends...

I've been running lots of skew-t plots recently for the skewmap project. Most of them plot exactly as expected. I must have generated thousands of skew-t's covering all shapes and sizes. I keep getting this odd error on one of my scripts plotting for the TX panhandle region:

Traceback (most recent call last):
  File "get_soundingmaps.py", line 389, in <module>
    smap.plot_soundings(fig,ax1,prs_temps,prs_relh,sfc_pressure,centerlat,centerlon,domainsize,model,cape=cape_bool,wetbulb=wetbulb_bool)
  File "C:\Users\Jack\Documents\GitHub\plotting\dev\soundingmaps.py", line 144, in plot_soundings
    cape = mpcalc.cape_cin(pres_abv_ground,abv_sfc_temp.data*units.degC,abv_sfc_dewp,parcel_prof)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\xarray.py", line 1206, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\units.py", line 246, in wrapper
    return func(*args, **kwargs)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\calc\thermo.py", line 1786, in cape_cin
    lfc_pressure, _ = lfc(pressure, temperature, dewpoint,
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\xarray.py", line 1206, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\units.py", line 246, in wrapper
    return func(*args, **kwargs)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\calc\thermo.py", line 524, in lfc
    if np.min(el_pres) > this_lcl[0]:
  File "<__array_function__ internals>", line 5, in amin
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\quantity.py", line 1641, in __array_function__
    return numpy_wrap("function", func, args, kwargs, types)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\numpy_func.py", line 904, in numpy_wrap
    return handled[name](*args, **kwargs)
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\pint\numpy_func.py", line 736, in implementation
    ret = func(*bound_args.args, **bound_args.kwargs)
  File "<__array_function__ internals>", line 5, in amin
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\numpy\core\fromnumeric.py", line 2830, in amin
    return _wrapreduction(a, np.minimum, 'min', axis, None, out,
  File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\numpy\core\fromnumeric.py", line 87, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
ValueError: zero-size array to reduction operation minimum which has no identity

The same run of the same script produced many nice maps, including several with CAPE/CIN shading:

(Plotting 18z fails).

This seems to suggest some data may be missing. But all the points are good, and I have the same script running over a larger domain inclusive of this area that happily plots 18z and beyond:

So I'm not totally sure what's going on. I do have a theory though. I only select model grids at and below 250mb so I don't have to deal with NaNs in the stratosphere and I'm not plotting that data anyway (only at and below 300mb) so I don't need it. Does MetPy require you have environmental data at/below the EL to shade CAPE/CIN? Seems weird that it would have taken me this long to find that out given how many convective soundings I've plotted but it's possible.

Any thoughts would be much appreciated!

[DanielAdriaansen]

Based on your Traceback (thanks for posting!) what I would do (since I don't know the answer either), is open this file:
File "C:\Users\Jack\anaconda3\envs\metpynew\lib\site-packages\metpy\calc\thermo.py"

And on Line 523 (one line before the error of 524), I would print the value of el_pres and this_lcl to see what they were. Then based on what came out (I suspect one or both may not have anything useful in them), I'd read the code around line 524 because there's a bunch of comments about assumptions of the LFC/LCL and EL in there.

One thing I find interesting, and I am not a thermo expert so there may be good reasons for this, is that cape_cin() uses the el() function to find EL, but the lfc() function called by cape_cin() does not use el(), rather it uses find_intersection() to find the EL(s). And your error is coming from the lfc() function so perhaps this is related to find_intersection()?

cape_cin():

MetPy/src/metpy/calc/thermo.py

Line 1711 in ada4ab2

def cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc='bottom',

lfc():

MetPy/src/metpy/calc/thermo.py

Line 419 in ada4ab2

def lfc(pressure, temperature, dewpoint, parcel_temperature_profile=None, dewpoint_start=None,

el():

MetPy/src/metpy/calc/thermo.py

Line 609 in ada4ab2

def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'):

find_intersection():

MetPy/src/metpy/calc/tools.py

Line 96 in ada4ab2

def find_intersections(x, a, b, direction='all', log_x=False):

I'll be following along and I'm curious what el_pres and this_lcl are in this case!

[dopplershift]

I wonder if this is related to #1735 . If you can reduce that down to a specific sounding we can use to reproduce, that would be helpful.

[jsillin]

Different model run so different sounding, but still the same error. Here's one sounding for which this error happens:

temps = [  9.34383798,   8.65136494,   7.01436763,   4.92487862,
         2.70710265,   0.47084482,  -1.8323332 ,  -4.39423443,
        -7.31241262, -10.34799136, -13.16191452, -16.77846247,
       -20.07967071, -23.53162645, -27.54043772, -32.19886153,
       -35.76403417, -39.5579374 , -43.41365981] *units.degC

pressures = [725.0, 700.0, 675.0, 650.0, 625.0, 600.0, 575.0, 550.0, 
        525.0, 500.0, 475.0, 450.0, 425.0, 400.0, 375.0, 350.0, 325.0, 300.0, 275.0]*units.hPa

dews = [-2.021389825590112, -15.228942071818679, -12.952620103490942, -13.132880753606077, -16.67285296152284, 
        -23.66782105106919, -27.75826410494708, -25.54051936892227, -21.55679924139037, -19.35074135957369, 
        -20.31607214051937, -19.146652748913436, -20.84407397974722, -24.44277893446173, -29.397610745430622, 
       -32.53159420035021, -43.00330760511325, -54.216772837434945, -59.706861547892]*units.degC

This seems to be a problem with the next line after this which is:

mu_parcel_prof = mpcalc.parcel_profile(mupressures,muparcel[1],muparcel[2])

returning the error:

----------------------------------------------------------
ValueError               Traceback (most recent call last)
<ipython-input-7-1c0d609674cd> in <module>
----> 1 mu_parcel_prof = mpcalc.parcel_profile(pressures,temps,dews)
      2 print(mu_parcel_prof)

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\xarray.py in wrapper(*args, **kwargs)
   1204 
   1205             # Evaluate inner calculation
-> 1206             result = func(*bound_args.args, **bound_args.kwargs)
   1207 
   1208             # Wrap output based on match and match_unit

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\units.py in wrapper(*args, **kwargs)
    244                                'that the function is being called properly.\n') + msg
    245                 raise ValueError(msg)
--> 246             return func(*args, **kwargs)
    247 
    248         return wrapper

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\calc\thermo.py in parcel_profile(pressure, temperature, dewpoint)
    724 
    725     """
--> 726     _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint)
    727     return concatenate((t_l, t_u))
    728 

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\calc\thermo.py in _parcel_profile_helper(pressure, temperature, dewpoint)
    869     """
    870     # Find the LCL
--> 871     press_lcl, temp_lcl = lcl(pressure[0], temperature, dewpoint)
    872     press_lcl = press_lcl.to(pressure.units)
    873 

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\xarray.py in wrapper(*args, **kwargs)
   1204 
   1205             # Evaluate inner calculation
-> 1206             result = func(*bound_args.args, **bound_args.kwargs)
   1207 
   1208             # Wrap output based on match and match_unit

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\units.py in wrapper(*args, **kwargs)
    244                                'that the function is being called properly.\n') + msg
    245                 raise ValueError(msg)
--> 246             return func(*args, **kwargs)
    247 
    248         return wrapper

~\anaconda3\envs\jupyterdev\lib\site-packages\metpy\calc\thermo.py in lcl(pressure, temperature, dewpoint, max_iters, eps)
    404 
    405     w = mixing_ratio(saturation_vapor_pressure(dewpoint), pressure)
--> 406     lcl_p = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature),
    407                            xtol=eps, maxiter=max_iters)
    408 

~\anaconda3\envs\jupyterdev\lib\site-packages\scipy\optimize\minpack.py in fixed_point(func, x0, args, xtol, maxiter, method)
    935     use_accel = {'del2': True, 'iteration': False}[method]
    936     x0 = _asarray_validated(x0, as_inexact=True)
--> 937     return _fixed_point_helper(func, x0, args, xtol, maxiter, use_accel)

~\anaconda3\envs\jupyterdev\lib\site-packages\scipy\optimize\minpack.py in _fixed_point_helper(func, x0, args, xtol, maxiter, use_accel)
    884         else:
    885             p = p1
--> 886         relerr = _lazywhere(p0 != 0, (p, p0), f=_relerr, fillvalue=p)
    887         if np.all(np.abs(relerr) < xtol):
    888             return p

~\anaconda3\envs\jupyterdev\lib\site-packages\scipy\_lib\_util.py in _lazywhere(cond, arrays, f, fillvalue, f2)
     48     tcode = np.mintypecode([a.dtype.char for a in arrays])
     49     out = np.full(np.shape(arrays[0]), fill_value=fillvalue, dtype=tcode)
---> 50     np.place(out, cond, f(*temp))
     51     if f2 is not None:
     52         temp = tuple(np.extract(~cond, arr) for arr in arrays)

<__array_function__ internals> in place(*args, **kwargs)

~\anaconda3\envs\jupyterdev\lib\site-packages\numpy\lib\function_base.py in place(arr, mask, vals)
   1746                         "not {name}".format(name=type(arr).__name__))
   1747 
-> 1748     return _insert(arr, mask, vals)
   1749 
   1750 

ValueError: place: mask and data must be the same size

yet printing the shape of pressures, temps, and dews gives

(19,)
(19,)
(19,)

so I'm not sure where the mismatch is.

I also printed out el_pres and this_lcl as per @DanielAdriaansen's suggestion. For this sounding, I get:

[] hectopascal
(<Quantity(591.126735, 'hectopascal')>, <Quantity(268.39331, 'kelvin')>)

so perhaps the lack of EL is the problem.

Let me know if any other information would be helpful!