diff --git a/python/lsst/pipe/tasks/background.py b/python/lsst/pipe/tasks/background.py
index 464e74132..060611843 100644
--- a/python/lsst/pipe/tasks/background.py
+++ b/python/lsst/pipe/tasks/background.py
@@ -30,22 +30,21 @@
     "SkyStatsConfig",
 ]
 
-import sys
-import numpy
 import importlib
 import itertools
-from scipy.ndimage import gaussian_filter
+import sys
 
-import lsst.afw.math as afwMath
-import lsst.afw.image as afwImage
-import lsst.afw.geom as afwGeom
 import lsst.afw.cameraGeom as afwCameraGeom
+import lsst.afw.geom as afwGeom
+import lsst.afw.image as afwImage
+import lsst.afw.math as afwMath
+import lsst.afw.table as afwTable
 import lsst.geom as geom
 import lsst.meas.algorithms as measAlg
-import lsst.afw.table as afwTable
-
-from lsst.pex.config import Config, Field, ListField, ChoiceField, ConfigField, RangeField, ConfigurableField
+import numpy
+from lsst.pex.config import ChoiceField, Config, ConfigField, ConfigurableField, Field, ListField, RangeField
 from lsst.pipe.base import Task
+from scipy.ndimage import gaussian_filter
 
 
 def robustMean(array, rej=3.0):
@@ -64,47 +63,62 @@ def robustMean(array, rej=3.0):
         Robust mean of `array`.
     """
     q1, median, q3 = numpy.percentile(array, [25.0, 50.0, 100.0])
-    good = numpy.abs(array - median) < rej*0.74*(q3 - q1)
+    good = numpy.abs(array - median) < rej * 0.74 * (q3 - q1)
     return array[good].mean()
 
 
 class BackgroundConfig(Config):
     """Configuration for background measurement"""
-    statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points",
-                            allowed={"MEANCLIP": "clipped mean",
-                                     "MEAN": "unclipped mean",
-                                     "MEDIAN": "median"})
+
+    statistic = ChoiceField(
+        dtype=str,
+        default="MEANCLIP",
+        doc="type of statistic to use for grid points",
+        allowed={"MEANCLIP": "clipped mean", "MEAN": "unclipped mean", "MEDIAN": "median"},
+    )
     xBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in x")
     yBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in y")
-    algorithm = ChoiceField(dtype=str, default="NATURAL_SPLINE", optional=True,
-                            doc="How to interpolate the background values. "
-                                "This maps to an enum; see afw::math::Background",
-                            allowed={
-                                "CONSTANT": "Use a single constant value",
-                                "LINEAR": "Use linear interpolation",
-                                "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
-                                "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust"
-                                                " to outliers",
-                                "NONE": "No background estimation is to be attempted",
-                            })
-    mask = ListField(dtype=str, default=["SAT", "BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA"],
-                     doc="Names of mask planes to ignore while estimating the background")
+    algorithm = ChoiceField(
+        dtype=str,
+        default="NATURAL_SPLINE",
+        optional=True,
+        doc="How to interpolate the background values. This maps to an enum; see afw::math::Background",
+        allowed={
+            "CONSTANT": "Use a single constant value",
+            "LINEAR": "Use linear interpolation",
+            "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
+            "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
+            "NONE": "No background estimation is to be attempted",
+        },
+    )
+    mask = ListField(
+        dtype=str,
+        default=["SAT", "BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA"],
+        doc="Names of mask planes to ignore while estimating the background",
+    )
 
 
 class SkyStatsConfig(Config):
     """Parameters controlling the measurement of sky statistics"""
-    statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points",
-                            allowed={"MEANCLIP": "clipped mean",
-                                     "MEAN": "unclipped mean",
-                                     "MEDIAN": "median"})
+
+    statistic = ChoiceField(
+        dtype=str,
+        default="MEANCLIP",
+        doc="type of statistic to use for grid points",
+        allowed={"MEANCLIP": "clipped mean", "MEAN": "unclipped mean", "MEDIAN": "median"},
+    )
     clip = Field(doc="Clipping threshold for background", dtype=float, default=3.0)
     nIter = Field(doc="Clipping iterations for background", dtype=int, default=3)
-    mask = ListField(doc="Mask planes to reject", dtype=str,
-                     default=["SAT", "DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA"])
+    mask = ListField(
+        doc="Mask planes to reject",
+        dtype=str,
+        default=["SAT", "DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA"],
+    )
 
 
 class SkyMeasurementConfig(Config):
     """Configuration for SkyMeasurementTask"""
+
     skyIter = Field(dtype=int, default=3, doc="k-sigma rejection iterations for sky scale")
     skyRej = Field(dtype=float, default=3.0, doc="k-sigma rejection threshold for sky scale")
     background = ConfigField(dtype=BackgroundConfig, doc="Background measurement")
@@ -122,6 +136,7 @@ class SkyMeasurementTask(Task):
     background model (a `lsst.afw.math.BackgroundMI`).  The sky frame represents
     the dominant response of the camera to the sky background.
     """
+
     ConfigClass = SkyMeasurementConfig
 
     @staticmethod
@@ -149,11 +164,16 @@ def exposureToBackground(bgExp):
         algorithm = header.getScalar("ALGORITHM")
         bbox = geom.Box2I(geom.Point2I(xMin, yMin), geom.Point2I(xMax, yMax))
         return afwMath.BackgroundList(
-            (afwMath.BackgroundMI(bbox, bgExp.getMaskedImage()),
-             afwMath.stringToInterpStyle(algorithm),
-             afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
-             afwMath.ApproximateControl.UNKNOWN,
-             0, 0, False))
+            (
+                afwMath.BackgroundMI(bbox, bgExp.getMaskedImage()),
+                afwMath.stringToInterpStyle(algorithm),
+                afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
+                afwMath.ApproximateControl.UNKNOWN,
+                0,
+                0,
+                False,
+            )
+        )
 
     def backgroundToExposure(self, statsImage, bbox):
         """Convert a background model to an exposure
@@ -207,22 +227,26 @@ def measureBackground(self, image):
         stats.setNanSafe(True)
         ctrl = afwMath.BackgroundControl(
             self.config.background.algorithm,
-            max(int(image.getWidth()/self.config.background.xBinSize + 0.5), 1),
-            max(int(image.getHeight()/self.config.background.yBinSize + 0.5), 1),
+            max(int(image.getWidth() / self.config.background.xBinSize + 0.5), 1),
+            max(int(image.getHeight() / self.config.background.yBinSize + 0.5), 1),
             "REDUCE_INTERP_ORDER",
             stats,
-            self.config.background.statistic
+            self.config.background.statistic,
         )
 
         bg = afwMath.makeBackground(image, ctrl)
 
-        return afwMath.BackgroundList((
-            bg,
-            afwMath.stringToInterpStyle(self.config.background.algorithm),
-            afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
-            afwMath.ApproximateControl.UNKNOWN,
-            0, 0, False
-        ))
+        return afwMath.BackgroundList(
+            (
+                bg,
+                afwMath.stringToInterpStyle(self.config.background.algorithm),
+                afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
+                afwMath.ApproximateControl.UNKNOWN,
+                0,
+                0,
+                False,
+            )
+        )
 
     def averageBackgrounds(self, bgList):
         """Average multiple background models
@@ -243,8 +267,9 @@ def averageBackgrounds(self, bgList):
         assert all(len(bg) == 1 for bg in bgList), "Mixed bgList: %s" % ([len(bg) for bg in bgList],)
         images = [bg[0][0].getStatsImage() for bg in bgList]
         boxes = [bg[0][0].getImageBBox() for bg in bgList]
-        assert len(set((box.getMinX(), box.getMinY(), box.getMaxX(), box.getMaxY()) for box in boxes)) == 1, \
-            "Bounding boxes not all equal"
+        assert (
+            len(set((box.getMinX(), box.getMinY(), box.getMaxX(), box.getMaxY()) for box in boxes)) == 1
+        ), "Bounding boxes not all equal"
         bbox = boxes.pop(0)
 
         # Ensure bad pixels are masked
@@ -343,18 +368,18 @@ def solveScales(self, scales):
         skySamples = numpy.array(skySamples)
 
         def solve(mask):
-            return afwMath.LeastSquares.fromDesignMatrix(skySamples[mask].reshape(mask.sum(), 1),
-                                                         imageSamples[mask],
-                                                         afwMath.LeastSquares.DIRECT_SVD).getSolution()
+            return afwMath.LeastSquares.fromDesignMatrix(
+                skySamples[mask].reshape(mask.sum(), 1), imageSamples[mask], afwMath.LeastSquares.DIRECT_SVD
+            ).getSolution()
 
         mask = numpy.isfinite(imageSamples) & numpy.isfinite(skySamples)
         for ii in range(self.config.skyIter):
             solution = solve(mask)
-            residuals = imageSamples - solution*skySamples
+            residuals = imageSamples - solution * skySamples
             lq, uq = numpy.percentile(residuals[mask], [25, 75])
-            stdev = 0.741*(uq - lq)  # Robust stdev from IQR
+            stdev = 0.741 * (uq - lq)  # Robust stdev from IQR
             with numpy.errstate(invalid="ignore"):  # suppress NAN warnings
-                bad = numpy.abs(residuals) > self.config.skyRej*stdev
+                bad = numpy.abs(residuals) > self.config.skyRej * stdev
             mask[bad] = False
 
         return solve(mask)
@@ -414,14 +439,15 @@ def interpolate1D(method, xSample, ySample, xInterp):
 
     """
     if len(xSample) == 0:
-        return numpy.ones_like(xInterp)*numpy.nan
+        return numpy.ones_like(xInterp) * numpy.nan
     try:
-        return afwMath.makeInterpolate(xSample.astype(float), ySample.astype(float),
-                                       method).interpolate(xInterp.astype(float))
+        return afwMath.makeInterpolate(xSample.astype(float), ySample.astype(float), method).interpolate(
+            xInterp.astype(float)
+        )
     except Exception:
         if method == afwMath.Interpolate.CONSTANT:
             # We've already tried the most basic interpolation and it failed
-            return numpy.ones_like(xInterp)*numpy.nan
+            return numpy.ones_like(xInterp) * numpy.nan
         newMethod = afwMath.lookupMaxInterpStyle(len(xSample))
         if newMethod == method:
             newMethod = afwMath.Interpolate.CONSTANT
@@ -453,15 +479,17 @@ def interpolateBadPixels(array, isBad, interpolationStyle):
     isGood = ~isBad
     for y in range(height):
         if numpy.any(isBad[y, :]) and numpy.any(isGood[y, :]):
-            array[y][isBad[y]] = interpolate1D(method, xIndices[isGood[y]], array[y][isGood[y]],
-                                               xIndices[isBad[y]])
+            array[y][isBad[y]] = interpolate1D(
+                method, xIndices[isGood[y]], array[y][isGood[y]], xIndices[isBad[y]]
+            )
 
     isBad = numpy.isnan(array)
     isGood = ~isBad
     for x in range(width):
         if numpy.any(isBad[:, x]) and numpy.any(isGood[:, x]):
-            array[:, x][isBad[:, x]] = interpolate1D(method, yIndices[isGood[:, x]],
-                                                     array[:, x][isGood[:, x]], yIndices[isBad[:, x]])
+            array[:, x][isBad[:, x]] = interpolate1D(
+                method, yIndices[isGood[:, x]], array[:, x][isGood[:, x]], yIndices[isBad[:, x]]
+            )
 
 
 class FocalPlaneBackgroundConfig(Config):
@@ -473,15 +501,21 @@ class FocalPlaneBackgroundConfig(Config):
     need to be revised according to each particular camera. For
     this reason, no defaults are set for those.
     """
+
     xSize = Field(dtype=float, doc="Bin size in x")
     ySize = Field(dtype=float, doc="Bin size in y")
     pixelSize = Field(dtype=float, default=1.0, doc="Pixel size in same units as xSize/ySize")
     minFrac = Field(dtype=float, default=0.1, doc="Minimum fraction of bin size for good measurement")
-    mask = ListField(dtype=str, doc="Mask planes to treat as bad",
-                     default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"])
+    mask = ListField(
+        dtype=str,
+        doc="Mask planes to treat as bad",
+        default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"],
+    )
     interpolation = ChoiceField(
         doc="how to interpolate the background values. This maps to an enum; see afw::math::Background",
-        dtype=str, default="AKIMA_SPLINE", optional=True,
+        dtype=str,
+        default="AKIMA_SPLINE",
+        optional=True,
         allowed={
             "CONSTANT": "Use a single constant value",
             "LINEAR": "Use linear interpolation",
@@ -520,6 +554,7 @@ class FocalPlaneBackground:
     Once you've built the background model, you can apply it to individual
     CCDs with the `toCcdBackground` method.
     """
+
     @classmethod
     def fromCamera(cls, config, camera):
         """Construct from a camera object
@@ -538,14 +573,17 @@ def fromCamera(cls, config, camera):
 
         width, height = cameraBox.getDimensions()
         # Offset so that we run from zero
-        offset = geom.Extent2D(cameraBox.getMin())*-1
+        offset = geom.Extent2D(cameraBox.getMin()) * -1
         # Add an extra pixel buffer on either side
-        dims = geom.Extent2I(int(numpy.ceil(width/config.xSize)) + 2,
-                             int(numpy.ceil(height/config.ySize)) + 2)
+        dims = geom.Extent2I(
+            int(numpy.ceil(width / config.xSize)) + 2, int(numpy.ceil(height / config.ySize)) + 2
+        )
         # Transform takes us from focal plane coordinates --> sample coordinates
-        transform = (geom.AffineTransform.makeTranslation(geom.Extent2D(1, 1))
-                     * geom.AffineTransform.makeScaling(1.0/config.xSize, 1.0/config.ySize)
-                     * geom.AffineTransform.makeTranslation(offset))
+        transform = (
+            geom.AffineTransform.makeTranslation(geom.Extent2D(1, 1))
+            * geom.AffineTransform.makeScaling(1.0 / config.xSize, 1.0 / config.ySize)
+            * geom.AffineTransform.makeTranslation(offset)
+        )
 
         return cls(config, dims, afwGeom.makeTransform(transform))
 
@@ -626,8 +664,9 @@ def addCcd(self, exposure):
             CCD exposure to measure
         """
         detector = exposure.getDetector()
-        transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS),
-                                                            detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE))
+        transform = detector.getTransformMap().getTransform(
+            detector.makeCameraSys(afwCameraGeom.PIXELS), detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE)
+        )
         image = exposure.getMaskedImage()
         maskVal = image.getMask().getPlaneBitMask(self.config.mask)
 
@@ -656,7 +695,7 @@ def addCcd(self, exposure):
             num = result.getValue(afwMath.NPOINT)
             if not numpy.isfinite(mean) or not numpy.isfinite(num):
                 continue
-            warped[xx, yy, afwImage.LOCAL] = mean*num
+            warped[xx, yy, afwImage.LOCAL] = mean * num
             warpedCounts[xx, yy, afwImage.LOCAL] = num
 
         self._values += warped
@@ -680,10 +719,12 @@ def toCcdBackground(self, detector, bbox):
         bg : `lsst.afw.math.BackgroundList`
             Background model for CCD.
         """
-        transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS),
-                                                            detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE))
-        binTransform = (geom.AffineTransform.makeScaling(self.config.binning)
-                        * geom.AffineTransform.makeTranslation(geom.Extent2D(0.5, 0.5)))
+        transform = detector.getTransformMap().getTransform(
+            detector.makeCameraSys(afwCameraGeom.PIXELS), detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE)
+        )
+        binTransform = geom.AffineTransform.makeScaling(
+            self.config.binning
+        ) * geom.AffineTransform.makeTranslation(geom.Extent2D(0.5, 0.5))
 
         # Binned image on CCD --> unbinned image on CCD --> focal plane --> binned focal plane
         toSample = afwGeom.makeTransform(binTransform).then(transform).then(self.transform)
@@ -692,7 +733,7 @@ def toCcdBackground(self, detector, bbox):
         fpNorm = afwImage.ImageF(focalPlane.getBBox())
         fpNorm.set(1.0)
 
-        image = afwImage.ImageF(bbox.getDimensions()//self.config.binning)
+        image = afwImage.ImageF(bbox.getDimensions() // self.config.binning)
         norm = afwImage.ImageF(image.getBBox())
         ctrl = afwMath.WarpingControl("bilinear")
         afwMath.warpImage(image, focalPlane, toSample.inverted(), ctrl)
@@ -705,11 +746,15 @@ def toCcdBackground(self, detector, bbox):
         image.getArray()[isBad] = image.getArray()[~isBad].mean()
 
         return afwMath.BackgroundList(
-            (afwMath.BackgroundMI(bbox, afwImage.makeMaskedImage(image, mask)),
-             afwMath.stringToInterpStyle(self.config.interpolation),
-             afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
-             afwMath.ApproximateControl.UNKNOWN,
-             0, 0, False)
+            (
+                afwMath.BackgroundMI(bbox, afwImage.makeMaskedImage(image, mask)),
+                afwMath.stringToInterpStyle(self.config.interpolation),
+                afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
+                afwMath.ApproximateControl.UNKNOWN,
+                0,
+                0,
+                False,
+            )
         )
 
     def merge(self, other):
@@ -730,8 +775,10 @@ def merge(self, other):
             The merged background model.
         """
         if (self.config.xSize, self.config.ySize) != (other.config.xSize, other.config.ySize):
-            raise RuntimeError("Size mismatch: %s vs %s" % ((self.config.xSize, self.config.ySize),
-                                                            (other.config.xSize, other.config.ySize)))
+            raise RuntimeError(
+                "Size mismatch: %s vs %s"
+                % ((self.config.xSize, self.config.ySize), (other.config.xSize, other.config.ySize))
+            )
         if self.dims != other.dims:
             raise RuntimeError("Dimensions mismatch: %s vs %s" % (self.dims, other.dims))
         self._values += other._values
@@ -760,8 +807,11 @@ def getStatsImage(self):
         """
         values = self._values.clone()
         values /= self._numbers
-        thresh = (self.config.minFrac
-                  * (self.config.xSize/self.config.pixelSize)*(self.config.ySize/self.config.pixelSize))
+        thresh = (
+            self.config.minFrac
+            * (self.config.xSize / self.config.pixelSize)
+            * (self.config.ySize / self.config.pixelSize)
+        )
         isBad = self._numbers.getArray() < thresh
         if self.config.doSmooth:
             array = values.getArray()
@@ -774,9 +824,11 @@ def getStatsImage(self):
 
 class MaskObjectsConfig(Config):
     """Configuration for MaskObjectsTask"""
+
     nIter = Field(dtype=int, default=3, doc="Number of iterations")
-    subtractBackground = ConfigurableField(target=measAlg.SubtractBackgroundTask,
-                                           doc="Background subtraction")
+    subtractBackground = ConfigurableField(
+        target=measAlg.SubtractBackgroundTask, doc="Background subtraction"
+    )
     detection = ConfigurableField(target=measAlg.SourceDetectionTask, doc="Source detection")
     detectSigma = Field(dtype=float, default=5.0, doc="Detection threshold (standard deviations)")
     doInterpolate = Field(dtype=bool, default=True, doc="Interpolate when removing objects?")
@@ -793,11 +845,15 @@ def setDefaults(self):
         self.interpolate.useApprox = False
 
     def validate(self):
-        if (self.detection.reEstimateBackground
-                or self.detection.doTempLocalBackground
-                or self.detection.doTempWideBackground):
-            raise RuntimeError("Incorrect settings for object masking: reEstimateBackground, "
-                               "doTempLocalBackground and doTempWideBackground must be False")
+        if (
+            self.detection.reEstimateBackground
+            or self.detection.doTempLocalBackground
+            or self.detection.doTempWideBackground
+        ):
+            raise RuntimeError(
+                "Incorrect settings for object masking: reEstimateBackground, "
+                "doTempLocalBackground and doTempWideBackground must be False"
+            )
 
 
 class MaskObjectsTask(Task):
@@ -811,6 +867,7 @@ class MaskObjectsTask(Task):
     We deliberately use the specified ``detectSigma`` instead of the PSF,
     in order to better pick up the faint wings of objects.
     """
+
     ConfigClass = MaskObjectsConfig
 
     def __init__(self, *args, **kwargs):
@@ -902,7 +959,7 @@ def smoothArray(array, bad, sigma):
     """
     convolved = gaussian_filter(numpy.where(bad, 0.0, array), sigma, mode="constant", cval=0.0)
     denominator = gaussian_filter(numpy.where(bad, 0.0, 1.0), sigma, mode="constant", cval=0.0)
-    return convolved/denominator
+    return convolved / denominator
 
 
 def _create_module_child(name):
diff --git a/python/lsst/pipe/tasks/skyCorrection.py b/python/lsst/pipe/tasks/skyCorrection.py
index d1d2a8fc4..2045f04ba 100644
--- a/python/lsst/pipe/tasks/skyCorrection.py
+++ b/python/lsst/pipe/tasks/skyCorrection.py
@@ -256,6 +256,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
         inputs = butlerQC.get(inputRefs)
         inputs.pop("rawLinker", None)
         outputs = self.run(**inputs)
+        outputs.pop("skyFrameScale", None)
         butlerQC.put(outputs, outputRefs)
 
     def run(self, calExps, calBkgs, skyFrames, camera):
@@ -300,6 +301,8 @@ def run(self, calExps, calBkgs, skyFrames, camera):
         Returns
         -------
         results : `Struct` containing:
+            skyFrameScale : `float`
+                Scale factor applied to the sky frame.
             skyCorr : `list` [`lsst.afw.math.BackgroundList`]
                 Detector-level sky correction background lists.
             calExpMosaic : `lsst.afw.image.ExposureF`
@@ -317,8 +320,9 @@ def run(self, calExps, calBkgs, skyFrames, camera):
         _ = self._subtractVisitBackground(calExps, calBkgs, camera, self.config.bgModel1)
 
         # Subtract a scaled sky frame from all input exposures
+        skyFrameScale = None
         if self.config.doSky:
-            self._subtractSkyFrame(calExps, skyFrames, calBkgs)
+            skyFrameScale = self._subtractSkyFrame(calExps, skyFrames, calBkgs)
 
         # Adds full-fp bg back onto exposures, removes it from list
         if self.config.undoBgModel1:
@@ -340,7 +344,9 @@ def run(self, calExps, calBkgs, skyFrames, camera):
             skyCorrExtras, camera, self.config.binning, ids=calExpIds, refExps=calExps
         )
 
-        return Struct(skyCorr=calBkgs, calExpMosaic=calExpMosaic, calBkgMosaic=calBkgMosaic)
+        return Struct(
+            skyFrameScale=skyFrameScale, skyCorr=calBkgs, calExpMosaic=calExpMosaic, calBkgMosaic=calBkgMosaic
+        )
 
     def _restoreOriginalBackgroundRefineMask(self, calExps, calBkgs):
         """Restore original background to each calexp and invert the related
@@ -559,6 +565,11 @@ def _subtractSkyFrame(self, calExps, skyFrames, calBkgs):
             Sky frame calibration data for the input detectors.
         calBkgs : `list` [`lsst.afw.math.BackgroundList`]
             Background lists associated with the input calibrated exposures.
+
+        Returns
+        -------
+        scale : `float`
+            Scale factor applied to the sky frame.
         """
         skyFrameBgModels = []
         scales = []
@@ -574,6 +585,7 @@ def _subtractSkyFrame(self, calExps, skyFrames, calBkgs):
             # also updating the calBkg list in-place
             self.sky.subtractSkyFrame(calExp.getMaskedImage(), skyFrameBgModel, scale, calBkg)
         self.log.info("Sky frame subtracted with a scale factor of %.5f", scale)
+        return scale
 
     def _binAndMosaic(self, exposures, camera, binning, ids=None, refExps=None):
         """Bin input exposures and mosaic across the entire focal plane.
diff --git a/tests/test_skyCorrection.py b/tests/test_skyCorrection.py
new file mode 100644
index 000000000..c6fc9118f
--- /dev/null
+++ b/tests/test_skyCorrection.py
@@ -0,0 +1,165 @@
+# This file is part of pipe_tasks.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import unittest
+from copy import deepcopy
+
+import lsst.utils.tests
+import numpy as np
+from lsst.afw.image import ExposureF
+from lsst.afw.math import BackgroundMI
+from lsst.obs.subaru import HyperSuprimeCam
+from lsst.pipe.tasks.skyCorrection import SkyCorrectionConfig, SkyCorrectionTask
+
+
+class SkyCorrectionTestCase(lsst.utils.tests.TestCase):
+
+    def setUp(self):
+        hsc = HyperSuprimeCam()
+        self.camera = hsc.getCamera()
+        bbox = self.camera[0].getBBox()
+
+        self.skyCorrectionConfig = SkyCorrectionConfig()
+        self.skyCorrectionConfig.bgModel1.xSize = 8192 * 0.015
+        self.skyCorrectionConfig.bgModel1.ySize = 8192 * 0.015
+        self.skyCorrectionConfig.bgModel1.pixelSize = 0.015
+        self.skyCorrectionConfig.bgModel2.xSize = 256 * 0.015
+        self.skyCorrectionConfig.bgModel2.ySize = 256 * 0.015
+        self.skyCorrectionConfig.bgModel2.pixelSize = 0.015
+
+        def _createGaussian(flux, fwhm, size):
+            """Create a 2D Gaussian image.
+
+            Parameters
+            ----------
+            flux : `float`
+                Total flux of the Gaussian.
+            fwhm : `float`
+                Full width at half maximum of the Gaussian.
+            size : `int`
+                Size of the square image in pixels.
+            """
+            sigma = fwhm / 2.355
+            x0, y0 = size // 2, size // 2  # Center of the Gaussian
+
+            # Create a 2D grid of (x, y) coordinates
+            y, x = np.mgrid[0:size, 0:size]
+            gaussian = np.exp(-((x - x0) ** 2 + (y - y0) ** 2) / (2 * sigma**2))
+
+            # Normalize to get the total flux of 100
+            gaussian *= flux / gaussian.sum()
+            return gaussian
+
+        # Generate calexp/calexpBackground/sky for the central 6 HSC detectors
+        self.calExps = []
+        self.calBkgs = []
+        self.skyFrames = []
+        sky_level = 10
+        for detector in [41, 42, 49, 50, 57, 58]:
+            rng = np.random.default_rng(detector)
+
+            # Science image
+            calexp = ExposureF(bbox)
+            calexp.maskedImage.set(0.0, 0x0, 1.0)
+            calexp.setDetector(self.camera[detector])
+            # Add a sky frame signature to a subregion of the image
+            calexp.image.array[:, (16 * 32) : (32 * 32)] += sky_level
+            calexp.image.array[(32 * 32) : (48 * 32)] += sky_level
+            # Add random noise
+            calexp.image.array += rng.normal(0.0, 1.0, (bbox.getDimensions().y, bbox.getDimensions().x))
+            # Add some fake Gaussian sources to the image
+            gaussian_size = 100
+            for _ in range(10):
+                x = round(rng.uniform(gaussian_size // 2, bbox.getDimensions().x - gaussian_size // 2))
+                y = round(rng.uniform(gaussian_size // 2, bbox.getDimensions().y - gaussian_size // 2))
+                calexp.image.array[
+                    y - gaussian_size // 2 : y + gaussian_size // 2,
+                    x - gaussian_size // 2 : x + gaussian_size // 2,
+                ] += _createGaussian(1e6, 10, gaussian_size)
+            self.calExps.append(calexp)
+
+            # Background image
+            backgroundImage = ExposureF(bbox)
+            backgroundImage.maskedImage.set(0.0, 0x0, 1.0)
+            backgroundImage.setDetector(self.camera[detector])
+            backgroundImage.image.array += 3000
+            background = BackgroundMI(bbox, backgroundImage.getMaskedImage())
+            calexpBackground = lsst.afw.math.BackgroundList(
+                (
+                    background,
+                    lsst.afw.math.Interpolate.AKIMA_SPLINE,
+                    lsst.afw.math.UndersampleStyle.REDUCE_INTERP_ORDER,
+                    lsst.afw.math.ApproximateControl.CHEBYSHEV,
+                    6,
+                    6,
+                    1,
+                )
+            )
+            self.calBkgs.append(calexpBackground)
+
+            # Sky frame
+            sky = ExposureF(64, 131)
+            sky.maskedImage.set(0.0, 0x0, 1.0)
+            sky.setDetector(self.camera[detector])
+            header = sky.getMetadata()
+            header.set("BOX.MINX", bbox.getMinX())
+            header.set("BOX.MINY", bbox.getMinY())
+            header.set("BOX.MAXX", bbox.getMaxX())
+            header.set("BOX.MAXY", bbox.getMaxY())
+            header.set("ALGORITHM", "NATURAL_SPLINE")
+            sky.image.array[:, 16:32] += 1
+            sky.image.array[32:48, :] += 1
+            # Add random noise
+            sky.image.array += rng.normal(0.0, 1.0, (131, 64))
+            # sky.image.array -= np.sum(sky.image.array) / sky.getBBox().getArea()
+            self.skyFrames.append(sky)
+
+    def tearDown(self):
+        del self.camera
+        del self.calExps
+        del self.calBkgs
+        del self.skyFrames
+
+    def testSkyCorrection(self):
+        """TEMP."""
+
+        skyCorrectionTask = SkyCorrectionTask(config=self.skyCorrectionConfig)
+        # Pass in deep copies, as the task modifies the input data
+        results = skyCorrectionTask.run(
+            deepcopy(self.calExps), deepcopy(self.calBkgs), self.skyFrames, self.camera
+        )
+        skyFrameScale = results.skyFrameScale
+        skyCorr = results.skyCorr
+        self.assertEqual(len(skyCorr), len(self.calExps))
+        self.assertAlmostEqual(skyFrameScale, 8.35748841, delta=1e-8)
+
+
+class MemoryTester(lsst.utils.tests.MemoryTestCase):
+    pass
+
+
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()