Value error calculating CAPE for a layer

[dopplershift]

This example from a StackOverflow question fails:

from datetime import datetime, timedelta   
from siphon.catalog import TDSCatalog
from siphon.ncss import NCSS
import numpy as np
from metpy.units import units
import cartopy.crs as ccrs
from metpy.calc import (dewpoint_from_relative_humidity, mixed_parcel, most_unstable_parcel, parcel_profile, pressure_to_height_std, lcl, height_to_pressure_std, get_layer, cape_cin)

year=2019
month=5
day=23
hour=0

cenlat = 34.91269
cenlon = -98.21048

time_start = datetime(year, month, day, hour, 0) #specified time
hour = time_start.hour
if hour < 10:
    hour = '0'+str(hour)
day = time_start.day
if day < 10:
    day = '0'+str(day)
month = time_start.month
if month < 10:
    month = '0'+str(month)

cat = TDSCatalog('https://www.ncdc.noaa.gov/thredds/catalog/model-rap130-old/'+str(time_start.year)+str(month)+'/'+str(time_start.year)+str(month)+str(day)+'/catalog.html?dataset=rap130-old/'+str(time_start.year)+str(month)+'/'+str(time_start.year)+str(month)+str(day)+'/rap_130_'+str(time_start.year)+str(month)+str(day)+'_'+str(hour)+'00_000.grb2')
latest_ds = list(cat.datasets.values())[0]
print(latest_ds.access_urls)
ncss = NCSS(latest_ds.access_urls['NetcdfSubset'])

query = ncss.query()
query.variables('Pressure_surface').variables('Geopotential_height_isobaric').variables('Geopotential_height_surface').variables('Relative_humidity_isobaric').variables('Temperature_isobaric').variables('Dewpoint_temperature_height_above_ground').variables('Temperature_height_above_ground').variables
query.add_lonlat().lonlat_box(cenlon-2.1, cenlon +2.1, cenlat-2.1, cenlat+2.1)
data1 = ncss.get_data(query)
dlev = data1.variables['Geopotential_height_isobaric'].dimensions[1]
dlat = data1.variables['Geopotential_height_isobaric'].dimensions[2]
dlon = data1.variables['Geopotential_height_isobaric'].dimensions[3]

SFCP = (np.asarray(data1.variables['Pressure_surface'][:])/100.) * units('hPa')
hgt = np.asarray(data1.variables['Geopotential_height_isobaric'][:]) * units('meter')
sfc_hgt = np.asarray(data1.variables['Geopotential_height_surface'][:]) * units('meter')
Temp_up = np.asarray(data1.variables['Temperature_isobaric'][:]) * units('kelvin')
RH_up = np.asarray(data1.variables['Relative_humidity_isobaric'][:])
Td = (np.asarray(data1.variables['Dewpoint_temperature_height_above_ground'][:]) * units('kelvin')).to('degC')
T = np.asarray(data1.variables['Temperature_height_above_ground'][:]) * units('kelvin')

# Get the dimension data
lats_r = data1.variables[dlat][:]
lons_r= data1.variables[dlon][:]
lev = (np.asarray(data1.variables[dlev][:])/100.) * units('hPa')

# Set up our array of latitude and longitude values and transform to the desired projection.
flon = float(cenlon)
flat = float(cenlat)
crs = ccrs.PlateCarree()
crlons, crlats = np.meshgrid(lons_r[:]*1000, lats_r[:]*1000)
trlatlons = crs.transform_points(ccrs.LambertConformal(central_longitude=265, central_latitude=25, standard_parallels=(25.,25.)),crlons,crlats)
trlons = trlatlons[:,:,0]
trlats = trlatlons[:,:,1]
dlon = np.abs(trlons - cenlon)
dlat = np.abs(trlats - cenlat)
ilon = np.where(dlon == np.min(dlon)) #position in the dlon array with minimal difference between gridpoint lon and input lon
ilat = np.where(dlat == np.min(dlat)) #position in the dlat array with minimal difference between gridpoint lat and input lat

Tdc_up = dewpoint_from_relative_humidity(Temp_up[0,:,ilat[0][0], ilon[1][0]],RH_up[0,:,ilat[0][0], ilon[1][0]]/100)

p_sounding = np.sort(np.append(lev, SFCP[0,ilat[0][0], ilon[1][0]]))
ind = np.where(p_sounding >= SFCP[0,ilat[0][0], ilon[1][0]])[0][0]
hgt_sounding = np.insert(hgt[0,:,ilat[0][0], ilon[1][0]].magnitude, ind, sfc_hgt[0,ilat[0][0], ilon[1][0]].magnitude) * hgt.units
T_sounding = (np.insert(Temp_up[0,:,ilat[0][0], ilon[1][0]].magnitude, ind, T[0,0,ilat[0][0], ilon[1][0]].magnitude) * T.units).to(Tdc_up.units)
Td_sounding = np.insert(Tdc_up.magnitude, ind, Td[0,0,ilat[0][0], ilon[1][0]].magnitude) * Tdc_up.units
p_skewt = p_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
hgt_skewt = hgt_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
T_skewt = T_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]
Td_skewt = Td_sounding[p_sounding <= SFCP[0,ilat[0][0], ilon[1][0]]]

AGLhgts = hgt_skewt[::-1]-hgt_skewt[-1]

ml_p, ml_T, ml_Td = mixed_parcel(np.flip(p_skewt), np.flip(T_skewt), np.flip(Td_skewt))
ml_profile = parcel_profile(p_skewt[::-1], ml_T, ml_Td)
ml_profile = (ml_profile - 273.15*units('kelvin')).magnitude*units('degC')

mu_p, mu_T, mu_Td, mu_index = most_unstable_parcel(np.flip(p_skewt), np.flip(T_skewt), np.flip(Td_skewt))
mu_profile = parcel_profile(p_skewt[::-1], mu_T, mu_Td)
mu_profile = (mu_profile - 273.15*units('kelvin')).magnitude*units('degC')

#Note: sbpcl_profile will have the exact same values of p_skewt, T_skewt, and Td_skewt in pprof below:
pprof = parcel_profile(p_skewt[::-1], T_skewt[-1], Td_skewt[-1])
pprof = (pprof - 273.15*units('kelvin')).magnitude*units('degC')

mllcl = lcl(ml_p, ml_T, ml_Td)
mllcl_h = pressure_to_height_std(mllcl[0]) - hgt_skewt[-1]

mulcl = lcl(mu_p, mu_T, mu_Td)
mulcl_h = pressure_to_height_std(mulcl[0]) - hgt_skewt[-1]

sblcl = lcl(p_skewt[-1], T_skewt[-1], Td_skewt[-1])
sblcl_h = pressure_to_height_std(sblcl[0]) - hgt_skewt[-1]

mllcl2000 = mllcl_h + 2*units('kilometer')
mulcl2000 = mulcl_h + 2*units('kilometer')
sblcl2000 = sblcl_h + 2*units('kilometer')
mllcl2000_p = height_to_pressure_std(mllcl2000)
mulcl2000_p = height_to_pressure_std(mulcl2000)
sblcl2000_p = height_to_pressure_std(sblcl2000)

ml_LCL_CAPE_layer = get_layer(p_skewt, T_skewt, Td_skewt, ml_profile[::-1], bottom = mllcl[0], depth = mllcl[0] - mllcl2000_p)
mu_LCL_CAPE_layer = get_layer(p_skewt, T_skewt, Td_skewt, mu_profile[::-1], bottom = mulcl[0], depth = mulcl[0] - mulcl2000_p)
sb_LCL_CAPE_layer = get_layer(p_skewt, T_skewt, Td_skewt, pprof[::-1], bottom = sblcl[0], depth = sblcl[0] - sblcl2000_p)

mlLCLCAPE = cape_cin(ml_LCL_CAPE_layer[0], ml_LCL_CAPE_layer[1], ml_LCL_CAPE_layer[2], ml_LCL_CAPE_layer[3])
muLCLCAPE = cape_cin(mu_LCL_CAPE_layer[0], mu_LCL_CAPE_layer[1], mu_LCL_CAPE_layer[2], mu_LCL_CAPE_layer[3])
sbLCLCAPE = cape_cin(sb_LCL_CAPE_layer[0], sb_LCL_CAPE_layer[1], sb_LCL_CAPE_layer[2], sb_LCL_CAPE_layer[3])

with the following:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/var/folders/r8/bhrhn4m56lgdckj491rm6_w40000gp/T/ipykernel_18027/3613122800.py in <module>
      4 
      5 mlLCLCAPE = cape_cin(ml_LCL_CAPE_layer[0], ml_LCL_CAPE_layer[1], ml_LCL_CAPE_layer[2], ml_LCL_CAPE_layer[3])
----> 6 muLCLCAPE = cape_cin(mu_LCL_CAPE_layer[0], mu_LCL_CAPE_layer[1], mu_LCL_CAPE_layer[2], mu_LCL_CAPE_layer[3])
      7 sbLCLCAPE = cape_cin(sb_LCL_CAPE_layer[0], sb_LCL_CAPE_layer[1], sb_LCL_CAPE_layer[2], sb_LCL_CAPE_layer[3])

~/repos/metpy/src/metpy/xarray.py in wrapper(*args, **kwargs)
   1228 
   1229             # Evaluate inner calculation
-> 1230             result = func(*bound_args.args, **bound_args.kwargs)
   1231 
   1232             # Wrap output based on match and match_unit

~/repos/metpy/src/metpy/units.py in wrapper(*args, **kwargs)
    286         def wrapper(*args, **kwargs):
    287             _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs)
--> 288             return func(*args, **kwargs)
    289 
    290         return wrapper

~/repos/metpy/src/metpy/calc/thermo.py in cape_cin(pressure, temperature, dewpoint, parcel_profile, which_lfc, which_el)
   1856                                                                    dewpoint, parcel_profile)
   1857     # Calculate LFC limit of integration
-> 1858     lfc_pressure, _ = lfc(pressure, temperature, dewpoint,
   1859                           parcel_temperature_profile=parcel_profile, which=which_lfc)
   1860 

~/repos/metpy/src/metpy/xarray.py in wrapper(*args, **kwargs)
   1228 
   1229             # Evaluate inner calculation
-> 1230             result = func(*bound_args.args, **bound_args.kwargs)
   1231 
   1232             # Wrap output based on match and match_unit

~/repos/metpy/src/metpy/units.py in wrapper(*args, **kwargs)
    286         def wrapper(*args, **kwargs):
    287             _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs)
--> 288             return func(*args, **kwargs)
    289 
    290         return wrapper

~/repos/metpy/src/metpy/calc/thermo.py in lfc(pressure, temperature, dewpoint, parcel_temperature_profile, dewpoint_start, which)
    556                                                 temperature[1:], direction='decreasing',
    557                                                 log_x=True)
--> 558             if np.min(el_pressure) > this_lcl[0]:
    559                 x = units.Quantity(np.nan, pressure.units)
    560                 y = units.Quantity(np.nan, temperature.units)

~/miniconda3/envs/py310/lib/python3.10/site-packages/numpy/core/overrides.py in amin(*args, **kwargs)

~/miniconda3/envs/py310/lib/python3.10/site-packages/pint/quantity.py in __array_function__(self, func, types, args, kwargs)
   1724 
   1725     def __array_function__(self, func, types, args, kwargs):
-> 1726         return numpy_wrap("function", func, args, kwargs, types)
   1727 
   1728     _wrapped_numpy_methods = ["flatten", "astype", "item"]

~/miniconda3/envs/py310/lib/python3.10/site-packages/pint/numpy_func.py in numpy_wrap(func_type, func, args, kwargs, types)
    919     if name not in handled or any(is_upcast_type(t) for t in types):
    920         return NotImplemented
--> 921     return handled[name](*args, **kwargs)

~/miniconda3/envs/py310/lib/python3.10/site-packages/pint/numpy_func.py in implementation(*args, **kwargs)
    751         for i, unwrapped_unit_arg in enumerate(unwrapped_unit_args):
    752             bound_args.arguments[valid_unit_arguments[i]] = unwrapped_unit_arg
--> 753         ret = func(*bound_args.args, **bound_args.kwargs)
    754 
    755         # Conditionally wrap output

~/miniconda3/envs/py310/lib/python3.10/site-packages/numpy/core/overrides.py in amin(*args, **kwargs)

~/miniconda3/envs/py310/lib/python3.10/site-packages/numpy/core/fromnumeric.py in amin(a, axis, out, keepdims, initial, where)
   2914     6
   2915     """
-> 2916     return _wrapreduction(a, np.minimum, 'min', axis, None, out,
   2917                           keepdims=keepdims, initial=initial, where=where)
   2918 

~/miniconda3/envs/py310/lib/python3.10/site-packages/numpy/core/fromnumeric.py in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
     84                 return reduction(axis=axis, out=out, **passkwargs)
     85 
---> 86     return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
     87 
     88 

ValueError: zero-size array to reduction operation minimum which has no identity

[dopplershift]

I think the failure has something to do with it taking the layer at the LCL, but that's merely a guess. Regardless, it shouldn't fail.

[Alireza-Amani]

Im giving this a try (as a first timer in the open source contribution).

when cape_cin function in calc/thermo.py is called in the last two lines in the snippet/example:

muLCLCAPE = cape_cin(mu_LCL_CAPE_layer[0], mu_LCL_CAPE_layer[1], mu_LCL_CAPE_layer[2], mu_LCL_CAPE_layer[3])
sbLCLCAPE = cape_cin(sb_LCL_CAPE_layer[0], sb_LCL_CAPE_layer[1], sb_LCL_CAPE_layer[2], sb_LCL_CAPE_layer[3])

in turn and inside it there is call to lfc function (lfc: Calculate the level of free convection (LFC)). Thats where the code breaks, unsuccessful call to lfc:

MetPy/src/metpy/calc/thermo.py

Lines 1858 to 1859 in b9d41bf

	lfc_pressure, _ = lfc(pressure, temperature, dewpoint,
	parcel_temperature_profile=parcel_profile, which=which_lfc)

The problem might be the fact that no return statement is invoked inside lfc. And the reason for nothing being returned is here inside lfc:

MetPy/src/metpy/calc/thermo.py

Lines 550 to 563 in b9d41bf

	# LFC exists. Make sure it is no lower than the LCL
	else:
	idx = x < this_lcl[0]
	# LFC height < LCL height, so set LFC = LCL
	if not any(idx):
	el_pressure, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:],
	temperature[1:], direction='decreasing',
	log_x=True)
	if np.min(el_pressure) > this_lcl[0]:
	x = units.Quantity(np.nan, pressure.units)
	y = units.Quantity(np.nan, temperature.units)
	else:
	x, y = this_lcl
	return x, y

When I print the el_pressure immediately after the call to find_intersections here is what I get: "[] hectopascal".
So when the condition np.min(el_pressure) is being evaluated, thats where it breaks.

So if what I pointed out makes any sense at all, I better be looking inside find_intersections function and see why it returns no pressure value. Or maybe change the if statement so that it can handle that case?

Pardon my bad markdown and probably bad communication here.

[Alireza-Amani]

I am now looking inside the find_intersection function when it is called by the last two lines in the OP snippet.
The reason why it returns empty objects is here:

MetPy/src/metpy/calc/tools.py

Lines 172 to 185 in b9d41bf

	# Check for duplicates
	duplicate_mask = (np.ediff1d(intersect_x, to_end=1) != 0)
	# Make a mask based on the direction of sign change desired
	if direction == 'increasing':
	mask = sign_change > 0
	elif direction == 'decreasing':
	mask = sign_change < 0
	elif direction == 'all':
	return intersect_x[duplicate_mask], intersect_y[duplicate_mask]
	else:
	raise ValueError(f'Unknown option for direction: {direction}')
	return intersect_x[mask & duplicate_mask], intersect_y[mask & duplicate_mask]

, where duplicate_mask is evaluated to be True but the mask variable is False. Therefore, the return in the last line returns an empty object.

The problematic call to find_intersection is when its called with direction == 'decreasing' and:

a = [22.72272666 21.78789749 20.82122583 19.82036577 18.78268829 17.70523254
 16.58466605 16.30972975]
b = b:[23.52751994 21.64722226 19.61214924 17.80605228 15.97829835 14.31277485
 12.74920203 12.42728438]

which renders the:

MetPy/src/metpy/calc/tools.py

Lines 137 to 138 in b9d41bf

	# Determine the sign of the change
	sign_change = np.sign(a[next_idx] - b[next_idx])

positive (next_idx == 1 here). And in turn makes the:

MetPy/src/metpy/calc/tools.py

Lines 178 to 179 in b9d41bf

	elif direction == 'decreasing':
	mask = sign_change < 0

eval to False.

• One thing about the direction argument: I havent read the scientific context/doc much, but here a doesnt seem to be decreasing relative to b. So, maybe find_intersection should not be called with direction == 'decreasing' inside the lfc function to begin with. Just a wild ignorant guess.

---

still not sure if its the right place to look at, or if I am making any sense here anyway.

[dopplershift]

@Alireza-Amani Welcome! Thanks a ton for the detailed analysis, it's very well structured and thought out. It does seem like get_intersection is the problem. I remember bumping up against it the last time I looked (very briefly) at this.

A wild guess would be to change to 'increasing' and see if anything else breaks...

[Alireza-Amani]

@dopplershift, thank you for taking the time to review the comments and for the suggestion.

In short, it doesn't break at any other place when I try that. Here are the last three variables of the OP snippet, printed:

mlLCLCAPE: (<Quantity(99.2639974, 'joule / kilogram')>, <Quantity(0.0, 'joule / kilogram')>)
muLCLCAPE: (<Quantity(0, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)
sbLCLCAPE: (<Quantity(0, 'joule / kilogram')>, <Quantity(0, 'joule / kilogram')>)

Details:

Changing direction to increasing inside lfc function at line 556:

MetPy/src/metpy/calc/thermo.py

Lines 553 to 563 in b9d41bf

	# LFC height < LCL height, so set LFC = LCL
	if not any(idx):
	el_pressure, _ = find_intersections(pressure[1:], parcel_temperature_profile[1:],
	temperature[1:], direction='decreasing',
	log_x=True)
	if np.min(el_pressure) > this_lcl[0]:
	x = units.Quantity(np.nan, pressure.units)
	y = units.Quantity(np.nan, temperature.units)
	else:
	x, y = this_lcl
	return x, y

will result in the return statement at line 563 being executed, i.e. returning nan lfc pressure to cape_cin function:

MetPy/src/metpy/calc/thermo.py

Lines 1858 to 1863 in b9d41bf

	lfc_pressure, _ = lfc(pressure, temperature, dewpoint,
	parcel_temperature_profile=parcel_profile, which=which_lfc)
	# If there is no LFC, no need to proceed.
	if np.isnan(lfc_pressure):
	return units.Quantity(0, 'J/kg'), units.Quantity(0, 'J/kg')

which will return at line 1863.

[dcamron]

Thanks for your help here @Alireza-Amani. Testing that out on my end, this change does break at least some other "do we still have an EL if LFC is lower than LCL" scenarios. We'll be diving into this ourselves for the next release, so if you do any more exploration, be sure to let us know!