lcl fails with NaNs #2191
Description
What went wrong?
Not sure if this is a bug or a feature request. My goal was to compute the condensation pressure on the 300-K surface, which is partially underground. I was expecting the lcl function to return NaNs at the gridpoints where the 300-K surface is underground, but instead, it crashed:
/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/interpolate/one_dimension.py:147: UserWarning: Interpolation point out of data bounds encountered
warnings.warn('Interpolation point out of data bounds encountered')
Traceback (most recent call last):
File "/home/decker/metpy/isentropic.py", line 55, in <module>
p, T = mpcalc.lcl(nam_isen['pressure'], nam_isen['temperature'], dwpk)
File "/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/xarray.py", line 1216, in wrapper
result = func(*bound_args.args, **bound_args.kwargs)
File "/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/units.py", line 246, in wrapper
return func(*args, **kwargs)
File "/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/calc/thermo.py", line 420, in lcl
lcl_p = so.fixed_point(_lcl_iter, pressure.m, args=(pressure.m, w, temperature),
File "/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/scipy/optimize/minpack.py", line 941, in fixed_point
x0 = _asarray_validated(x0, as_inexact=True)
File "/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/scipy/_lib/_util.py", line 293, in _asarray_validated
a = toarray(a)
File "/home/decker/miniconda3/envs/met433/lib/python3.9/site-packages/numpy/lib/function_base.py", line 488, in asarray_chkfinite
raise ValueError(
ValueError: array must not contain infs or NaNs
I think I can manually write a nested loop to calculate condensation pressure gridpoint by gridpoint, checking for NaNs before calling lcl(), but I don't think this workaround should be necessary.
Operating System
Linux
Version
1.1.0
Python Version
3.9.7
Code to Reproduce
from datetime import datetime, timedelta
import xarray as xr
from xarray.backends import NetCDF4DataStore
from metpy.units import units
import metpy.calc as mpcalc
from siphon.catalog import TDSCatalog
def get_nam212(init_time, valid_time):
ymd = init_time.strftime('%Y%m%d')
hr = init_time.strftime('%H')
filename = f'{ymd}_{hr}00.grib2'
ds_name = 'NAM_CONUS_40km_conduit_' + filename
cat_name = ('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NAM/'
'CONUS_40km/conduit/' + ds_name + '/catalog.xml')
cat = TDSCatalog(cat_name)
ds = cat.datasets[ds_name]
ncss = ds.subset()
query = ncss.query()
query.time(valid_time).variables('all')
nc = ncss.get_data(query)
data = xr.open_dataset(NetCDF4DataStore(nc)).metpy.parse_cf()
return data
init_time = datetime(2021, 11, 7, 12)
plot_time = init_time + timedelta(hours=12)
nam = get_nam212(init_time, plot_time)
hght = nam['Geopotential_height_isobaric'].squeeze().metpy.quantify()
tmpk = nam['Temperature_isobaric'].squeeze().metpy.quantify()
urel = nam['u-component_of_wind_isobaric'].squeeze().metpy.quantify()
vrel = nam['v-component_of_wind_isobaric'].squeeze().metpy.quantify()
spfh = nam['Specific_humidity_isobaric'].squeeze().metpy.quantify()
hght = mpcalc.smooth_gaussian(hght, 16).rename('hght')
tmpk = mpcalc.smooth_gaussian(tmpk, 16).rename('tmpk')
urel = mpcalc.smooth_gaussian(urel, 16).rename('urel')
vrel = mpcalc.smooth_gaussian(vrel, 16).rename('vrel')
spfh = mpcalc.smooth_gaussian(spfh, 16).rename('spfh')
isen_lev = [300] * units('K')
nam_isen = mpcalc.isentropic_interpolation_as_dataset(isen_lev, tmpk, hght,
urel, vrel, spfh)
dwpk = mpcalc.dewpoint_from_specific_humidity(nam_isen['pressure'],
nam_isen['temperature'],
nam_isen['spfh'])
p, T = mpcalc.lcl(nam_isen['pressure'], nam_isen['temperature'], dwpk)