Add 'first' and 'last' moasicking modes to reproject_and_coadd

reproject/mosaicking/coadd.py

@@ -73,9 +73,12 @@
         `~astropy.io.fits.HDUList` instance, specifies the HDU to use.
     reproject_function : callable
         The function to use for the reprojection
-    combine_function : { 'mean', 'sum', 'median' }
+    combine_function : { 'mean', 'sum', 'median', 'first', 'last' }
         The type of function to use for combining the values into the final
-        image.
+        image. For 'first' and 'last', respectively, the reprojected images are
+        simply overlaid on top of each other. With respect to the order of the
+        input images in ``input_data``, either the first or the last image to
+        cover a region of overlap determines the output data for that region.
     match_background : bool
         Whether to match the backgrounds of the images.
     background_reference : `None` or `int`
@@ -94,8 +97,8 @@
 
     # Validate inputs
 
-    if combine_function not in ("mean", "sum", "median"):
-        raise ValueError("combine_function should be one of mean/sum/median")
+    if combine_function not in ("mean", "sum", "median", "first", "last"):
+        raise ValueError("combine_function should be one of mean/sum/median/first/last")
 
     if reproject_function is None:
         raise ValueError(
@@ -233,7 +236,23 @@
         if combine_function == "mean":
             with np.errstate(invalid="ignore"):
                 final_array /= final_footprint
-
+    elif combine_function == "first":
+        for array in arrays:
+            mask = final_footprint[array.view_in_original_array] == 0
+            final_footprint[array.view_in_original_array] = np.where(
+                mask, array.footprint, final_footprint[array.view_in_original_array]
+            )
+            final_array[array.view_in_original_array] = np.where(
+                mask, array.array, final_array[array.view_in_original_array]
+            )
+    elif combine_function == "last":
+        for array in arrays:
+            final_footprint[array.view_in_original_array] = np.where(
+                array.footprint, array.footprint, final_footprint[array.view_in_original_array]
+            )
+            final_array[array.view_in_original_array] = np.where(
+                array.footprint > 0, array.array, final_array[array.view_in_original_array]
+            )
     elif combine_function == "median":
         # Here we need to operate in chunks since we could otherwise run
         # into memory issues

reproject/mosaicking/tests/test_coadd.py

@@ -8,7 +8,7 @@
 from astropy.io import fits
 from astropy.io.fits import Header
 from astropy.wcs import WCS
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_equal
 
 from reproject import reproject_exact, reproject_interp
 from reproject.mosaicking.coadd import reproject_and_coadd
@@ -41,11 +41,10 @@ def _get_tiles(self, views):
 
         input_data = []
 
-        for jmin, jmax, imin, imax in views:
-            array = self.array[jmin:jmax, imin:imax].copy()
+        for view in views:
+            array = self.array[view].copy()
             wcs = self.wcs.deepcopy()
-            wcs.wcs.crpix[0] -= imin
-            wcs.wcs.crpix[1] -= jmin
+            wcs = wcs[view]
             input_data.append((array, wcs))
 
         return input_data
@@ -58,7 +57,7 @@ def _nonoverlapping_views(self):
         views = []
         for i in range(4):
             for j in range(5):
-                views.append((je[j], je[j + 1], ie[i], ie[i + 1]))
+                views.append(np.s_[je[j] : je[j + 1], ie[i] : ie[i + 1]])
 
         return views
 
@@ -70,11 +69,11 @@ def _overlapping_views(self):
         views = []
         for i in range(4):
             for j in range(5):
-                views.append((je[j], je[j + 1] + 10, ie[i], ie[i + 1] + 10))
+                views.append(np.s_[je[j] : je[j + 1] + 10, ie[i] : ie[i + 1] + 10])
 
         return views
 
-    @pytest.mark.parametrize("combine_function", ["mean", "sum"])
+    @pytest.mark.parametrize("combine_function", ["mean", "sum", "first", "last"])
     def test_coadd_no_overlap(self, combine_function, reproject_function):
         # Make sure that if all tiles are exactly non-overlapping, and
         # we use 'sum' or 'mean', we get the exact input array back.
@@ -109,6 +108,38 @@ def test_coadd_with_overlap(self, reproject_function):
 
         assert_allclose(array, self.array, atol=ATOL)
 
+    @pytest.mark.parametrize("combine_function", ["first", "last"])
+    def test_coadd_with_overlap_first_last(self, reproject_function, combine_function):
+        views = self._overlapping_views
+        input_data = self._get_tiles(views)
+
+        # Make each of the overlapping tiles different
+        for i, (array, wcs) in enumerate(input_data):
+            input_data[i] = (np.full_like(array, i), wcs)
+
+        array, footprint = reproject_and_coadd(
+            input_data,
+            self.wcs,
+            shape_out=self.array.shape,
+            combine_function=combine_function,
+            reproject_function=reproject_function,
+        )
+
+        # Test that either the correct tile sets the output value in the overlap regions
+        test_sequence = list(enumerate(views))
+        if combine_function == "last":
+            test_sequence = test_sequence[::-1]
+        for i, view in test_sequence:
+            # Each tile in test_sequence should overwrite teh following tiles
+            # in the overlap regions. We'll use nans to mark pixels in the
+            # output array that have already been set by a preceeding tile, so
+            # we'll go through, check that each tile matches the non-nan pixels
+            # in its region, and then set that whole region to nan.
+            output_tile = array[view]
+            output_values = output_tile[np.isfinite(output_tile)]
+            assert_equal(output_values, i)
+            array[view] = np.nan
+
     def test_coadd_background_matching(self, reproject_function):
         # Test out the background matching