Errors when using metpy to calculate specific humidity

Question

I would like to use metpy to calculate the near-surface (i.e. 2m) specific humidity using ERA5 hourly reanalysis data. I just installed metpy locally yesterday via pip, so I'm assuming that my code is up-to-date. My problem is that I keep running into the errors below.

Here is my code at present:

# import modules
import numpy as np
import xarray as xr
from metpy.units import units
import metpy.calc as mpcalc

# read data
d2m = xr.open_dataset('netcdf/ERA5_dewpt2m_1992.nc')
sp = xr.open_dataset('netcdf/ERA5_pres_1992.nc')

# assign units (approach 1)
#d2m = d2m*units.kelvin
#sp = sp*units.pascal

# assign units (approach 2)
d2m = units.Quantity(d2m, "kelvin")
sp = units.Quantity(sp, "pascal")

# calculate specific humidity
aqh2m = mpcalc.specific_humidity_from_dewpoint(sp, d2m)

If I ignore the "units" step, then the function of course complains about the lack of units. Note that I've tried two different approaches to handling units above, but neither of them work.

If I try this approach:

# assign units (approach 1)
d2m = d2m*units.kelvin
sp = sp*units.pascal

# calculate specific humidity
aqh2m = mpcalc.specific_humidity_from_dewpoint(sp, d2m)

then I get the following error:

ValueError: This function changed in 1.0--double check that the function is being called properly.
`specific_humidity_from_dewpoint` given arguments with incorrect units: `pressure` requires "[pressure]" but given "none", `dewpoint` requires "[temperature]" but given "none"
Any variable `x` can be assigned a unit as follows:
    from metpy.units import units
    x = units.Quantity(x, "m/s")

However, if I take the unit assignment approach suggested in the error message, i.e.:

# assign units (approach 2)
d2m = units.Quantity(d2m, "kelvin")
sp = units.Quantity(sp, "pascal")

# calculate specific humidity
aqh2m = mpcalc.specific_humidity_from_dewpoint(sp, d2m)

then I simply get a different error message:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-86573b0b8029> in <module>()
----> 1 d2m = units.Quantity(d2m, "kelvin")
      2 sp = units.Quantity(sp, "pascal")
      3 
      4 # calculate specific humidity
      5 aqh2m = mpcalc.specific_humidity_from_dewpoint(sp, d2m)

~/.local/lib/python3.6/site-packages/pint/quantity.py in __new__(cls, value, units)
    201     def __new__(cls, value, units=None):
    202         if is_upcast_type(type(value)):
--> 203             raise TypeError(f"Quantity cannot wrap upcast type {type(value)}")
    204         elif units is None:
    205             if isinstance(value, str):

TypeError: Quantity cannot wrap upcast type xarray.core.dataset.Dataset

What should I do? I don't know where to go from here. Any advice/guidance you can provide would be appreciated.

FYI, here are the versions of my various modules. I'll list them all, since I'm not sure what metpy uses internally for unit support:

python 3.6.3
numpy 1.19.5
xarray 0.16.2
pandas 1.1.5
pint 0.17
pooch 1.4.0
pyproj 3.0.1
matplotlib 2.1.2
traitlets 4.3.2
scipy 1.0.0

Thanks in advance.

score 5 · Answer 1 · answered Jun 24 '21 at 13:06

It appears that units.Quantity cannot process Xarray objects. Instead xarray.DataArray objects have the Xarray.DataArray.metpy.quantify method to turn the data into a metpy pint.Quantity (note that this loads non-dask arrays into memory). This will move the unit from an attribute to the data.

You can go as:

# import modules
import numpy as np
import xarray as xr
from metpy.units import units
import metpy.calc as mpcalc

# read data
d2m = xr.open_dataset('netcdf/ERA5_dewpt2m_1992.nc')
sp = xr.open_dataset('netcdf/ERA5_pres_1992.nc')

# assign units
d2m = d2m.metpy.quantify()
sp = sp.metpy.quantify()

# calculate specific humidity
aqh2m = mpcalc.specific_humidity_from_dewpoint(sp, d2m)

I expect your ERA5 data to be formatted CF-compliant so this should work. If not, maybe the section on "Non CF-compliant Datasets" at the page for Xarray-specific Metpy functionality might help: https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html

score 5 · Answer 2 · answered Jun 24 '21 at 20:03

One of the problems you're having is that the result of xr.open_dataset() is a Dataset, not a DataArray, which is what you really want to be working with for MetPy's calculation functions. Once you have those, MetPy should be able to automatically handle DataArrays with unit metadata like:

import xarray as xr
import metpy.calc as mpcalc

# read data
d2m_ds = xr.open_dataset('netcdf/ERA5_dewpt2m_1992.nc')
sp_ds = xr.open_dataset('netcdf/ERA5_pres_1992.nc')

# Pull out DataArrays and ask MetPy to parse CF metadata
d2m_arr = d2m_ds.metpy.parse_cf('dewpt2m')
sp_arr = sp_ds.metpy.parse_cf('press')

# calculate specific humidity
aqh2m = mpcalc.specific_humidity_from_dewpoint(sp_arr, d2m_arr)

Errors when using metpy to calculate specific humidity

2 Answers2

Linked