Unit/Pint Troubles

[jsillin]

Hello!

I have been trying to plot frontogenesis using GFS data pulled from NOMADS:

ds = xr.open_dataset(url)
init_hr = dt.datetime(int(year),int(month),int(day),int(init_hour))
times = ds['tmp2m'].metpy.time
init_time = ds['time'][0]

lats = np.arange(15,70,0.25)
lons = np.arange(220,310,0.25)

forecast_hour = times[0].values

data = ds.metpy.parse_cf()
data = data.isel(time=6)
#Rename variables to useful things
data = data.rename({
    'absvprs':'avort',
    'hgtprs':'gph',
    'rhprs':'rh',
    'tmpprs':'temp',
    'ugrdprs':'u',
    'vgrdprs': 'v',
})

vertical, = data['temp'].metpy.coordinates('vertical')
time = data['temp'].metpy.time
zH5_crs = data['temp'].metpy.cartopy_crs

### Pull dataarrays for each parameter of interest
t7 = data['temp'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['u'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()*1.94384449
v7 = data['v'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()*1.94384449
h7 = data['gph'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()

x, y = t7.metpy.coordinates('x', 'y')
dx, dy = mpcalc.grid_deltas_from_dataarray(t7)
lat, lon = xr.broadcast(y, x)
wind_slice = slice(5,-5,5)

t7 = t7 *units.K
u7 = u7*units.knots
v7 = v7*units.knots
h7 = h7*units.m

h7_fgen = mpcalc.frontogenesis(t7,u7,v7,dx,dy)

despite explicitly tagging each array with units, I am still getting the following error:
ValueError: frontogenesisgiven arguments with incorrect units:thtarequires "[temperature]" but given "dimensionless",urequires "[speed]" but given "dimensionless",v requires "[speed]" but given "dimensionless".

I am currently running:
numpy 1.18.5
xarray 0.15.1
pint 0.14
metpy 0.12.1

I tried printing out one of these arrays to see if the unit attribution actually worked, and it seemed to look good:
<xarray.DataArray 'v' (lat: 221, lon: 361)> <Quantity([[ -1.3879315 -1.4656852 -1.7961389 ... 5.065632 5.0850706 4.5213556] [ -0.8242166 -1.3101777 -1.6600697 ... 5.2405777 6.912284 5.7071004] [ -2.1265924 -1.8350157 -1.6989466 ... 3.977079 6.173623 6.484638 ] ... [-12.623352 -13.031559 -13.420328 ... 28.625025 29.363686 30.005156 ] [-13.7119055 -14.158989 -14.528319 ... 29.44144 28.411203 29.110987 ] [-14.858773 -15.26698 -15.616873 ... 28.041874 28.08075 27.167143 ]], 'knot')>

Any ideas what might be causing the trouble and how I might get metpy to see the units?

[jthielen]

In this most recent example, it looks like there are issues with using xarray DataArrays with Quantity arrays inside them (what is created when multiplying by DataArrays by units). While this a well-supported data structure in MetPy v1.0 and later, it can cause issues with earlier MetPy versions. Instead, for MetPy v0.12.1, we would need units to be on the xarray 'units' attribute. And so, I'd recommend seeing if this adaptation resolves your issue:

ds = xr.open_dataset(url)
init_hr = dt.datetime(int(year),int(month),int(day),int(init_hour))
times = ds['tmp2m'].metpy.time
init_time = ds['time'][0]

lats = np.arange(15,70,0.25)
lons = np.arange(220,310,0.25)

forecast_hour = times[0].values

data = ds.metpy.parse_cf()
data = data.isel(time=6)
#Rename variables to useful things
data = data.rename({
    'absvprs':'avort',
    'hgtprs':'gph',
    'rhprs':'rh',
    'tmpprs':'temp',
    'ugrdprs':'u',
    'vgrdprs': 'v',
})

vertical, = data['temp'].metpy.coordinates('vertical')
time = data['temp'].metpy.time
zH5_crs = data['temp'].metpy.cartopy_crs

### Pull dataarrays for each parameter of interest
t7 = data['temp'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310))
u7 = data['u'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()*1.94384449
v7 = data['v'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()*1.94384449
h7 = data['gph'].sel(lev=700.0,lat=slice(15, 70),lon=slice(220, 310)).squeeze()

x, y = t7.metpy.coordinates('x', 'y')
dx, dy = mpcalc.grid_deltas_from_dataarray(t7)
lat, lon = xr.broadcast(y, x)
wind_slice = slice(5,-5,5)

t7.attrs['units'] = 'K'
u7.attrs['units'] = v7.attrs['units'] = 'knots'
h7.attrs['units'] = 'm'

h7_fgen = mpcalc.frontogenesis(t7,u7,v7,dx,dy)

[jsillin]

yep same message

[jthielen]

Alright, that's good info, thanks. For some reason the change to the units attribute is not taking effect. What does print(t7.attrs['units'], t7.metpy.units) placed right before the frontogenesis call return?

[jsillin]

This appears to give:
K 1 dimensionless

[jthielen]

Thanks for that! Based on that and reexamining the MetPy v0.12.1 source, I think I finally got this traced down to a cache issue in pre v1.0 versions of MetPy (it caches the unit attribute on the accessor in a hidden attribute, and so doesn't respond well to manual changes to that attribute). Not sure the best way to fix it with this pathway, so it may just be best to use the other solution suggested below.

[jthielen]

@dopplershift Would this be worth a patch to the 0.12.x branch?

[DanielAdriaansen]

Since you're using 0.12.1, try this:

h7_fgen = mpcalc.frontogenesis(t7.values,u7.values,v7.values,dx,dy)

I'm not sure 0.12.1 can handle a pint quantity wrapped by a DataArray but I think 1.0 can. By using values, you'll pull out just the pint quantity object from the DataArray object for frontogenesis() to work with.

[jsillin]

Thanks so much for the suggestion! When I tried this, I ended up back with metpy not seeming to see any units:
'''
Traceback (most recent call last):
File "gfs_fgen_test.py", line 117, in
h7_fgen = mpcalc.frontogenesis(t7.values,u7.values,v7.values,dx,dy)
File "C:\Users\Jack\anaconda3\envs\metpy\lib\site-packages\metpy\xarray.py", line 677, in wrapper
return func(*args, **kwargs)
File "C:\Users\Jack\anaconda3\envs\metpy\lib\site-packages\metpy\calc\kinematics.py", line 58, in wrapper
ret = func(*args, **kwargs)
File "C:\Users\Jack\anaconda3\envs\metpy\lib\site-packages\metpy\units.py", line 319, in wrapper
raise ValueError(msg)
ValueError: frontogenesis given arguments with incorrect units: thta requires "[temperature]" but given "none", u requires "[speed]" but given "none", v requires "[speed]" but given "none".
Any variable x can be assigned a unit as follows:
from metpy.units import units
x = x * units.meter / units.second
'''

[jthielen]

.values actually always downcasts to a NumPy array...I think what you'd want here is .data to get the Pint Quantity contained within.

[jsillin]

aha!

t7 = t7 *units.K
u7 = u7*units.knots
v7 = v7*units.knots
h7 = h7*units.m

h7_fgen = mpcalc.frontogenesis(t7.data,u7.data,v7.data,dx,dy)

works! hooray! thanks so much to you both!