Coverage for python / lsst / pipe / tasks / skyCorrection.py: 17%
233 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-06 08:52 +0000
1# This file is part of pipe_tasks.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = ["SkyCorrectionTask", "SkyCorrectionConfig"]
26import warnings
28import numpy as np
30from lsst.afw.image import ExposureF, makeMaskedImage
31from lsst.afw.math import BackgroundMI, binImage
32from lsst.daf.butler import DeferredDatasetHandle
33from lsst.pex.config import Config, ConfigField, ConfigurableField, Field
34from lsst.pipe.base import (
35 AlgorithmError,
36 AnnotatedPartialOutputsError,
37 PipelineTask,
38 PipelineTaskConfig,
39 PipelineTaskConnections,
40 Struct,
41)
42from lsst.pipe.base.connectionTypes import Input, Output, PrerequisiteInput
43from lsst.pipe.tasks.background import (
44 FocalPlaneBackground,
45 FocalPlaneBackgroundConfig,
46 MaskObjectsTask,
47 SkyMeasurementTask,
48)
49from lsst.pipe.tasks.visualizeVisit import VisualizeMosaicExpConfig, VisualizeMosaicExpTask
52def _reorderAndPadList(inputList, inputKeys, outputKeys, padWith=None):
53 """Match the order of one list to another, padding if necessary.
55 Parameters
56 ----------
57 inputList : `list`
58 List to be reordered and padded. Elements can be any type.
59 inputKeys : iterable
60 Iterable of values to be compared with outputKeys.
61 Length must match `inputList`.
62 outputKeys : iterable
63 Iterable of values to be compared with inputKeys.
64 padWith :
65 Any value to be inserted where one of inputKeys is not in outputKeys.
67 Returns
68 -------
69 outputList : `list`
70 Copy of inputList reordered per outputKeys and padded with `padWith`
71 so that the length matches length of outputKeys.
72 """
73 outputList = []
74 for outputKey in outputKeys:
75 if outputKey in inputKeys:
76 outputList.append(inputList[inputKeys.index(outputKey)])
77 else:
78 outputList.append(padWith)
79 return outputList
82class SkyCorrectionConnections(PipelineTaskConnections, dimensions=("instrument", "visit")):
83 calExps = Input(
84 doc="Background-subtracted calibrated exposures.",
85 name="calexp",
86 multiple=True,
87 deferLoad=True,
88 storageClass="ExposureF",
89 dimensions=["instrument", "visit", "detector"],
90 )
91 calBkgs = Input(
92 doc="Subtracted backgrounds for input calibrated exposures.",
93 name="calexpBackground",
94 multiple=True,
95 storageClass="Background",
96 dimensions=["instrument", "visit", "detector"],
97 )
98 backgroundToPhotometricRatioHandles = Input(
99 doc="Ratio of a background-flattened image to a photometric-flattened image. "
100 "Only used if doApplyFlatBackgroundRatio is True.",
101 multiple=True,
102 name="background_to_photometric_ratio",
103 storageClass="Image",
104 dimensions=["instrument", "visit", "detector"],
105 deferLoad=True,
106 )
107 skyFrames = PrerequisiteInput(
108 doc="Calibration sky frames.",
109 name="sky",
110 multiple=True,
111 deferLoad=True,
112 storageClass="ExposureF",
113 dimensions=["instrument", "physical_filter", "detector"],
114 isCalibration=True,
115 )
116 camera = PrerequisiteInput(
117 doc="Input camera.",
118 name="camera",
119 storageClass="Camera",
120 dimensions=["instrument"],
121 isCalibration=True,
122 )
123 skyCorr = Output(
124 doc="Sky correction data, to be subtracted from the calibrated exposures.",
125 name="skyCorr",
126 multiple=True,
127 storageClass="Background",
128 dimensions=["instrument", "visit", "detector"],
129 )
130 calExpMosaic = Output(
131 doc="Full focal plane mosaicked image of the sky corrected calibrated exposures.",
132 name="calexp_skyCorr_visit_mosaic",
133 storageClass="ImageF",
134 dimensions=["instrument", "visit"],
135 )
136 calBkgMosaic = Output(
137 doc="Full focal plane mosaicked image of the sky corrected calibrated exposure backgrounds.",
138 name="calexpBackground_skyCorr_visit_mosaic",
139 storageClass="ImageF",
140 dimensions=["instrument", "visit"],
141 )
143 def __init__(self, *, config: "SkyCorrectionConfig | None" = None):
144 super().__init__(config=config)
145 assert config is not None
146 if not config.doSky:
147 del self.skyFrames
148 if not config.doApplyFlatBackgroundRatio:
149 del self.backgroundToPhotometricRatioHandles
152class SkyCorrectionConfig(PipelineTaskConfig, pipelineConnections=SkyCorrectionConnections):
153 doApplyFlatBackgroundRatio = Field(
154 dtype=bool,
155 default=False,
156 doc="This should be True if the input image was processed with an illumination correction.",
157 )
158 maskObjects = ConfigurableField(
159 target=MaskObjectsTask,
160 doc="Mask Objects",
161 )
162 doMaskObjects = Field(
163 dtype=bool,
164 default=True,
165 doc="Iteratively mask objects to find good sky?",
166 )
167 bgModel1 = ConfigField(
168 dtype=FocalPlaneBackgroundConfig,
169 doc="Initial background model, prior to sky frame subtraction",
170 )
171 undoBgModel1 = Field(
172 dtype=bool,
173 default=False,
174 doc="If True, adds back initial background model after sky and removes bgModel1 from the list",
175 )
176 sky = ConfigurableField(
177 target=SkyMeasurementTask,
178 doc="Sky measurement",
179 )
180 doSky = Field(
181 dtype=bool,
182 default=True,
183 doc="Do sky frame subtraction?",
184 )
185 bgModel2 = ConfigField(
186 dtype=FocalPlaneBackgroundConfig,
187 doc="Final (cleanup) background model, after sky frame subtraction",
188 )
189 doBgModel2 = Field(
190 dtype=bool,
191 default=True,
192 doc="Do final (cleanup) background model subtraction, after sky frame subtraction?",
193 )
194 binning = Field(
195 dtype=int,
196 default=8,
197 doc="Binning factor for constructing full focal plane '*_camera' output datasets",
198 )
200 def setDefaults(self):
201 Config.setDefaults(self)
202 self.bgModel2.doSmooth = True
203 self.bgModel2.minFrac = 0.3
204 self.bgModel2.xSize = 256
205 self.bgModel2.ySize = 256
206 self.bgModel2.smoothScale = 1.0
208 def validate(self):
209 super().validate()
210 if self.undoBgModel1 and not self.doSky and not self.doBgModel2:
211 raise ValueError("If undoBgModel1 is True, task requires at least one of doSky or doBgModel2.")
214class SkyCorrectionTask(PipelineTask):
215 """Perform a full focal plane sky correction."""
217 ConfigClass = SkyCorrectionConfig
218 _DefaultName = "skyCorr"
220 def __init__(self, *args, **kwargs):
221 super().__init__(**kwargs)
222 self.makeSubtask("sky")
223 self.makeSubtask("maskObjects")
225 def runQuantum(self, butlerQC, inputRefs, outputRefs):
226 # Sort input/output connections by detector ID, padding where
227 # necessary, to ensure that all detectors are processed consistently.
228 # Detector IDs are defined from the intersection of calExps, calBkgs,
229 # and optionally skyFrames and backgroundToPhotometricRatioHandles.
230 # This resolves potential missing data issues when processing a visit
231 # that contains only partial inputs.
232 calExpOrder = {ref.dataId["detector"] for ref in inputRefs.calExps}
233 calBkgOrder = {ref.dataId["detector"] for ref in inputRefs.calBkgs}
234 detectorOrder = calExpOrder & calBkgOrder
235 if self.config.doApplyFlatBackgroundRatio:
236 ratioOrder = {ref.dataId["detector"] for ref in inputRefs.backgroundToPhotometricRatioHandles}
237 detectorOrder &= ratioOrder
238 if self.config.doSky:
239 skyFrameOrder = {ref.dataId["detector"] for ref in inputRefs.skyFrames}
240 detectorOrder &= skyFrameOrder
241 detectorOrder = sorted(detectorOrder)
242 inputRefs.calExps = _reorderAndPadList(
243 inputRefs.calExps, [ref.dataId["detector"] for ref in inputRefs.calExps], detectorOrder
244 )
245 inputRefs.calBkgs = _reorderAndPadList(
246 inputRefs.calBkgs, [ref.dataId["detector"] for ref in inputRefs.calBkgs], detectorOrder
247 )
248 # Only attempt to fetch sky frames if they are going to be applied.
249 if self.config.doSky:
250 inputRefs.skyFrames = _reorderAndPadList(
251 inputRefs.skyFrames, [ref.dataId["detector"] for ref in inputRefs.skyFrames], detectorOrder
252 )
253 else:
254 inputRefs.skyFrames = []
255 # Only attempt to fetch flat ratios if they are going to be applied.
256 if self.config.doApplyFlatBackgroundRatio:
257 inputRefs.backgroundToPhotometricRatioHandles = _reorderAndPadList(
258 inputRefs.backgroundToPhotometricRatioHandles,
259 [ref.dataId["detector"] for ref in inputRefs.backgroundToPhotometricRatioHandles],
260 detectorOrder,
261 )
262 else:
263 inputRefs.backgroundToPhotometricRatioHandles = []
264 outputRefs.skyCorr = _reorderAndPadList(
265 outputRefs.skyCorr, [ref.dataId["detector"] for ref in outputRefs.skyCorr], detectorOrder
266 )
267 inputs = butlerQC.get(inputRefs)
268 try:
269 outputs = self.run(**inputs)
270 except AlgorithmError as e:
271 error = AnnotatedPartialOutputsError.annotate(
272 e,
273 self, # type: ignore
274 log=self.log,
275 )
276 raise error from e
277 butlerQC.put(outputs, outputRefs)
279 def run(self, calExps, calBkgs, skyFrames, camera, backgroundToPhotometricRatioHandles=[]):
280 """Perform sky correction on a visit.
282 The original visit-level background is first restored to the calibrated
283 exposure and the existing background model is inverted in-place. If
284 doMaskObjects is True, the mask map associated with this exposure will
285 be iteratively updated (over nIter loops) by re-estimating the
286 background each iteration and redetecting footprints.
288 An initial full focal plane sky subtraction (bgModel1) will take place
289 prior to scaling and subtracting the sky frame.
291 If doSky is True, the sky frame will be scaled to the flux in the input
292 visit.
294 If doBgModel2 is True, a final full focal plane sky subtraction will
295 take place after the sky frame has been subtracted.
297 The first N elements of the returned skyCorr will consist of inverted
298 elements of the calexpBackground model (i.e., subtractive). All
299 subsequent elements appended to skyCorr thereafter will be additive
300 such that, when skyCorr is subtracted from a calexp, the net result
301 will be to undo the initial per-detector background solution and then
302 apply the skyCorr model thereafter. Adding skyCorr to a
303 calexpBackground will effectively negate the calexpBackground,
304 returning only the additive background components of the skyCorr
305 background model.
307 Parameters
308 ----------
309 calExps : `list` [`lsst.afw.image.ExposureF`]
310 Detector calibrated exposure images for the visit.
311 calBkgs : `list` [`lsst.afw.math.BackgroundList`]
312 Detector background lists matching the calibrated exposures.
313 skyFrames : `list` [`lsst.afw.image.ExposureF`]
314 Sky frame calibration data for the input detectors.
315 camera : `lsst.afw.cameraGeom.Camera`
316 Camera matching the input data to process.
317 backgroundToPhotometricRatioHandles :
318 `list` [`lsst.daf.butler.DeferredDatasetHandle`], optional
319 Deferred dataset handles pointing to the Background to photometric
320 ratio images for the input detectors.
322 Returns
323 -------
324 results : `Struct` containing:
325 skyFrameScale : `float`
326 Scale factor applied to the sky frame.
327 skyCorr : `list` [`lsst.afw.math.BackgroundList`]
328 Detector-level sky correction background lists.
329 calExpMosaic : `lsst.afw.image.ExposureF`
330 Visit-level mosaic of the sky corrected data, binned.
331 Analogous to `calexp - skyCorr`.
332 calBkgMosaic : `lsst.afw.image.ExposureF`
333 Visit-level mosaic of the sky correction background, binned.
334 Analogous to `calexpBackground + skyCorr`.
335 """
336 # Process each detector separately and merge into bgModel1.
337 # This avoids storing every full-res image in-memory at once.
338 bgModel1 = FocalPlaneBackground.fromCamera(self.config.bgModel1, camera)
339 detectors = []
340 masks = []
341 skyCorrs = []
342 bgModel1Indices = []
343 if not self.config.doApplyFlatBackgroundRatio:
344 backgroundToPhotometricRatioHandles = [None] * len(calExps)
345 for calExpHandle, calBkg, backgroundToPhotometricRatioHandle in zip(
346 calExps, calBkgs, backgroundToPhotometricRatioHandles
347 ):
348 calExp = self._getCalExp(
349 calExpHandle, backgroundToPhotometricRatioHandle=backgroundToPhotometricRatioHandle
350 )
351 detectors.append(calExp.getDetector())
353 # Restore original background in-place; optionally refine mask maps
354 _ = self._restoreOriginalBackgroundRefineMask(calExp, calBkg)
355 masks.append(calExp.mask)
356 skyCorrs.append(calBkg) # Contains only the inverted original background elements at this stage
357 bgModel1Indices.append(len(calBkg)) # Index of the original background element
359 # Make a background model for the image, using bgModel1 configs
360 bgModel1Detector = FocalPlaneBackground.fromCamera(self.config.bgModel1, camera)
361 bgModel1Detector.addCcd(calExp)
362 bgModel1.merge(bgModel1Detector)
363 self.log.info(
364 "Detector %d: Merged %d unmasked pixels (%.1f%s of detector area) into initial BG model",
365 calExp.getDetector().getId(),
366 bgModel1Detector._numbers.getArray().sum(),
367 100 * bgModel1Detector._numbers.getArray().sum() / calExp.getBBox().getArea(),
368 "%",
369 )
371 # Validate bgModel1
372 self._validateBgModel("bgModel1", bgModel1, self.config.bgModel1)
374 # Update skyCorrs with new bgModel1 background elements
375 for detector, skyCorr in zip(detectors, skyCorrs):
376 with warnings.catch_warnings():
377 warnings.filterwarnings("ignore", "invalid value encountered")
378 calBkgElement = bgModel1.toCcdBackground(detector, detector.getBBox())
379 skyCorr.append(calBkgElement[0])
381 # Fit a scaled sky frame to all input exposures
382 skyFrameScale = None
383 if self.config.doSky:
384 skyFrameScale = self._fitSkyFrame(
385 calExps, masks, skyCorrs, skyFrames, backgroundToPhotometricRatioHandles
386 )
388 # Remove the initial background model (bgModel1) from every skyCorr
389 if self.config.undoBgModel1:
390 for skyCorr, bgModel1Index in zip(skyCorrs, bgModel1Indices):
391 skyCorr._backgrounds.pop(bgModel1Index)
392 self.log.info(
393 "Initial background models (bgModel1s) have been removed from all skyCorr background lists",
394 )
396 # Bin exposures, generate full-fp bg, map to CCDs and subtract in-place
397 if self.config.doBgModel2:
398 bgModel2 = FocalPlaneBackground.fromCamera(self.config.bgModel2, camera)
399 for calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle in zip(
400 calExps, masks, skyCorrs, backgroundToPhotometricRatioHandles
401 ):
402 calExp = self._getCalExp(calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle)
404 # Make a background model for the image, using bgModel2 configs
405 bgModel2Detector = FocalPlaneBackground.fromCamera(self.config.bgModel2, camera)
406 bgModel2Detector.addCcd(calExp)
407 bgModel2.merge(bgModel2Detector)
408 self.log.info(
409 "Detector %d: Merged %d unmasked pixels (%.1f%s of detector area) into final BG model",
410 calExp.getDetector().getId(),
411 bgModel2Detector._numbers.getArray().sum(),
412 100 * bgModel2Detector._numbers.getArray().sum() / calExp.getBBox().getArea(),
413 "%",
414 )
416 # Validate bgModel2
417 self._validateBgModel("bgModel2", bgModel2, self.config.bgModel2)
419 # Update skyCorrs with new bgModel2 background elements
420 for detector, skyCorr in zip(detectors, skyCorrs):
421 with warnings.catch_warnings():
422 warnings.filterwarnings("ignore", "invalid value encountered")
423 calBkgElement = bgModel2.toCcdBackground(detector, detector.getBBox())
424 skyCorr.append(calBkgElement[0])
426 # Make camera-level mosaics of bg subtracted calexps and subtracted bgs
427 calExpsBinned = []
428 calBkgsBinned = []
429 for calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle, bgModel1Index in zip(
430 calExps, masks, skyCorrs, backgroundToPhotometricRatioHandles, bgModel1Indices
431 ):
432 calExp = self._getCalExp(calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle)
434 skyCorrExtra = skyCorr.clone() # new skyCorr elements created in this task
435 skyCorrExtra._backgrounds = skyCorrExtra._backgrounds[bgModel1Index:]
436 skyCorrExtraMI = makeMaskedImage(skyCorrExtra.getImage())
437 skyCorrExtraMI.setMask(calExp.getMask())
439 calExpsBinned.append(binImage(calExp.getMaskedImage(), self.config.binning))
440 calBkgsBinned.append(binImage(skyCorrExtraMI, self.config.binning))
442 mosConfig = VisualizeMosaicExpConfig()
443 mosConfig.binning = self.config.binning
444 mosTask = VisualizeMosaicExpTask(config=mosConfig)
445 detectorIds = [detector.getId() for detector in detectors]
446 calExpMosaic = mosTask.run(calExpsBinned, camera, inputIds=detectorIds).outputData
447 calBkgMosaic = mosTask.run(calBkgsBinned, camera, inputIds=detectorIds).outputData
449 return Struct(
450 skyFrameScale=skyFrameScale,
451 skyCorr=skyCorrs,
452 calExpMosaic=calExpMosaic,
453 calBkgMosaic=calBkgMosaic,
454 )
456 def _getCalExp(self, calExpHandle, mask=None, skyCorr=None, backgroundToPhotometricRatioHandle=None):
457 """Get a calexp from a DeferredDatasetHandle, and optionally apply an
458 updated mask and skyCorr.
460 Parameters
461 ----------
462 calExpHandle : `~lsst.afw.image.ExposureF`
463 | `lsst.daf.butler.DeferredDatasetHandle`
464 Either the image exposure data or a handle to the calexp dataset.
465 mask : `lsst.afw.image.Mask`, optional
466 Mask to apply to the calexp.
467 skyCorr : `lsst.afw.math.BackgroundList`, optional
468 Background list to subtract from the calexp.
470 Returns
471 -------
472 calExp : `lsst.afw.image.ExposureF`
473 The calexp with the mask and skyCorr applied.
474 """
475 if isinstance(calExpHandle, DeferredDatasetHandle):
476 calExp: ExposureF = calExpHandle.get()
477 else:
478 # Here we clone the imaging data to avoid modifying data which is
479 # used in downstream processing.
480 calExp: ExposureF = calExpHandle.clone()
482 # Convert from background-flattened to photometric-flattened images
483 # Note: remember to convert back to background-flattened images
484 # if science images are to be output by this task.
485 if self.config.doApplyFlatBackgroundRatio:
486 if not backgroundToPhotometricRatioHandle:
487 raise ValueError(
488 "A list of backgroundToPhotometricRatioHandles must be supplied if "
489 "config.doApplyFlatBackgroundRatio=True.",
490 )
491 ratioImage = backgroundToPhotometricRatioHandle.get()
492 calExp.maskedImage *= ratioImage
493 self.log.info(
494 "Detector %d: Converted background-flattened image to a photometric-flattened image",
495 calExp.getDetector().getId(),
496 )
498 if mask is not None:
499 calExp.setMask(mask)
500 if skyCorr is not None:
501 image = calExp.getMaskedImage()
502 image -= skyCorr.getImage()
504 return calExp
506 def _getSkyFrame(self, skyFrameHandle):
507 """Get a calexp from a DeferredDatasetHandle, and optionally apply an
508 updated mask and skyCorr.
510 Parameters
511 ----------
512 skyFrameHandle : `lsst.daf.butler.DeferredDatasetHandle`
513 Either the sky frame data or a handle to the sky frame dataset.
515 Returns
516 -------
517 skyFrame : `lsst.afw.image.ExposureF`
518 The calexp with the mask and skyCorr applied.
519 """
520 if isinstance(skyFrameHandle, DeferredDatasetHandle):
521 skyFrame: ExposureF = skyFrameHandle.get()
522 else:
523 skyFrame: ExposureF = skyFrameHandle
524 return skyFrame
526 def _restoreOriginalBackgroundRefineMask(self, calExp, calBkg):
527 """Restore original background to a calexp and invert the related
528 background model; optionally refine the mask plane.
530 The original visit-level background is restored to the calibrated
531 exposure and the existing background model is inverted in-place. If
532 doMaskObjects is True, the mask map associated with the exposure will
533 be iteratively updated (over nIter loops) by re-estimating the
534 background each iteration and redetecting footprints.
536 The background model modified in-place in this method will comprise the
537 first N elements of the skyCorr dataset type, i.e., these N elements
538 are the inverse of the calexpBackground model. All subsequent elements
539 appended to skyCorr will be additive such that, when skyCorr is
540 subtracted from a calexp, the net result will be to undo the initial
541 per-detector background solution and then apply the skyCorr model
542 thereafter. Adding skyCorr to a calexpBackground will effectively
543 negate the calexpBackground, returning only the additive background
544 components of the skyCorr background model.
546 Parameters
547 ----------
548 calExp : `lsst.afw.image.ExposureF`
549 Detector level calexp image.
550 calBkg : `lsst.afw.math.BackgroundList`
551 Detector level background lists associated with the calexp.
553 Returns
554 -------
555 calExp : `lsst.afw.image.ExposureF`
556 The calexp with the originally subtracted background restored.
557 skyCorrBase : `lsst.afw.math.BackgroundList`
558 The inverted original background models; the genesis for skyCorr.
559 """
560 image = calExp.getMaskedImage()
562 # Invert all elements of the existing bg model; restore in calexp
563 for calBkgElement in calBkg:
564 statsImage = calBkgElement[0].getStatsImage()
565 statsImage *= -1
566 skyCorrBase = calBkg.getImage()
567 image -= skyCorrBase
569 stats = np.nanpercentile(skyCorrBase.array, [50, 75, 25])
570 self.log.info(
571 "Detector %d: Original background restored (BG median = %.1f counts, BG IQR = %.1f counts)",
572 calExp.getDetector().getId(),
573 -stats[0],
574 np.subtract(*stats[1:]),
575 )
577 # Iteratively subtract bg, re-detect sources, and add bg back on
578 if self.config.doMaskObjects:
579 maskFrac0 = 1 - np.sum(calExp.mask.array == 0) / calExp.mask.array.size
580 self.maskObjects.findObjects(calExp)
581 maskFrac1 = 1 - np.sum(calExp.mask.array == 0) / calExp.mask.array.size
583 self.log.info(
584 "Detector %d: Iterative source detection and mask growth has increased masked area by %.1f%%",
585 calExp.getDetector().getId(),
586 (100 * (maskFrac1 - maskFrac0)),
587 )
589 return calExp, skyCorrBase
591 def _validateBgModel(self, bgModelID, bgModel, config):
592 """Check that the background model contains enough valid superpixels,
593 and raise a useful error if not.
595 Parameters
596 ----------
597 bgModelID : `str`
598 Identifier for the background model.
599 bgModel : `~lsst.pipe.tasks.background.FocalPlaneBackground`
600 Background model to check.
601 config : `~lsst.pipe.tasks.background.FocalPlaneBackgroundConfig`
602 Configuration used to create the background model.
603 """
604 bgModelArray = bgModel._numbers.getArray()
605 spArea = (config.xSize / config.pixelSize) * (config.ySize / config.pixelSize)
606 self.log.info(
607 "%s: FP background model constructed using %.2f x %.2f mm (%d x %d pixel) superpixels",
608 bgModelID,
609 config.xSize,
610 config.ySize,
611 int(config.xSize / config.pixelSize),
612 int(config.ySize / config.pixelSize),
613 )
614 self.log.info(
615 "%s: Pixel data exists in %d of %d superpixels; the most populated superpixel is %.1f%% filled",
616 bgModelID,
617 np.sum(bgModelArray > 0),
618 bgModelArray.size,
619 100 * np.max(bgModelArray) / spArea,
620 )
622 thresh = config.minFrac * spArea
623 if np.all(bgModelArray < thresh):
624 raise NoUsableBackgroundSuperpixelsError(bgModelID, config.minFrac)
626 with warnings.catch_warnings():
627 warnings.filterwarnings("ignore", r"invalid value encountered")
628 stats = np.nanpercentile(bgModel.getStatsImage().array, [50, 75, 25])
629 self.log.info(
630 "%s: FP BG median = %.1f counts, FP BG IQR = %.1f counts",
631 bgModelID,
632 stats[0],
633 np.subtract(*stats[1:]),
634 )
636 def _fitSkyFrame(self, calExps, masks, skyCorrs, skyFrames, backgroundToPhotometricRatioHandles):
637 """Determine the full focal plane sky frame scale factor relative to
638 an input list of calibrated exposures.
640 This method measures the sky frame scale on all inputs, resulting in
641 values equal to the background method solveScales().
643 Input skyCorrs are updated in-place.
645 Parameters
646 ----------
647 calExps : `list` [`lsst.afw.image.ExposureF`]
648 Calibrated exposures to be background subtracted.
649 masks : `list` [`lsst.afw.image.Mask`]
650 Masks associated with the input calibrated exposures.
651 skyCorrs : `list` [`lsst.afw.math.BackgroundList`]
652 Background lists associated with the input calibrated exposures.
653 skyFrames : `list` [`lsst.afw.image.ExposureF`]
654 Sky frame calibration data for the input detectors.
656 Returns
657 -------
658 scale : `float`
659 Fitted scale factor applied to the sky frame.
660 """
661 skyBkgs = []
662 scales = []
663 for calExpHandle, mask, skyCorr, skyFrameHandle, backgroundToPhotometricRatioHandle in zip(
664 calExps, masks, skyCorrs, skyFrames, backgroundToPhotometricRatioHandles
665 ):
666 calExp = self._getCalExp(calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle)
667 skyFrame = self._getSkyFrame(skyFrameHandle)
668 skyBkg = self.sky.exposureToBackground(skyFrame)
669 del skyFrame # Free up memory
670 skyBkgs.append(skyBkg)
671 # Return a tuple of gridded image and sky frame clipped means
672 samples = self.sky.measureScale(calExp.getMaskedImage(), skyBkg)
673 scales.append(samples)
674 scale = self.sky.solveScales(scales)
675 for skyCorr, skyBkg in zip(skyCorrs, skyBkgs):
676 bgData = list(skyBkg[0])
677 bg = bgData[0]
678 statsImage = bg.getStatsImage().clone()
679 statsImage *= scale
680 newBg = BackgroundMI(bg.getImageBBox(), statsImage)
681 newBgData = [newBg] + bgData[1:]
682 skyCorr.append(newBgData)
683 self.log.info("Sky frame subtracted with a scale factor of %.5f", scale)
684 return scale
687class NoUsableBackgroundSuperpixelsError(AlgorithmError):
688 def __init__(self, bgModelID: str, minFrac: float):
689 self.bgModelID = bgModelID
690 self.minFrac = float(minFrac)
691 message = (
692 f"No background model superpixels in {bgModelID} are more than {100 * self.minFrac:.1f}% filled. "
693 "Try decreasing the minFrac configuration parameter, optimizing the subset of detectors "
694 "being processed, or increasing the number of detectors being processed."
695 )
696 super().__init__(message)
698 @property
699 def metadata(self) -> dict:
700 return {
701 "bgModelID": self.bgModelID,
702 }