Error in absolute_momentum (bug in normal_component?)

[ShunsukeHoshino]

I want to get absolute momentum with xarray.DataArray like:

ds = xr.DataArray(...)  # set data, including 'u' and 'v'
cross = cross_section(data, ...)  # create cross section
m = absolute_momentum(cross['u'], cross['v'])

But DimensionalityError occurs:

> DimensionalityError                       Traceback (most recent call last)
> <ipython-input-269-605840055d4d> in <module>
> ----> 1 m = mp_calc.absolute_momentum(cross['u'], cross['v'])
> 
> ~/myproject-env/lib/python3.8/site-packages/metpy/xarray.py in wrapper(*args, **kwargs)
>    1300                 if not first.metpy.coordinates_identical(other):
>    1301                     raise ValueError('Input DataArray arguments must be on same coordinates.')
> -> 1302         return func(*args, **kwargs)
>    1303     return wrapper
>    1304 
> 
> ~/myproject-env/lib/python3.8/site-packages/metpy/calc/cross_sections.py in absolute_momentum(u, v, index)
>     293     """
>     294     # Get the normal component of the wind
> --> 295     norm_wind = normal_component(u, v, index=index).metpy.convert_units('m/s')
>     296 
>     297     # Get other pieces of calculation (all as ndarrays matching shape of norm_wind)
> 
> ~/myproject-env/lib/python3.8/site-packages/metpy/xarray.py in convert_units(self, units)
>     170         convert_coordinate_units
>     171         """
> --> 172         return self.quantify().copy(data=self.unit_array.to(units))
>     173 
>     174     def convert_coordinate_units(self, coord, units):
> 
> ~/myproject-env/lib/python3.8/site-packages/pint/quantity.py in to(self, other, *contexts, **ctx_kwargs)
>     658         other = to_units_container(other, self._REGISTRY)
>     659 
> --> 660         magnitude = self._convert_magnitude_not_inplace(other, *contexts, **ctx_kwargs)
>     661 
>     662         return self.__class__(magnitude, other)
> 
> ~/myproject-env/lib/python3.8/site-packages/pint/quantity.py in _convert_magnitude_not_inplace(self, other, *contexts, **ctx_kwargs)
>     607                 return self._REGISTRY.convert(self._magnitude, self._units, other)
>     608 
> --> 609         return self._REGISTRY.convert(self._magnitude, self._units, other)
>     610 
>     611     def _convert_magnitude(self, other, *contexts, **ctx_kwargs):
> 
> ~/myproject-env/lib/python3.8/site-packages/pint/registry.py in convert(self, value, src, dst, inplace)
>     945             return value
>     946 
> --> 947         return self._convert(value, src, dst, inplace)
>     948 
>     949     def _convert(self, value, src, dst, inplace=False, check_dimensionality=True):
> 
> ~/myproject-env/lib/python3.8/site-packages/pint/registry.py in _convert(self, value, src, dst, inplace)
>    1828                 value, src = src._magnitude, src._units
>    1829 
> -> 1830         return super()._convert(value, src, dst, inplace)
>    1831 
>    1832     def _get_compatible_units(self, input_units, group_or_system):
> 
> ~/myproject-env/lib/python3.8/site-packages/pint/registry.py in _convert(self, value, src, dst, inplace)
>    1431 
>    1432         if not (src_offset_unit or dst_offset_unit):
> -> 1433             return super()._convert(value, src, dst, inplace)
>    1434 
>    1435         src_dim = self._get_dimensionality(src)
> 
> ~/myproject-env/lib/python3.8/site-packages/pint/registry.py in _convert(self, value, src, dst, inplace, check_dimensionality)
>     978             # then the conversion cannot be performed.
>     979             if src_dim != dst_dim:
> --> 980                 raise DimensionalityError(src, dst, src_dim, dst_dim)
>     981 
>     982         # Here src and dst have only multiplicative units left. Thus we can
> 
> DimensionalityError: Cannot convert from 'dimensionless' (dimensionless) to 'meter / second' ([length] / [time])

So I try:
nc = normal_component(cross['u'], cross['v'])

but result (nc) has no attribute.

How can I use absolute_momentum ?

My environment is follows:

• Platform: Ubuntu Linux
• python version: 3.8.10
• metpy version: 1.0.1

[23ccozad]

The likely problem is that your ds['u_wind'].data and ds['v_wind'].data are arrays, when they should be pint Quantity objects. If you add the following lines of code right after you read in ds, that should fix the issue. Of course, your data may not be in knots, so feel free to change that to the appropriate unit.

ds['u_wind'].data = units.Quantity(ds['u_wind'].data, 'knot')
ds['v_wind'].data = units.Quantity(ds['v_wind'].data, 'knot')

I also think we should adjust our documentation to make it clear that the data in the xarray DataArray needs to be a pint Quantity.

[jthielen]

@23ccozad This is intended to work without having Pint Quantities inside your DataArrays, as MetPy's xarray accessor is supposed to flexibly handle both Quantity-in-DataArray and unit-on-attribute cases for input data.

@ShunsukeHoshino Would you be able to report what the units attribute is on your u and v fields in both your original dataset and interpolated cross section dataset? Perhaps the units attribute is missing in your original dataset (in which case you will likely need to attach units manually either to the attribute or using a Quantity array as @23ccozad suggests), or it got stripped in the cross section interpolation (which would be a bug to be fixed in the cross_section or interpolate_to_slice functions).

[23ccozad]

@jthielen Good to know that this should work using the units attribute in the DataArray. I'm guessing that @ShunsukeHoshino has appropriate units on the original dataset and interpolated cross-section dataset and that maybe the units attribute is being lost in normal_component() as part of absolute_momentum(). If that's the case, normal_component() could be modified to ensure that the units attribute is part of the returned DataArray.

[jthielen]

@23ccozad That's a good possibility as well, since the multiplications in normal_component can strip attributes depending on xarray version and set options (xref pydata/xarray#3891). Working around that, or perhaps changing

MetPy/src/metpy/calc/cross_sections.py

Line 295 in edb9828

norm_wind = normal_component(u, v, index=index).metpy.convert_units('m/s')

to

norm_wind = normal_component(u.metpy.convert_units('m/s'), v.metpy.convert_units('m/s'), index=index)

may solve that kind of issue.

[23ccozad]

@jthielen I was originally thinking about putting the units back on the DataArray by adding

if 'units' in data_x.attrs:
    component_norm.attrs['units'] = data_x.attrs['units']

following

MetPy/src/metpy/calc/cross_sections.py

Line 205 in 15d0c80

component_norm = data_x * unit_norm[0] + data_y * unit_norm[1]

But I think both of our possible solutions produce the same result since the data ends up in meters per second one way or another.

[dopplershift]

Well, I had a minor preference for @23ccozad's solution since it mirrored the handling of grid_mapping and had fewer side effects.

However, @jthielen's suggestion is more robust to the potential case of u and v components having different units (however unlikely that is). So I think that's the way we should handle it.

[dopplershift]

@jthielen What's the rationale for explicitly converting to "m/s" rather than just using quantify()?

[jthielen]

@dopplershift I don't think there was a true rationale, just that it was the easiest-to-apply carry-over from the original implementation where units had to be stripped. Only reason to keep it now is that momentum in another speed unit like knots or mph just feels weird (since, given the arithmetic order, I believe the units of the normal wind component are what the returned units will be in absence of another change).

Also, fun little note on the side: unit attribute handling was part of the original implementation, but I believe it got cleaned away somewhere in the "great xarray refactor."

[dopplershift]

The refactor making a bunch of that superfluous is what I was thinking too. I think if we want to return momentum in "m/s" that's fine, but the code might be better structured to convert as the last step, making the purpose explicit.

[ShunsukeHoshino]

@23ccozad @jthielen

Thank you for suggestion.
The units attribute of ds['u'] and ds['v'] are both 'm/s' and these attributes are passed to cross['u'] and cross['v'] correctly, so I tried

m = absolute_momentum(cross['u'].metpy.quantify(), cross['v'].metpy.quantify()) .

This works well.

Thank you very much!

[23ccozad]

Glad you were able to find a solution @ShunsukeHoshino. Thanks for bringing up this issue so that absolute_momentum can be adjusted to better handle units from xarray DataArrays in the future.