MetPy
MetPy copied to clipboard
Unidata
Value error calculating CAPE for a layer
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
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.
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:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/thermo.py#L1858-L1859
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:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/thermo.py#L550-L563
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.
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:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/tools.py#L172-L185
, 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:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/tools.py#L137-L138
positive (next_idx == 1 here). And in turn makes the:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/tools.py#L178-L179
eval to False.
- One thing about the
directionargument: I havent read the scientific context/doc much, but hereadoesnt seem to bedecreasingrelative tob. So, maybefind_intersectionshould not be called withdirection == 'decreasing'inside thelfcfunction 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.
@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...
@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:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/thermo.py#L553-L563
will result in the return statement at line 563 being executed, i.e. returning nan lfc pressure to cape_cin function:
https://github.com/Unidata/MetPy/blob/b9d41bfe28044ead56ba56f5dd510203548a3d7a/src/metpy/calc/thermo.py#L1858-L1863
which will return at line 1863.
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!
Another stackoverflow question with no EL on either end hit this. We'll prioritize this for the November milestone.