Issue with mpcalc.equivalent_potential_temperature

[SPiltz01]

I noticed that mpcalc.potential_temperature(pmsl, T_Sfc) and mpcalc.equivalent_potential_temperature(pmsl, T_Sfc, Td_sfc) return two different array types. It looks like mpcalc.potential_temperature is returning a 2D array with units attached, while mpcalc.equivalent_potential_temperature is returning a 2D array of individual unit-ed quantities.

#Sfc theta
sfc_Theta = mpcalc.potential_temperature(pmsl, T_Sfc)

print('Theta:')
print(sfc_Theta)

Theta:
[[294.52630615234375 294.7657165527344 297.162841796875 ... 295.77618408203125 294.7249450683594 294.43292236328125] [294.4823303222656 294.2461853027344 298.4425048828125 ... 295.5538635253906 294.7125549316406 294.36090087890625] [298.4830017089844 299.8385009765625 299.524169921875 ... 295.846923828125 295.6717834472656 294.9158020019531] ... [297.68463134765625 297.5202331542969 297.4456787109375 ... 296.65972900390625 296.6708984375 296.67218017578125] [297.57684326171875 297.4230651855469 297.5305480957031 ... 296.80523681640625 296.8392333984375 296.8185729980469] [297.8320007324219 297.47601318359375 297.5194396972656 ... 297.0497741699219 297.0862121582031 297.0396423339844]] kelvin

#Sfc thetae
sfc_Thetae = mpcalc.equivalent_potential_temperature(pmsl, T_Sfc, Td_sfc)

print('Theta-E:')
print(sfc_Thetae)

Theta-E:
[[<Quantity(303.69500732421875, 'kelvin')>
<Quantity(303.3419494628906, 'kelvin')>
<Quantity(305.7081604003906, 'kelvin')> ...
<Quantity(329.5704040527344, 'kelvin')>
<Quantity(329.3699645996094, 'kelvin')>
<Quantity(329.0014953613281, 'kelvin')>]
[<Quantity(303.3021545410156, 'kelvin')>
<Quantity(302.5059814453125, 'kelvin')>
<Quantity(306.9647521972656, 'kelvin')> ...
<Quantity(331.018798828125, 'kelvin')>
<Quantity(329.92242431640625, 'kelvin')>
<Quantity(329.8907470703125, 'kelvin')>]
[<Quantity(307.45062255859375, 'kelvin')>
<Quantity(308.7452697753906, 'kelvin')>
<Quantity(308.14593505859375, 'kelvin')> ...
<Quantity(330.7307434082031, 'kelvin')>
<Quantity(331.5306701660156, 'kelvin')>
<Quantity(331.3540954589844, 'kelvin')>]
...
[<Quantity(332.9311218261719, 'kelvin')>
<Quantity(333.3431701660156, 'kelvin')>
<Quantity(333.7388000488281, 'kelvin')> ...
<Quantity(345.5285949707031, 'kelvin')>
<Quantity(345.1958923339844, 'kelvin')>
<Quantity(344.6487731933594, 'kelvin')>]
[<Quantity(333.20770263671875, 'kelvin')>
<Quantity(333.5695495605469, 'kelvin')>
<Quantity(334.2283935546875, 'kelvin')> ...
<Quantity(346.56817626953125, 'kelvin')>
<Quantity(346.26495361328125, 'kelvin')>
<Quantity(345.5577087402344, 'kelvin')>]
[<Quantity(332.3034973144531, 'kelvin')>
<Quantity(333.88092041015625, 'kelvin')>
<Quantity(334.7017517089844, 'kelvin')> ...
<Quantity(347.37225341796875, 'kelvin')>
<Quantity(347.0073547363281, 'kelvin')>
<Quantity(346.72216796875, 'kelvin')>]]

Thanks!

[jthielen]

Thank you for the report! The second case (array of scalar quantities) is indeed unintended behavior.

To aid in debugging, would you be able to let us know what the type of your inputs pmsl, T_Sfc, Td_sfc are (e.g., xarray.DataArray, pint.Quantity, numpy.ndarray)? Also, because much of the unit compatibility features have been in flux recently, could you also let us know what versions of MetPy, NumPy, Pint, and xarray you are using?

[SPiltz01]

MetPy - 0.12.0
NumPy - 1.18.1
Pint - 0.11
xarray - 0.15.0

Input was from the UCAR Thredds server...

ncss = NCSS('http://thredds.ucar.edu'
'/thredds/ncss/grib/NCEP/GFS/Global_0p25deg/GFS_Global_0p25deg_'+m_date+'_'+m_hour+'00.grib2')

Query for required variables

data = ncss.query().time(modtime)

data.variables('MSLP_Eta_model_reduction_msl',
'u-component_of_wind_height_above_ground',
'v-component_of_wind_height_above_ground',
'Temperature_height_above_ground',
'Dewpoint_temperature_height_above_ground',
'v-component_of_wind_isobaric',
).add_lonlat()

Set the lat/lon box for the data to pull in.

data.lonlat_box(wlon, elon, slat, nlat)

############################

Actually getting the data

data = ncss.get_data(data)
#print (data.variables)

Assign variable names to collected data

dtime = data.variables['Temperature_height_above_ground'].dimensions[0]
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
times = data.variables[dtime]
vtimes = num2date(times[:], times.units)

T_Sfc = (units(data.variables['Temperature_height_above_ground'].units) *
data.variables['Temperature_height_above_ground'][0, 0, :])

Td_sfc = (units(data.variables['Dewpoint_temperature_height_above_ground'].units) *
data.variables['Dewpoint_temperature_height_above_ground'][0, 0, :])

uwnd_10m = units('m/s') * data.variables['u-component_of_wind_height_above_ground'][0, 0, :]
vwnd_10m = units('m/s') * data.variables['v-component_of_wind_height_above_ground'][0, 0, :]

pmsl = (units(data.variables['MSLP_Eta_model_reduction_msl'].units) *
data.variables['MSLP_Eta_model_reduction_msl'][0, :])

[jthielen]

@SPiltz01 Thank you for that additional information. Based on this, it looks like a masked array issue within equivalent_potential_temperature:

https://github.com/Unidata/MetPy/blob/v0.12.1/src/metpy/calc/thermo.py#L992-L1002

As your own code correctly accounts for, multiplying masked arrays by units on the right does not work for masked arrays as of recent versions of Pint. If you wanted a workaround for this bug for now, you can either wrap your NetCDF4 variables with np.array() before adding units to unmask your arrays, or utilize xarray instead of the NetCDF4 library directly (which will also use fill values instead of masked arrays). I had coincidentally helped someone with a very similar issue a couple days ago, so the example notebooks I prepared may be helpful to you too: array-wrapping fix, xarray fix.

That being said, this is still a bug in MetPy (that was missed in #983 ... it will be nice when #1214 is in place). I believe changing the concluding line to instead be

return units.Quantity(th_e, 'kelvin')

will resolve this. saturation_equivalent_potential_temperature will need to be updated as well.

[SPiltz01]

Thanks Jon. I appreciate the help. - Steve