Change default to not include local background in PSF model/residual images

CHANGES.rst

@@ -58,6 +58,11 @@ API Changes
   - The ``GridddedPSFModel`` string representations now include the
     model ``flux``, ``x_0``, and ``y_0`` parameters. [#1680]
 
+  - The PSF photometry ``make_model_image`` and ``make_residual_image``
+    methods no longer include the local background by default. This is a
+    backwards-incompatible change. If the previous behavior is desired,
+    set ``include_localbkg=True``. [#1703]
+
 - ``photutils.segmentation``
 
   - The ``SourceCatalog`` ``get_label`` and ``get_labels`` methods now

photutils/psf/photometry.py

@@ -1027,9 +1027,10 @@ def __call__(self, data, *, mask=None, error=None, init_params=None):
 
         return source_tbl
 
-    def make_model_image(self, shape, psf_shape, include_localbkg=True):
+    def make_model_image(self, shape, psf_shape, *, include_localbkg=False):
         """
-        Create a 2D image from the fit PSF models and local background.
+        Create a 2D image from the fit PSF models and optional local
+        background.
 
         Parameters
         ----------
@@ -1042,8 +1043,11 @@ def make_model_image(self, shape, psf_shape, include_localbkg=True):
 
         include_localbkg : bool, optional
             Whether to include the local background in the rendered
-            output image.
-            Default is True.
+            output image. Note that the local background level is
+            included around each source over the region defined by
+            ``psf_shape``. Thus, regions where the ``psf_shape`` of
+            sources overlap will have the local background added
+            multiple times.
 
         Returns
         -------
@@ -1078,7 +1082,7 @@ def make_model_image(self, shape, psf_shape, include_localbkg=True):
 
         return data
 
-    def make_residual_image(self, data, psf_shape, include_localbkg=True):
+    def make_residual_image(self, data, psf_shape, *, include_localbkg=False):
         """
         Create a 2D residual image from the fit PSF models and local
         background.
@@ -1096,8 +1100,10 @@ def make_residual_image(self, data, psf_shape, include_localbkg=True):
 
         include_localbkg : bool, optional
             Whether to include the local background in the subtracted
-            model.
-            Default is True.
+            model. Note that the local background level is subtracted
+            around each source over the region defined by ``psf_shape``.
+            Thus, regions where the ``psf_shape`` of sources overlap
+            will have the local background subtracted multiple times.
 
         Returns
         -------
@@ -1107,8 +1113,8 @@ def make_residual_image(self, data, psf_shape, include_localbkg=True):
         """
         if isinstance(data, NDData):
             residual = deepcopy(data)
-            residual.data[:] = self.make_residual_image(data.data, psf_shape,
-                                                        include_localbkg=include_localbkg)
+            residual.data[:] = self.make_residual_image(
+                data.data, psf_shape, include_localbkg=include_localbkg)
         else:
             unit = None
             if isinstance(data, u.Quantity):

photutils/psf/tests/test_photometry.py

@@ -205,6 +205,44 @@ def test_psf_photometry(test_data):
     assert resid_data3.unit == unit
 
 
+@pytest.mark.skipif(not HAS_SCIPY, reason='scipy is required')
+def test_model_residual_image(test_data):
+    data, error, _ = test_data
+
+    data = data + 10
+    psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35)
+    fit_shape = (5, 5)
+    finder = DAOStarFinder(16.0, 2.0)
+    bkgstat = MMMBackground()
+    localbkg_estimator = LocalBackground(5, 10, bkgstat)
+    psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder,
+                            aperture_radius=4,
+                            localbkg_estimator=localbkg_estimator)
+    psfphot(data, error=error)
+
+    psf_shape = (25, 25)
+    model1 = psfphot.make_model_image(data.shape, psf_shape,
+                                      include_localbkg=False)
+    model2 = psfphot.make_model_image(data.shape, psf_shape,
+                                      include_localbkg=True)
+    resid1 = psfphot.make_residual_image(data, psf_shape,
+                                         include_localbkg=False)
+    resid2 = psfphot.make_residual_image(data, psf_shape,
+                                         include_localbkg=True)
+
+    x, y = 0, 100
+    assert model1[y, x] < 0.1
+    assert model2[y, x] > 9
+    assert resid1[y, x] > 9
+    assert resid2[y, x] < 0
+
+    x, y = 0, 80
+    assert model1[y, x] < 0.1
+    assert model2[y, x] > 18
+    assert resid1[y, x] > 9
+    assert resid2[y, x] < -9
+
+
 @pytest.mark.skipif(not HAS_SCIPY, reason='scipy is required')
 def test_psf_photometry_compound(test_data):
     """