metpy.calc.el on entire grid points.

[krishnaap]

I am trying to calculate the equilibrium level on the entire grid points of a grib file.

ds
Out[253]: 
<xarray.Dataset>
Dimensions:        (step: 2, isobaricInhPa: 24, latitude: 501, longitude: 751)
Coordinates:
    time           datetime64[ns] ...
  * step           (step) timedelta64[ns] 00:00:00 03:00:00
  * isobaricInhPa  (isobaricInhPa) float64 1e+03 975.0 950.0 ... 30.0 20.0 10.0
  * latitude       (latitude) float64 -15.0 -14.88 -14.76 ... 44.76 44.88 45.0
  * longitude      (longitude) float64 30.0 30.12 30.24 ... 119.8 119.9 120.0
    valid_time     (step) datetime64[ns] ...
Data variables:
    t              (step, isobaricInhPa, latitude, longitude) float32 ...
    gh             (step, isobaricInhPa, latitude, longitude) float32 ...
    r              (step, isobaricInhPa, latitude, longitude) float32 ...
    u              (step, isobaricInhPa, latitude, longitude) float32 ...
    v              (step, isobaricInhPa, latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             dems
    GRIB_centreDescription:  New Delhi 
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             New Delhi 
    history:                 2022-03-16T15:56 GRIB to CDM+CF via cfgrib-0.9.9...

Where the variables used in el calculations, t and td are in 3D (ignored the step), whereas the pressure 'p' is a dimension, which I want to use in the calculation as a variable (presently, p is in 1 dimension). In this case how to calculate 'el' in entire grid points? Can I use the data array method to calculate it in the entire grid points instead of a loop?

ds=xr.open_dataset(fname,engine='cfgrib')
varlist=list(ds.variables)
print(varlist)
temp = ds.t[0,:,:,:]
rel = ds.r[[0,:,:,:]
hei= ds.gh[0,:,:,:]
p = ds.isobaricInhPa

t = (temp-273.15)*units.degC
td = mpcalc.dewpoint_from_relative_humidity(t, rel)
el = mpcalc.el(p,t,td,which='top')   

Then I am getting an error

ValueError: operands could not be broadcast together with shapes (24,) (24,501,751) (24,)

Traceback (most recent call last):

File "/media/Dataset/EL.py", line 59, in
el = mpcalc.el(p,t,td,which='top')

File "/home/anaconda3/envs/ide/lib/python3.9/site-packages/metpy/xarray.py", line 1216, in wrapper
result = func(*bound_args.args, **bound_args.kwargs)

File "/home/anaconda3/envs/ide/lib/python3.9/site-packages/metpy/units.py", line 246, in wrapper
return func(*args, **kwargs)

File "/home/anaconda3/envs/ide/lib/python3.9/site-packages/metpy/calc/thermo.py", line 673, in el
pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint)

File "/home/anaconda3/envs/ide/lib/python3.9/site-packages/metpy/calc/tools.py", line 1475, in _remove_nans
mask |= np.isnan(v)
ValueError: operands could not be broadcast together with shapes (24,) (24,501,751) (24,)

[DanielAdriaansen]

Pressure has to have the same shape as the other variables in order for this to work, so you have to take the values for the pressure coordinate (1D) and expand them to 3D. You can do this with Xarray using broadcast_like():
https://docs.xarray.dev/en/latest/generated/xarray.DataArray.broadcast_like.html

ds=xr.open_dataset(fname,engine='cfgrib')
varlist=list(ds.variables)
print(varlist)
temp = ds.t[0,:,:,:]
rel = ds.r[[0,:,:,:]
hei= ds.gh[0,:,:,:]
p = ds.isobaricInhPa

# Create 3D pressure to use with el() function
p3d = xr.DataArray(p,dims=['isobaricInhPa'],coords={'isobaricInhPa':p},attrs={'units':'hectopascals'}).broadcast_like(temp)

t = (temp-273.15)*units.degC
td = mpcalc.dewpoint_from_relative_humidity(t, rel)
el = mpcalc.el(**p3d**,t,td,which='top')

[DanielAdriaansen]

Unfortunately, no. Per the MetPy 1.2 documentation for el() here: https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.el.html

Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).

So, you will need to loop across the x-y dimensions and call el() for each grid cell and store the result in a 2D array.

The ability to be able to do this without looping is on the MetPy radar (see here: #1533) but still being worked on. There are some suggested workarounds to speed this up on that issue, and I also created another workaround here: #2259 that may be useful for you.

[DanielAdriaansen]

Also, I apologize for not recognizing this earlier, but you probably don't need to create a 3D pressure since you'll need to loop. You can just use the 1D pressure variable at each x-y point.

[dopplershift]

Solving problems like this is a high priority item for us, but won't be solved in the short term.

[krishnaap]

Also, I apologize for not recognizing this earlier, but you probably don't need to create a 3D pressure since you'll need to loop. You can just use the 1D pressure variable at each x-y point.

Looping is working fine, but I don't think it's a good idea if there are too many grid points just like in my case. Is there anyways to use the function just like xarray, that is like calculating using data arrays or data frames to speed up the process?

dt=xr.open_dataset(fname,engine='cfgrib')
varlist=list(dt.variables)
print(varlist)
ds=dt.sel(latitude=slice(0,40),longitude=slice(60,100)) 

# p = ds.isobaricInhPa                                 #dummy variables to copy attrs
temp = ds.t[0,:,:,:]
p3d = xr.DataArray(p,dims=['isobaricInhPa'],coords={'isobaricInhPa':p},attrs={'units':'hPa'}).broadcast_like(temp)

lat = ds.latitude
lon = ds.longitude
for i in range(0,len(lat)):
    for j in range(0,len(lon)):
        temp = ds.t[0,:,i,j]
        rel = ds.r[0,:,i,j]
        hei= ds.gh[0,:,i,j]
        p = ds.isobaricInhPa
        t = (temp-273.15)*units.degC
        td = mpcalc.dewpoint_from_relative_humidity(t, rel)  
        pres =p3d[:,i,j]
        # pres = p3d[:,i,:]
        # print(pres)
        # print(pres,t,td)
        el = mpcalc.el(pres,t,td,which='top') 
        print(el[0])

I think this is the function from thermo.py (anaconda3/envs/ide/lib/python3.9/site-packages/metpy/calc/thermo.py)

def el(pressure, temperature, dewpoint, parcel_temperature_profile=None, which='top'):
r"""Calculate the equilibrium level.

This works by finding the last intersection of the ideal parcel path and
the measured environmental temperature. If there is one or fewer intersections, there is
no equilibrium level.

Parameters
----------
pressure : `pint.Quantity`
    Atmospheric pressure profile

temperature : `pint.Quantity`
    Temperature at the levels given by `pressure`

dewpoint : `pint.Quantity`
    Dewpoint at the levels given by `pressure`

parcel_temperature_profile: `pint.Quantity`, optional
    The parcel's temperature profile from which to calculate the EL. Defaults to the
    surface parcel profile.

which: str, optional
    Pick which LFC to return. Options are 'top', 'bottom', 'wide', 'most_cape', and 'all'.
    'top' returns the lowest-pressure EL, default.
    'bottom' returns the highest-pressure EL.
    'wide' returns the EL whose corresponding LFC is farthest away.
    'most_cape' returns the EL that results in the most CAPE in the profile.

Returns
-------
`pint.Quantity`
    EL pressure, or array of same if which='all'

`pint.Quantity`
    EL temperature, or array of same if which='all'

See Also
--------
parcel_profile

Notes
-----
Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
Since this function returns scalar values when given a profile, this will return Pint
Quantities even when given xarray DataArray profiles.

.. versionchanged:: 1.0
   Renamed ``dewpt`` parameter to ``dewpoint``

"""
pressure, temperature, dewpoint = _remove_nans(pressure, temperature, dewpoint)
# Default to surface parcel if no profile or starting pressure level is given
if parcel_temperature_profile is None:
    new_stuff = parcel_profile_with_lcl(pressure, temperature, dewpoint)
    pressure, temperature, dewpoint, parcel_temperature_profile = new_stuff
    parcel_temperature_profile = parcel_temperature_profile.to(temperature.units)

# If the top of the sounding parcel is warmer than the environment, there is no EL
if parcel_temperature_profile[-1] > temperature[-1]:
    return (units.Quantity(np.nan, pressure.units),
            units.Quantity(np.nan, temperature.units))

# Interpolate in log space to find the appropriate pressure - units have to be stripped
# and reassigned to allow np.log() to function properly.
x, y = find_intersections(pressure[1:], parcel_temperature_profile[1:], temperature[1:],
                          direction='decreasing', log_x=True)
lcl_p, _ = lcl(pressure[0], temperature[0], dewpoint[0])
idx = x < lcl_p
if len(x) > 0 and x[-1] < lcl_p:
    return _multiple_el_lfc_options(x, y, idx, which, pressure,
                                    parcel_temperature_profile, temperature, dewpoint,
                                    intersect_type='EL')
else:
    return (units.Quantity(np.nan, pressure.units),
            units.Quantity(np.nan, temperature.units))

Is there anyway to change it into a xarray array and calculate in a faster method using it?

[dopplershift]

Looping is working fine, but I don't think it's a good idea if there are too many grid points just like in my case. Is there anyways to use the function just like xarray, that is like calculating using data arrays or data frames to speed up the process?

Speed is definitely the problem with looping. What you're asking for here is the work we still need to do to make MetPy faster. As far as I know, it's not trivial, but I'd be happy to be proven wrong. I'd love to have you or anyone else submit a PR to make it work.