Math between metpy return dataarrays and xarray dataarrays

[Marston]

Doing some math between some data returned from Metpy and other purely xarray dataarrays. I thought this would be a simple equation of the arrays but I got his error involving the units, which I believed is caused by the Metpy generated/returned dataarrays which includes Pint units:

    # get pressure in Pa
    pa = grb.data_to_da(pkey)
    pa.values = pa.values * 1e-2 # change to hPa
    pa.attrs['units'] = 'hPa'
    # get mixing ratio
    mr = grb.data_to_da(mrkey)
    mr.values = mr.values * 1e3 #Changed units to g/kg according to metpy example
    mr.attrs['units'] = 'g/kg'  
    print('calculate the water vapor pressure...') 
    pw = vapor_pressure(pa, mr)
    print('Extracting the temperature...')
    t = grb.data_to_da(tkey)
    print('Calculate the refractivity...')
    ref = ( b1 * ((pa.data-pw.data)/t.data) ) + \
             ( b2 * ((pw.data) /t) ) + \
             ( b3 * (pw.data/(t.data*t.data)) )

Error:

Cannot convert from 'dimensionless' to 'hectopascal'
  File "/home/Projects/gl17tozarr/myfile.py", line 52, in <module>
    ref = ( b1 * ((pa.data-pw.data)/t.data) ) + ( b2 * ((pw.data) /t) ) + ( b3 * (pw.data/(t.data*t.data)) )
pint.errors.DimensionalityError: Cannot convert from 'dimensionless' to 'hectopascal'

How do I remove the effects of pint so that I can do the calculations in the form:
ref = ( b1 * ((pa-pw)/t) ) + ( b2 * ((pw) /t) ) + ( b3 * (pw/(t*t)) )?

This solution works but is the effects of Pint supposed follow the returned dataarray?
ref = ( b1 * (( pa.values - pw.values )/t.values ) ) + ( b2 * (( pw.values ) / t ) ) + ( b3 * ( pw.values/( t.values*t.values )) )

[dcamron]

There's nothing in Pint keeping you from doing appropriate math with units. Your issue here is that your pa DataArray has a units attribute but is not using Pint under the hood, meanwhile your pw DataArray is. So pa - pw (eg 3 - 100 hPa) is then an invalid operation as far as Pint is concerned.

The xarray with MetPy Tutorial in the docs has a units section for working with Pint-in-xarray. You can DataArray.metpy.quantify() or DataArray.metpy.dequantify() across the board to standardize working with or without Pint, respectively.

[dopplershift]

So I just wanted to follow up and show how I would do this so that you're not fighting the unit framework, but using it to make sure you've got everything right and in proper unit agreement:

from metpy.units import units

# Get pressure and mixing ratio and have as Pint underneath--relies on built-in unit metadata
pa = grb.data_to_da(pkey).metpy.quantify()
mr = grb.data_to_da(mrkey).metpy.quantify()
pw = vapor_pressure(pa, mr)

t = grb.data_to_da(tkey).metpy.quantify()

# Fill in your values here, but notice how I use the units
b1 = 1 * units('K / hPa')
b2 = 2 * units('K / hPa')
b3 = 3 * units('K^2 / hPa')

ref = b1 * ((pa - pw) / t) + b2 * (pw / t) + b3 * (pw / t**2)

By using the built-in metadata from the data file, and specifying the constants with unit information, you avoid any need to manually deal with units and it makes the code robust about underlying data changes. If your data doesn't actually have built-in metadata, I can modify the above code to deal with that.

[Marston]

This is an interesting solution. I am not sure how to add the method .metpy.quantify() to my custom function grb.data_to_da. Can you explain a bit how that is intended to work? From how I understand it, I would need to add this method to my underlying function. I simply tried to chain it, but it didn't work. All the dataArrays have units assigned as attributes.

[Marston]

I noticed that .metpy.quantify() removes the pa.attrs['units'] from the dataArray.

[dopplershift]

Right, .metpy.quantify() removes the units attribute in lieu of having that information tracked by Pint--instead of there being multiple locations of information that need synchronization, we just rely on Pint having that information. You can still get it with .metpy.units.

You shouldn't have to modify your function so long as it's returning an Xarray DataArray instance. MetPy registers what xarray calls an "accessor" to add the custom namespace (.metpy) of attributes and methods, like (.quantify()). So as long as you import metpy (from metpy.units import units counts) before you try to do .metpy.quantify(), it should work.

[dcamron]

And if you decide to .metpy.dequantify() to return the underlying values to a NumPy array (eg to save the file to NetCDF), MetPy will reattach a new DataArray.attrs['units'] with the correct units coming from Pint.