search

Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.25k stars 415 forks source link

lcl fails with NaNs #2191

Open sgdecker opened 2 years ago

sgdecker commented 2 years ago

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)
dopplershift commented 2 years ago

If you look through the traceback, you can see the error is coming from so.fixed_point, which is the function we're using from scipy.optimize to solve for the LCL. So not a bug, but would be nice to have support for NaNs. Our options in this case are: