Coverage for python/lsst/pipe/tasks/skyCorrection.py: 17%
229 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 02:06 -0700
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 PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
35from lsst.pipe.base.connectionTypes import Input, Output, PrerequisiteInput
36from lsst.pipe.tasks.background import (
37 FocalPlaneBackground,
38 FocalPlaneBackgroundConfig,
39 MaskObjectsTask,
40 SkyMeasurementTask,
41)
42from lsst.pipe.tasks.visualizeVisit import VisualizeMosaicExpConfig, VisualizeMosaicExpTask
45def _skyFrameLookup(datasetType, registry, quantumDataId, collections):
46 """Lookup function to identify sky frames.
48 Parameters
49 ----------
50 datasetType : `lsst.daf.butler.DatasetType`
51 Dataset to lookup.
52 registry : `lsst.daf.butler.Registry`
53 Butler registry to query.
54 quantumDataId : `lsst.daf.butler.DataCoordinate`
55 Data id to transform to find sky frames.
56 The ``detector`` entry will be stripped.
57 collections : `lsst.daf.butler.CollectionSearch`
58 Collections to search through.
60 Returns
61 -------
62 results : `list` [`lsst.daf.butler.DatasetRef`]
63 List of datasets that will be used as sky calibration frames.
64 """
65 newDataId = quantumDataId.subset(registry.dimensions.conform(["instrument", "visit"]))
66 skyFrames = []
67 for dataId in registry.queryDataIds(["visit", "detector"], dataId=newDataId).expanded():
68 skyFrame = registry.findDataset(
69 datasetType, dataId, collections=collections, timespan=dataId.timespan
70 )
71 skyFrames.append(skyFrame)
72 return skyFrames
75def _reorderAndPadList(inputList, inputKeys, outputKeys, padWith=None):
76 """Match the order of one list to another, padding if necessary.
78 Parameters
79 ----------
80 inputList : `list`
81 List to be reordered and padded. Elements can be any type.
82 inputKeys : iterable
83 Iterable of values to be compared with outputKeys.
84 Length must match `inputList`.
85 outputKeys : iterable
86 Iterable of values to be compared with inputKeys.
87 padWith :
88 Any value to be inserted where one of inputKeys is not in outputKeys.
90 Returns
91 -------
92 outputList : `list`
93 Copy of inputList reordered per outputKeys and padded with `padWith`
94 so that the length matches length of outputKeys.
95 """
96 outputList = []
97 for outputKey in outputKeys:
98 if outputKey in inputKeys:
99 outputList.append(inputList[inputKeys.index(outputKey)])
100 else:
101 outputList.append(padWith)
102 return outputList
105class SkyCorrectionConnections(PipelineTaskConnections, dimensions=("instrument", "visit")):
106 rawLinker = Input(
107 doc="Raw data to provide exp-visit linkage to connect calExp inputs to camera/sky calibs.",
108 name="raw",
109 multiple=True,
110 deferLoad=True,
111 storageClass="Exposure",
112 dimensions=["instrument", "exposure", "detector"],
113 )
114 calExps = Input(
115 doc="Background-subtracted calibrated exposures.",
116 name="calexp",
117 multiple=True,
118 deferLoad=True,
119 storageClass="ExposureF",
120 dimensions=["instrument", "visit", "detector"],
121 )
122 calBkgs = Input(
123 doc="Subtracted backgrounds for input calibrated exposures.",
124 name="calexpBackground",
125 multiple=True,
126 storageClass="Background",
127 dimensions=["instrument", "visit", "detector"],
128 )
129 backgroundToPhotometricRatioHandles = Input(
130 doc="Ratio of a background-flattened image to a photometric-flattened image. "
131 "Only used if doApplyFlatBackgroundRatio is True.",
132 multiple=True,
133 name="background_to_photometric_ratio",
134 storageClass="Image",
135 dimensions=["instrument", "visit", "detector"],
136 deferLoad=True,
137 )
138 skyFrames = PrerequisiteInput(
139 doc="Calibration sky frames.",
140 name="sky",
141 multiple=True,
142 deferLoad=True,
143 storageClass="ExposureF",
144 dimensions=["instrument", "physical_filter", "detector"],
145 isCalibration=True,
146 lookupFunction=_skyFrameLookup,
147 )
148 camera = PrerequisiteInput(
149 doc="Input camera.",
150 name="camera",
151 storageClass="Camera",
152 dimensions=["instrument"],
153 isCalibration=True,
154 )
155 skyCorr = Output(
156 doc="Sky correction data, to be subtracted from the calibrated exposures.",
157 name="skyCorr",
158 multiple=True,
159 storageClass="Background",
160 dimensions=["instrument", "visit", "detector"],
161 )
162 calExpMosaic = Output(
163 doc="Full focal plane mosaicked image of the sky corrected calibrated exposures.",
164 name="calexp_skyCorr_visit_mosaic",
165 storageClass="ImageF",
166 dimensions=["instrument", "visit"],
167 )
168 calBkgMosaic = Output(
169 doc="Full focal plane mosaicked image of the sky corrected calibrated exposure backgrounds.",
170 name="calexpBackground_skyCorr_visit_mosaic",
171 storageClass="ImageF",
172 dimensions=["instrument", "visit"],
173 )
175 def __init__(self, *, config: "SkyCorrectionConfig | None" = None):
176 super().__init__(config=config)
177 assert config is not None
178 if not config.doSky:
179 del self.skyFrames
180 if not config.doApplyFlatBackgroundRatio:
181 del self.backgroundToPhotometricRatioHandles
184class SkyCorrectionConfig(PipelineTaskConfig, pipelineConnections=SkyCorrectionConnections):
185 doApplyFlatBackgroundRatio = Field(
186 dtype=bool,
187 default=False,
188 doc="This should be True if the input image was processed with an illumination correction.",
189 )
190 maskObjects = ConfigurableField(
191 target=MaskObjectsTask,
192 doc="Mask Objects",
193 )
194 doMaskObjects = Field(
195 dtype=bool,
196 default=True,
197 doc="Iteratively mask objects to find good sky?",
198 )
199 bgModel1 = ConfigField(
200 dtype=FocalPlaneBackgroundConfig,
201 doc="Initial background model, prior to sky frame subtraction",
202 )
203 undoBgModel1 = Field(
204 dtype=bool,
205 default=False,
206 doc="If True, adds back initial background model after sky and removes bgModel1 from the list",
207 )
208 sky = ConfigurableField(
209 target=SkyMeasurementTask,
210 doc="Sky measurement",
211 )
212 doSky = Field(
213 dtype=bool,
214 default=True,
215 doc="Do sky frame subtraction?",
216 )
217 bgModel2 = ConfigField(
218 dtype=FocalPlaneBackgroundConfig,
219 doc="Final (cleanup) background model, after sky frame subtraction",
220 )
221 doBgModel2 = Field(
222 dtype=bool,
223 default=True,
224 doc="Do final (cleanup) background model subtraction, after sky frame subtraction?",
225 )
226 binning = Field(
227 dtype=int,
228 default=8,
229 doc="Binning factor for constructing full focal plane '*_camera' output datasets",
230 )
232 def setDefaults(self):
233 Config.setDefaults(self)
234 self.bgModel2.doSmooth = True
235 self.bgModel2.minFrac = 0.3
236 self.bgModel2.xSize = 256
237 self.bgModel2.ySize = 256
238 self.bgModel2.smoothScale = 1.0
240 def validate(self):
241 super().validate()
242 if self.undoBgModel1 and not self.doSky and not self.doBgModel2:
243 raise ValueError("If undoBgModel1 is True, task requires at least one of doSky or doBgModel2.")
246class SkyCorrectionTask(PipelineTask):
247 """Perform a full focal plane sky correction."""
249 ConfigClass = SkyCorrectionConfig
250 _DefaultName = "skyCorr"
252 def __init__(self, *args, **kwargs):
253 super().__init__(**kwargs)
254 self.makeSubtask("sky")
255 self.makeSubtask("maskObjects")
257 def runQuantum(self, butlerQC, inputRefs, outputRefs):
258 # Sort input/output connections by detector ID, padding where
259 # necessary, to ensure that all detectors are processed consistently.
260 # Detector IDs are defined from the intersection of calExps, calBkgs,
261 # and optionally skyFrames and backgroundToPhotometricRatioHandles.
262 # This resolves potential missing data issues when processing a visit
263 # that contains only partial inputs.
264 calExpOrder = {ref.dataId["detector"] for ref in inputRefs.calExps}
265 calBkgOrder = {ref.dataId["detector"] for ref in inputRefs.calBkgs}
266 detectorOrder = calExpOrder & calBkgOrder
267 if self.config.doApplyFlatBackgroundRatio:
268 ratioOrder = {ref.dataId["detector"] for ref in inputRefs.backgroundToPhotometricRatioHandles}
269 detectorOrder &= ratioOrder
270 if self.config.doSky:
271 skyFrameOrder = {ref.dataId["detector"] for ref in inputRefs.skyFrames}
272 detectorOrder &= skyFrameOrder
273 detectorOrder = sorted(detectorOrder)
274 inputRefs.calExps = _reorderAndPadList(
275 inputRefs.calExps, [ref.dataId["detector"] for ref in inputRefs.calExps], detectorOrder
276 )
277 inputRefs.calBkgs = _reorderAndPadList(
278 inputRefs.calBkgs, [ref.dataId["detector"] for ref in inputRefs.calBkgs], detectorOrder
279 )
280 # Only attempt to fetch sky frames if they are going to be applied.
281 if self.config.doSky:
282 inputRefs.skyFrames = _reorderAndPadList(
283 inputRefs.skyFrames, [ref.dataId["detector"] for ref in inputRefs.skyFrames], detectorOrder
284 )
285 else:
286 inputRefs.skyFrames = []
287 # Only attempt to fetch flat ratios if they are going to be applied.
288 if self.config.doApplyFlatBackgroundRatio:
289 inputRefs.backgroundToPhotometricRatioHandles = _reorderAndPadList(
290 inputRefs.backgroundToPhotometricRatioHandles,
291 [ref.dataId["detector"] for ref in inputRefs.backgroundToPhotometricRatioHandles],
292 detectorOrder,
293 )
294 else:
295 inputRefs.backgroundToPhotometricRatioHandles = []
296 outputRefs.skyCorr = _reorderAndPadList(
297 outputRefs.skyCorr, [ref.dataId["detector"] for ref in outputRefs.skyCorr], detectorOrder
298 )
299 inputs = butlerQC.get(inputRefs)
300 inputs.pop("rawLinker", None)
301 outputs = self.run(**inputs)
302 butlerQC.put(outputs, outputRefs)
304 def run(self, calExps, calBkgs, skyFrames, camera, backgroundToPhotometricRatioHandles=[]):
305 """Perform sky correction on a visit.
307 The original visit-level background is first restored to the calibrated
308 exposure and the existing background model is inverted in-place. If
309 doMaskObjects is True, the mask map associated with this exposure will
310 be iteratively updated (over nIter loops) by re-estimating the
311 background each iteration and redetecting footprints.
313 An initial full focal plane sky subtraction (bgModel1) will take place
314 prior to scaling and subtracting the sky frame.
316 If doSky is True, the sky frame will be scaled to the flux in the input
317 visit.
319 If doBgModel2 is True, a final full focal plane sky subtraction will
320 take place after the sky frame has been subtracted.
322 The first N elements of the returned skyCorr will consist of inverted
323 elements of the calexpBackground model (i.e., subtractive). All
324 subsequent elements appended to skyCorr thereafter will be additive
325 such that, when skyCorr is subtracted from a calexp, the net result
326 will be to undo the initial per-detector background solution and then
327 apply the skyCorr model thereafter. Adding skyCorr to a
328 calexpBackground will effectively negate the calexpBackground,
329 returning only the additive background components of the skyCorr
330 background model.
332 Parameters
333 ----------
334 calExps : `list` [`lsst.afw.image.ExposureF`]
335 Detector calibrated exposure images for the visit.
336 calBkgs : `list` [`lsst.afw.math.BackgroundList`]
337 Detector background lists matching the calibrated exposures.
338 skyFrames : `list` [`lsst.afw.image.ExposureF`]
339 Sky frame calibration data for the input detectors.
340 camera : `lsst.afw.cameraGeom.Camera`
341 Camera matching the input data to process.
342 backgroundToPhotometricRatioHandles :
343 `list` [`lsst.daf.butler.DeferredDatasetHandle`], optional
344 Deferred dataset handles pointing to the Background to photometric
345 ratio images for the input detectors.
347 Returns
348 -------
349 results : `Struct` containing:
350 skyFrameScale : `float`
351 Scale factor applied to the sky frame.
352 skyCorr : `list` [`lsst.afw.math.BackgroundList`]
353 Detector-level sky correction background lists.
354 calExpMosaic : `lsst.afw.image.ExposureF`
355 Visit-level mosaic of the sky corrected data, binned.
356 Analogous to `calexp - skyCorr`.
357 calBkgMosaic : `lsst.afw.image.ExposureF`
358 Visit-level mosaic of the sky correction background, binned.
359 Analogous to `calexpBackground + skyCorr`.
360 """
361 # Process each detector separately and merge into bgModel1.
362 # This avoids storing every full-res image in-memory at once.
363 bgModel1 = FocalPlaneBackground.fromCamera(self.config.bgModel1, camera)
364 detectors = []
365 masks = []
366 skyCorrs = []
367 bgModel1Indices = []
368 if not self.config.doApplyFlatBackgroundRatio:
369 backgroundToPhotometricRatioHandles = [None] * len(calExps)
370 for calExpHandle, calBkg, backgroundToPhotometricRatioHandle in zip(
371 calExps, calBkgs, backgroundToPhotometricRatioHandles
372 ):
373 calExp = self._getCalExp(
374 calExpHandle, backgroundToPhotometricRatioHandle=backgroundToPhotometricRatioHandle
375 )
376 detectors.append(calExp.getDetector())
378 # Restore original background in-place; optionally refine mask maps
379 _ = self._restoreOriginalBackgroundRefineMask(calExp, calBkg)
380 masks.append(calExp.mask)
381 skyCorrs.append(calBkg) # Contains only the inverted original background elements at this stage
382 bgModel1Indices.append(len(calBkg)) # Index of the original background element
384 # Make a background model for the image, using bgModel1 configs
385 bgModel1Detector = FocalPlaneBackground.fromCamera(self.config.bgModel1, camera)
386 bgModel1Detector.addCcd(calExp)
387 bgModel1.merge(bgModel1Detector)
388 self.log.info(
389 "Detector %d: Merged %d unmasked pixels (%.1f%s of detector area) into initial BG model",
390 calExp.getDetector().getId(),
391 bgModel1Detector._numbers.getArray().sum(),
392 100 * bgModel1Detector._numbers.getArray().sum() / calExp.getBBox().getArea(),
393 "%",
394 )
396 # Validate bgModel1
397 self._validateBgModel("bgModel1", bgModel1, self.config.bgModel1)
399 # Update skyCorrs with new bgModel1 background elements
400 for detector, skyCorr in zip(detectors, skyCorrs):
401 with warnings.catch_warnings():
402 warnings.filterwarnings("ignore", "invalid value encountered")
403 calBkgElement = bgModel1.toCcdBackground(detector, detector.getBBox())
404 skyCorr.append(calBkgElement[0])
406 # Fit a scaled sky frame to all input exposures
407 skyFrameScale = None
408 if self.config.doSky:
409 skyFrameScale = self._fitSkyFrame(
410 calExps, masks, skyCorrs, skyFrames, backgroundToPhotometricRatioHandles
411 )
413 # Remove the initial background model (bgModel1) from every skyCorr
414 if self.config.undoBgModel1:
415 for skyCorr, bgModel1Index in zip(skyCorrs, bgModel1Indices):
416 skyCorr._backgrounds.pop(bgModel1Index)
417 self.log.info(
418 "Initial background models (bgModel1s) have been removed from all skyCorr background lists",
419 )
421 # Bin exposures, generate full-fp bg, map to CCDs and subtract in-place
422 if self.config.doBgModel2:
423 bgModel2 = FocalPlaneBackground.fromCamera(self.config.bgModel2, camera)
424 for calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle in zip(
425 calExps, masks, skyCorrs, backgroundToPhotometricRatioHandles
426 ):
427 calExp = self._getCalExp(calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle)
429 # Make a background model for the image, using bgModel2 configs
430 bgModel2Detector = FocalPlaneBackground.fromCamera(self.config.bgModel2, camera)
431 bgModel2Detector.addCcd(calExp)
432 bgModel2.merge(bgModel2Detector)
433 self.log.info(
434 "Detector %d: Merged %d unmasked pixels (%.1f%s of detector area) into final BG model",
435 calExp.getDetector().getId(),
436 bgModel2Detector._numbers.getArray().sum(),
437 100 * bgModel2Detector._numbers.getArray().sum() / calExp.getBBox().getArea(),
438 "%",
439 )
441 # Validate bgModel2
442 self._validateBgModel("bgModel2", bgModel2, self.config.bgModel2)
444 # Update skyCorrs with new bgModel2 background elements
445 for detector, skyCorr in zip(detectors, skyCorrs):
446 with warnings.catch_warnings():
447 warnings.filterwarnings("ignore", "invalid value encountered")
448 calBkgElement = bgModel2.toCcdBackground(detector, detector.getBBox())
449 skyCorr.append(calBkgElement[0])
451 # Make camera-level mosaics of bg subtracted calexps and subtracted bgs
452 calExpsBinned = []
453 calBkgsBinned = []
454 for calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle, bgModel1Index in zip(
455 calExps, masks, skyCorrs, backgroundToPhotometricRatioHandles, bgModel1Indices
456 ):
457 calExp = self._getCalExp(calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle)
459 skyCorrExtra = skyCorr.clone() # new skyCorr elements created in this task
460 skyCorrExtra._backgrounds = skyCorrExtra._backgrounds[bgModel1Index:]
461 skyCorrExtraMI = makeMaskedImage(skyCorrExtra.getImage())
462 skyCorrExtraMI.setMask(calExp.getMask())
464 calExpsBinned.append(binImage(calExp.getMaskedImage(), self.config.binning))
465 calBkgsBinned.append(binImage(skyCorrExtraMI, self.config.binning))
467 mosConfig = VisualizeMosaicExpConfig()
468 mosConfig.binning = self.config.binning
469 mosTask = VisualizeMosaicExpTask(config=mosConfig)
470 detectorIds = [detector.getId() for detector in detectors]
471 calExpMosaic = mosTask.run(calExpsBinned, camera, inputIds=detectorIds).outputData
472 calBkgMosaic = mosTask.run(calBkgsBinned, camera, inputIds=detectorIds).outputData
474 return Struct(
475 skyFrameScale=skyFrameScale,
476 skyCorr=skyCorrs,
477 calExpMosaic=calExpMosaic,
478 calBkgMosaic=calBkgMosaic,
479 )
481 def _getCalExp(self, calExpHandle, mask=None, skyCorr=None, backgroundToPhotometricRatioHandle=None):
482 """Get a calexp from a DeferredDatasetHandle, and optionally apply an
483 updated mask and skyCorr.
485 Parameters
486 ----------
487 calExpHandle : `~lsst.afw.image.ExposureF`
488 | `lsst.daf.butler.DeferredDatasetHandle`
489 Either the image exposure data or a handle to the calexp dataset.
490 mask : `lsst.afw.image.Mask`, optional
491 Mask to apply to the calexp.
492 skyCorr : `lsst.afw.math.BackgroundList`, optional
493 Background list to subtract from the calexp.
495 Returns
496 -------
497 calExp : `lsst.afw.image.ExposureF`
498 The calexp with the mask and skyCorr applied.
499 """
500 if isinstance(calExpHandle, DeferredDatasetHandle):
501 calExp: ExposureF = calExpHandle.get()
502 else:
503 # Here we clone the imaging data to avoid modifying data which is
504 # used in downstream processing.
505 calExp: ExposureF = calExpHandle.clone()
507 # Convert from background-flattened to photometric-flattened images
508 # Note: remember to convert back to background-flattened images
509 # if science images are to be output by this task.
510 if self.config.doApplyFlatBackgroundRatio:
511 if not backgroundToPhotometricRatioHandle:
512 raise ValueError(
513 "A list of backgroundToPhotometricRatioHandles must be supplied if "
514 "config.doApplyFlatBackgroundRatio=True.",
515 )
516 ratioImage = backgroundToPhotometricRatioHandle.get()
517 calExp.maskedImage *= ratioImage
518 self.log.info(
519 "Detector %d: Converted background-flattened image to a photometric-flattened image",
520 calExp.getDetector().getId(),
521 )
523 if mask is not None:
524 calExp.setMask(mask)
525 if skyCorr is not None:
526 image = calExp.getMaskedImage()
527 image -= skyCorr.getImage()
529 return calExp
531 def _getSkyFrame(self, skyFrameHandle):
532 """Get a calexp from a DeferredDatasetHandle, and optionally apply an
533 updated mask and skyCorr.
535 Parameters
536 ----------
537 skyFrameHandle : `lsst.daf.butler.DeferredDatasetHandle`
538 Either the sky frame data or a handle to the sky frame dataset.
540 Returns
541 -------
542 skyFrame : `lsst.afw.image.ExposureF`
543 The calexp with the mask and skyCorr applied.
544 """
545 if isinstance(skyFrameHandle, DeferredDatasetHandle):
546 skyFrame: ExposureF = skyFrameHandle.get()
547 else:
548 skyFrame: ExposureF = skyFrameHandle
549 return skyFrame
551 def _restoreOriginalBackgroundRefineMask(self, calExp, calBkg):
552 """Restore original background to a calexp and invert the related
553 background model; optionally refine the mask plane.
555 The original visit-level background is restored to the calibrated
556 exposure and the existing background model is inverted in-place. If
557 doMaskObjects is True, the mask map associated with the exposure will
558 be iteratively updated (over nIter loops) by re-estimating the
559 background each iteration and redetecting footprints.
561 The background model modified in-place in this method will comprise the
562 first N elements of the skyCorr dataset type, i.e., these N elements
563 are the inverse of the calexpBackground model. All subsequent elements
564 appended to skyCorr will be additive such that, when skyCorr is
565 subtracted from a calexp, the net result will be to undo the initial
566 per-detector background solution and then apply the skyCorr model
567 thereafter. Adding skyCorr to a calexpBackground will effectively
568 negate the calexpBackground, returning only the additive background
569 components of the skyCorr background model.
571 Parameters
572 ----------
573 calExp : `lsst.afw.image.ExposureF`
574 Detector level calexp image.
575 calBkg : `lsst.afw.math.BackgroundList`
576 Detector level background lists associated with the calexp.
578 Returns
579 -------
580 calExp : `lsst.afw.image.ExposureF`
581 The calexp with the originally subtracted background restored.
582 skyCorrBase : `lsst.afw.math.BackgroundList`
583 The inverted original background models; the genesis for skyCorr.
584 """
585 image = calExp.getMaskedImage()
587 # Invert all elements of the existing bg model; restore in calexp
588 for calBkgElement in calBkg:
589 statsImage = calBkgElement[0].getStatsImage()
590 statsImage *= -1
591 skyCorrBase = calBkg.getImage()
592 image -= skyCorrBase
594 stats = np.nanpercentile(skyCorrBase.array, [50, 75, 25])
595 self.log.info(
596 "Detector %d: Original background restored (BG median = %.1f counts, BG IQR = %.1f counts)",
597 calExp.getDetector().getId(),
598 -stats[0],
599 np.subtract(*stats[1:]),
600 )
602 # Iteratively subtract bg, re-detect sources, and add bg back on
603 if self.config.doMaskObjects:
604 maskFrac0 = 1 - np.sum(calExp.mask.array == 0) / calExp.mask.array.size
605 self.maskObjects.findObjects(calExp)
606 maskFrac1 = 1 - np.sum(calExp.mask.array == 0) / calExp.mask.array.size
608 self.log.info(
609 "Detector %d: Iterative source detection and mask growth has increased masked area by %.1f%%",
610 calExp.getDetector().getId(),
611 (100 * (maskFrac1 - maskFrac0)),
612 )
614 return calExp, skyCorrBase
616 def _validateBgModel(self, bgModelID, bgModel, config):
617 """Check that the background model contains enough valid superpixels,
618 and raise a useful error if not.
620 Parameters
621 ----------
622 bgModelID : `str`
623 Identifier for the background model.
624 bgModel : `~lsst.pipe.tasks.background.FocalPlaneBackground`
625 Background model to check.
626 config : `~lsst.pipe.tasks.background.FocalPlaneBackgroundConfig`
627 Configuration used to create the background model.
628 """
629 bgModelArray = bgModel._numbers.getArray()
630 spArea = (config.xSize / config.pixelSize) * (config.ySize / config.pixelSize)
631 self.log.info(
632 "%s: FP background model constructed using %.2f x %.2f mm (%d x %d pixel) superpixels",
633 bgModelID,
634 config.xSize,
635 config.ySize,
636 int(config.xSize / config.pixelSize),
637 int(config.ySize / config.pixelSize),
638 )
639 self.log.info(
640 "%s: Pixel data exists in %d of %d superpixels; the most populated superpixel is %.1f%% filled",
641 bgModelID,
642 np.sum(bgModelArray > 0),
643 bgModelArray.size,
644 100 * np.max(bgModelArray) / spArea,
645 )
647 thresh = config.minFrac * spArea
648 if np.all(bgModelArray < thresh):
649 raise RuntimeError(
650 f"No background model superpixels are more than {100*config.minFrac}% filled. "
651 "Try decreasing the minFrac configuration parameter, optimizing the subset of detectors "
652 "being processed, or increasing the number of detectors being processed."
653 )
655 with warnings.catch_warnings():
656 warnings.filterwarnings("ignore", r"invalid value encountered")
657 stats = np.nanpercentile(bgModel.getStatsImage().array, [50, 75, 25])
658 self.log.info(
659 "%s: FP BG median = %.1f counts, FP BG IQR = %.1f counts",
660 bgModelID,
661 stats[0],
662 np.subtract(*stats[1:]),
663 )
665 def _fitSkyFrame(self, calExps, masks, skyCorrs, skyFrames, backgroundToPhotometricRatioHandles):
666 """Determine the full focal plane sky frame scale factor relative to
667 an input list of calibrated exposures.
669 This method measures the sky frame scale on all inputs, resulting in
670 values equal to the background method solveScales().
672 Input skyCorrs are updated in-place.
674 Parameters
675 ----------
676 calExps : `list` [`lsst.afw.image.ExposureF`]
677 Calibrated exposures to be background subtracted.
678 masks : `list` [`lsst.afw.image.Mask`]
679 Masks associated with the input calibrated exposures.
680 skyCorrs : `list` [`lsst.afw.math.BackgroundList`]
681 Background lists associated with the input calibrated exposures.
682 skyFrames : `list` [`lsst.afw.image.ExposureF`]
683 Sky frame calibration data for the input detectors.
685 Returns
686 -------
687 scale : `float`
688 Fitted scale factor applied to the sky frame.
689 """
690 skyBkgs = []
691 scales = []
692 for calExpHandle, mask, skyCorr, skyFrameHandle, backgroundToPhotometricRatioHandle in zip(
693 calExps, masks, skyCorrs, skyFrames, backgroundToPhotometricRatioHandles
694 ):
695 calExp = self._getCalExp(calExpHandle, mask, skyCorr, backgroundToPhotometricRatioHandle)
696 skyFrame = self._getSkyFrame(skyFrameHandle)
697 skyBkg = self.sky.exposureToBackground(skyFrame)
698 del skyFrame # Free up memory
699 skyBkgs.append(skyBkg)
700 # Return a tuple of gridded image and sky frame clipped means
701 samples = self.sky.measureScale(calExp.getMaskedImage(), skyBkg)
702 scales.append(samples)
703 scale = self.sky.solveScales(scales)
704 for skyCorr, skyBkg in zip(skyCorrs, skyBkgs):
705 bgData = list(skyBkg[0])
706 bg = bgData[0]
707 statsImage = bg.getStatsImage().clone()
708 statsImage *= scale
709 newBg = BackgroundMI(bg.getImageBBox(), statsImage)
710 newBgData = [newBg] + bgData[1:]
711 skyCorr.append(newBgData)
712 self.log.info("Sky frame subtracted with a scale factor of %.5f", scale)
713 return scale