diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py
index 535114d92..1e68e2ce2 100644
--- a/reproject/mosaicking/coadd.py
+++ b/reproject/mosaicking/coadd.py
@@ -276,11 +276,12 @@ def reproject_and_coadd(
         if combine_function == "mean":
             with np.errstate(invalid="ignore"):
                 output_array /= output_footprint
+                output_array[output_footprint == 0] = 0
 
     elif combine_function in ("first", "last", "min", "max"):
         for array in arrays:
             if combine_function == "first":
-                mask = output_footprint[array.view_in_original_array] == 0
+                mask = (output_footprint[array.view_in_original_array] == 0) & (array.footprint > 0)
             elif combine_function == "last":
                 mask = array.footprint > 0
             elif combine_function == "min":
diff --git a/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits b/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits
index 1ec779a00..7d8747d26 100644
Binary files a/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits and b/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits differ
diff --git a/reproject/mosaicking/tests/test_coadd.py b/reproject/mosaicking/tests/test_coadd.py
index 9738cdf2f..b8ff25e0b 100644
--- a/reproject/mosaicking/tests/test_coadd.py
+++ b/reproject/mosaicking/tests/test_coadd.py
@@ -236,6 +236,39 @@ def test_coadd_background_matching_one_array(self, reproject_function):
         np.testing.assert_allclose(array, array_matched)
         np.testing.assert_allclose(footprint, footprint_matched)
 
+    @pytest.mark.parametrize("combine_function", ["first", "last", "min", "max", "sum", "mean"])
+    @pytest.mark.parametrize("match_background", [True, False])
+    def test_footprint_correct(self, reproject_function, combine_function, match_background):
+        # Test that the output array is zero outside the returned footprint
+        # We're running this test over a somewhat large grid of parameters, so
+        # cut down the array size to avoid increasing the total test runtime
+        # too much.
+        slice = np.s_[::3, ::3]
+        wcs1 = self.wcs[slice]
+        wcs2 = self.wcs.deepcopy()
+        # Add a 45-degree rotation
+        wcs2.wcs.pc = np.array([[0.5, -0.5], [0.5, 0.5]])
+
+        wcs_out = wcs1.deepcopy()
+        # Expand the output WCS to go beyond the input images
+        wcs_out.wcs.cdelt = 2 * wcs1.wcs.cdelt[0], 2 * wcs1.wcs.cdelt[1]
+
+        # Ensure the input data is fully non-zero, so we can tell where data
+        # got projected to in the output image.
+        array1 = np.full_like(self.array[slice], 2)
+        array2 = array1 + 5
+
+        array, footprint = reproject_and_coadd(
+            [(array1, wcs1), (array2, wcs2)],
+            wcs_out,
+            shape_out=self.array.shape,
+            combine_function=combine_function,
+            reproject_function=reproject_function,
+            match_background=match_background,
+        )
+
+        assert np.all((array != 0) == (footprint > 0))
+
     def test_coadd_background_matching_with_nan(self, reproject_function):
         # Test out the background matching when NaN values are present. We do
         # this by using three arrays with the same footprint but with different