Calculating temperature advection with ERA5 data and metpy 0.12

[Z-Richard]

For quick reference I'm using metpy 0.12.2, and I am trying to compute temperature advection:
temp_advection = mpcalc.advection(temperature, [u, v], (dx, dy), dim_order='yx')
Here, temperature, u, v have a dimension of 181 * 341, and are DataArray; dx 181 * 340, dy, 180 * 341, both of which are computed using mpcalc.lat_lon_grid_deltas, and are Pint.Quantity. However, an error was thrown:

TypeError                                 Traceback (most recent call last)
<ipython-input-134-ce0f612590b3> in <module>
     38     return dy.shape
     39 
---> 40 slp_plots(t)

<ipython-input-134-ce0f612590b3> in slp_plots(time)
     32     dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
     33 
---> 34     temp_advection = mpcalc.advection(temperature, [u, v], (dx, dy), dim_order='yx')
     35 
     36     plt.show()

/opt/anaconda/envs/adm/lib/python3.7/site-packages/metpy/xarray.py in wrapper(*args, **kwargs)
    675         kwargs = {name: (v.metpy.unit_array if isinstance(v, xr.DataArray) else v)
    676                   for name, v in kwargs.items()}
--> 677         return func(*args, **kwargs)
    678     return wrapper
    679 

/opt/anaconda/envs/adm/lib/python3.7/site-packages/metpy/calc/kinematics.py in wrapper(dim_order, *args, **kwargs)
     56                 kwargs[k] = _check_and_flip(v)
     57 
---> 58         ret = func(*args, **kwargs)
     59 
     60         # If we flipped on the way in, need to flip on the way out so that output array(s)

/opt/anaconda/envs/adm/lib/python3.7/site-packages/metpy/calc/kinematics.py in advection(scalar, wind, deltas)
    315     """
    316     # This allows passing in a list of wind components or an array.
--> 317     wind = _stack(wind)
    318 
    319     # If we have more than one component, we need to reverse the order along the first

/opt/anaconda/envs/adm/lib/python3.7/site-packages/metpy/calc/kinematics.py in _stack(arrs)
     20 
     21 def _stack(arrs):
---> 22     return concatenate([a[np.newaxis] if iterable(a) else a for a in arrs], axis=0)
     23 
     24 

/opt/anaconda/envs/adm/lib/python3.7/site-packages/metpy/calc/kinematics.py in <listcomp>(.0)
     20 
     21 def _stack(arrs):
---> 22     return concatenate([a[np.newaxis] if iterable(a) else a for a in arrs], axis=0)
     23 
     24 

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/dataarray.py in __getitem__(self, key)
    640         else:
    641             # xarray-style array indexing
--> 642             return self.isel(indexers=self._item_key_to_dict(key))
    643 
    644     def __setitem__(self, key: Any, value: Any) -> None:

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/dataarray.py in isel(self, indexers, drop, missing_dims, **indexers_kwargs)
   1038         if any(is_fancy_indexer(idx) for idx in indexers.values()):
   1039             ds = self._to_temp_dataset()._isel_fancy(
-> 1040                 indexers, drop=drop, missing_dims=missing_dims
   1041             )
   1042             return self._from_temp_dataset(ds)

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/dataset.py in _isel_fancy(self, indexers, drop, missing_dims)
   2012             if name in self.indexes:
   2013                 new_var, new_index = isel_variable_and_index(
-> 2014                     name, var, self.indexes[name], var_indexers
   2015                 )
   2016                 if new_index is not None:

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/indexes.py in isel_variable_and_index(name, variable, index, indexers)
    104         )
    105 
--> 106     new_variable = variable.isel(indexers)
    107 
    108     if new_variable.dims != (name,):

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/variable.py in isel(self, indexers, missing_dims, **indexers_kwargs)
   1116 
   1117         key = tuple(indexers.get(dim, slice(None)) for dim in self.dims)
-> 1118         return self[key]
   1119 
   1120     def squeeze(self, dim=None):

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/variable.py in __getitem__(self, key)
    764         array `x.values` directly.
    765         """
--> 766         dims, indexer, new_order = self._broadcast_indexes(key)
    767         data = as_indexable(self._data)[indexer]
    768         if new_order:

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/variable.py in _broadcast_indexes(self, key)
    610         # key can be mapped as an OuterIndexer.
    611         if all(not isinstance(k, Variable) for k in key):
--> 612             return self._broadcast_indexes_outer(key)
    613 
    614         # If all key is 1-dimensional and there are no duplicate labels,

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/variable.py in _broadcast_indexes_outer(self, key)
    686             new_key.append(k)
    687 
--> 688         return dims, OuterIndexer(tuple(new_key)), None
    689 
    690     def _nonzero(self):

/opt/anaconda/envs/adm/lib/python3.7/site-packages/xarray/core/indexing.py in __init__(self, key)
    408                 if not np.issubdtype(k.dtype, np.integer):
    409                     raise TypeError(
--> 410                         f"invalid indexer array, does not have integer dtype: {k!r}"
    411                     )
    412                 if k.ndim != 1:

TypeError: invalid indexer array, does not have integer dtype: array(None, dtype=object)

I am not sure what's going on here after looking through the 0.12 documentation. Any insights are appreciated!

[dcamron]

If you need to work with MetPy<1.0, you'll have to

• explicitly pull out U, V numpy arrays from your data arrays: u = ds['u'].values; use these as-is or manually create Pint quantities with units, or
• use one of the MetPy xarray accessors that was available in 0.12: u = ds['u'].metpy.unit_array to create Pint quantities with units quickly.

You'll have to specifically do this for U, V for this to work. As the return will be a Pint quantity with units, I'd suggest a consistent application of above to Temperature as well, so the result is physically meaningful. Either way, I'd strongly suggest upgrading to a newer MetPy if at all possible. Thanks!

[Z-Richard]

Thanks for the reply! I agree that I should upgrade to a new MetPy (and a later version of Python).