Multiplying two cubes with different units

[will-henney]

Hi,

I am new to MPDAF and have run into a problem when trying to multiply together two Cube objects that have different units. Here is a minimal example:

import numpy as np
from mpdaf.obj import Cube
import astropy.units as u
a = Cube(data=np.ones((1, 1, 1)), unit=u.s)
b = Cube(data=np.ones((1, 1, 1)), unit=u.s)
c = Cube(data=np.ones((1, 1, 1)), unit=u.cm)
# Multiplying cubes with the same units is OK
a * b
# Multiplying cubes with different units yields an error
a * c

This throws an error UnitConversionError: 'cm' (length) and 's' (time) are not convertible (full traceback below).

I can see why adding cubes with different units should yield an error, but multiplying should be OK, right? I would expect to get a result with units u.cm * u.s.

I can work around the issue by using the .data and .unit attributes directly, but I was hoping to take advantage of the automatic variance propagation.

Thanks for any advice you might have

Will

Full error traceback:

---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
<ipython-input-2-7d2666f7101b> in <module>
      8 a * b
      9 # Multiplying cubes with different units yields an error
---> 10 a * c

~/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/arithmetic.py in __mul__(self, other)
    198             return res
    199         else:
--> 200             return _arithmetic(ma.multiply, self, other)
    201 
    202     def __div__(self, other):

~/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/arithmetic.py in _arithmetic(operation, a, b)
    162     return a.__class__.new_from_obj(
    163         a, copy=False, unit=unit,
--> 164         data=_arithmetic_data(operation, a, b, newshape=newshape),
    165         var=_arithmetic_var(operation, a, b, newshape=newshape)
    166     )

~/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/arithmetic.py in _arithmetic_data(operation, a, b, newshape)
     84 def _arithmetic_data(operation, a, b, newshape=None):
     85     if a.unit != b.unit:
---> 86         data = UnitMaskedArray(b.data, b.unit, a.unit)
     87     else:
     88         data = b.data

~/miniconda3/envs/py39/lib/python3.9/site-packages/mpdaf/obj/objs.py in UnitMaskedArray(mask_array, old_unit, new_unit)
    253         return mask_array
    254     else:
--> 255         return np.ma.array(Quantity(mask_array.data, old_unit, copy=False)
    256                            .to(new_unit).value, mask=mask_array.mask)

~/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/quantity.py in to(self, unit, equivalencies, copy)
    696             # Avoid using to_value to ensure that we make a copy. We also
    697             # don't want to slow down this method (esp. the scalar case).
--> 698             value = self._to_value(unit, equivalencies)
    699         else:
    700             # to_value only copies if necessary

~/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/quantity.py in _to_value(self, unit, equivalencies)
    660         if equivalencies == []:
    661             equivalencies = self._equivalencies
--> 662         return self.unit.to(unit, self.view(np.ndarray),
    663                             equivalencies=equivalencies)
    664 

~/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/core.py in to(self, other, value, equivalencies)
    985             return UNITY
    986         else:
--> 987             return self._get_converter(other, equivalencies=equivalencies)(value)
    988 
    989     def in_units(self, other, value=1.0, equivalencies=[]):

~/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/core.py in _get_converter(self, other, equivalencies)
    916                             pass
    917 
--> 918             raise exc
    919 
    920     def _to(self, other):

~/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/core.py in _get_converter(self, other, equivalencies)
    901         # if that doesn't work, maybe we can do it with equivalencies?
    902         try:
--> 903             return self._apply_equivalencies(
    904                 self, other, self._normalize_equivalencies(equivalencies))
    905         except UnitsError as exc:

~/miniconda3/envs/py39/lib/python3.9/site-packages/astropy/units/core.py in _apply_equivalencies(self, unit, other, equivalencies)
    884         other_str = get_err_str(other)
    885 
--> 886         raise UnitConversionError(
    887             "{} and {} are not convertible".format(
    888                 unit_str, other_str))

UnitConversionError: 'cm' (length) and 's' (time) are not convertible