Merge pull request #508 from hletrd/sum-combiner

Added sum combiner (#500)

AUTHORS.rst

@@ -43,6 +43,7 @@ Alphabetical list of contributors
 * Zè Vinícius (@mirca)
 * Josh Walawender (@joshwalawender)
 * Nathan Walker (@walkerna22)
+* Jiyong Youn (@hletrd)
 
 (If you have contributed to the ccdproc project and your name is missing,
 please send an email to the coordinators, or

CHANGES.rst

@@ -26,6 +26,8 @@ New Features
 
 - ``combine`` now accepts ``numpy.ndarray`` as the input ``img_list``. [#493, #503]
 
+- Added ``sum`` option in method for ``combime``. [#500, #508]
+
 
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^

ccdproc/combiner.py

@@ -447,6 +447,74 @@ def average_combine(self, scale_func=ma.average, scale_to=None,
         # return the combined image
         return combined_image
 
+    def sum_combine(self, sum_func=ma.sum, scale_to=None,
+                        uncertainty_func=ma.std):
+        """
+        Sum combine together a set of arrays.
+
+        A `~ccdproc.CCDData` object is returned with the data property
+        set to the sum of the arrays. If the data was masked or any
+        data have been rejected, those pixels will not be included in the
+        sum. A mask will be returned, and if a pixel has been
+        rejected in all images, it will be masked. The uncertainty of
+        the combined image is set by the multiplication of summation of
+        standard deviation of the input by square root of number of images.
+        Because sum_combine returns 'pure sum' with masked pixels ignored, if
+        re-scaled sum is needed, average_combine have to be used with
+        multiplication by number of images combined.
+
+        Parameters
+        ----------
+        scale_func : function, optional
+            Function to calculate the sum. Defaults to
+            `numpy.ma.sum`.
+
+        scale_to : float or None, optional
+            Scaling factor used in the sum combined image. If given,
+            it overrides `scaling`. Defaults to ``None``.
+
+        uncertainty_func : function, optional
+            Function to calculate uncertainty. Defaults to `numpy.ma.std`.
+
+        Returns
+        -------
+        combined_image: `~ccdproc.CCDData`
+            CCDData object based on the combined input of CCDData objects.
+        """
+        if scale_to is not None:
+            scalings = scale_to
+        elif self.scaling is not None:
+            scalings = self.scaling
+        else:
+            scalings = 1.0
+
+        # set up the data
+        data = sum_func(scalings * self.data_arr, axis=0)
+
+        # set up the mask
+        masked_values = self.data_arr.mask.sum(axis=0)
+        mask = (masked_values == len(self.data_arr))
+
+        # set up the deviation
+        uncertainty = uncertainty_func(self.data_arr, axis=0)
+        # Divide uncertainty by the number of pixel (#309)
+        uncertainty /= np.sqrt(len(self.data_arr) - masked_values)
+        # Convert uncertainty to plain numpy array (#351)
+        uncertainty = np.asarray(uncertainty)
+        # Multiply uncertainty by square root of the number of images
+        uncertainty *= len(self.data_arr) - masked_values
+
+        # create the combined image with a dtype that matches the combiner
+        combined_image = CCDData(np.asarray(data.data, dtype=self.dtype),
+                                 mask=mask, unit=self.unit,
+                                 uncertainty=StdDevUncertainty(uncertainty))
+
+        # update the meta data
+        combined_image.meta['NCOMBINE'] = len(self.data_arr)
+
+        # return the combined image
+        return combined_image
+
 
 def combine(img_list, output_file=None,
             method='average', weights=None, scale=None, mem_limit=16e9,
@@ -476,6 +544,7 @@ def combine(img_list, output_file=None,
 
         - ``'average'`` : To combine by calculating the average.
         - ``'median'`` : To combine by calculating the median.
+        - ``'sum'`` : To combine by calculating the sum.
 
         Default is ``'average'``.
 
@@ -547,8 +616,8 @@ def combine(img_list, output_file=None,
         Default is ``None``.
 
     combine_uncertainty_function : callable, None, optional
-        If ``None`` use the default uncertainty func when using average or
-        median combine, otherwise use the function provided.
+        If ``None`` use the default uncertainty func when using average, median or
+        sum combine, otherwise use the function provided.
         Default is ``None``.
 
     ccdkwargs : Other keyword arguments for `ccdproc.fits_ccddata_reader`.
@@ -574,6 +643,8 @@ def combine(img_list, output_file=None,
         combine_function = 'average_combine'
     elif method == 'median':
         combine_function = 'median_combine'
+    elif method == 'sum':
+        combine_function = 'sum_combine'
     else:
         raise ValueError("unrecognised combine method : {0}.".format(method))
 

ccdproc/tests/test_combiner.py

@@ -73,6 +73,9 @@ def test_combiner_dtype(ccd_data):
     med = c.median_combine()
     # dtype of median should match dtype of input
     assert med.dtype == c.dtype
+    result_sum = c.sum_combine()
+    # dtype of sum should match dtype of input
+    assert result_sum.dtype == c.dtype
 
 
 #test mask is created from ccd.data
@@ -205,6 +208,16 @@ def test_combiner_average(ccd_data):
     assert ccd.unit == u.adu
     assert ccd.meta['NCOMBINE'] == len(ccd_list)
 
+#test that the sum combination works and returns a ccddata object
+def test_combiner_sum(ccd_data):
+    ccd_list = [ccd_data, ccd_data, ccd_data]
+    c = Combiner(ccd_list)
+    ccd = c.sum_combine()
+    assert isinstance(ccd, CCDData)
+    assert ccd.shape == (100, 100)
+    assert ccd.unit == u.adu
+    assert ccd.meta['NCOMBINE'] == len(ccd_list)
+
 
 #test data combined with mask is created correctly
 def test_combiner_mask_average(ccd_data):
@@ -270,6 +283,22 @@ def test_combiner_mask_median(ccd_data):
     assert ccd.mask[0, 0]
     assert not ccd.mask[5, 5]
 
+
+#test data combined with mask is created correctly
+def test_combiner_mask_sum(ccd_data):
+    data = np.zeros((10, 10))
+    data[5, 5] = 1
+    mask = (data == 0)
+    ccd = CCDData(data, unit=u.adu, mask=mask)
+    ccd_list = [ccd, ccd, ccd]
+    c = Combiner(ccd_list)
+    ccd = c.sum_combine()
+    assert ccd.data[0, 0] == 0
+    assert ccd.data[5, 5] == 3
+    assert ccd.mask[0, 0]
+    assert not ccd.mask[5, 5]
+
+
 #test combiner convenience function reads fits file and combine as expected
 def test_combine_average_fitsimages():
     fitsfile = get_pkg_data_filename('data/a8280271.fits')
@@ -395,6 +424,19 @@ def test_median_combine_uncertainty(ccd_data):
     np.testing.assert_array_equal(ccd.data, ccd2.data)
     np.testing.assert_array_equal(ccd.uncertainty.array, ccd2.uncertainty.array)
 
+#test the optional uncertainty function in sum_combine
+def test_sum_combine_uncertainty(ccd_data):
+    ccd_list = [ccd_data, ccd_data, ccd_data]
+    c = Combiner(ccd_list)
+    ccd = c.sum_combine(uncertainty_func=np.sum)
+    uncert_ref = np.sum(c.data_arr, 0) * np.sqrt(3)
+    np.testing.assert_almost_equal(ccd.uncertainty.array, uncert_ref)
+
+    # Compare this also to the "combine" call
+    ccd2 = combine(ccd_list, method='sum', combine_uncertainty_function=np.sum)
+    np.testing.assert_array_equal(ccd.data, ccd2.data)
+    np.testing.assert_array_equal(ccd.uncertainty.array, ccd2.uncertainty.array)
+
 
 # test resulting uncertainty is corrected for the number of images
 def test_combiner_uncertainty_average():
@@ -428,6 +470,23 @@ def test_combiner_uncertainty_average_mask():
     np.testing.assert_array_almost_equal(ccd.uncertainty.array,
                                          ref_uncertainty)
 
+# test resulting uncertainty is corrected for the number of images (with mask)
+def test_combiner_uncertainty_sum_mask():
+    mask = np.zeros((10, 10), dtype=np.bool_)
+    mask[5, 5] = True
+    ccd_with_mask = CCDData(np.ones((10, 10)), unit=u.adu, mask=mask)
+    ccd_list = [ccd_with_mask,
+                CCDData(np.ones((10, 10))*2, unit=u.adu),
+                CCDData(np.ones((10, 10))*3, unit=u.adu)]
+    c = Combiner(ccd_list)
+    ccd = c.sum_combine()
+    # Just the standard deviation of ccd data.
+    ref_uncertainty = np.ones((10, 10)) * np.std([1, 2, 3])
+    ref_uncertainty *= np.sqrt(3)
+    ref_uncertainty[5, 5] = np.std([2, 3]) * np.sqrt(2)
+    np.testing.assert_array_almost_equal(ccd.uncertainty.array,
+                                         ref_uncertainty)
+
 def test_combiner_3d():
     data1 = CCDData(3 * np.ones((5,5,5)), unit=u.adu)
     data2 = CCDData(2 * np.ones((5,5,5)), unit=u.adu)
@@ -487,7 +546,7 @@ def test_clip_extrema_3d():
     np.testing.assert_array_equal(result, expected)
 
 
-@pytest.mark.parametrize('comb_func', ['average_combine', 'median_combine'])
+@pytest.mark.parametrize('comb_func', ['average_combine', 'median_combine', 'sum_combine'])
 def test_writeable_after_combine(ccd_data, tmpdir, comb_func):
     tmp_file = tmpdir.join('tmp.fits')
     from ..combiner import Combiner