diff --git a/python/lsst/pipe/tasks/skyCorrection.py b/python/lsst/pipe/tasks/skyCorrection.py index 6292dd23a..cabe67d43 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` @@ -318,8 +321,9 @@ def run(self, calExps, calBkgs, skyFrames, camera): initialBackgroundIndex = len(calBkgs[0]._backgrounds) - 1 # 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: @@ -341,7 +345,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 @@ -562,6 +568,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 = [] @@ -577,6 +588,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..cb81e1228 --- /dev/null +++ b/tests/test_skyCorrection.py @@ -0,0 +1,154 @@ +# 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 . + +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.doMaskObjects = False + # Set bgModel1 size to a single bin for the whole plane (aka constant) + 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 + + # Generate calexp/calexpBackground/sky for the central 6 HSC detectors + self.calExps = [] + self.calBkgs = [] + self.skyFrames = [] + self.background_level = 3000 + self.sky_level = 5 + for detector in [42, 50]: + rng = np.random.default_rng(detector) + + # Science image + calexp = ExposureF(bbox) + calexp.maskedImage.set(0.0, 0x0, 650.0) + calexp.setDetector(self.camera[detector]) + # Add a sky frame signature to a subregion of the image + calexp.image.array[:, (16 * 32) : (32 * 32)] += self.sky_level + calexp.image.array[(32 * 32) : (64 * 32)] += self.sky_level + # Add random noise + calexp.image.array += rng.normal(0.0, 25.0, (bbox.getDimensions().y, bbox.getDimensions().x)) + 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 += self.background_level + background = BackgroundMI(bbox, backgroundImage.getMaskedImage()) + calexpBackground = lsst.afw.math.BackgroundList( + ( + background, + lsst.afw.math.Interpolate.CONSTANT, + lsst.afw.math.UndersampleStyle.REDUCE_INTERP_ORDER, + lsst.afw.math.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) + 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:64, :] += 1 + # Add random noise + sky.image.array += rng.normal(0.0, 0.1, (131, 64)) + sky.image.array -= np.sum(sky.image.array) / (131 * 64) + self.skyFrames.append(sky) + + def tearDown(self): + del self.camera + del self.calExps + del self.calBkgs + del self.skyFrames + + def testSkyCorrectionDefault(self): + """Test SkyCorrectionTask with mostly default configuration values.""" + + 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, self.sky_level, delta=2e-1) + self.assertAlmostEqual(np.nanmean(results.calExpMosaic.array), 0, delta=1e-2) + + def testSkyCorrectionSkyFrameOnly(self): + """Test SkyCorrectionTask with the config undoBgModel1 set to True.""" + + skyCorrectionConfig = deepcopy(self.skyCorrectionConfig) + skyCorrectionConfig.undoBgModel1 = True + skyCorrectionConfig.doBgModel2 = False + skyCorrectionTask = SkyCorrectionTask(config=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 + ) + self.assertAlmostEqual( + np.nanmean(results.calExpMosaic.array), + np.nanmean(self.calExps[0].image.array) + self.background_level, + delta=1e-2, + ) + + +class MemoryTester(lsst.utils.tests.MemoryTestCase): + pass + + +def setup_module(module): + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main()