Coverage for python/lsst/ip/isr/isrTask.py: 15%
1005 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-07 02:17 -0800
1# This file is part of ip_isr.
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/>.
22__all__ = ["IsrTask", "IsrTaskConfig"]
24import math
25import numpy
26import numbers
28import lsst.geom
29import lsst.afw.image as afwImage
30import lsst.afw.math as afwMath
31import lsst.pex.config as pexConfig
32import lsst.pipe.base as pipeBase
33import lsst.pipe.base.connectionTypes as cT
35from contextlib import contextmanager
36from deprecated.sphinx import deprecated
37from lsstDebug import getDebugFrame
39from lsst.afw.cameraGeom import NullLinearityType
40from lsst.afw.display import getDisplay
41from lsst.meas.algorithms.detection import SourceDetectionTask
42from lsst.utils.timer import timeMethod
44from . import isrFunctions
45from . import isrQa
46from . import linearize
47from .defects import Defects
49from .assembleCcdTask import AssembleCcdTask
50from .binImageDataTask import BinImageDataTask
51from .crosstalk import CrosstalkTask, CrosstalkCalib
52from .fringe import FringeTask
53from .isr import maskNans
54from .masking import MaskingTask
55from .overscan import OverscanCorrectionTask
56from .straylight import StrayLightTask
57from .vignette import VignetteTask
58from .ampOffset import AmpOffsetTask
59from .deferredCharge import DeferredChargeTask
60from .isrStatistics import IsrStatisticsTask
61from .ptcDataset import PhotonTransferCurveDataset
62from .isrFunctions import compareCameraKeywords
65def crosstalkSourceLookup(datasetType, registry, quantumDataId, collections):
66 """Lookup function to identify crosstalkSource entries.
68 This should return an empty list under most circumstances. Only
69 when inter-chip crosstalk has been identified should this be
70 populated.
72 Parameters
73 ----------
74 datasetType : `str`
75 Dataset to lookup.
76 registry : `lsst.daf.butler.Registry`
77 Butler registry to query.
78 quantumDataId : `lsst.daf.butler.DataCoordinate`
79 Expanded data id to transform to identify crosstalkSources. The
80 ``detector`` entry will be stripped.
81 collections : `lsst.daf.butler.CollectionSearch`
82 Collections to search through.
84 Returns
85 -------
86 results : `list` [`lsst.daf.butler.DatasetRef`]
87 List of datasets that match the query that will be used as
88 crosstalkSources.
89 """
90 newDataId = quantumDataId.subset(registry.dimensions.conform(["instrument", "exposure"]))
91 results = set(registry.queryDatasets(datasetType, collections=collections, dataId=newDataId,
92 findFirst=True))
93 # In some contexts, calling `.expanded()` to expand all data IDs in the
94 # query results can be a lot faster because it vectorizes lookups. But in
95 # this case, expandDataId shouldn't need to hit the database at all in the
96 # steady state, because only the detector record is unknown and those are
97 # cached in the registry.
98 records = {k: newDataId.records[k] for k in newDataId.dimensions.elements}
99 return [ref.expanded(registry.expandDataId(ref.dataId, records=records)) for ref in results]
102class IsrTaskConnections(pipeBase.PipelineTaskConnections,
103 dimensions={"instrument", "exposure", "detector"},
104 defaultTemplates={}):
105 ccdExposure = cT.Input(
106 name="raw",
107 doc="Input exposure to process.",
108 storageClass="Exposure",
109 dimensions=["instrument", "exposure", "detector"],
110 )
111 camera = cT.PrerequisiteInput(
112 name="camera",
113 storageClass="Camera",
114 doc="Input camera to construct complete exposures.",
115 dimensions=["instrument"],
116 isCalibration=True,
117 )
119 crosstalk = cT.PrerequisiteInput(
120 name="crosstalk",
121 doc="Input crosstalk object",
122 storageClass="CrosstalkCalib",
123 dimensions=["instrument", "detector"],
124 isCalibration=True,
125 minimum=0, # can fall back to cameraGeom
126 )
127 crosstalkSources = cT.PrerequisiteInput(
128 name="isrOverscanCorrected",
129 doc="Overscan corrected input images.",
130 storageClass="Exposure",
131 dimensions=["instrument", "exposure", "detector"],
132 deferLoad=True,
133 multiple=True,
134 lookupFunction=crosstalkSourceLookup,
135 minimum=0, # not needed for all instruments, no config to control this
136 )
137 bias = cT.PrerequisiteInput(
138 name="bias",
139 doc="Input bias calibration.",
140 storageClass="ExposureF",
141 dimensions=["instrument", "detector"],
142 isCalibration=True,
143 )
144 dark = cT.PrerequisiteInput(
145 name='dark',
146 doc="Input dark calibration.",
147 storageClass="ExposureF",
148 dimensions=["instrument", "detector"],
149 isCalibration=True,
150 )
151 flat = cT.PrerequisiteInput(
152 name="flat",
153 doc="Input flat calibration.",
154 storageClass="ExposureF",
155 dimensions=["instrument", "physical_filter", "detector"],
156 isCalibration=True,
157 )
158 ptc = cT.PrerequisiteInput(
159 name="ptc",
160 doc="Input Photon Transfer Curve dataset",
161 storageClass="PhotonTransferCurveDataset",
162 dimensions=["instrument", "detector"],
163 isCalibration=True,
164 )
165 fringes = cT.PrerequisiteInput(
166 name="fringe",
167 doc="Input fringe calibration.",
168 storageClass="ExposureF",
169 dimensions=["instrument", "physical_filter", "detector"],
170 isCalibration=True,
171 minimum=0, # only needed for some bands, even when enabled
172 )
173 strayLightData = cT.PrerequisiteInput(
174 name='yBackground',
175 doc="Input stray light calibration.",
176 storageClass="StrayLightData",
177 dimensions=["instrument", "physical_filter", "detector"],
178 deferLoad=True,
179 isCalibration=True,
180 minimum=0, # only needed for some bands, even when enabled
181 )
182 bfKernel = cT.PrerequisiteInput(
183 name='bfKernel',
184 doc="Input brighter-fatter kernel.",
185 storageClass="NumpyArray",
186 dimensions=["instrument"],
187 isCalibration=True,
188 minimum=0, # can use either bfKernel or newBFKernel
189 )
190 newBFKernel = cT.PrerequisiteInput(
191 name='brighterFatterKernel',
192 doc="Newer complete kernel + gain solutions.",
193 storageClass="BrighterFatterKernel",
194 dimensions=["instrument", "detector"],
195 isCalibration=True,
196 minimum=0, # can use either bfKernel or newBFKernel
197 )
198 defects = cT.PrerequisiteInput(
199 name='defects',
200 doc="Input defect tables.",
201 storageClass="Defects",
202 dimensions=["instrument", "detector"],
203 isCalibration=True,
204 )
205 linearizer = cT.PrerequisiteInput(
206 name='linearizer',
207 storageClass="Linearizer",
208 doc="Linearity correction calibration.",
209 dimensions=["instrument", "detector"],
210 isCalibration=True,
211 minimum=0, # can fall back to cameraGeom
212 )
213 opticsTransmission = cT.PrerequisiteInput(
214 name="transmission_optics",
215 storageClass="TransmissionCurve",
216 doc="Transmission curve due to the optics.",
217 dimensions=["instrument"],
218 isCalibration=True,
219 )
220 filterTransmission = cT.PrerequisiteInput(
221 name="transmission_filter",
222 storageClass="TransmissionCurve",
223 doc="Transmission curve due to the filter.",
224 dimensions=["instrument", "physical_filter"],
225 isCalibration=True,
226 )
227 sensorTransmission = cT.PrerequisiteInput(
228 name="transmission_sensor",
229 storageClass="TransmissionCurve",
230 doc="Transmission curve due to the sensor.",
231 dimensions=["instrument", "detector"],
232 isCalibration=True,
233 )
234 atmosphereTransmission = cT.PrerequisiteInput(
235 name="transmission_atmosphere",
236 storageClass="TransmissionCurve",
237 doc="Transmission curve due to the atmosphere.",
238 dimensions=["instrument"],
239 isCalibration=True,
240 )
241 illumMaskedImage = cT.PrerequisiteInput(
242 name="illum",
243 doc="Input illumination correction.",
244 storageClass="MaskedImageF",
245 dimensions=["instrument", "physical_filter", "detector"],
246 isCalibration=True,
247 )
248 deferredChargeCalib = cT.PrerequisiteInput(
249 name="cpCtiCalib",
250 doc="Deferred charge/CTI correction dataset.",
251 storageClass="IsrCalib",
252 dimensions=["instrument", "detector"],
253 isCalibration=True,
254 )
256 outputExposure = cT.Output(
257 name='postISRCCD',
258 doc="Output ISR processed exposure.",
259 storageClass="Exposure",
260 dimensions=["instrument", "exposure", "detector"],
261 )
262 preInterpExposure = cT.Output(
263 name='preInterpISRCCD',
264 doc="Output ISR processed exposure, with pixels left uninterpolated.",
265 storageClass="ExposureF",
266 dimensions=["instrument", "exposure", "detector"],
267 )
268 outputBin1Exposure = cT.Output(
269 name="postIsrBin1",
270 doc="First binned image.",
271 storageClass="ExposureF",
272 dimensions=["instrument", "exposure", "detector"],
273 )
274 outputBin2Exposure = cT.Output(
275 name="postIsrBin2",
276 doc="Second binned image.",
277 storageClass="ExposureF",
278 dimensions=["instrument", "exposure", "detector"],
279 )
281 outputOssThumbnail = cT.Output(
282 name="OssThumb",
283 doc="Output Overscan-subtracted thumbnail image.",
284 storageClass="Thumbnail",
285 dimensions=["instrument", "exposure", "detector"],
286 )
287 outputFlattenedThumbnail = cT.Output(
288 name="FlattenedThumb",
289 doc="Output flat-corrected thumbnail image.",
290 storageClass="Thumbnail",
291 dimensions=["instrument", "exposure", "detector"],
292 )
293 outputStatistics = cT.Output(
294 name="isrStatistics",
295 doc="Output of additional statistics table.",
296 storageClass="StructuredDataDict",
297 dimensions=["instrument", "exposure", "detector"],
298 )
300 def __init__(self, *, config=None):
301 super().__init__(config=config)
303 if config.doBias is not True:
304 self.prerequisiteInputs.remove("bias")
305 if config.doLinearize is not True:
306 self.prerequisiteInputs.remove("linearizer")
307 if config.doCrosstalk is not True:
308 self.prerequisiteInputs.remove("crosstalkSources")
309 self.prerequisiteInputs.remove("crosstalk")
310 if config.doBrighterFatter is not True:
311 self.prerequisiteInputs.remove("bfKernel")
312 self.prerequisiteInputs.remove("newBFKernel")
313 if config.doDefect is not True:
314 self.prerequisiteInputs.remove("defects")
315 if config.doDark is not True:
316 self.prerequisiteInputs.remove("dark")
317 if config.doFlat is not True:
318 self.prerequisiteInputs.remove("flat")
319 if config.doFringe is not True:
320 self.prerequisiteInputs.remove("fringes")
321 if config.doStrayLight is not True:
322 self.prerequisiteInputs.remove("strayLightData")
323 if config.usePtcGains is not True and config.usePtcReadNoise is not True:
324 self.prerequisiteInputs.remove("ptc")
325 if config.doAttachTransmissionCurve is not True:
326 self.prerequisiteInputs.remove("opticsTransmission")
327 self.prerequisiteInputs.remove("filterTransmission")
328 self.prerequisiteInputs.remove("sensorTransmission")
329 self.prerequisiteInputs.remove("atmosphereTransmission")
330 else:
331 if config.doUseOpticsTransmission is not True:
332 self.prerequisiteInputs.remove("opticsTransmission")
333 if config.doUseFilterTransmission is not True:
334 self.prerequisiteInputs.remove("filterTransmission")
335 if config.doUseSensorTransmission is not True:
336 self.prerequisiteInputs.remove("sensorTransmission")
337 if config.doUseAtmosphereTransmission is not True:
338 self.prerequisiteInputs.remove("atmosphereTransmission")
339 if config.doIlluminationCorrection is not True:
340 self.prerequisiteInputs.remove("illumMaskedImage")
341 if config.doDeferredCharge is not True:
342 self.prerequisiteInputs.remove("deferredChargeCalib")
344 if config.doWrite is not True:
345 self.outputs.remove("outputExposure")
346 self.outputs.remove("preInterpExposure")
347 self.outputs.remove("outputFlattenedThumbnail")
348 self.outputs.remove("outputOssThumbnail")
349 self.outputs.remove("outputStatistics")
350 self.outputs.remove("outputBin1Exposure")
351 self.outputs.remove("outputBin2Exposure")
352 else:
353 if config.doBinnedExposures is not True:
354 self.outputs.remove("outputBin1Exposure")
355 self.outputs.remove("outputBin2Exposure")
356 if config.doSaveInterpPixels is not True:
357 self.outputs.remove("preInterpExposure")
358 if config.qa.doThumbnailOss is not True:
359 self.outputs.remove("outputOssThumbnail")
360 if config.qa.doThumbnailFlattened is not True:
361 self.outputs.remove("outputFlattenedThumbnail")
362 if config.doCalculateStatistics is not True:
363 self.outputs.remove("outputStatistics")
366class IsrTaskConfig(pipeBase.PipelineTaskConfig,
367 pipelineConnections=IsrTaskConnections):
368 """Configuration parameters for IsrTask.
370 Items are grouped in the order in which they are executed by the task.
371 """
372 datasetType = pexConfig.Field(
373 dtype=str,
374 doc="Dataset type for input data; users will typically leave this alone, "
375 "but camera-specific ISR tasks will override it",
376 default="raw",
377 )
379 fallbackFilterName = pexConfig.Field(
380 dtype=str,
381 doc="Fallback default filter name for calibrations.",
382 optional=True
383 )
384 useFallbackDate = pexConfig.Field(
385 dtype=bool,
386 doc="Pass observation date when using fallback filter.",
387 default=False,
388 )
389 expectWcs = pexConfig.Field(
390 dtype=bool,
391 default=True,
392 doc="Expect input science images to have a WCS (set False for e.g. spectrographs)."
393 )
394 fwhm = pexConfig.Field(
395 dtype=float,
396 doc="FWHM of PSF in arcseconds (currently unused).",
397 default=1.0,
398 )
399 qa = pexConfig.ConfigField(
400 dtype=isrQa.IsrQaConfig,
401 doc="QA related configuration options.",
402 )
403 doHeaderProvenance = pexConfig.Field(
404 dtype=bool,
405 default=True,
406 doc="Write calibration identifiers into output exposure header?",
407 )
409 # Calib checking configuration:
410 doRaiseOnCalibMismatch = pexConfig.Field(
411 dtype=bool,
412 default=False,
413 doc="Should IsrTask halt if exposure and calibration header values do not match?",
414 )
415 cameraKeywordsToCompare = pexConfig.ListField(
416 dtype=str,
417 doc="List of header keywords to compare between exposure and calibrations.",
418 default=[],
419 )
421 # Image conversion configuration
422 doConvertIntToFloat = pexConfig.Field(
423 dtype=bool,
424 doc="Convert integer raw images to floating point values?",
425 default=True,
426 )
428 # Saturated pixel handling.
429 doSaturation = pexConfig.Field(
430 dtype=bool,
431 doc="Mask saturated pixels? NB: this is totally independent of the"
432 " interpolation option - this is ONLY setting the bits in the mask."
433 " To have them interpolated make sure doSaturationInterpolation=True",
434 default=True,
435 )
436 saturatedMaskName = pexConfig.Field(
437 dtype=str,
438 doc="Name of mask plane to use in saturation detection and interpolation",
439 default="SAT",
440 )
441 saturation = pexConfig.Field(
442 dtype=float,
443 doc="The saturation level to use if no Detector is present in the Exposure (ignored if NaN)",
444 default=float("NaN"),
445 )
446 growSaturationFootprintSize = pexConfig.Field(
447 dtype=int,
448 doc="Number of pixels by which to grow the saturation footprints",
449 default=1,
450 )
452 # Suspect pixel handling.
453 doSuspect = pexConfig.Field(
454 dtype=bool,
455 doc="Mask suspect pixels?",
456 default=False,
457 )
458 suspectMaskName = pexConfig.Field(
459 dtype=str,
460 doc="Name of mask plane to use for suspect pixels",
461 default="SUSPECT",
462 )
463 numEdgeSuspect = pexConfig.Field(
464 dtype=int,
465 doc="Number of edge pixels to be flagged as untrustworthy.",
466 default=0,
467 )
468 edgeMaskLevel = pexConfig.ChoiceField(
469 dtype=str,
470 doc="Mask edge pixels in which coordinate frame: DETECTOR or AMP?",
471 default="DETECTOR",
472 allowed={
473 'DETECTOR': 'Mask only the edges of the full detector.',
474 'AMP': 'Mask edges of each amplifier.',
475 },
476 )
478 # Initial masking options.
479 doSetBadRegions = pexConfig.Field(
480 dtype=bool,
481 doc="Should we set the level of all BAD patches of the chip to the chip's average value?",
482 default=True,
483 )
484 badStatistic = pexConfig.ChoiceField(
485 dtype=str,
486 doc="How to estimate the average value for BAD regions.",
487 default='MEANCLIP',
488 allowed={
489 "MEANCLIP": "Correct using the (clipped) mean of good data",
490 "MEDIAN": "Correct using the median of the good data",
491 },
492 )
494 # Overscan subtraction configuration.
495 doOverscan = pexConfig.Field(
496 dtype=bool,
497 doc="Do overscan subtraction?",
498 default=True,
499 )
500 overscan = pexConfig.ConfigurableField(
501 target=OverscanCorrectionTask,
502 doc="Overscan subtraction task for image segments.",
503 )
505 # Amplifier to CCD assembly configuration
506 doAssembleCcd = pexConfig.Field(
507 dtype=bool,
508 default=True,
509 doc="Assemble amp-level exposures into a ccd-level exposure?"
510 )
511 assembleCcd = pexConfig.ConfigurableField(
512 target=AssembleCcdTask,
513 doc="CCD assembly task",
514 )
516 # General calibration configuration.
517 doAssembleIsrExposures = pexConfig.Field(
518 dtype=bool,
519 default=False,
520 doc="Assemble amp-level calibration exposures into ccd-level exposure?"
521 )
522 doTrimToMatchCalib = pexConfig.Field(
523 dtype=bool,
524 default=False,
525 doc="Trim raw data to match calibration bounding boxes?"
526 )
528 # Bias subtraction.
529 doBias = pexConfig.Field(
530 dtype=bool,
531 doc="Apply bias frame correction?",
532 default=True,
533 )
534 biasDataProductName = pexConfig.Field(
535 dtype=str,
536 doc="Name of the bias data product",
537 default="bias",
538 )
539 doBiasBeforeOverscan = pexConfig.Field(
540 dtype=bool,
541 doc="Reverse order of overscan and bias correction.",
542 default=False
543 )
545 # Deferred charge correction.
546 doDeferredCharge = pexConfig.Field(
547 dtype=bool,
548 doc="Apply deferred charge correction?",
549 default=False,
550 )
551 deferredChargeCorrection = pexConfig.ConfigurableField(
552 target=DeferredChargeTask,
553 doc="Deferred charge correction task.",
554 )
556 # Variance construction
557 doVariance = pexConfig.Field(
558 dtype=bool,
559 doc="Calculate variance?",
560 default=True
561 )
562 gain = pexConfig.Field(
563 dtype=float,
564 doc="The gain to use if no Detector is present in the Exposure (ignored if NaN)",
565 default=float("NaN"),
566 )
567 readNoise = pexConfig.Field(
568 dtype=float,
569 doc="The read noise to use if no Detector is present in the Exposure",
570 default=0.0,
571 )
572 doEmpiricalReadNoise = pexConfig.Field(
573 dtype=bool,
574 default=False,
575 doc="Calculate empirical read noise instead of value from AmpInfo data?"
576 )
577 usePtcReadNoise = pexConfig.Field(
578 dtype=bool,
579 default=False,
580 doc="Use read noise values from the Photon Transfer Curve?"
581 )
582 maskNegativeVariance = pexConfig.Field(
583 dtype=bool,
584 default=True,
585 doc="Mask pixels that claim a negative variance? This likely indicates a failure "
586 "in the measurement of the overscan at an edge due to the data falling off faster "
587 "than the overscan model can account for it."
588 )
589 negativeVarianceMaskName = pexConfig.Field(
590 dtype=str,
591 default="BAD",
592 doc="Mask plane to use to mark pixels with negative variance, if `maskNegativeVariance` is True.",
593 )
594 # Linearization.
595 doLinearize = pexConfig.Field(
596 dtype=bool,
597 doc="Correct for nonlinearity of the detector's response?",
598 default=True,
599 )
601 # Crosstalk.
602 doCrosstalk = pexConfig.Field(
603 dtype=bool,
604 doc="Apply intra-CCD crosstalk correction?",
605 default=False,
606 )
607 doCrosstalkBeforeAssemble = pexConfig.Field(
608 dtype=bool,
609 doc="Apply crosstalk correction before CCD assembly, and before trimming?",
610 default=False,
611 )
612 crosstalk = pexConfig.ConfigurableField(
613 target=CrosstalkTask,
614 doc="Intra-CCD crosstalk correction",
615 )
617 # Masking options.
618 doDefect = pexConfig.Field(
619 dtype=bool,
620 doc="Apply correction for CCD defects, e.g. hot pixels?",
621 default=True,
622 )
623 doNanMasking = pexConfig.Field(
624 dtype=bool,
625 doc="Mask non-finite (NAN, inf) pixels?",
626 default=True,
627 )
628 doWidenSaturationTrails = pexConfig.Field(
629 dtype=bool,
630 doc="Widen bleed trails based on their width?",
631 default=True
632 )
634 # Brighter-Fatter correction.
635 doBrighterFatter = pexConfig.Field(
636 dtype=bool,
637 default=False,
638 doc="Apply the brighter-fatter correction?"
639 )
640 doFluxConservingBrighterFatterCorrection = pexConfig.Field(
641 dtype=bool,
642 default=False,
643 doc="Apply the flux-conserving BFE correction by Miller et al.?"
644 )
645 brighterFatterLevel = pexConfig.ChoiceField(
646 dtype=str,
647 default="DETECTOR",
648 doc="The level at which to correct for brighter-fatter.",
649 allowed={
650 "AMP": "Every amplifier treated separately.",
651 "DETECTOR": "One kernel per detector",
652 }
653 )
654 brighterFatterMaxIter = pexConfig.Field(
655 dtype=int,
656 default=10,
657 doc="Maximum number of iterations for the brighter-fatter correction"
658 )
659 brighterFatterThreshold = pexConfig.Field(
660 dtype=float,
661 default=1000,
662 doc="Threshold used to stop iterating the brighter-fatter correction. It is the "
663 "absolute value of the difference between the current corrected image and the one "
664 "from the previous iteration summed over all the pixels."
665 )
666 brighterFatterApplyGain = pexConfig.Field(
667 dtype=bool,
668 default=True,
669 doc="Should the gain be applied when applying the brighter-fatter correction?"
670 )
671 brighterFatterMaskListToInterpolate = pexConfig.ListField(
672 dtype=str,
673 doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
674 "correction.",
675 default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
676 )
677 brighterFatterMaskGrowSize = pexConfig.Field(
678 dtype=int,
679 default=0,
680 doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
681 "when brighter-fatter correction is applied."
682 )
684 # Dark subtraction.
685 doDark = pexConfig.Field(
686 dtype=bool,
687 doc="Apply dark frame correction?",
688 default=True,
689 )
690 darkDataProductName = pexConfig.Field(
691 dtype=str,
692 doc="Name of the dark data product",
693 default="dark",
694 )
696 # Camera-specific stray light removal.
697 doStrayLight = pexConfig.Field(
698 dtype=bool,
699 doc="Subtract stray light in the y-band (due to encoder LEDs)?",
700 default=False,
701 )
702 strayLight = pexConfig.ConfigurableField(
703 target=StrayLightTask,
704 doc="y-band stray light correction"
705 )
707 # Flat correction.
708 doFlat = pexConfig.Field(
709 dtype=bool,
710 doc="Apply flat field correction?",
711 default=True,
712 )
713 flatDataProductName = pexConfig.Field(
714 dtype=str,
715 doc="Name of the flat data product",
716 default="flat",
717 )
718 flatScalingType = pexConfig.ChoiceField(
719 dtype=str,
720 doc="The method for scaling the flat on the fly.",
721 default='USER',
722 allowed={
723 "USER": "Scale by flatUserScale",
724 "MEAN": "Scale by the inverse of the mean",
725 "MEDIAN": "Scale by the inverse of the median",
726 },
727 )
728 flatUserScale = pexConfig.Field(
729 dtype=float,
730 doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise",
731 default=1.0,
732 )
733 doTweakFlat = pexConfig.Field(
734 dtype=bool,
735 doc="Tweak flats to match observed amplifier ratios?",
736 default=False
737 )
739 # Amplifier normalization based on gains instead of using flats
740 # configuration.
741 doApplyGains = pexConfig.Field(
742 dtype=bool,
743 doc="Correct the amplifiers for their gains instead of applying flat correction",
744 default=False,
745 )
746 usePtcGains = pexConfig.Field(
747 dtype=bool,
748 doc="Use the gain values from the input Photon Transfer Curve?",
749 default=False,
750 )
751 normalizeGains = pexConfig.Field(
752 dtype=bool,
753 doc="Normalize all the amplifiers in each CCD to have the same median value.",
754 default=False,
755 )
757 # Fringe correction.
758 doFringe = pexConfig.Field(
759 dtype=bool,
760 doc="Apply fringe correction?",
761 default=True,
762 )
763 fringe = pexConfig.ConfigurableField(
764 target=FringeTask,
765 doc="Fringe subtraction task",
766 )
767 fringeAfterFlat = pexConfig.Field(
768 dtype=bool,
769 doc="Do fringe subtraction after flat-fielding?",
770 default=True,
771 )
773 # Amp offset correction.
774 doAmpOffset = pexConfig.Field(
775 doc="Calculate amp offset corrections?",
776 dtype=bool,
777 default=False,
778 )
779 ampOffset = pexConfig.ConfigurableField(
780 doc="Amp offset correction task.",
781 target=AmpOffsetTask,
782 )
784 # Initial CCD-level background statistics options.
785 doMeasureBackground = pexConfig.Field(
786 dtype=bool,
787 doc="Measure the background level on the reduced image?",
788 default=False,
789 )
791 # Camera-specific masking configuration.
792 doCameraSpecificMasking = pexConfig.Field(
793 dtype=bool,
794 doc="Mask camera-specific bad regions?",
795 default=False,
796 )
797 masking = pexConfig.ConfigurableField(
798 target=MaskingTask,
799 doc="Masking task."
800 )
802 # Interpolation options.
803 doInterpolate = pexConfig.Field(
804 dtype=bool,
805 doc="Interpolate masked pixels?",
806 default=True,
807 )
808 doSaturationInterpolation = pexConfig.Field(
809 dtype=bool,
810 doc="Perform interpolation over pixels masked as saturated?"
811 " NB: This is independent of doSaturation; if that is False this plane"
812 " will likely be blank, resulting in a no-op here.",
813 default=True,
814 )
815 doNanInterpolation = pexConfig.Field(
816 dtype=bool,
817 doc="Perform interpolation over pixels masked as NaN?"
818 " NB: This is independent of doNanMasking; if that is False this plane"
819 " will likely be blank, resulting in a no-op here.",
820 default=True,
821 )
822 doNanInterpAfterFlat = pexConfig.Field(
823 dtype=bool,
824 doc=("If True, ensure we interpolate NaNs after flat-fielding, even if we "
825 "also have to interpolate them before flat-fielding."),
826 default=False,
827 )
828 maskListToInterpolate = pexConfig.ListField(
829 dtype=str,
830 doc="List of mask planes that should be interpolated.",
831 default=['SAT', 'BAD'],
832 )
833 doSaveInterpPixels = pexConfig.Field(
834 dtype=bool,
835 doc="Save a copy of the pre-interpolated pixel values?",
836 default=False,
837 )
838 useLegacyInterp = pexConfig.Field(
839 dtype=bool,
840 doc="Use the legacy interpolation algorithm. If False use Gaussian Process.",
841 default=True,
842 )
844 # Default photometric calibration options.
845 fluxMag0T1 = pexConfig.DictField(
846 keytype=str,
847 itemtype=float,
848 doc="The approximate flux of a zero-magnitude object in a one-second exposure, per filter.",
849 default=dict((f, pow(10.0, 0.4*m)) for f, m in (("Unknown", 28.0),
850 ))
851 )
852 defaultFluxMag0T1 = pexConfig.Field(
853 dtype=float,
854 doc="Default value for fluxMag0T1 (for an unrecognized filter).",
855 default=pow(10.0, 0.4*28.0)
856 )
858 # Vignette correction configuration.
859 doVignette = pexConfig.Field(
860 dtype=bool,
861 doc=("Compute and attach the validPolygon defining the unvignetted region to the exposure "
862 "according to vignetting parameters?"),
863 default=False,
864 )
865 doMaskVignettePolygon = pexConfig.Field(
866 dtype=bool,
867 doc=("Add a mask bit for pixels within the vignetted region. Ignored if doVignette "
868 "is False"),
869 default=True,
870 )
871 vignetteValue = pexConfig.Field(
872 dtype=float,
873 doc="Value to replace image array pixels with in the vignetted region? Ignored if None.",
874 optional=True,
875 default=None,
876 )
877 vignette = pexConfig.ConfigurableField(
878 target=VignetteTask,
879 doc="Vignetting task.",
880 )
882 # Transmission curve configuration.
883 doAttachTransmissionCurve = pexConfig.Field(
884 dtype=bool,
885 default=False,
886 doc="Construct and attach a wavelength-dependent throughput curve for this CCD image?"
887 )
888 doUseOpticsTransmission = pexConfig.Field(
889 dtype=bool,
890 default=True,
891 doc="Load and use transmission_optics (if doAttachTransmissionCurve is True)?"
892 )
893 doUseFilterTransmission = pexConfig.Field(
894 dtype=bool,
895 default=True,
896 doc="Load and use transmission_filter (if doAttachTransmissionCurve is True)?"
897 )
898 doUseSensorTransmission = pexConfig.Field(
899 dtype=bool,
900 default=True,
901 doc="Load and use transmission_sensor (if doAttachTransmissionCurve is True)?"
902 )
903 doUseAtmosphereTransmission = pexConfig.Field(
904 dtype=bool,
905 default=True,
906 doc="Load and use transmission_atmosphere (if doAttachTransmissionCurve is True)?"
907 )
909 # Illumination correction.
910 doIlluminationCorrection = pexConfig.Field(
911 dtype=bool,
912 default=False,
913 doc="Perform illumination correction?"
914 )
915 illuminationCorrectionDataProductName = pexConfig.Field(
916 dtype=str,
917 doc="Name of the illumination correction data product.",
918 default="illumcor",
919 )
920 illumScale = pexConfig.Field(
921 dtype=float,
922 doc="Scale factor for the illumination correction.",
923 default=1.0,
924 )
925 illumFilters = pexConfig.ListField(
926 dtype=str,
927 default=[],
928 doc="Only perform illumination correction for these filters."
929 )
931 # Calculate image quality statistics?
932 doStandardStatistics = pexConfig.Field(
933 dtype=bool,
934 doc="Should standard image quality statistics be calculated?",
935 default=True,
936 )
937 # Calculate additional statistics?
938 doCalculateStatistics = pexConfig.Field(
939 dtype=bool,
940 doc="Should additional ISR statistics be calculated?",
941 default=False,
942 )
943 isrStats = pexConfig.ConfigurableField(
944 target=IsrStatisticsTask,
945 doc="Task to calculate additional statistics.",
946 )
948 # Make binned images?
949 doBinnedExposures = pexConfig.Field(
950 dtype=bool,
951 doc="Should binned exposures be calculated?",
952 default=False,
953 )
954 binning = pexConfig.ConfigurableField(
955 target=BinImageDataTask,
956 doc="Task to bin the exposure.",
957 )
958 binFactor1 = pexConfig.Field(
959 dtype=int,
960 doc="Binning factor for first binned exposure. This is intended for a finely binned output.",
961 default=8,
962 check=lambda x: x > 1,
963 )
964 binFactor2 = pexConfig.Field(
965 dtype=int,
966 doc="Binning factor for second binned exposure. This is intended for a coarsely binned output.",
967 default=64,
968 check=lambda x: x > 1,
969 )
971 # Write the outputs to disk. If ISR is run as a subtask, this may not
972 # be needed.
973 doWrite = pexConfig.Field(
974 dtype=bool,
975 doc="Persist postISRCCD?",
976 default=True,
977 )
979 def validate(self):
980 super().validate()
981 if self.doFlat and self.doApplyGains:
982 raise ValueError("You may not specify both doFlat and doApplyGains")
983 if self.doBiasBeforeOverscan and self.doTrimToMatchCalib:
984 raise ValueError("You may not specify both doBiasBeforeOverscan and doTrimToMatchCalib")
985 if self.doSaturationInterpolation and self.saturatedMaskName not in self.maskListToInterpolate:
986 self.maskListToInterpolate.append(self.saturatedMaskName)
987 if not self.doSaturationInterpolation and self.saturatedMaskName in self.maskListToInterpolate:
988 self.maskListToInterpolate.remove(self.saturatedMaskName)
989 if self.doNanInterpolation and "UNMASKEDNAN" not in self.maskListToInterpolate:
990 self.maskListToInterpolate.append("UNMASKEDNAN")
991 if self.ampOffset.doApplyAmpOffset and not self.doAmpOffset:
992 raise ValueError("ampOffset.doApplyAmpOffset requires doAmpOffset to be True.")
995class IsrTask(pipeBase.PipelineTask):
996 """Apply common instrument signature correction algorithms to a raw frame.
998 The process for correcting imaging data is very similar from
999 camera to camera. This task provides a vanilla implementation of
1000 doing these corrections, including the ability to turn certain
1001 corrections off if they are not needed. The inputs to the primary
1002 method, `run()`, are a raw exposure to be corrected and the
1003 calibration data products. The raw input is a single chip sized
1004 mosaic of all amps including overscans and other non-science
1005 pixels.
1007 The __init__ method sets up the subtasks for ISR processing, using
1008 the defaults from `lsst.ip.isr`.
1010 Parameters
1011 ----------
1012 args : `list`
1013 Positional arguments passed to the Task constructor.
1014 None used at this time.
1015 kwargs : `dict`, optional
1016 Keyword arguments passed on to the Task constructor.
1017 None used at this time.
1018 """
1019 ConfigClass = IsrTaskConfig
1020 _DefaultName = "isr"
1022 def __init__(self, **kwargs):
1023 super().__init__(**kwargs)
1024 self.makeSubtask("assembleCcd")
1025 self.makeSubtask("crosstalk")
1026 self.makeSubtask("strayLight")
1027 self.makeSubtask("fringe")
1028 self.makeSubtask("masking")
1029 self.makeSubtask("overscan")
1030 self.makeSubtask("vignette")
1031 self.makeSubtask("ampOffset")
1032 self.makeSubtask("deferredChargeCorrection")
1033 self.makeSubtask("isrStats")
1034 self.makeSubtask("binning")
1036 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1037 inputs = butlerQC.get(inputRefs)
1039 try:
1040 inputs['detectorNum'] = inputRefs.ccdExposure.dataId['detector']
1041 except Exception as e:
1042 raise ValueError("Failure to find valid detectorNum value for Dataset %s: %s." %
1043 (inputRefs, e))
1045 detector = inputs['ccdExposure'].getDetector()
1047 # This is use for header provenance.
1048 additionalInputDates = {}
1050 if self.config.doCrosstalk is True:
1051 # Crosstalk sources need to be defined by the pipeline
1052 # yaml if they exist.
1053 if 'crosstalk' in inputs and inputs['crosstalk'] is not None:
1054 if not isinstance(inputs['crosstalk'], CrosstalkCalib):
1055 inputs['crosstalk'] = CrosstalkCalib.fromTable(inputs['crosstalk'])
1056 else:
1057 coeffVector = (self.config.crosstalk.crosstalkValues
1058 if self.config.crosstalk.useConfigCoefficients else None)
1059 crosstalkCalib = CrosstalkCalib().fromDetector(detector, coeffVector=coeffVector)
1060 inputs['crosstalk'] = crosstalkCalib
1061 if inputs['crosstalk'].interChip and len(inputs['crosstalk'].interChip) > 0:
1062 if 'crosstalkSources' not in inputs:
1063 self.log.warning("No crosstalkSources found for chip with interChip terms!")
1065 if self.doLinearize(detector) is True:
1066 if 'linearizer' in inputs:
1067 if isinstance(inputs['linearizer'], dict):
1068 linearizer = linearize.Linearizer(detector=detector, log=self.log)
1069 linearizer.fromYaml(inputs['linearizer'])
1070 self.log.warning("Dictionary linearizers will be deprecated in DM-28741.")
1071 elif isinstance(inputs['linearizer'], numpy.ndarray):
1072 linearizer = linearize.Linearizer(table=inputs.get('linearizer', None),
1073 detector=detector,
1074 log=self.log)
1075 self.log.warning("Bare lookup table linearizers will be deprecated in DM-28741.")
1076 else:
1077 linearizer = inputs['linearizer']
1078 self.log.info("Loading linearizer from the Butler.")
1079 linearizer.log = self.log
1080 inputs['linearizer'] = linearizer
1081 else:
1082 inputs['linearizer'] = linearize.Linearizer(detector=detector, log=self.log)
1083 self.log.info("Constructing linearizer from cameraGeom information.")
1085 if self.config.doDefect is True:
1086 if "defects" in inputs and inputs['defects'] is not None:
1087 # defects is loaded as a BaseCatalog with columns
1088 # x0, y0, width, height. Masking expects a list of defects
1089 # defined by their bounding box
1090 if not isinstance(inputs["defects"], Defects):
1091 inputs["defects"] = Defects.fromTable(inputs["defects"])
1093 # Load the correct style of brighter-fatter kernel, and repack
1094 # the information as a numpy array.
1095 brighterFatterSource = None
1096 if self.config.doBrighterFatter:
1097 brighterFatterKernel = inputs.pop('newBFKernel', None)
1098 if brighterFatterKernel is None:
1099 # This type of kernel must be in (y, x) index
1100 # ordering, as it used directly as the .array
1101 # component of the afwImage kernel.
1102 brighterFatterKernel = inputs.get('bfKernel', None)
1103 brighterFatterSource = 'bfKernel'
1104 additionalInputDates[brighterFatterSource] = self.extractCalibDate(brighterFatterKernel)
1106 if brighterFatterKernel is None:
1107 # This was requested by the config, but none were found.
1108 raise RuntimeError("No brighter-fatter kernel was supplied.")
1109 elif not isinstance(brighterFatterKernel, numpy.ndarray):
1110 # This is a ISR calib kernel. These kernels are
1111 # generated in (x, y) index ordering, and need to be
1112 # transposed to be used directly as the .array
1113 # component of the afwImage kernel. This is done
1114 # explicitly below when setting the ``bfKernel``
1115 # input.
1116 brighterFatterSource = 'newBFKernel'
1117 additionalInputDates[brighterFatterSource] = self.extractCalibDate(brighterFatterKernel)
1119 detName = detector.getName()
1120 level = brighterFatterKernel.level
1122 # This is expected to be a dictionary of amp-wise gains.
1123 inputs['bfGains'] = brighterFatterKernel.gain
1124 if self.config.brighterFatterLevel == 'DETECTOR':
1125 kernel = None
1126 if level == 'DETECTOR':
1127 if detName in brighterFatterKernel.detKernels:
1128 kernel = brighterFatterKernel.detKernels[detName]
1129 else:
1130 raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1131 elif level == 'AMP':
1132 self.log.warning("Making DETECTOR level kernel from AMP based brighter "
1133 "fatter kernels.")
1134 brighterFatterKernel.makeDetectorKernelFromAmpwiseKernels(detName)
1135 kernel = brighterFatterKernel.detKernels[detName]
1136 if kernel is None:
1137 raise RuntimeError("Could not identify brighter-fatter kernel!")
1138 # Do the one single transpose here so the kernel
1139 # can be directly loaded into the afwImage .array
1140 # component.
1141 inputs['bfKernel'] = numpy.transpose(kernel)
1142 elif self.config.brighterFatterLevel == 'AMP':
1143 raise NotImplementedError("Per-amplifier brighter-fatter correction not implemented")
1145 if self.config.doFringe is True and self.fringe.checkFilter(inputs['ccdExposure']):
1146 expId = inputs['ccdExposure'].info.id
1147 inputs['fringes'] = self.fringe.loadFringes(inputs['fringes'],
1148 expId=expId,
1149 assembler=self.assembleCcd
1150 if self.config.doAssembleIsrExposures else None)
1151 else:
1152 inputs['fringes'] = pipeBase.Struct(fringes=None)
1154 if self.config.doStrayLight is True and self.strayLight.checkFilter(inputs['ccdExposure']):
1155 if 'strayLightData' not in inputs:
1156 inputs['strayLightData'] = None
1158 if self.config.doHeaderProvenance:
1159 # Add calibration provenanace info to header.
1160 exposureMetadata = inputs['ccdExposure'].getMetadata()
1162 # These inputs change name during this step. These should
1163 # have matching entries in the additionalInputDates dict.
1164 additionalInputs = []
1165 if self.config.doBrighterFatter:
1166 additionalInputs.append(brighterFatterSource)
1168 for inputName in sorted(list(inputs.keys()) + additionalInputs):
1169 reference = getattr(inputRefs, inputName, None)
1170 if reference is not None and hasattr(reference, "run"):
1171 runKey = f"LSST CALIB RUN {inputName.upper()}"
1172 runValue = reference.run
1173 idKey = f"LSST CALIB UUID {inputName.upper()}"
1174 idValue = str(reference.id)
1175 dateKey = f"LSST CALIB DATE {inputName.upper()}"
1177 if inputName in additionalInputDates:
1178 dateValue = additionalInputDates[inputName]
1179 else:
1180 dateValue = self.extractCalibDate(inputs[inputName])
1182 exposureMetadata[runKey] = runValue
1183 exposureMetadata[idKey] = idValue
1184 exposureMetadata[dateKey] = dateValue
1186 outputs = self.run(**inputs)
1187 butlerQC.put(outputs, outputRefs)
1189 @timeMethod
1190 def run(self, ccdExposure, *, camera=None, bias=None, linearizer=None,
1191 crosstalk=None, crosstalkSources=None,
1192 dark=None, flat=None, ptc=None, bfKernel=None, bfGains=None, defects=None,
1193 fringes=pipeBase.Struct(fringes=None), opticsTransmission=None, filterTransmission=None,
1194 sensorTransmission=None, atmosphereTransmission=None,
1195 detectorNum=None, strayLightData=None, illumMaskedImage=None,
1196 deferredChargeCalib=None,
1197 ):
1198 """Perform instrument signature removal on an exposure.
1200 Steps included in the ISR processing, in order performed, are:
1202 - saturation and suspect pixel masking
1203 - overscan subtraction
1204 - CCD assembly of individual amplifiers
1205 - bias subtraction
1206 - variance image construction
1207 - linearization of non-linear response
1208 - crosstalk masking
1209 - brighter-fatter correction
1210 - dark subtraction
1211 - fringe correction
1212 - stray light subtraction
1213 - flat correction
1214 - masking of known defects and camera specific features
1215 - vignette calculation
1216 - appending transmission curve and distortion model
1218 Parameters
1219 ----------
1220 ccdExposure : `lsst.afw.image.Exposure`
1221 The raw exposure that is to be run through ISR. The
1222 exposure is modified by this method.
1223 camera : `lsst.afw.cameraGeom.Camera`, optional
1224 The camera geometry for this exposure. Required if
1225 one or more of ``ccdExposure``, ``bias``, ``dark``, or
1226 ``flat`` does not have an associated detector.
1227 bias : `lsst.afw.image.Exposure`, optional
1228 Bias calibration frame.
1229 linearizer : `lsst.ip.isr.linearize.LinearizeBase`, optional
1230 Functor for linearization.
1231 crosstalk : `lsst.ip.isr.crosstalk.CrosstalkCalib`, optional
1232 Calibration for crosstalk.
1233 crosstalkSources : `list`, optional
1234 List of possible crosstalk sources.
1235 dark : `lsst.afw.image.Exposure`, optional
1236 Dark calibration frame.
1237 flat : `lsst.afw.image.Exposure`, optional
1238 Flat calibration frame.
1239 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
1240 Photon transfer curve dataset, with, e.g., gains
1241 and read noise.
1242 bfKernel : `numpy.ndarray`, optional
1243 Brighter-fatter kernel.
1244 bfGains : `dict` of `float`, optional
1245 Gains used to override the detector's nominal gains for the
1246 brighter-fatter correction. A dict keyed by amplifier name for
1247 the detector in question.
1248 defects : `lsst.ip.isr.Defects`, optional
1249 List of defects.
1250 fringes : `lsst.pipe.base.Struct`, optional
1251 Struct containing the fringe correction data, with
1252 elements:
1254 ``fringes``
1255 fringe calibration frame (`lsst.afw.image.Exposure`)
1256 ``seed``
1257 random seed derived from the ``ccdExposureId`` for random
1258 number generator (`numpy.uint32`)
1259 opticsTransmission: `lsst.afw.image.TransmissionCurve`, optional
1260 A ``TransmissionCurve`` that represents the throughput of the,
1261 optics, to be evaluated in focal-plane coordinates.
1262 filterTransmission : `lsst.afw.image.TransmissionCurve`
1263 A ``TransmissionCurve`` that represents the throughput of the
1264 filter itself, to be evaluated in focal-plane coordinates.
1265 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1266 A ``TransmissionCurve`` that represents the throughput of the
1267 sensor itself, to be evaluated in post-assembly trimmed detector
1268 coordinates.
1269 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1270 A ``TransmissionCurve`` that represents the throughput of the
1271 atmosphere, assumed to be spatially constant.
1272 detectorNum : `int`, optional
1273 The integer number for the detector to process.
1274 strayLightData : `object`, optional
1275 Opaque object containing calibration information for stray-light
1276 correction. If `None`, no correction will be performed.
1277 illumMaskedImage : `lsst.afw.image.MaskedImage`, optional
1278 Illumination correction image.
1280 Returns
1281 -------
1282 result : `lsst.pipe.base.Struct`
1283 Result struct with component:
1285 ``exposure``
1286 The fully ISR corrected exposure.
1287 (`lsst.afw.image.Exposure`)
1288 ``outputExposure``
1289 An alias for ``exposure``. (`lsst.afw.image.Exposure`)
1290 ``ossThumb``
1291 Thumbnail image of the exposure after overscan subtraction.
1292 (`numpy.ndarray`)
1293 ``flattenedThumb``
1294 Thumbnail image of the exposure after flat-field correction.
1295 (`numpy.ndarray`)
1296 ``outputStatistics``
1297 Values of the additional statistics calculated.
1299 Raises
1300 ------
1301 RuntimeError
1302 Raised if a configuration option is set to `True`, but the
1303 required calibration data has not been specified.
1305 Notes
1306 -----
1307 The current processed exposure can be viewed by setting the
1308 appropriate `lsstDebug` entries in the ``debug.display``
1309 dictionary. The names of these entries correspond to some of
1310 the `IsrTaskConfig` Boolean options, with the value denoting the
1311 frame to use. The exposure is shown inside the matching
1312 option check and after the processing of that step has
1313 finished. The steps with debug points are:
1315 * doAssembleCcd
1316 * doBias
1317 * doCrosstalk
1318 * doBrighterFatter
1319 * doDark
1320 * doFringe
1321 * doStrayLight
1322 * doFlat
1324 In addition, setting the ``postISRCCD`` entry displays the
1325 exposure after all ISR processing has finished.
1326 """
1328 ccdExposure = self.ensureExposure(ccdExposure, camera, detectorNum)
1329 bias = self.ensureExposure(bias, camera, detectorNum)
1330 dark = self.ensureExposure(dark, camera, detectorNum)
1331 flat = self.ensureExposure(flat, camera, detectorNum)
1333 ccd = ccdExposure.getDetector()
1334 filterLabel = ccdExposure.getFilter()
1335 physicalFilter = isrFunctions.getPhysicalFilter(filterLabel, self.log)
1337 if not ccd:
1338 assert not self.config.doAssembleCcd, "You need a Detector to run assembleCcd."
1339 ccd = [FakeAmp(ccdExposure, self.config)]
1341 # Validate Input
1342 if self.config.doBias and bias is None:
1343 raise RuntimeError("Must supply a bias exposure if config.doBias=True.")
1344 if self.doLinearize(ccd) and linearizer is None:
1345 raise RuntimeError("Must supply a linearizer if config.doLinearize=True for this detector.")
1346 if self.config.doBrighterFatter and bfKernel is None:
1347 raise RuntimeError("Must supply a kernel if config.doBrighterFatter=True.")
1348 if self.config.doDark and dark is None:
1349 raise RuntimeError("Must supply a dark exposure if config.doDark=True.")
1350 if self.config.doFlat and flat is None:
1351 raise RuntimeError("Must supply a flat exposure if config.doFlat=True.")
1352 if self.config.doDefect and defects is None:
1353 raise RuntimeError("Must supply defects if config.doDefect=True.")
1354 if (self.config.doFringe and physicalFilter in self.fringe.config.filters
1355 and fringes.fringes is None):
1356 # The `fringes` object needs to be a pipeBase.Struct, as
1357 # we use it as a `dict` for the parameters of
1358 # `FringeTask.run()`. The `fringes.fringes` `list` may
1359 # not be `None` if `doFringe=True`. Otherwise, raise.
1360 raise RuntimeError("Must supply fringe exposure as a pipeBase.Struct.")
1361 if (self.config.doIlluminationCorrection and physicalFilter in self.config.illumFilters
1362 and illumMaskedImage is None):
1363 raise RuntimeError("Must supply an illumcor if config.doIlluminationCorrection=True.")
1364 if (self.config.doDeferredCharge and deferredChargeCalib is None):
1365 raise RuntimeError("Must supply a deferred charge calibration if config.doDeferredCharge=True.")
1366 if (self.config.usePtcGains and ptc is None):
1367 raise RuntimeError("No ptcDataset provided to use PTC gains.")
1368 if (self.config.usePtcReadNoise and ptc is None):
1369 raise RuntimeError("No ptcDataset provided to use PTC read noise.")
1371 # Validate that the inputs match the exposure configuration.
1372 exposureMetadata = ccdExposure.getMetadata()
1373 doRaise = self.config.doRaiseOnCalibMismatch
1374 keywords = self.config.cameraKeywordsToCompare
1375 if self.config.doBias:
1376 compareCameraKeywords(doRaise, keywords, exposureMetadata, bias, "bias", log=self.log)
1377 self.compareUnits(bias.metadata, "bias")
1378 if self.config.doBrighterFatter:
1379 compareCameraKeywords(doRaise, keywords, exposureMetadata, bfKernel, "bf", log=self.log)
1380 if self.config.doCrosstalk:
1381 compareCameraKeywords(doRaise, keywords, exposureMetadata, crosstalk, "crosstalk", log=self.log)
1382 if self.config.doDark:
1383 compareCameraKeywords(doRaise, keywords, exposureMetadata, dark, "dark", log=self.log)
1384 self.compareUnits(dark.metadata, "dark")
1385 if self.config.doDefect:
1386 compareCameraKeywords(doRaise, keywords, exposureMetadata, defects, "defects", log=self.log)
1387 if self.config.doDeferredCharge:
1388 compareCameraKeywords(
1389 doRaise,
1390 keywords,
1391 exposureMetadata,
1392 deferredChargeCalib,
1393 "CTI",
1394 log=self.log,
1395 )
1396 if self.config.doFlat:
1397 compareCameraKeywords(doRaise, keywords, exposureMetadata, flat, "flat", log=self.log)
1398 self.compareUnits(flat.metadata, "flat")
1399 if (self.config.doFringe and physicalFilter in self.fringe.config.filters):
1400 compareCameraKeywords(
1401 doRaise,
1402 keywords,
1403 exposureMetadata,
1404 fringes.fringes,
1405 "fringe",
1406 log=self.log,
1407 )
1408 if (self.config.doIlluminationCorrection and physicalFilter in self.config.illumFilters):
1409 compareCameraKeywords(
1410 doRaise,
1411 keywords,
1412 exposureMetadata,
1413 illumMaskedImage,
1414 "illumination",
1415 log=self.log,
1416 )
1417 if self.doLinearize(ccd):
1418 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer, "linearizer", log=self.log)
1419 if self.config.usePtcGains or self.config.usePtcReadNoise:
1420 compareCameraKeywords(doRaise, keywords, exposureMetadata, ptc, "PTC", log=self.log)
1421 if self.config.doStrayLight:
1422 compareCameraKeywords(
1423 doRaise,
1424 keywords,
1425 exposureMetadata,
1426 strayLightData,
1427 "straylight",
1428 log=self.log,
1429 )
1431 # Start in adu. Update units to electrons when gain is applied:
1432 # updateVariance, applyGains
1433 # Check if needed during/after BFE correction, CTI correction.
1434 exposureMetadata["LSST ISR UNITS"] = "adu"
1435 exposureMetadata["LSST ISR CROSSTALK APPLIED"] = False
1436 exposureMetadata["LSST ISR LINEARIZER APPLIED"] = False
1437 exposureMetadata["LSST ISR CTI APPLIED"] = False
1438 exposureMetadata["LSST ISR BIAS APPLIED"] = False
1439 exposureMetadata["LSST ISR DARK APPLIED"] = False
1440 exposureMetadata["LSST ISR BF APPLIED"] = False
1441 exposureMetadata["LSST ISR FLAT APPLIED"] = False
1442 exposureMetadata["LSST ISR DEFECTS APPLIED"] = False
1444 # Begin ISR processing.
1445 if self.config.doConvertIntToFloat:
1446 self.log.info("Converting exposure to floating point values.")
1447 ccdExposure = self.convertIntToFloat(ccdExposure)
1449 if self.config.doBias and self.config.doBiasBeforeOverscan:
1450 self.log.info("Applying bias correction.")
1451 isrFunctions.biasCorrection(ccdExposure.getMaskedImage(), bias.getMaskedImage(),
1452 trimToFit=self.config.doTrimToMatchCalib)
1453 self.debugView(ccdExposure, "doBias")
1454 ccdExposure.metadata["LSST ISR BIAS APPLIED"] = True
1456 # Amplifier level processing.
1457 overscans = []
1459 if self.config.doOverscan and self.config.overscan.doParallelOverscan:
1460 # This will attempt to mask bleed pixels across all amplifiers.
1461 self.overscan.maskParallelOverscan(ccdExposure, ccd)
1463 for amp in ccd:
1464 # if ccdExposure is one amp,
1465 # check for coverage to prevent performing ops multiple times
1466 if ccdExposure.getBBox().contains(amp.getBBox()):
1467 # Check for fully masked bad amplifiers,
1468 # and generate masks for SUSPECT and SATURATED values.
1469 badAmp = self.maskAmplifier(ccdExposure, amp, defects)
1471 if self.config.doOverscan and not badAmp:
1472 # Overscan correction on amp-by-amp basis.
1473 overscanResults = self.overscanCorrection(ccdExposure, amp)
1474 self.log.debug("Corrected overscan for amplifier %s.", amp.getName())
1475 if overscanResults is not None and \
1476 self.config.qa is not None and self.config.qa.saveStats is True:
1477 if isinstance(overscanResults.overscanMean, float):
1478 # Only serial overscan was run
1479 mean = overscanResults.overscanMean
1480 median = overscanResults.overscanMedian
1481 sigma = overscanResults.overscanSigma
1482 residMean = overscanResults.residualMean
1483 residMedian = overscanResults.residualMedian
1484 residSigma = overscanResults.residualSigma
1485 else:
1486 # Both serial and parallel overscan were
1487 # run. Only report serial here.
1488 mean = overscanResults.overscanMean[0]
1489 median = overscanResults.overscanMedian[0]
1490 sigma = overscanResults.overscanSigma[0]
1491 residMean = overscanResults.residualMean[0]
1492 residMedian = overscanResults.residualMedian[0]
1493 residSigma = overscanResults.residualSigma[0]
1495 self.metadata[f"FIT MEDIAN {amp.getName()}"] = median
1496 self.metadata[f"FIT MEAN {amp.getName()}"] = mean
1497 self.metadata[f"FIT STDEV {amp.getName()}"] = sigma
1498 self.log.debug(" Overscan stats for amplifer %s: %f +/- %f",
1499 amp.getName(), mean, sigma)
1501 self.metadata[f"RESIDUAL MEDIAN {amp.getName()}"] = residMedian
1502 self.metadata[f"RESIDUAL MEAN {amp.getName()}"] = residMean
1503 self.metadata[f"RESIDUAL STDEV {amp.getName()}"] = residSigma
1504 self.log.debug(" Overscan stats for amplifer %s after correction: %f +/- %f",
1505 amp.getName(), residMean, residSigma)
1507 ccdExposure.getMetadata().set('OVERSCAN', "Overscan corrected")
1508 else:
1509 if badAmp:
1510 self.log.warning("Amplifier %s is bad.", amp.getName())
1511 overscanResults = None
1513 overscans.append(overscanResults if overscanResults is not None else None)
1514 else:
1515 self.log.info("Skipped OSCAN for %s.", amp.getName())
1517 # Define an effective PTC that will contain the gain and readout
1518 # noise to be used throughout the ISR task.
1519 ptc = self.defineEffectivePtc(ptc, ccd, bfGains, overscans, exposureMetadata)
1521 if self.config.doDeferredCharge:
1522 self.log.info("Applying deferred charge/CTI correction.")
1523 self.deferredChargeCorrection.run(
1524 ccdExposure,
1525 deferredChargeCalib,
1526 gains=ptc.gain,
1527 )
1528 self.debugView(ccdExposure, "doDeferredCharge")
1529 ccdExposure.metadata["LSST ISR CTI APPLIED"] = True
1531 if self.config.doCrosstalk and self.config.doCrosstalkBeforeAssemble:
1532 self.log.info("Applying crosstalk correction.")
1533 self.crosstalk.run(ccdExposure, crosstalk=crosstalk,
1534 crosstalkSources=crosstalkSources, camera=camera)
1535 self.debugView(ccdExposure, "doCrosstalk")
1536 ccdExposure.metadata["LSST ISR CROSSTALK APPLIED"] = True
1538 if self.config.doAssembleCcd:
1539 self.log.info("Assembling CCD from amplifiers.")
1540 ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
1542 if self.config.expectWcs and not ccdExposure.getWcs():
1543 self.log.warning("No WCS found in input exposure.")
1544 self.debugView(ccdExposure, "doAssembleCcd")
1546 ossThumb = None
1547 if self.config.qa.doThumbnailOss:
1548 ossThumb = isrQa.makeThumbnail(ccdExposure, isrQaConfig=self.config.qa)
1550 if self.config.doBias and not self.config.doBiasBeforeOverscan:
1551 self.log.info("Applying bias correction.")
1552 isrFunctions.biasCorrection(ccdExposure.getMaskedImage(), bias.getMaskedImage(),
1553 trimToFit=self.config.doTrimToMatchCalib)
1554 self.debugView(ccdExposure, "doBias")
1555 ccdExposure.metadata["LSST ISR BIAS APPLIED"] = True
1557 if self.config.doVariance:
1558 for amp in ccd:
1559 if ccdExposure.getBBox().contains(amp.getBBox()):
1560 self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1561 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1562 self.updateVariance(ampExposure, amp, ptc)
1564 if self.config.qa is not None and self.config.qa.saveStats is True:
1565 qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1566 afwMath.MEDIAN | afwMath.STDEVCLIP)
1567 self.metadata[f"ISR VARIANCE {amp.getName()} MEDIAN"] = \
1568 qaStats.getValue(afwMath.MEDIAN)
1569 self.metadata[f"ISR VARIANCE {amp.getName()} STDEV"] = \
1570 qaStats.getValue(afwMath.STDEVCLIP)
1571 self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1572 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1573 qaStats.getValue(afwMath.STDEVCLIP))
1574 if self.config.maskNegativeVariance:
1575 self.maskNegativeVariance(ccdExposure)
1577 if self.doLinearize(ccd):
1578 self.log.info("Applying linearizer.")
1579 linearizer.applyLinearity(image=ccdExposure.getMaskedImage().getImage(),
1580 detector=ccd, log=self.log)
1581 ccdExposure.metadata["LSST ISR LINEARIZER APPLIED"] = True
1583 if self.config.doCrosstalk and not self.config.doCrosstalkBeforeAssemble:
1584 self.log.info("Applying crosstalk correction.")
1585 self.crosstalk.run(ccdExposure, crosstalk=crosstalk,
1586 crosstalkSources=crosstalkSources)
1587 self.debugView(ccdExposure, "doCrosstalk")
1588 ccdExposure.metadata["LSST ISR CROSSTALK APPLIED"] = True
1590 # Masking block. Optionally mask known defects,NaN/inf pixels,
1591 # widen trails, and do anything else the camera needs. Saturated and
1592 # suspect pixels have already been masked.
1593 if self.config.doDefect:
1594 self.log.info("Masking defects.")
1595 self.maskDefect(ccdExposure, defects)
1596 ccdExposure.metadata["LSST ISR DEFECTS APPLIED"] = True
1598 if self.config.numEdgeSuspect > 0:
1599 self.log.info("Masking edges as SUSPECT.")
1600 self.maskEdges(ccdExposure, numEdgePixels=self.config.numEdgeSuspect,
1601 maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
1603 if self.config.doNanMasking:
1604 self.log.info("Masking non-finite (NAN, inf) value pixels.")
1605 self.maskNan(ccdExposure)
1607 if self.config.doWidenSaturationTrails:
1608 self.log.info("Widening saturation trails.")
1609 isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
1611 if self.config.doCameraSpecificMasking:
1612 self.log.info("Masking regions for camera specific reasons.")
1613 self.masking.run(ccdExposure)
1615 if self.config.doBrighterFatter:
1616 # We need to apply flats and darks before we can interpolate, and
1617 # we need to interpolate before we do B-F, but we do B-F without
1618 # the flats and darks applied so we can work in units of electrons
1619 # or holes. This context manager applies and then removes the darks
1620 # and flats.
1621 #
1622 # We also do not want to interpolate values here, so operate on
1623 # temporary images so we can apply only the BF-correction and roll
1624 # back the interpolation.
1625 interpExp = ccdExposure.clone()
1626 with self.flatContext(interpExp, flat, dark):
1627 isrFunctions.interpolateFromMask(
1628 maskedImage=interpExp.getMaskedImage(),
1629 fwhm=self.config.fwhm,
1630 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1631 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1632 useLegacyInterp=self.config.useLegacyInterp,
1633 )
1634 bfExp = interpExp.clone()
1636 self.log.info("Applying brighter-fatter correction using kernel type %s / gains %s.",
1637 type(bfKernel), type(bfGains))
1638 if self.config.doFluxConservingBrighterFatterCorrection:
1639 bfResults = isrFunctions.fluxConservingBrighterFatterCorrection(
1640 bfExp,
1641 bfKernel,
1642 self.config.brighterFatterMaxIter,
1643 self.config.brighterFatterThreshold,
1644 self.config.brighterFatterApplyGain,
1645 bfGains
1646 )
1647 else:
1648 bfResults = isrFunctions.brighterFatterCorrection(
1649 bfExp,
1650 bfKernel,
1651 self.config.brighterFatterMaxIter,
1652 self.config.brighterFatterThreshold,
1653 self.config.brighterFatterApplyGain,
1654 bfGains
1655 )
1656 bfCorrIters = bfResults[1]
1657 self.metadata["LSST ISR BF ITERS"] = bfCorrIters
1658 if bfCorrIters == self.config.brighterFatterMaxIter - 1:
1659 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1660 bfResults[0])
1661 else:
1662 self.log.info("Finished brighter-fatter correction in %d iterations.",
1663 bfResults[1])
1664 image = ccdExposure.getMaskedImage().getImage()
1665 bfCorr = bfExp.getMaskedImage().getImage()
1666 bfCorr -= interpExp.getMaskedImage().getImage()
1667 image += bfCorr
1669 # Applying the brighter-fatter correction applies a
1670 # convolution to the science image. At the edges this
1671 # convolution may not have sufficient valid pixels to
1672 # produce a valid correction. Mark pixels within the size
1673 # of the brighter-fatter kernel as EDGE to warn of this
1674 # fact.
1675 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1676 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1677 maskPlane="EDGE")
1679 if self.config.brighterFatterMaskGrowSize > 0:
1680 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1681 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1682 isrFunctions.growMasks(ccdExposure.getMask(),
1683 radius=self.config.brighterFatterMaskGrowSize,
1684 maskNameList=maskPlane,
1685 maskValue=maskPlane)
1687 self.debugView(ccdExposure, "doBrighterFatter")
1688 ccdExposure.metadata["LSST ISR BF APPLIED"] = True
1690 if self.config.doDark:
1691 self.log.info("Applying dark correction.")
1692 self.darkCorrection(ccdExposure, dark)
1693 self.debugView(ccdExposure, "doDark")
1694 ccdExposure.metadata["LSST ISR DARK APPLIED"] = True
1696 if self.config.doFringe and not self.config.fringeAfterFlat:
1697 self.log.info("Applying fringe correction before flat.")
1698 self.fringe.run(ccdExposure, **fringes.getDict())
1699 self.debugView(ccdExposure, "doFringe")
1701 if self.config.doStrayLight and self.strayLight.check(ccdExposure):
1702 self.log.info("Checking strayLight correction.")
1703 self.strayLight.run(ccdExposure, strayLightData)
1704 self.debugView(ccdExposure, "doStrayLight")
1706 if self.config.doFlat:
1707 self.log.info("Applying flat correction.")
1708 self.flatCorrection(ccdExposure, flat)
1709 self.debugView(ccdExposure, "doFlat")
1710 ccdExposure.metadata["LSST ISR FLAT APPLIED"] = True
1711 # TODO: DM-49159
1712 # Add metadata re: type of flat.
1714 if self.config.doApplyGains:
1715 self.log.info("Applying gain correction instead of flat.")
1716 isrFunctions.applyGains(ccdExposure, self.config.normalizeGains,
1717 ptcGains=ptc.gain)
1718 exposureMetadata["LSST ISR UNITS"] = "electron"
1720 if self.config.doFringe and self.config.fringeAfterFlat:
1721 self.log.info("Applying fringe correction after flat.")
1722 self.fringe.run(ccdExposure, **fringes.getDict())
1724 if self.config.doVignette:
1725 if self.config.doMaskVignettePolygon:
1726 self.log.info("Constructing, attaching, and masking vignette polygon.")
1727 else:
1728 self.log.info("Constructing and attaching vignette polygon.")
1729 self.vignettePolygon = self.vignette.run(
1730 exposure=ccdExposure, doUpdateMask=self.config.doMaskVignettePolygon,
1731 vignetteValue=self.config.vignetteValue, log=self.log)
1733 if self.config.doAttachTransmissionCurve:
1734 self.log.info("Adding transmission curves.")
1735 isrFunctions.attachTransmissionCurve(ccdExposure, opticsTransmission=opticsTransmission,
1736 filterTransmission=filterTransmission,
1737 sensorTransmission=sensorTransmission,
1738 atmosphereTransmission=atmosphereTransmission)
1740 flattenedThumb = None
1741 if self.config.qa.doThumbnailFlattened:
1742 flattenedThumb = isrQa.makeThumbnail(ccdExposure, isrQaConfig=self.config.qa)
1744 if self.config.doIlluminationCorrection and physicalFilter in self.config.illumFilters:
1745 self.log.info("Performing illumination correction.")
1746 isrFunctions.illuminationCorrection(ccdExposure.getMaskedImage(),
1747 illumMaskedImage, illumScale=self.config.illumScale,
1748 trimToFit=self.config.doTrimToMatchCalib)
1750 preInterpExp = None
1751 if self.config.doSaveInterpPixels:
1752 preInterpExp = ccdExposure.clone()
1754 # Reset and interpolate bad pixels.
1755 #
1756 # Large contiguous bad regions (which should have the BAD mask
1757 # bit set) should have their values set to the image median.
1758 # This group should include defects and bad amplifiers. As the
1759 # area covered by these defects are large, there's little
1760 # reason to expect that interpolation would provide a more
1761 # useful value.
1762 #
1763 # Smaller defects can be safely interpolated after the larger
1764 # regions have had their pixel values reset. This ensures
1765 # that the remaining defects adjacent to bad amplifiers (as an
1766 # example) do not attempt to interpolate extreme values.
1767 if self.config.doSetBadRegions:
1768 badPixelCount, badPixelValue = isrFunctions.setBadRegions(ccdExposure)
1769 if badPixelCount > 0:
1770 self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1772 if self.config.doInterpolate:
1773 self.log.info("Interpolating masked pixels.")
1774 isrFunctions.interpolateFromMask(
1775 maskedImage=ccdExposure.getMaskedImage(),
1776 fwhm=self.config.fwhm,
1777 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1778 maskNameList=list(self.config.maskListToInterpolate),
1779 useLegacyInterp=self.config.useLegacyInterp,
1780 )
1782 self.roughZeroPoint(ccdExposure)
1784 # Calculate amp offset corrections within the CCD.
1785 if self.config.doAmpOffset:
1786 if self.config.ampOffset.doApplyAmpOffset:
1787 self.log.info("Measuring and applying amp offset corrections.")
1788 else:
1789 self.log.info("Measuring amp offset corrections only, without applying them.")
1790 self.ampOffset.run(ccdExposure)
1792 if self.config.doMeasureBackground:
1793 self.log.info("Measuring background level.")
1794 self.measureBackground(ccdExposure, self.config.qa)
1796 if self.config.qa is not None and self.config.qa.saveStats is True:
1797 for amp in ccd:
1798 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1799 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
1800 afwMath.MEDIAN | afwMath.STDEVCLIP)
1801 self.metadata[f"ISR BACKGROUND {amp.getName()} MEDIAN"] = qaStats.getValue(afwMath.MEDIAN)
1802 self.metadata[f"ISR BACKGROUND {amp.getName()} STDEV"] = \
1803 qaStats.getValue(afwMath.STDEVCLIP)
1804 self.log.debug(" Background stats for amplifer %s: %f +/- %f",
1805 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1806 qaStats.getValue(afwMath.STDEVCLIP))
1808 # Calculate standard image quality statistics
1809 if self.config.doStandardStatistics:
1810 metadata = ccdExposure.getMetadata()
1811 for amp in ccd:
1812 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1813 ampName = amp.getName()
1814 metadata[f"LSST ISR MASK SAT {ampName}"] = isrFunctions.countMaskedPixels(
1815 ampExposure.getMaskedImage(),
1816 [self.config.saturatedMaskName]
1817 )
1818 metadata[f"LSST ISR MASK BAD {ampName}"] = isrFunctions.countMaskedPixels(
1819 ampExposure.getMaskedImage(),
1820 ["BAD"]
1821 )
1822 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
1823 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP)
1825 metadata[f"LSST ISR FINAL MEAN {ampName}"] = qaStats.getValue(afwMath.MEAN)
1826 metadata[f"LSST ISR FINAL MEDIAN {ampName}"] = qaStats.getValue(afwMath.MEDIAN)
1827 metadata[f"LSST ISR FINAL STDEV {ampName}"] = qaStats.getValue(afwMath.STDEVCLIP)
1829 k1 = f"LSST ISR FINAL MEDIAN {ampName}"
1830 k2 = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
1831 if self.config.doOverscan and k1 in metadata and k2 in metadata:
1832 metadata[f"LSST ISR LEVEL {ampName}"] = metadata[k1] - metadata[k2]
1833 else:
1834 metadata[f"LSST ISR LEVEL {ampName}"] = numpy.nan
1836 # calculate additional statistics.
1837 outputStatistics = None
1838 if self.config.doCalculateStatistics:
1839 outputStatistics = self.isrStats.run(ccdExposure, serialOverscanResults=overscans,
1840 parallelOverscanResults=[None for _ in overscans],
1841 bias=bias, dark=dark, flat=flat, ptc=ptc,
1842 defects=defects,
1843 doLegacyCtiStatistics=True).results
1845 # do any binning.
1846 outputBin1Exposure = None
1847 outputBin2Exposure = None
1848 if self.config.doBinnedExposures:
1849 self.log.info("Creating binned exposures.")
1850 outputBin1Exposure = self.binning.run(
1851 ccdExposure,
1852 binFactor=self.config.binFactor1,
1853 ).outputData
1854 outputBin2Exposure = self.binning.run(
1855 ccdExposure,
1856 binFactor=self.config.binFactor2,
1857 ).outputData
1859 self.debugView(ccdExposure, "postISRCCD")
1861 return pipeBase.Struct(
1862 exposure=ccdExposure,
1863 ossThumb=ossThumb,
1864 flattenedThumb=flattenedThumb,
1866 outputBin1Exposure=outputBin1Exposure,
1867 outputBin2Exposure=outputBin2Exposure,
1869 preInterpExposure=preInterpExp,
1870 outputExposure=ccdExposure,
1871 outputOssThumbnail=ossThumb,
1872 outputFlattenedThumbnail=flattenedThumb,
1873 outputStatistics=outputStatistics,
1874 )
1876 def defineEffectivePtc(self, ptcDataset, detector, bfGains, overScans, metadata):
1877 """Define an effective Photon Transfer Curve dataset
1878 with nominal gains and noise.
1880 Parameters
1881 ----------
1882 ptcDataset : `lsst.ip.isr.PhotonTransferCurveDataset`
1883 Input Photon Transfer Curve dataset.
1884 detector : `lsst.afw.cameraGeom.Detector`
1885 Detector object.
1886 bfGains : `dict`
1887 Gains from running the brighter-fatter code.
1888 A dict keyed by amplifier name for the detector
1889 in question.
1890 overScans : `list` [`lsst.pipe.base.Struct`]
1891 List of overscanResults structures
1892 metadata : `lsst.daf.base.PropertyList`
1893 Exposure metadata to update gain and read noise
1894 provenance.
1896 Returns
1897 -------
1898 effectivePtc : `lsst.ip.isr.PhotonTransferCurveDataset`
1899 PTC dataset containing gains and readout noise
1900 values to be used throughout
1901 Instrument Signature Removal.
1902 """
1903 amps = detector.getAmplifiers()
1904 ampNames = [amp.getName() for amp in amps]
1905 detName = detector.getName()
1906 effectivePtc = PhotonTransferCurveDataset(ampNames, 'EFFECTIVE_PTC', 1)
1907 boolGainMismatch = False
1908 doWarningPtcValidation = True
1910 for amp, overscanResults in zip(amps, overScans):
1911 ampName = amp.getName()
1912 # Gain:
1913 # Try first with the PTC gains.
1914 gainProvenanceString = "amp"
1915 if self.config.usePtcGains:
1916 gain = ptcDataset.gain[ampName]
1917 gainProvenanceString = "ptc"
1918 self.log.debug("Using gain from Photon Transfer Curve.")
1919 else:
1920 # Try then with the amplifier gain.
1921 # We already have a detector at this point. If there was no
1922 # detector to begin with, one would have been created with
1923 # self.config.gain and self.config.noise. Same comment
1924 # applies for the noise block below.
1925 gain = amp.getGain()
1927 # Check if the gain up to this point differs from the
1928 # gain in bfGains. If so, raise or warn, accordingly.
1929 if not boolGainMismatch and bfGains is not None and ampName in bfGains:
1930 bfGain = bfGains[ampName]
1931 if not math.isclose(gain, bfGain, rel_tol=1e-4):
1932 if self.config.doRaiseOnCalibMismatch:
1933 raise RuntimeError("Gain mismatch for det %s amp %s: "
1934 "(gain (%s): %s, bfGain: %s)",
1935 detName, ampName, gainProvenanceString,
1936 gain, bfGain)
1937 else:
1938 self.log.warning("Gain mismatch for det %s amp %s: "
1939 "(gain (%s): %s, bfGain: %s)",
1940 detName, ampName, gainProvenanceString,
1941 gain, bfGain)
1942 boolGainMismatch = True
1944 # Gain:
1945 if math.isnan(gain):
1946 gain = 1.0
1947 self.log.warning("Gain for amp %s set to NaN! Updating to"
1948 " 1.0 to generate Poisson variance.", ampName)
1949 elif gain <= 0:
1950 patchedGain = 1.0
1951 self.log.warning("Gain for amp %s == %g <= 0; setting to %f.",
1952 ampName, gain, patchedGain)
1953 gain = patchedGain
1955 # Noise:
1956 # Try first with the empirical noise from the overscan.
1957 noiseProvenanceString = "amp"
1958 if self.config.doEmpiricalReadNoise and overscanResults is not None:
1959 noiseProvenanceString = "serial overscan"
1960 if isinstance(overscanResults.residualSigma, float):
1961 # Only serial overscan was run
1962 noise = overscanResults.residualSigma
1963 else:
1964 # Both serial and parallel overscan were
1965 # run. Only report noise from serial here.
1966 noise = overscanResults.residualSigma[0]
1968 # Overscan noise is always in adu; we must standardize
1969 # the read noise attributes of all
1970 # PhotonTransferCurveDataset objects to units of
1971 # electrons.
1972 noise *= gain
1973 elif self.config.usePtcReadNoise:
1974 # Try then with the PTC noise.
1975 self.log.debug("Using noise from Photon Transfer Curve.")
1976 noise = ptcDataset.noise[ampName]
1977 noiseProvenanceString = "ptc"
1978 else:
1979 # Finally, try with the amplifier noise.
1980 # We already have a detector at this point. If there was no
1981 # detector to begin with, one would have been created with
1982 # self.config.gain and self.config.noise.
1983 # Amplifier object read noise is always in electron
1984 noise = amp.getReadNoise()
1986 # PTC Turnoff:
1987 # Copy it over from the input PTC if it's positive. If it's a nan
1988 # set it to a high value.
1989 if ptcDataset is not None:
1990 ptcTurnoff = ptcDataset.ptcTurnoff[ampName]
1991 else:
1992 ptcTurnoff = 2e19
1994 if (isinstance(ptcTurnoff, numbers.Real) and ptcTurnoff > 0):
1995 effectivePtc.ptcTurnoff[ampName] = ptcTurnoff
1996 elif math.isnan(ptcTurnoff):
1997 effectivePtc.ptcTurnoff[ampName] = 2e19
1999 effectivePtc.gain[ampName] = gain
2000 effectivePtc.noise[ampName] = noise
2001 # Make sure read noise, turnoff, and gain make sense
2002 effectivePtc.validateGainNoiseTurnoffValues(ampName, doWarn=doWarningPtcValidation)
2003 doWarningPtcValidation = False
2005 # These keys are duplicated for compatability with isrTaskLSST
2006 metadata[f"LSST GAIN {amp.getName()}"] = effectivePtc.gain[ampName]
2007 metadata[f"LSST ISR GAIN {amp.getName()}"] = effectivePtc.gain[ampName]
2008 metadata[f"LSST READNOISE {amp.getName()}"] = effectivePtc.noise[ampName]
2009 metadata[f"LSST ISR READNOISE {amp.getName()}"] = effectivePtc.noise[ampName]
2011 self.log.info("Det: %s - Noise provenance: %s, Gain provenance: %s",
2012 detName,
2013 noiseProvenanceString,
2014 gainProvenanceString)
2015 metadata["LSST ISR GAIN SOURCE"] = gainProvenanceString
2016 metadata["LSST ISR NOISE SOURCE"] = noiseProvenanceString
2017 metadata["LSST ISR READNOISE UNITS"] = "electron"
2019 return effectivePtc
2021 def ensureExposure(self, inputExp, camera=None, detectorNum=None):
2022 """Ensure that the data returned by Butler is a fully constructed exp.
2024 ISR requires exposure-level image data for historical reasons, so if we
2025 did not recieve that from Butler, construct it from what we have,
2026 modifying the input in place.
2028 Parameters
2029 ----------
2030 inputExp : `lsst.afw.image` image-type.
2031 The input data structure obtained from Butler.
2032 Can be `lsst.afw.image.Exposure`,
2033 `lsst.afw.image.DecoratedImageU`,
2034 or `lsst.afw.image.ImageF`
2035 camera : `lsst.afw.cameraGeom.camera`, optional
2036 The camera associated with the image. Used to find the appropriate
2037 detector if detector is not already set.
2038 detectorNum : `int`, optional
2039 The detector in the camera to attach, if the detector is not
2040 already set.
2042 Returns
2043 -------
2044 inputExp : `lsst.afw.image.Exposure`
2045 The re-constructed exposure, with appropriate detector parameters.
2047 Raises
2048 ------
2049 TypeError
2050 Raised if the input data cannot be used to construct an exposure.
2051 """
2052 if isinstance(inputExp, afwImage.DecoratedImageU):
2053 inputExp = afwImage.makeExposure(afwImage.makeMaskedImage(inputExp))
2054 elif isinstance(inputExp, afwImage.ImageF):
2055 inputExp = afwImage.makeExposure(afwImage.makeMaskedImage(inputExp))
2056 elif isinstance(inputExp, afwImage.MaskedImageF):
2057 inputExp = afwImage.makeExposure(inputExp)
2058 elif isinstance(inputExp, afwImage.Exposure):
2059 pass
2060 elif inputExp is None:
2061 # Assume this will be caught by the setup if it is a problem.
2062 return inputExp
2063 else:
2064 raise TypeError("Input Exposure is not known type in isrTask.ensureExposure: %s." %
2065 (type(inputExp), ))
2067 if inputExp.getDetector() is None:
2068 if camera is None or detectorNum is None:
2069 raise RuntimeError('Must supply both a camera and detector number when using exposures '
2070 'without a detector set.')
2071 inputExp.setDetector(camera[detectorNum])
2073 return inputExp
2075 @staticmethod
2076 def extractCalibDate(calib):
2077 """Extract common calibration metadata values that will be written to
2078 output header.
2080 Parameters
2081 ----------
2082 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
2083 Calibration to pull date information from.
2085 Returns
2086 -------
2087 dateString : `str`
2088 Calibration creation date string to add to header.
2089 """
2090 if hasattr(calib, "getMetadata"):
2091 if 'CALIB_CREATION_DATE' in calib.getMetadata():
2092 return " ".join((calib.getMetadata().get("CALIB_CREATION_DATE", "Unknown"),
2093 calib.getMetadata().get("CALIB_CREATION_TIME", "Unknown")))
2094 else:
2095 return " ".join((calib.getMetadata().get("CALIB_CREATE_DATE", "Unknown"),
2096 calib.getMetadata().get("CALIB_CREATE_TIME", "Unknown")))
2097 else:
2098 return "Unknown Unknown"
2100 def compareUnits(self, calibMetadata, calibName):
2101 """Compare units from calibration to ISR units.
2103 For the regular IsrTask this is used to confirm that calibs
2104 suitable for IsrTaskLSST are not used with the old IsrTask.
2106 Parameters
2107 ----------
2108 calibMetadata : `lsst.daf.base.PropertyList`
2109 Calibration metadata from header.
2110 calibName : `str`
2111 Calibration name for log message.
2112 """
2113 calibUnits = calibMetadata.get("LSST ISR UNITS", "adu")
2114 isrUnits = "adu"
2115 if calibUnits != isrUnits:
2116 if self.config.doRaiseOnCalibMismatch:
2117 raise RuntimeError(
2118 "Unit mismatch: isr has %s units but %s has %s units",
2119 isrUnits,
2120 calibName,
2121 calibUnits,
2122 )
2123 else:
2124 self.log.warning(
2125 "Unit mismatch: isr has %s units but %s has %s units",
2126 isrUnits,
2127 calibName,
2128 calibUnits,
2129 )
2131 def convertIntToFloat(self, exposure):
2132 """Convert exposure image from uint16 to float.
2134 If the exposure does not need to be converted, the input is
2135 immediately returned. For exposures that are converted to use
2136 floating point pixels, the variance is set to unity and the
2137 mask to zero.
2139 Parameters
2140 ----------
2141 exposure : `lsst.afw.image.Exposure`
2142 The raw exposure to be converted.
2144 Returns
2145 -------
2146 newexposure : `lsst.afw.image.Exposure`
2147 The input ``exposure``, converted to floating point pixels.
2149 Raises
2150 ------
2151 RuntimeError
2152 Raised if the exposure type cannot be converted to float.
2154 """
2155 if isinstance(exposure, afwImage.ExposureF):
2156 # Nothing to be done
2157 self.log.debug("Exposure already of type float.")
2158 return exposure
2159 if not hasattr(exposure, "convertF"):
2160 raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
2162 newexposure = exposure.convertF()
2163 newexposure.variance[:] = 1
2164 newexposure.mask[:] = 0x0
2166 return newexposure
2168 def maskAmplifier(self, ccdExposure, amp, defects):
2169 """Identify bad amplifiers, saturated and suspect pixels.
2171 Parameters
2172 ----------
2173 ccdExposure : `lsst.afw.image.Exposure`
2174 Input exposure to be masked.
2175 amp : `lsst.afw.cameraGeom.Amplifier`
2176 Catalog of parameters defining the amplifier on this
2177 exposure to mask.
2178 defects : `lsst.ip.isr.Defects`
2179 List of defects. Used to determine if the entire
2180 amplifier is bad.
2182 Returns
2183 -------
2184 badAmp : `Bool`
2185 If this is true, the entire amplifier area is covered by
2186 defects and unusable.
2188 """
2189 maskedImage = ccdExposure.getMaskedImage()
2191 badAmp = False
2193 # Check if entire amp region is defined as a defect
2194 # NB: need to use amp.getBBox() for correct comparison with current
2195 # defects definition.
2196 if defects is not None:
2197 badAmp = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
2199 # In the case of a bad amp, we will set mask to "BAD"
2200 # (here use amp.getRawBBox() for correct association with pixels in
2201 # current ccdExposure).
2202 if badAmp:
2203 dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
2204 afwImage.PARENT)
2205 maskView = dataView.getMask()
2206 maskView |= maskView.getPlaneBitMask("BAD")
2207 del maskView
2208 return badAmp
2210 # Mask remaining defects after assembleCcd() to allow for defects that
2211 # cross amplifier boundaries. Saturation and suspect pixels can be
2212 # masked now, though.
2213 limits = dict()
2214 if self.config.doSaturation and not badAmp:
2215 # Set to the default from the camera model.
2216 limits.update({self.config.saturatedMaskName: amp.getSaturation()})
2217 # And update if it is set in the config.
2218 if math.isfinite(self.config.saturation):
2219 limits.update({self.config.saturatedMaskName: self.config.saturation})
2220 if self.config.doSuspect and not badAmp:
2221 limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
2223 for maskName, maskThreshold in limits.items():
2224 if not math.isnan(maskThreshold):
2225 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2226 isrFunctions.makeThresholdMask(
2227 maskedImage=dataView,
2228 threshold=maskThreshold,
2229 growFootprints=0,
2230 maskName=maskName
2231 )
2233 # Determine if we've fully masked this amplifier with SUSPECT and
2234 # SAT pixels.
2235 maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
2236 afwImage.PARENT)
2237 maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
2238 self.config.suspectMaskName])
2239 if numpy.all(maskView.getArray() & maskVal > 0):
2240 badAmp = True
2241 maskView |= maskView.getPlaneBitMask("BAD")
2243 return badAmp
2245 def overscanCorrection(self, ccdExposure, amp):
2246 """Apply overscan correction in place.
2248 This method does initial pixel rejection of the overscan
2249 region. The overscan can also be optionally segmented to
2250 allow for discontinuous overscan responses to be fit
2251 separately. The actual overscan subtraction is performed by
2252 the `lsst.ip.isr.overscan.OverscanTask`, which is called here
2253 after the amplifier is preprocessed.
2255 Parameters
2256 ----------
2257 ccdExposure : `lsst.afw.image.Exposure`
2258 Exposure to have overscan correction performed.
2259 amp : `lsst.afw.cameraGeom.Amplifer`
2260 The amplifier to consider while correcting the overscan.
2262 Returns
2263 -------
2264 overscanResults : `lsst.pipe.base.Struct`
2265 Result struct with components:
2267 ``imageFit``
2268 Value or fit subtracted from the amplifier image data.
2269 (scalar or `lsst.afw.image.Image`)
2270 ``overscanFit``
2271 Value or fit subtracted from the overscan image data.
2272 (scalar or `lsst.afw.image.Image`)
2273 ``overscanImage``
2274 Image of the overscan region with the overscan
2275 correction applied. This quantity is used to estimate
2276 the amplifier read noise empirically.
2277 (`lsst.afw.image.Image`)
2278 ``edgeMask``
2279 Mask of the suspect pixels. (`lsst.afw.image.Mask`)
2280 ``overscanMean``
2281 Median overscan fit value. (`float`)
2282 ``overscanSigma``
2283 Clipped standard deviation of the overscan after
2284 correction. (`float`)
2286 Raises
2287 ------
2288 RuntimeError
2289 Raised if the ``amp`` does not contain raw pixel information.
2291 See Also
2292 --------
2293 lsst.ip.isr.overscan.OverscanTask
2294 """
2295 if amp.getRawHorizontalOverscanBBox().isEmpty():
2296 self.log.info("ISR_OSCAN: No overscan region. Not performing overscan correction.")
2297 return None
2299 # Perform overscan correction on subregions.
2300 overscanResults = self.overscan.run(ccdExposure, amp)
2302 metadata = ccdExposure.getMetadata()
2303 ampName = amp.getName()
2305 keyBase = "LSST ISR OVERSCAN"
2306 # The overscan statistics units will always match the units of
2307 # the image at the point they are calculated.
2308 metadata[f"{keyBase} SERIAL UNITS"] = metadata.get("LSST ISR UNITS")
2310 # Updated quantities
2311 if isinstance(overscanResults.overscanMean, float):
2312 # Serial overscan correction only:
2313 metadata[f"{keyBase} SERIAL MEAN {ampName}"] = overscanResults.overscanMean
2314 metadata[f"{keyBase} SERIAL MEDIAN {ampName}"] = overscanResults.overscanMedian
2315 metadata[f"{keyBase} SERIAL STDEV {ampName}"] = overscanResults.overscanSigma
2317 metadata[f"{keyBase} RESIDUAL SERIAL MEAN {ampName}"] = overscanResults.residualMean
2318 metadata[f"{keyBase} RESIDUAL SERIAL MEDIAN {ampName}"] = overscanResults.residualMedian
2319 metadata[f"{keyBase} RESIDUAL SERIAL STDEV {ampName}"] = overscanResults.residualSigma
2320 elif isinstance(overscanResults.overscanMean, tuple):
2321 # Both serial and parallel overscan have run:
2322 metadata[f"{keyBase} PARALLEL UNITS"] = metadata.get("LSST ISR UNITS")
2323 metadata[f"{keyBase} SERIAL MEAN {ampName}"] = overscanResults.overscanMean[0]
2324 metadata[f"{keyBase} SERIAL MEDIAN {ampName}"] = overscanResults.overscanMedian[0]
2325 metadata[f"{keyBase} SERIAL STDEV {ampName}"] = overscanResults.overscanSigma[0]
2327 metadata[f"{keyBase} PARALLEL MEAN {ampName}"] = overscanResults.overscanMean[1]
2328 metadata[f"{keyBase} PARALLEL MEDIAN {ampName}"] = overscanResults.overscanMedian[1]
2329 metadata[f"{keyBase} PARALLEL STDEV {ampName}"] = overscanResults.overscanSigma[1]
2331 metadata[f"{keyBase} RESIDUAL SERIAL MEAN {ampName}"] = overscanResults.residualMean[0]
2332 metadata[f"{keyBase} RESIDUAL SERIAL MEDIAN {ampName}"] = overscanResults.residualMedian[0]
2333 metadata[f"{keyBase} RESIDUAL SERIAL STDEV {ampName}"] = overscanResults.residualSigma[0]
2335 metadata[f"{keyBase} RESIDUAL PARALLEL MEAN {ampName}"] = overscanResults.residualMean[1]
2336 metadata[f"{keyBase} RESIDUAL PARALLEL MEDIAN {ampName}"] = overscanResults.residualMedian[1]
2337 metadata[f"{keyBase} RESIDUAL PARALLEL STDEV {ampName}"] = overscanResults.residualSigma[1]
2338 else:
2339 self.log.warning("Unexpected type for overscan values; none added to header.")
2341 return overscanResults
2343 def updateVariance(self, ampExposure, amp, ptcDataset):
2344 """Set the variance plane using the gain and read noise
2346 The read noise is calculated from the ``overscanImage`` if the
2347 ``doEmpiricalReadNoise`` option is set in the configuration; otherwise
2348 the value from the amplifier data is used.
2350 Parameters
2351 ----------
2352 ampExposure : `lsst.afw.image.Exposure`
2353 Exposure to process.
2354 amp : `lsst.afw.cameraGeom.Amplifier` or `FakeAmp`
2355 Amplifier detector data.
2356 ptcDataset : `lsst.ip.isr.PhotonTransferCurveDataset`
2357 Effective PTC dataset containing the gains and read noise.
2359 See also
2360 --------
2361 lsst.ip.isr.isrFunctions.updateVariance
2362 """
2363 ampName = amp.getName()
2364 # At this point, the effective PTC should have
2365 # gain and read noise values.
2366 gain = ptcDataset.gain[ampName]
2367 readNoise = ptcDataset.noise[ampName]
2369 # The image will always be in adu and the noise
2370 # will always be in electrons, and we will output
2371 # the variance plane in the image units (adu^2).
2372 isrFunctions.updateVariance(
2373 maskedImage=ampExposure.getMaskedImage(),
2374 gain=gain,
2375 readNoise=readNoise,
2376 )
2378 def maskNegativeVariance(self, exposure):
2379 """Identify and mask pixels with negative variance values.
2381 Parameters
2382 ----------
2383 exposure : `lsst.afw.image.Exposure`
2384 Exposure to process.
2386 See Also
2387 --------
2388 lsst.ip.isr.isrFunctions.updateVariance
2389 """
2390 maskPlane = exposure.getMask().getPlaneBitMask(self.config.negativeVarianceMaskName)
2391 bad = numpy.where(exposure.getVariance().getArray() <= 0.0)
2392 exposure.mask.array[bad] |= maskPlane
2394 def darkCorrection(self, exposure, darkExposure, invert=False):
2395 """Apply dark correction in place.
2397 Parameters
2398 ----------
2399 exposure : `lsst.afw.image.Exposure`
2400 Exposure to process.
2401 darkExposure : `lsst.afw.image.Exposure`
2402 Dark exposure of the same size as ``exposure``.
2403 invert : `Bool`, optional
2404 If True, re-add the dark to an already corrected image.
2406 Raises
2407 ------
2408 RuntimeError
2409 Raised if either ``exposure`` or ``darkExposure`` do not
2410 have their dark time defined.
2412 See Also
2413 --------
2414 lsst.ip.isr.isrFunctions.darkCorrection
2415 """
2416 expScale = exposure.getInfo().getVisitInfo().getDarkTime()
2417 if math.isnan(expScale):
2418 raise RuntimeError("Exposure darktime is NaN.")
2419 if darkExposure.getInfo().getVisitInfo() is not None \
2420 and not math.isnan(darkExposure.getInfo().getVisitInfo().getDarkTime()):
2421 darkScale = darkExposure.getInfo().getVisitInfo().getDarkTime()
2422 else:
2423 # DM-17444: darkExposure.getInfo.getVisitInfo() is None
2424 # so getDarkTime() does not exist.
2425 self.log.warning("darkExposure.getInfo().getVisitInfo() does not exist. Using darkScale = 1.0.")
2426 darkScale = 1.0
2428 isrFunctions.darkCorrection(
2429 maskedImage=exposure.getMaskedImage(),
2430 darkMaskedImage=darkExposure.getMaskedImage(),
2431 expScale=expScale,
2432 darkScale=darkScale,
2433 invert=invert,
2434 trimToFit=self.config.doTrimToMatchCalib
2435 )
2437 def doLinearize(self, detector):
2438 """Check if linearization is needed for the detector cameraGeom.
2440 Checks config.doLinearize and the linearity type of the first
2441 amplifier.
2443 Parameters
2444 ----------
2445 detector : `lsst.afw.cameraGeom.Detector`
2446 Detector to get linearity type from.
2448 Returns
2449 -------
2450 doLinearize : `Bool`
2451 If True, linearization should be performed.
2452 """
2453 return self.config.doLinearize and \
2454 detector.getAmplifiers()[0].getLinearityType() != NullLinearityType
2456 def flatCorrection(self, exposure, flatExposure, invert=False):
2457 """Apply flat correction in place.
2459 Parameters
2460 ----------
2461 exposure : `lsst.afw.image.Exposure`
2462 Exposure to process.
2463 flatExposure : `lsst.afw.image.Exposure`
2464 Flat exposure of the same size as ``exposure``.
2465 invert : `Bool`, optional
2466 If True, unflatten an already flattened image.
2468 See Also
2469 --------
2470 lsst.ip.isr.isrFunctions.flatCorrection
2471 """
2472 isrFunctions.flatCorrection(
2473 maskedImage=exposure.getMaskedImage(),
2474 flatMaskedImage=flatExposure.getMaskedImage(),
2475 scalingType=self.config.flatScalingType,
2476 userScale=self.config.flatUserScale,
2477 invert=invert,
2478 trimToFit=self.config.doTrimToMatchCalib
2479 )
2481 def saturationDetection(self, exposure, amp):
2482 """Detect and mask saturated pixels in config.saturatedMaskName.
2484 Parameters
2485 ----------
2486 exposure : `lsst.afw.image.Exposure`
2487 Exposure to process. Only the amplifier DataSec is processed.
2488 amp : `lsst.afw.cameraGeom.Amplifier`
2489 Amplifier detector data.
2491 See Also
2492 --------
2493 lsst.ip.isr.isrFunctions.makeThresholdMask
2494 """
2495 if not math.isnan(amp.getSaturation()):
2496 maskedImage = exposure.getMaskedImage()
2497 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2498 isrFunctions.makeThresholdMask(
2499 maskedImage=dataView,
2500 threshold=amp.getSaturation(),
2501 growFootprints=0,
2502 maskName=self.config.saturatedMaskName,
2503 )
2505 def saturationInterpolation(self, exposure):
2506 """Interpolate over saturated pixels, in place.
2508 This method should be called after `saturationDetection`, to
2509 ensure that the saturated pixels have been identified in the
2510 SAT mask. It should also be called after `assembleCcd`, since
2511 saturated regions may cross amplifier boundaries.
2513 Parameters
2514 ----------
2515 exposure : `lsst.afw.image.Exposure`
2516 Exposure to process.
2518 See Also
2519 --------
2520 lsst.ip.isr.isrTask.saturationDetection
2521 lsst.ip.isr.isrFunctions.interpolateFromMask
2522 """
2523 isrFunctions.interpolateFromMask(
2524 maskedImage=exposure.getMaskedImage(),
2525 fwhm=self.config.fwhm,
2526 growSaturatedFootprints=self.config.growSaturationFootprintSize,
2527 maskNameList=list(self.config.saturatedMaskName),
2528 useLegacyInterp=self.config.useLegacyInterp,
2529 )
2531 def suspectDetection(self, exposure, amp):
2532 """Detect and mask suspect pixels in config.suspectMaskName.
2534 Parameters
2535 ----------
2536 exposure : `lsst.afw.image.Exposure`
2537 Exposure to process. Only the amplifier DataSec is processed.
2538 amp : `lsst.afw.cameraGeom.Amplifier`
2539 Amplifier detector data.
2541 See Also
2542 --------
2543 lsst.ip.isr.isrFunctions.makeThresholdMask
2545 Notes
2546 -----
2547 Suspect pixels are pixels whose value is greater than
2548 amp.getSuspectLevel(). This is intended to indicate pixels that may be
2549 affected by unknown systematics; for example if non-linearity
2550 corrections above a certain level are unstable then that would be a
2551 useful value for suspectLevel. A value of `nan` indicates that no such
2552 level exists and no pixels are to be masked as suspicious.
2553 """
2554 suspectLevel = amp.getSuspectLevel()
2555 if math.isnan(suspectLevel):
2556 return
2558 maskedImage = exposure.getMaskedImage()
2559 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2560 isrFunctions.makeThresholdMask(
2561 maskedImage=dataView,
2562 threshold=suspectLevel,
2563 growFootprints=0,
2564 maskName=self.config.suspectMaskName,
2565 )
2567 def maskDefect(self, exposure, defectBaseList):
2568 """Mask defects using mask plane "BAD", in place.
2570 Parameters
2571 ----------
2572 exposure : `lsst.afw.image.Exposure`
2573 Exposure to process.
2574 defectBaseList : defect-type
2575 List of defects to mask. Can be of type `lsst.ip.isr.Defects`
2576 or `list` of `lsst.afw.image.DefectBase`.
2578 Notes
2579 -----
2580 Call this after CCD assembly, since defects may cross amplifier
2581 boundaries.
2582 """
2583 maskedImage = exposure.getMaskedImage()
2584 if not isinstance(defectBaseList, Defects):
2585 # Promotes DefectBase to Defect
2586 defectList = Defects(defectBaseList)
2587 else:
2588 defectList = defectBaseList
2589 defectList.maskPixels(maskedImage, maskName="BAD")
2591 def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
2592 """Mask edge pixels with applicable mask plane.
2594 Parameters
2595 ----------
2596 exposure : `lsst.afw.image.Exposure`
2597 Exposure to process.
2598 numEdgePixels : `int`, optional
2599 Number of edge pixels to mask.
2600 maskPlane : `str`, optional
2601 Mask plane name to use.
2602 level : `str`, optional
2603 Level at which to mask edges.
2604 """
2605 maskedImage = exposure.getMaskedImage()
2606 maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
2608 if numEdgePixels > 0:
2609 if level == 'DETECTOR':
2610 boxes = [maskedImage.getBBox()]
2611 elif level == 'AMP':
2612 boxes = [amp.getBBox() for amp in exposure.getDetector()]
2614 for box in boxes:
2615 # This makes a bbox numEdgeSuspect pixels smaller than the
2616 # image on each side
2617 subImage = maskedImage[box]
2618 box.grow(-numEdgePixels)
2619 # Mask pixels outside box
2620 SourceDetectionTask.setEdgeBits(
2621 subImage,
2622 box,
2623 maskBitMask)
2625 def maskAndInterpolateDefects(self, exposure, defectBaseList):
2626 """Mask and interpolate defects using mask plane "BAD", in place.
2628 Parameters
2629 ----------
2630 exposure : `lsst.afw.image.Exposure`
2631 Exposure to process.
2632 defectBaseList : defects-like
2633 List of defects to mask and interpolate. Can be
2634 `lsst.ip.isr.Defects` or `list` of `lsst.afw.image.DefectBase`.
2636 See Also
2637 --------
2638 lsst.ip.isr.isrTask.maskDefect
2639 """
2640 self.maskDefect(exposure, defectBaseList)
2641 self.maskEdges(exposure, numEdgePixels=self.config.numEdgeSuspect,
2642 maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
2643 isrFunctions.interpolateFromMask(
2644 maskedImage=exposure.getMaskedImage(),
2645 fwhm=self.config.fwhm,
2646 growSaturatedFootprints=0,
2647 maskNameList=["BAD"],
2648 useLegacyInterp=self.config.useLegacyInterp,
2649 )
2651 def maskNan(self, exposure):
2652 """Mask NaNs using mask plane "UNMASKEDNAN", in place.
2654 Parameters
2655 ----------
2656 exposure : `lsst.afw.image.Exposure`
2657 Exposure to process.
2659 Notes
2660 -----
2661 We mask over all non-finite values (NaN, inf), including those
2662 that are masked with other bits (because those may or may not be
2663 interpolated over later, and we want to remove all NaN/infs).
2664 Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
2665 preserve the historical name.
2666 """
2667 maskedImage = exposure.getMaskedImage()
2669 # Find and mask NaNs
2670 maskedImage.getMask().addMaskPlane("UNMASKEDNAN")
2671 maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
2672 numNans = maskNans(maskedImage, maskVal)
2673 self.metadata["NUMNANS"] = numNans
2674 if numNans > 0:
2675 self.log.warning("There were %d unmasked NaNs.", numNans)
2677 def maskAndInterpolateNan(self, exposure):
2678 """"Mask and interpolate NaN/infs using mask plane "UNMASKEDNAN",
2679 in place.
2681 Parameters
2682 ----------
2683 exposure : `lsst.afw.image.Exposure`
2684 Exposure to process.
2686 See Also
2687 --------
2688 lsst.ip.isr.isrTask.maskNan
2689 """
2690 self.maskNan(exposure)
2691 isrFunctions.interpolateFromMask(
2692 maskedImage=exposure.getMaskedImage(),
2693 fwhm=self.config.fwhm,
2694 growSaturatedFootprints=0,
2695 maskNameList=["UNMASKEDNAN"],
2696 useLegacyInterp=self.config.useLegacyInterp,
2697 )
2699 def measureBackground(self, exposure, IsrQaConfig=None):
2700 """Measure the image background in subgrids, for quality control.
2702 Parameters
2703 ----------
2704 exposure : `lsst.afw.image.Exposure`
2705 Exposure to process.
2706 IsrQaConfig : `lsst.ip.isr.isrQa.IsrQaConfig`
2707 Configuration object containing parameters on which background
2708 statistics and subgrids to use.
2709 """
2710 if IsrQaConfig is not None:
2711 statsControl = afwMath.StatisticsControl(IsrQaConfig.flatness.clipSigma,
2712 IsrQaConfig.flatness.nIter)
2713 maskVal = exposure.getMaskedImage().getMask().getPlaneBitMask(["BAD", "SAT", "DETECTED"])
2714 statsControl.setAndMask(maskVal)
2715 maskedImage = exposure.getMaskedImage()
2716 stats = afwMath.makeStatistics(maskedImage, afwMath.MEDIAN | afwMath.STDEVCLIP, statsControl)
2717 skyLevel = stats.getValue(afwMath.MEDIAN)
2718 skySigma = stats.getValue(afwMath.STDEVCLIP)
2719 self.log.info("Flattened sky level: %f +/- %f.", skyLevel, skySigma)
2720 metadata = exposure.getMetadata()
2721 metadata["SKYLEVEL"] = skyLevel
2722 metadata["SKYSIGMA"] = skySigma
2724 # calcluating flatlevel over the subgrids
2725 stat = afwMath.MEANCLIP if IsrQaConfig.flatness.doClip else afwMath.MEAN
2726 meshXHalf = int(IsrQaConfig.flatness.meshX/2.)
2727 meshYHalf = int(IsrQaConfig.flatness.meshY/2.)
2728 nX = int((exposure.getWidth() + meshXHalf) / IsrQaConfig.flatness.meshX)
2729 nY = int((exposure.getHeight() + meshYHalf) / IsrQaConfig.flatness.meshY)
2730 skyLevels = numpy.zeros((nX, nY))
2732 for j in range(nY):
2733 yc = meshYHalf + j * IsrQaConfig.flatness.meshY
2734 for i in range(nX):
2735 xc = meshXHalf + i * IsrQaConfig.flatness.meshX
2737 xLLC = xc - meshXHalf
2738 yLLC = yc - meshYHalf
2739 xURC = xc + meshXHalf - 1
2740 yURC = yc + meshYHalf - 1
2742 bbox = lsst.geom.Box2I(lsst.geom.Point2I(xLLC, yLLC), lsst.geom.Point2I(xURC, yURC))
2743 miMesh = maskedImage.Factory(exposure.getMaskedImage(), bbox, afwImage.LOCAL)
2745 skyLevels[i, j] = afwMath.makeStatistics(miMesh, stat, statsControl).getValue()
2747 good = numpy.where(numpy.isfinite(skyLevels))
2748 if len(good[0]) == 0:
2749 # There are no good pixels.
2750 self.log.warning("No good pixels to measure sky levels.")
2751 skyMedian = numpy.nan
2752 flatness = numpy.nan
2753 flatness_rms = numpy.nan
2754 flatness_pp = numpy.nan
2755 else:
2756 skyMedian = numpy.median(skyLevels[good])
2757 flatness = (skyLevels[good] - skyMedian) / skyMedian
2758 flatness_rms = numpy.std(flatness)
2759 flatness_pp = flatness.max() - flatness.min() if len(flatness) > 0 else numpy.nan
2761 self.log.info("Measuring sky levels in %dx%d grids: %f.", nX, nY, skyMedian)
2762 self.log.info("Sky flatness in %dx%d grids - pp: %f rms: %f.",
2763 nX, nY, flatness_pp, flatness_rms)
2765 metadata["FLATNESS_PP"] = float(flatness_pp)
2766 metadata["FLATNESS_RMS"] = float(flatness_rms)
2767 metadata["FLATNESS_NGRIDS"] = '%dx%d' % (nX, nY)
2768 metadata["FLATNESS_MESHX"] = IsrQaConfig.flatness.meshX
2769 metadata["FLATNESS_MESHY"] = IsrQaConfig.flatness.meshY
2771 def roughZeroPoint(self, exposure):
2772 """Set an approximate magnitude zero point for the exposure.
2774 Parameters
2775 ----------
2776 exposure : `lsst.afw.image.Exposure`
2777 Exposure to process.
2778 """
2779 filterLabel = exposure.getFilter()
2780 physicalFilter = isrFunctions.getPhysicalFilter(filterLabel, self.log)
2782 if physicalFilter in self.config.fluxMag0T1:
2783 fluxMag0 = self.config.fluxMag0T1[physicalFilter]
2784 else:
2785 self.log.warning("No rough magnitude zero point defined for filter %s.", physicalFilter)
2786 fluxMag0 = self.config.defaultFluxMag0T1
2788 expTime = exposure.getInfo().getVisitInfo().getExposureTime()
2789 if expTime == 0.0:
2790 self.log.debug("Received exposure with 0.0 expTime; skipping rough zero point.")
2791 return
2792 elif not expTime > 0: # handle NaN as well as <= 0
2793 self.log.warning("Non-positive exposure time; skipping rough zero point.")
2794 return
2796 self.log.info("Setting rough magnitude zero point for filter %s: %f",
2797 physicalFilter, 2.5*math.log10(fluxMag0*expTime))
2798 exposure.setPhotoCalib(afwImage.makePhotoCalibFromCalibZeroPoint(fluxMag0*expTime, 0.0))
2800 @contextmanager
2801 def flatContext(self, exp, flat, dark=None):
2802 """Context manager that applies and removes flats and darks,
2803 if the task is configured to apply them.
2805 Parameters
2806 ----------
2807 exp : `lsst.afw.image.Exposure`
2808 Exposure to process.
2809 flat : `lsst.afw.image.Exposure`
2810 Flat exposure the same size as ``exp``.
2811 dark : `lsst.afw.image.Exposure`, optional
2812 Dark exposure the same size as ``exp``.
2814 Yields
2815 ------
2816 exp : `lsst.afw.image.Exposure`
2817 The flat and dark corrected exposure.
2818 """
2819 if self.config.doDark and dark is not None:
2820 self.darkCorrection(exp, dark)
2821 if self.config.doFlat:
2822 self.flatCorrection(exp, flat)
2823 try:
2824 yield exp
2825 finally:
2826 if self.config.doFlat:
2827 self.flatCorrection(exp, flat, invert=True)
2828 if self.config.doDark and dark is not None:
2829 self.darkCorrection(exp, dark, invert=True)
2831 @deprecated(
2832 reason=(
2833 "makeBinnedImages is no longer used. "
2834 "Please subtask lsst.ip.isr.BinImageDataTask instead."
2835 ),
2836 version="v28", category=FutureWarning
2837 )
2838 def makeBinnedImages(self, exposure):
2839 """Make visualizeVisit style binned exposures.
2841 Parameters
2842 ----------
2843 exposure : `lsst.afw.image.Exposure`
2844 Exposure to bin.
2846 Returns
2847 -------
2848 bin1 : `lsst.afw.image.Exposure`
2849 Binned exposure using binFactor1.
2850 bin2 : `lsst.afw.image.Exposure`
2851 Binned exposure using binFactor2.
2852 """
2853 mi = exposure.getMaskedImage()
2855 bin1 = afwMath.binImage(mi, self.config.binFactor1)
2856 bin2 = afwMath.binImage(mi, self.config.binFactor2)
2858 bin1 = afwImage.makeExposure(bin1)
2859 bin2 = afwImage.makeExposure(bin2)
2861 bin1.setInfo(exposure.getInfo())
2862 bin2.setInfo(exposure.getInfo())
2864 return bin1, bin2
2866 def debugView(self, exposure, stepname):
2867 """Utility function to examine ISR exposure at different stages.
2869 Parameters
2870 ----------
2871 exposure : `lsst.afw.image.Exposure`
2872 Exposure to view.
2873 stepname : `str`
2874 State of processing to view.
2875 """
2876 frame = getDebugFrame(self._display, stepname)
2877 if frame:
2878 display = getDisplay(frame)
2879 display.scale('asinh', 'zscale')
2880 display.mtv(exposure)
2881 prompt = "Press Enter to continue [c]... "
2882 while True:
2883 ans = input(prompt).lower()
2884 if ans in ("", "c",):
2885 break
2888class FakeAmp(object):
2889 """A Detector-like object that supports returning gain and saturation level
2891 This is used when the input exposure does not have a detector.
2893 Parameters
2894 ----------
2895 exposure : `lsst.afw.image.Exposure`
2896 Exposure to generate a fake amplifier for.
2897 config : `lsst.ip.isr.isrTaskConfig`
2898 Configuration to apply to the fake amplifier.
2899 """
2901 def __init__(self, exposure, config):
2902 self._bbox = exposure.getBBox(afwImage.LOCAL)
2903 self._RawHorizontalOverscanBBox = lsst.geom.Box2I()
2904 self._gain = config.gain
2905 self._readNoise = config.readNoise
2906 self._saturation = config.saturation
2908 def getBBox(self):
2909 return self._bbox
2911 def getRawBBox(self):
2912 return self._bbox
2914 def getRawHorizontalOverscanBBox(self):
2915 return self._RawHorizontalOverscanBBox
2917 def getGain(self):
2918 return self._gain
2920 def getReadNoise(self):
2921 return self._readNoise
2923 def getSaturation(self):
2924 return self._saturation
2926 def getSuspectLevel(self):
2927 return float("NaN")