Support Subtract and Division of NDCube by NDData

changelog/880.bugfix.1.rst

@@ -0,0 +1 @@
+Adding unitful `~ndcube.NDCube` and ``astropy.nddata.NDData`` objects backed by ``dask`` did not preserve underlying arrays as ``dask`` arrays. This is now fixed.

changelog/880.bugfix.2.rst

@@ -0,0 +1 @@
+Fix bug where error was returned rather than raised with trying to perform arithmetic operation between `~ndcube.NDCube` and an object whose WCS attribute is not ``None``.

changelog/880.feature.rst

@@ -0,0 +1 @@
+Enable subtraction and division of `~ndcube.NDCube` by an `~astropy.nddata.NDData` instance, including uncertainty, mask and unit support.

docs/explaining_ndcube/index.rst

@@ -14,3 +14,4 @@ Explaining ``ndcube``
    tabular_coordinates
    reproject
    visualization
+   arithmetic

ndcube/conftest.py

@@ -832,6 +832,13 @@ def ndcube_2d_dask(wcs_2d_lt_ln):
     return NDCube(da, wcs=wcs_2d_lt_ln, uncertainty=da_uncert, mask=da_mask, unit=u.J)
 
 
+@pytest.fixture
+def nddata_2d_dask(ndcube_2d_dask):
+    value = astropy.nddata.NDData(ndcube_2d_dask)
+    value._wcs = None
+    return value
+
+
 @pytest.fixture
 def ndcube_2d(request):
     """

ndcube/ndcube.py

@@ -980,34 +980,49 @@ def _arithmetic_handle_mask(self, self_mask, value_mask):
     def _arithmetic_operate_with_nddata(self, operation, value):
         handle_mask = self._arithmetic_handle_mask
         if value.wcs is not None:
-            return TypeError("Cannot add coordinate-aware NDCubes together.")
-
+            raise TypeError("Cannot add coordinate-aware objects to NDCubes.")
         kwargs = {}
         if operation == "add":
             # Handle units
             if self.unit is not None and value.unit is not None:
-                value_data = (value.data * value.unit).to_value(self.unit)
+                value_data = value.data * (value.unit / self.unit).to(u.dimensionless_unscaled)
             elif self.unit is None and value.unit is None:
                 value_data = value.data
             else:
                 raise TypeError("Adding objects requires that both have a unit or neither has a unit.") # change the test as well.
             # Handle data and uncertainty
             kwargs["data"] = self.data + value_data
             uncert_op = np.add
-        elif operation == "multiply":
+        elif operation in ("multiply", "true_divide"):
             # Handle units
             if self.unit is not None or value.unit is not None:
                 cube_unit = u.Unit('') if self.unit is None else self.unit
                 value_unit = u.Unit('') if value.unit is None else value.unit
-                kwargs["unit"] = cube_unit * value_unit
-            kwargs["data"] = self.data * value.data
-            uncert_op = np.multiply
+                kwargs["unit"] = (cube_unit * value_unit if operation == "multiply"
+                                  else cube_unit / value_unit)
+            if operation == "multiply":
+                kwargs["data"] = self.data * value.data
+                uncert_op = np.multiply
+            else:
+                kwargs["data"] = self.data / value.data
+                uncert_op = np.true_divide
         else:
             raise ValueError("Value of operation argument is not recognized.")
-        kwargs["uncertainty"] = self._combine_uncertainty(uncert_op, value, kwargs["data"])
+        # Calculate uncertainty.
+        new_uncert = self._combine_uncertainty(uncert_op, value, kwargs["data"])
+        if new_uncert:
+            # New uncertainty object must be decoupled from its original
+            # parent_nddata object. Set this to None here, and the parent_nddata
+            # will become the new cube on instantiation.
+            new_uncert.parent_nddata = None
+            uncert_unit = kwargs.get("unit", self.unit)
+            if uncert_unit:
+                # Give uncertainty object the same units as the new NDCube.
+                new_uncert.unit = uncert_unit
+        kwargs["uncertainty"] = new_uncert
         kwargs["mask"] = handle_mask(self.mask, value.mask)
 
-        return kwargs  # return the new NDCube instance
+        return kwargs
 
     def __add__(self, value):
         kwargs = {}
@@ -1052,7 +1067,12 @@ def __radd__(self, value):
         return self.__add__(value)
 
     def __sub__(self, value):
-        return self.__add__(-value)
+        if isinstance(value, NDData):
+            new_value = copy.copy(value)
+            new_value._data = -value.data
+        else:
+            new_value = -value
+        return self.__add__(new_value)
 
     def __rsub__(self, value):
         return self.__neg__().__add__(value)
@@ -1084,6 +1104,9 @@ def __rmul__(self, value):
         return self.__mul__(value)
 
     def __truediv__(self, value):
+        if isinstance(value, NDData):
+            kwargs = self._arithmetic_operate_with_nddata("true_divide", value)
+            return self._new_instance(**kwargs)
         return self.__mul__(1/value)
 
     def __rtruediv__(self, value):
@@ -1483,6 +1506,10 @@ def to_nddata(self,
         Attribute values can be altered on the output object by setting a kwarg with the new
         value, e.g. ``data=new_data``.
         Any attributes not supported by the new class (``nddata_type``), will be discarded.
+        A motivating use case for this method is in enabling arithmetic operations between
+        `~ndcube.NDCube` instances by removing coordinate-awareness.  See the section of the
+        ndcube documentation on
+        'Enabling Arithmetic Operations between NDCubes with NDCube.to_nddata'.
 
         Parameters
         ----------

ndcube/tests/test_ndcube_arithmetic.py

@@ -122,6 +122,73 @@ def test_cube_arithmetic_subtract(ndcube_2d_ln_lt_units, value):
     check_arithmetic_value_and_units(new_cube, cube_quantity - value)
 
 
+@pytest.mark.parametrize(("ndc", "value", "expected_kwargs"),
+                        [(
+                          "ndcube_2d_ln_lt_no_unit_no_unc_no_mask_2",
+                          NDData(np.ones((2, 3)), wcs=None),
+                          {"data": np.array([[-1, 0, 1], [2, 3, 4]])}
+                         ),
+                         (
+                          "ndcube_2d_ln_lt_no_unit_no_unc_no_mask_2",
+                          NDData(np.ones((2, 3)),
+                                 wcs=None,
+                                 uncertainty=StdDevUncertainty(np.ones((2, 3))*0.1),
+                                 mask=np.array([[True, False, False], [False, True, False]])),
+                          {"data": np.array([[-1, 0, 1], [2, 3, 4]]),
+                           "uncertainty": astropy.nddata.StdDevUncertainty(np.ones((2, 3))*0.1),
+                           "mask": np.array([[True, False, False], [False, True, False]])}
+                         ), # ndc has no mask no uncertainty no unit, but nddata has all.
+                         (
+                          "ndcube_2d_ln_lt_unit_unc_mask",
+                          NDData(np.ones((2, 3)), wcs=None, unit=u.ct),
+                          {"data": np.array([[-1, 0, 1], [2, 3, 4]]),
+                           "uncertainty": astropy.nddata.StdDevUncertainty(np.array([[0, 0.05, 0.1],
+                                                                                     [0.15, 0.2, 0.25]])),
+                           "mask": np.array([[False, True, True], [False, True, True]])}
+                         ), # ndc has mask and uncertainty unit, but nddata doesn't.
+                         (
+                          "ndcube_2d_ln_lt_unit_unc_mask",
+                          NDData(np.ones((2, 3)),
+                                 wcs=None,
+                                 uncertainty=StdDevUncertainty(np.ones((2, 3))*0.1),
+                                 mask=np.array([[True, False, False], [False, True, False]]),
+                                 unit=u.ct),
+                          {"unit": u.ct,
+                           "data": np.array([[-1, 0, 1], [2, 3, 4]]),
+                           "uncertainty": astropy.nddata.StdDevUncertainty(np.array([[0.1       , 0.1118034 , 0.14142136],
+                                                                                     [0.18027756, 0.2236068 , 0.26925824]])),
+                           "mask": np.array([[True, True, True], [False, True, True]])}
+                         ) # both of them have uncertainty and mask and unit.
+                        ],
+                        indirect=("ndc",))
+def test_cube_arithmetic_subtract_nddata(ndc, value, expected_kwargs, wcs_2d_lt_ln):
+    output_cube = ndc - value  # perform the subtraction
+
+    expected_kwargs["wcs"] = wcs_2d_lt_ln
+    expected_cube = NDCube(**expected_kwargs)
+
+    # Assert output cube is same as expected cube
+    assert_cubes_equal(output_cube, expected_cube, check_uncertainty_values=True)
+
+
+def test_cube_dask_arithmetic_subtract_nddata(ndcube_2d_dask):
+    ndc = ndcube_2d_dask
+    value = NDData(np.ones(ndc.data.shape), wcs=None, unit=ndc.unit)
+    output_cube = ndc - value
+    assert type(output_cube.data) is type(ndc.data)
+    assert type(output_cube.uncertainty.array) is type(ndc.uncertainty.array)
+    assert type(output_cube.mask) is type(ndc.mask)
+
+
+def test_cube_arithmetic_subtract_nddata_dask(wcs_2d_lt_ln, nddata_2d_dask):
+    value = nddata_2d_dask
+    ndc = NDCube(np.ones(value.data.shape), wcs=wcs_2d_lt_ln, unit=value.unit)
+    output_cube = ndc - value
+    assert type(output_cube.data) is type(value.data)
+    assert type(output_cube.uncertainty.array) is type(value.uncertainty.array)
+    assert type(output_cube.mask) is type(value.mask)
+
+
 @pytest.mark.parametrize('value', [
     10 * u.ct,
     u.Quantity([10], u.ct),
@@ -134,6 +201,24 @@ def test_cube_arithmetic_rsubtract(ndcube_2d_ln_lt_units, value):
     check_arithmetic_value_and_units(new_cube, value - cube_quantity)
 
 
+def test_cube_dask_arithmetic_rsubtract_nddata(ndcube_2d_dask):
+    ndc = ndcube_2d_dask
+    value = NDData(np.ones(ndc.data.shape), wcs=None, unit=ndc.unit)
+    output_cube = value - ndc
+    assert type(output_cube.data) is type(ndc.data)
+    assert type(output_cube.uncertainty.array) is type(ndc.uncertainty.array)
+    assert type(output_cube.mask) is type(ndc.mask)
+
+
+def test_cube_arithmetic_rsubtract_nddata_dask(wcs_2d_lt_ln, nddata_2d_dask):
+    value = nddata_2d_dask
+    ndc = NDCube(np.ones(value.data.shape), wcs=wcs_2d_lt_ln, unit=value.unit)
+    output_cube = value - ndc
+    assert type(output_cube.data) is type(value.data)
+    assert type(output_cube.uncertainty.array) is type(value.uncertainty.array)
+    assert type(output_cube.mask) is type(value.mask)
+
+
 @pytest.mark.parametrize('value', [
     10 * u.ct,
     u.Quantity([10], u.ct),
@@ -231,13 +316,84 @@ def test_cube_arithmetic_divide(ndcube_2d_ln_lt_units, value):
     new_cube = ndcube_2d_ln_lt_units / value
     check_arithmetic_value_and_units(new_cube, cube_quantity / value)
 
+
+@pytest.mark.parametrize(("ndc", "value", "expected_kwargs"),
+                        [(
+                          "ndcube_2d_ln_lt_no_unit_no_unc_no_mask_2",
+                          NDData(np.ones((2, 3)) + 1, wcs=None),
+                          {"data": np.array([[0, 0.5, 1], [1.5, 2, 2.5]])},
+                         ),
+                         (
+                          "ndcube_2d_ln_lt_no_unit_no_unc_no_mask_2",
+                          NDData(np.ones((2, 3)) + 1,
+                                 wcs=None,
+                                 uncertainty=StdDevUncertainty((np.ones((2, 3)) + 1) * 0.1),
+                                 mask=np.array([[True, False, False], [False, True, False]]),
+                                 unit=u.ct),
+                           {"data": np.array([[0, 0.5, 1], [1.5, 2, 2.5]]),
+                           "uncertainty": astropy.nddata.StdDevUncertainty((np.ones((2, 3)) + 1) * 0.1),
+                           "mask": np.array([[True, False, False], [False, True, False]]),
+                           "unit": u.dimensionless_unscaled / u.ct} # ndc has no mask no uncertainty no unit, but nddata has all.
+                         ),
+                         (
+                          "ndcube_2d_ln_lt_unit_unc_mask",
+                          NDData(np.ones((2, 3)) * 2, wcs=None),
+                           {"data": np.array([[0, 0.5, 1], [1.5, 2, 2.5]]),
+                           "uncertainty": astropy.nddata.StdDevUncertainty(np.array([[0, 0.05, 0.1],
+                                                                                     [0.15, 0.2, 0.25]])),
+                           "mask": np.array([[False, True, True], [False, True, True]])}
+                         ), # ndc has mask, uncertainty and unit, but nddata doesn't.
+                         (
+                         "ndcube_2d_ln_lt_unit_unc_mask",
+                          NDData(np.ones((2, 3)) + 1,
+                                 wcs=None,
+                                 uncertainty=StdDevUncertainty((np.ones((2, 3)) + 1) * 0.1),
+                                 mask=np.array([[True, False, False], [False, True, False]]),
+                                 unit=u.ct),
+                          {"unit": u.dimensionless_unscaled,
+                           "data": np.array([[0, 0.5, 1], [1.5, 2, 2.5]]),
+                           "uncertainty": astropy.nddata.StdDevUncertainty(np.array([[0.       , 0.0559017, 0.1118034],
+                                                                                     [0.1677051, 0.2236068, 0.2795085]])),
+                           "mask": np.array([[True, True, True], [False, True, True]])}
+                         ) # both of them have uncertainty and mask and unit.
+                        ],
+                        indirect=("ndc",))
+def test_cube_arithmetic_divide_nddata(ndc, value, expected_kwargs, wcs_2d_lt_ln):
+    output_cube = ndc / value  # perform the division
+
+    expected_kwargs["wcs"] = wcs_2d_lt_ln
+    expected_cube = NDCube(**expected_kwargs)
+
+    # Assert output cube is same as expected cube
+    assert_cubes_equal(output_cube, expected_cube, check_uncertainty_values=True)
+
+
+def test_cube_dask_arithmetic_divide_nddata(ndcube_2d_dask):
+    ndc = ndcube_2d_dask
+    value = NDData(np.ones(ndc.data.shape), wcs=None, unit=ndc.unit)
+    output_cube = ndc / value
+    assert type(output_cube.data) is type(ndc.data)
+    assert type(output_cube.uncertainty.array) is type(ndc.uncertainty.array)
+    assert type(output_cube.mask) is type(ndc.mask)
+
+
+def test_cube_arithmetic_divide_nddata_dask(wcs_2d_lt_ln, nddata_2d_dask):
+    value = nddata_2d_dask
+    ndc = NDCube(np.ones(value.data.shape), wcs=wcs_2d_lt_ln, unit=value.unit)
+    output_cube = ndc / value
+    assert type(output_cube.data) is type(value.data)
+    assert type(output_cube.uncertainty.array) is type(value.uncertainty.array)
+    assert type(output_cube.mask) is type(value.mask)
+
+
 @pytest.mark.parametrize('value', [1, 2, -1])
 def test_cube_arithmetic_rdivide(ndcube_2d_ln_lt_units, value):
     cube_quantity = u.Quantity(ndcube_2d_ln_lt_units.data, ndcube_2d_ln_lt_units.unit)
     with np.errstate(divide='ignore'):
         new_cube =  value / ndcube_2d_ln_lt_units
         check_arithmetic_value_and_units(new_cube,  value / cube_quantity)
 
+
 @pytest.mark.parametrize('value', [1, 2, -1])
 def test_cube_arithmetic_rdivide_uncertainty(ndcube_4d_unit_uncertainty, value):
     cube_quantity = u.Quantity(ndcube_4d_unit_uncertainty.data, ndcube_4d_unit_uncertainty.unit)
@@ -247,6 +403,27 @@ def test_cube_arithmetic_rdivide_uncertainty(ndcube_4d_unit_uncertainty, value):
             new_cube =  value / ndcube_4d_unit_uncertainty
             check_arithmetic_value_and_units(new_cube,  value / cube_quantity)
 
+
+def test_cube_dask_arithmetic_rdivide_nddata(ndcube_2d_dask):
+    ndc = ndcube_2d_dask
+    value = NDData(np.ones(ndc.data.shape), wcs=None, unit=ndc.unit)
+    match = "does not support propagation of uncertainties for power. Setting uncertainties to None."
+    with pytest.warns(NDCubeUserWarning, match=match):  # noqa: PT031
+        with np.errstate(divide='ignore'):
+            output_cube = value / ndc
+            assert type(output_cube.data) is type(ndc.data)
+            assert type(output_cube.mask) is type(ndc.mask)
+
+
+def test_cube_arithmetic_rdivide_nddata_dask(wcs_2d_lt_ln, nddata_2d_dask):
+    value = nddata_2d_dask
+    ndc = NDCube(np.ones(value.data.shape), wcs=wcs_2d_lt_ln, unit=value.unit)
+    output_cube = value / ndc
+    assert type(output_cube.data) is type(value.data)
+    assert type(output_cube.uncertainty.array) is type(value.uncertainty.array)
+    assert type(output_cube.mask) is type(value.mask)
+
+
 def test_cube_arithmetic_neg(ndcube_2d_ln_lt_units):
     check_arithmetic_value_and_units(
         -ndcube_2d_ln_lt_units,