Coverage for python/lsst/ip/isr/isrTaskLSST.py: 13%
812 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-07 02:17 -0800
1__all__ = ["IsrTaskLSST", "IsrTaskLSSTConfig"]
3import numpy
4import math
6from . import isrFunctions
7from . import isrQa
8from .defects import Defects
10from contextlib import contextmanager
11from deprecated.sphinx import deprecated
12import lsst.pex.config as pexConfig
13import lsst.afw.math as afwMath
14import lsst.pipe.base as pipeBase
15import lsst.afw.image as afwImage
16import lsst.pipe.base.connectionTypes as cT
17from lsst.pipe.base import UnprocessableDataError
18from lsst.meas.algorithms.detection import SourceDetectionTask
19import lsst.afw.detection as afwDetection
21from .ampOffset import AmpOffsetTask
22from .binImageDataTask import BinImageDataTask
23from .overscan import SerialOverscanCorrectionTask, ParallelOverscanCorrectionTask
24from .overscanAmpConfig import OverscanCameraConfig
25from .assembleCcdTask import AssembleCcdTask
26from .deferredCharge import DeferredChargeTask
27from .crosstalk import CrosstalkTask
28from .masking import MaskingTask
29from .isrStatistics import IsrStatisticsTask
30from .isr import maskNans
31from .ptcDataset import PhotonTransferCurveDataset
32from .isrFunctions import isTrimmedExposure, compareCameraKeywords
35class IsrTaskLSSTConnections(pipeBase.PipelineTaskConnections,
36 dimensions={"instrument", "exposure", "detector"},
37 defaultTemplates={}):
38 ccdExposure = cT.Input(
39 name="raw",
40 doc="Input exposure to process.",
41 storageClass="Exposure",
42 dimensions=["instrument", "exposure", "detector"],
43 )
44 camera = cT.PrerequisiteInput(
45 name="camera",
46 storageClass="Camera",
47 doc="Input camera to construct complete exposures.",
48 dimensions=["instrument"],
49 isCalibration=True,
50 )
51 dnlLUT = cT.PrerequisiteInput(
52 name="dnlLUT",
53 doc="Look-up table for differential non-linearity.",
54 storageClass="IsrCalib",
55 dimensions=["instrument", "exposure", "detector"],
56 isCalibration=True,
57 # TODO DM 36636
58 )
59 bias = cT.PrerequisiteInput(
60 name="bias",
61 doc="Input bias calibration.",
62 storageClass="ExposureF",
63 dimensions=["instrument", "detector"],
64 isCalibration=True,
65 )
66 deferredChargeCalib = cT.PrerequisiteInput(
67 name="cti",
68 doc="Deferred charge/CTI correction dataset.",
69 storageClass="IsrCalib",
70 dimensions=["instrument", "detector"],
71 isCalibration=True,
72 )
73 linearizer = cT.PrerequisiteInput(
74 name='linearizer',
75 storageClass="Linearizer",
76 doc="Linearity correction calibration.",
77 dimensions=["instrument", "detector"],
78 isCalibration=True,
79 )
80 ptc = cT.PrerequisiteInput(
81 name="ptc",
82 doc="Input Photon Transfer Curve dataset",
83 storageClass="PhotonTransferCurveDataset",
84 dimensions=["instrument", "detector"],
85 isCalibration=True,
86 )
87 crosstalk = cT.PrerequisiteInput(
88 name="crosstalk",
89 doc="Input crosstalk object",
90 storageClass="CrosstalkCalib",
91 dimensions=["instrument", "detector"],
92 isCalibration=True,
93 )
94 defects = cT.PrerequisiteInput(
95 name='defects',
96 doc="Input defect tables.",
97 storageClass="Defects",
98 dimensions=["instrument", "detector"],
99 isCalibration=True,
100 )
101 bfKernel = cT.PrerequisiteInput(
102 name="bfk",
103 doc="Complete kernel + gain solutions.",
104 storageClass="BrighterFatterKernel",
105 dimensions=["instrument", "detector"],
106 isCalibration=True,
107 )
108 dark = cT.PrerequisiteInput(
109 name='dark',
110 doc="Input dark calibration.",
111 storageClass="ExposureF",
112 dimensions=["instrument", "detector"],
113 isCalibration=True,
114 )
115 flat = cT.PrerequisiteInput(
116 name="flat",
117 doc="Input flat calibration.",
118 storageClass="ExposureF",
119 dimensions=["instrument", "detector", "physical_filter"],
120 isCalibration=True,
121 )
122 outputExposure = cT.Output(
123 name='postISRCCD',
124 doc="Output ISR processed exposure.",
125 storageClass="Exposure",
126 dimensions=["instrument", "exposure", "detector"],
127 )
128 preInterpExposure = cT.Output(
129 name='preInterpISRCCD',
130 doc="Output ISR processed exposure, with pixels left uninterpolated.",
131 storageClass="ExposureF",
132 dimensions=["instrument", "exposure", "detector"],
133 )
134 outputBin1Exposure = cT.Output(
135 name="postIsrBin1",
136 doc="First binned image.",
137 storageClass="ExposureF",
138 dimensions=["instrument", "exposure", "detector"],
139 )
140 outputBin2Exposure = cT.Output(
141 name="postIsrBin2",
142 doc="Second binned image.",
143 storageClass="ExposureF",
144 dimensions=["instrument", "exposure", "detector"],
145 )
147 outputStatistics = cT.Output(
148 name="isrStatistics",
149 doc="Output of additional statistics table.",
150 storageClass="StructuredDataDict",
151 dimensions=["instrument", "exposure", "detector"],
152 )
154 def __init__(self, *, config=None):
155 super().__init__(config=config)
157 if config.doBootstrap:
158 del self.ptc
159 if config.doDiffNonLinearCorrection is not True:
160 del self.dnlLUT
161 if config.doBias is not True:
162 del self.bias
163 if config.doDeferredCharge is not True:
164 del self.deferredChargeCalib
165 if config.doLinearize is not True:
166 del self.linearizer
167 if not config.doCrosstalk:
168 del self.crosstalk
169 if config.doDefect is not True:
170 del self.defects
171 if config.doBrighterFatter is not True:
172 del self.bfKernel
173 if config.doDark is not True:
174 del self.dark
175 if config.doFlat is not True:
176 del self.flat
178 if config.doBinnedExposures is not True:
179 del self.outputBin1Exposure
180 del self.outputBin2Exposure
181 if config.doSaveInterpPixels is not True:
182 del self.preInterpExposure
184 if config.doCalculateStatistics is not True:
185 del self.outputStatistics
188class IsrTaskLSSTConfig(pipeBase.PipelineTaskConfig,
189 pipelineConnections=IsrTaskLSSTConnections):
190 """Configuration parameters for IsrTaskLSST.
192 Items are grouped in the order in which they are executed by the task.
193 """
194 expectWcs = pexConfig.Field(
195 dtype=bool,
196 default=True,
197 doc="Expect input science images to have a WCS (set False for e.g. spectrographs)."
198 )
199 qa = pexConfig.ConfigField(
200 dtype=isrQa.IsrQaConfig,
201 doc="QA related configuration options.",
202 )
203 doHeaderProvenance = pexConfig.Field(
204 dtype=bool,
205 default=True,
206 doc="Write calibration identifiers into output exposure header.",
207 )
209 # Calib checking configuration:
210 doRaiseOnCalibMismatch = pexConfig.Field(
211 dtype=bool,
212 default=False,
213 doc="Should IsrTaskLSST halt if exposure and calibration header values do not match?",
214 )
215 cameraKeywordsToCompare = pexConfig.ListField(
216 dtype=str,
217 doc="List of header keywords to compare between exposure and calibrations.",
218 default=[],
219 )
221 # Differential non-linearity correction.
222 doDiffNonLinearCorrection = pexConfig.Field(
223 dtype=bool,
224 doc="Do differential non-linearity correction?",
225 default=False,
226 )
228 doBootstrap = pexConfig.Field(
229 dtype=bool,
230 default=False,
231 doc="Is this task to be run in a ``bootstrap`` fashion that does not require "
232 "a PTC or full calibrations?",
233 )
235 doCheckUnprocessableData = pexConfig.Field(
236 dtype=bool,
237 default=True,
238 doc="Check if this image is completely unprocessable due to all bad amps.",
239 )
241 overscanCamera = pexConfig.ConfigField(
242 dtype=OverscanCameraConfig,
243 doc="Per-detector and per-amplifier overscan configurations.",
244 )
246 serialOverscanMedianShiftSigmaThreshold = pexConfig.Field(
247 dtype=float,
248 default=numpy.inf,
249 doc="Number of sigma difference from per-amp overscan median (as compared to PTC) to "
250 "check if an amp is in a different state than the baseline PTC calib and should "
251 "be marked BAD. Set to np.inf/np.nan to turn off overscan median checking.",
252 )
254 ampNoiseThreshold = pexConfig.Field(
255 dtype=float,
256 default=25.0,
257 doc="Maximum amplifier noise (e-) that is allowed before an amp is masked as bad. "
258 "Set to np.inf/np.nan to turn off noise checking.",
259 )
261 bssVoltageMinimum = pexConfig.Field(
262 dtype=float,
263 default=5.0,
264 doc="Minimum back-side bias voltage. Below this the detector is ``off`` and an "
265 "UnprocessableDataError will be logged. Check will be skipped if doCheckUnprocessableData "
266 "is False or if value is less than or equal to 0.",
267 )
268 bssVoltageKeyword = pexConfig.Field(
269 dtype=str,
270 default="BSSVBS",
271 doc="Back-side bias voltage header keyword. Only checked if doCheckUnprocessableData is True "
272 "and bssVoltageMinimum is greater than 0.",
273 )
274 hvBiasKeyword = pexConfig.Field(
275 dtype=str,
276 default="HVBIAS",
277 doc="Back-side bias voltage on/off header keyword. Only checked if doCheckUnprocessableData is True "
278 "and bssVoltageMinimum is greater than 0.",
279 )
281 # Amplifier to CCD assembly configuration.
282 doAssembleCcd = pexConfig.Field(
283 dtype=bool,
284 default=True,
285 doc="Assemble amp-level exposures into a ccd-level exposure?"
286 )
287 assembleCcd = pexConfig.ConfigurableField(
288 target=AssembleCcdTask,
289 doc="CCD assembly task.",
290 )
292 # Bias subtraction.
293 doBias = pexConfig.Field(
294 dtype=bool,
295 doc="Apply bias frame correction?",
296 default=True,
297 )
299 # Deferred charge correction.
300 doDeferredCharge = pexConfig.Field(
301 dtype=bool,
302 doc="Apply deferred charge correction?",
303 default=True,
304 )
305 deferredChargeCorrection = pexConfig.ConfigurableField(
306 target=DeferredChargeTask,
307 doc="Deferred charge correction task.",
308 )
310 # Linearization.
311 doLinearize = pexConfig.Field(
312 dtype=bool,
313 doc="Correct for nonlinearity of the detector's response?",
314 default=True,
315 )
317 # Gains.
318 doCorrectGains = pexConfig.Field(
319 dtype=bool,
320 doc="Apply temperature correction to the gains?",
321 default=False,
322 )
323 doApplyGains = pexConfig.Field(
324 dtype=bool,
325 doc="Apply gains to the image?",
326 default=True,
327 )
329 # Variance construction.
330 doVariance = pexConfig.Field(
331 dtype=bool,
332 doc="Calculate variance?",
333 default=True
334 )
335 maskNegativeVariance = pexConfig.Field(
336 dtype=bool,
337 doc="Mask pixels that claim a negative variance. This likely indicates a failure "
338 "in the measurement of the overscan at an edge due to the data falling off faster "
339 "than the overscan model can account for it.",
340 default=True,
341 )
342 negativeVarianceMaskName = pexConfig.Field(
343 dtype=str,
344 doc="Mask plane to use to mark pixels with negative variance, if `maskNegativeVariance` is True.",
345 default="BAD",
346 )
347 doSaturation = pexConfig.Field(
348 dtype=bool,
349 doc="Mask saturated pixels? NB: this is totally independent of the"
350 " interpolation option - this is ONLY setting the bits in the mask."
351 " To have them interpolated make sure doInterpolate=True and"
352 " maskListToInterpolate includes SAT.",
353 default=True,
354 )
355 saturatedMaskName = pexConfig.Field(
356 dtype=str,
357 doc="Name of mask plane to use in saturation detection and interpolation.",
358 default="SAT",
359 )
360 defaultSaturationSource = pexConfig.ChoiceField(
361 dtype=str,
362 doc="Source to retrieve default amp-level saturation values.",
363 allowed={
364 "NONE": "No default saturation values; only config overrides will be used.",
365 "CAMERAMODEL": "Use the default from the camera model (old defaults).",
366 "PTCTURNOFF": "Use the ptcTurnoff value as the saturation level.",
367 },
368 default="PTCTURNOFF",
369 )
370 doSuspect = pexConfig.Field(
371 dtype=bool,
372 doc="Mask suspect pixels?",
373 default=True,
374 )
375 suspectMaskName = pexConfig.Field(
376 dtype=str,
377 doc="Name of mask plane to use for suspect pixels.",
378 default="SUSPECT",
379 )
380 defaultSuspectSource = pexConfig.ChoiceField(
381 dtype=str,
382 doc="Source to retrieve default amp-level suspect values.",
383 allowed={
384 "NONE": "No default suspect values; only config overrides will be used.",
385 "CAMERAMODEL": "Use the default from the camera model (old defaults).",
386 "PTCTURNOFF": "Use the ptcTurnoff value as the suspect level.",
387 },
388 default="PTCTURNOFF",
389 )
391 # Crosstalk.
392 doCrosstalk = pexConfig.Field(
393 dtype=bool,
394 doc="Apply intra-CCD crosstalk correction?",
395 default=True,
396 )
397 crosstalk = pexConfig.ConfigurableField(
398 target=CrosstalkTask,
399 doc="Intra-CCD crosstalk correction.",
400 )
402 # Masking options.
403 doITLDipMask = pexConfig.Field(
404 dtype=bool,
405 doc="Apply ``ITL dip`` masking. The ``itlDipMaskPlane`` mask plane "
406 "will be added even if this configuration is False.",
407 default=True,
408 )
409 itlDipMaskPlanes = pexConfig.ListField(
410 dtype=str,
411 doc="Mask plane to use for ITL dip pixels.",
412 default=["SUSPECT", "ITL_DIP"],
413 )
414 doDefect = pexConfig.Field(
415 dtype=bool,
416 doc="Apply correction for CCD defects, e.g. hot pixels?",
417 default=True,
418 )
419 badAmps = pexConfig.ListField(
420 dtype=str,
421 doc="List of bad amps that should be masked as BAD in the defect code. "
422 "Value should be of form {detector_name}_{amp_name}, e.g. ``R42_S21_C07``. "
423 "Only used if doDefect is True.",
424 default=[],
425 )
426 doNanMasking = pexConfig.Field(
427 dtype=bool,
428 doc="Mask non-finite (NAN, inf) pixels. The UNMASKEDNAN mask plane "
429 "will be added even if this configuration is False.",
430 default=True,
431 )
432 doWidenSaturationTrails = pexConfig.Field(
433 dtype=bool,
434 doc="Widen bleed trails based on their width.",
435 default=False,
436 )
437 masking = pexConfig.ConfigurableField(
438 target=MaskingTask,
439 doc="Masking task."
440 )
441 doE2VEdgeBleedMask = pexConfig.Field(
442 dtype=bool,
443 doc="Mask flag-like edge bleeds from saturated columns "
444 "in E2V amplifiers.",
445 default=True,
446 )
447 e2vEdgeBleedSatMinArea = pexConfig.Field(
448 dtype=int,
449 doc="Minimum limit of saturated cores footprint area to apply edge"
450 "bleed masking in E2V amplifiers.",
451 default=10000,
452 )
453 e2vEdgeBleedSatMaxArea = pexConfig.Field(
454 dtype=int,
455 doc="Maximum limit of saturated cores footprint area to apply edge"
456 "bleed masking in E2V amplifiers.",
457 default=100000,
458 )
459 e2vEdgeBleedYMax = pexConfig.Field(
460 dtype=int,
461 doc="Height in pixels of edge bleed masking in E2V amplifiers (width"
462 "is the width of the amplifier).",
463 default=350,
464 )
465 doITLEdgeBleedMask = pexConfig.Field(
466 dtype=bool,
467 doc="Mask edge bleeds from saturated columns in ITL amplifiers.",
468 default=True,
469 )
470 doITLSatSagMask = pexConfig.Field(
471 dtype=bool,
472 doc="Mask columns presenting saturation sag.",
473 default=True,
474 )
475 itlEdgeBleedSatMinArea = pexConfig.Field(
476 dtype=int,
477 doc="Minimum limit for saturated cores footprint area.",
478 default=10000,
479 )
480 itlEdgeBleedSatMaxArea = pexConfig.Field(
481 dtype=int,
482 doc="Maximum limit for saturated cores footprint area.",
483 default=100000,
484 )
485 itlEdgeBleedThreshold = pexConfig.Field(
486 dtype=float,
487 doc="Sky background threshold for edge bleed detection.",
488 default=5000.,
489 )
490 itlEdgeBleedModelConstant = pexConfig.Field(
491 dtype=float,
492 doc="Constant in the edge bleed exponential decay model.",
493 default=0.02,
494 )
496 # Interpolation options.
497 doInterpolate = pexConfig.Field(
498 dtype=bool,
499 doc="Interpolate masked pixels?",
500 default=True,
501 )
502 maskListToInterpolate = pexConfig.ListField(
503 dtype=str,
504 doc="List of mask planes that should be interpolated.",
505 default=["SAT", "BAD", "UNMASKEDNAN"],
506 )
507 doSaveInterpPixels = pexConfig.Field(
508 dtype=bool,
509 doc="Save a copy of the pre-interpolated pixel values?",
510 default=False,
511 )
512 useLegacyInterp = pexConfig.Field(
513 dtype=bool,
514 doc="Use the legacy interpolation algorithm. If False use Gaussian Process.",
515 default=True,
516 )
518 # Amp offset correction.
519 doAmpOffset = pexConfig.Field(
520 doc="Calculate amp offset corrections?",
521 dtype=bool,
522 default=False,
523 )
524 ampOffset = pexConfig.ConfigurableField(
525 doc="Amp offset correction task.",
526 target=AmpOffsetTask,
527 )
529 # Initial masking options.
530 doSetBadRegions = pexConfig.Field(
531 dtype=bool,
532 doc="Should we set the level of all BAD patches of the chip to the chip's average value?",
533 default=True,
534 )
536 # Brighter-Fatter correction.
537 doBrighterFatter = pexConfig.Field(
538 dtype=bool,
539 doc="Apply the brighter-fatter correction?",
540 default=True,
541 )
542 brighterFatterLevel = pexConfig.ChoiceField(
543 dtype=str,
544 doc="The level at which to correct for brighter-fatter.",
545 allowed={
546 "AMP": "Every amplifier treated separately.",
547 "DETECTOR": "One kernel per detector.",
548 },
549 default="DETECTOR",
550 )
551 brighterFatterMaxIter = pexConfig.Field(
552 dtype=int,
553 doc="Maximum number of iterations for the brighter-fatter correction.",
554 default=10,
555 )
556 brighterFatterThreshold = pexConfig.Field(
557 dtype=float,
558 doc="Threshold used to stop iterating the brighter-fatter correction. It is the "
559 "absolute value of the difference between the current corrected image and the one "
560 "from the previous iteration summed over all the pixels.",
561 default=1000,
562 )
563 brighterFatterMaskListToInterpolate = pexConfig.ListField(
564 dtype=str,
565 doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
566 "correction.",
567 default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
568 )
569 brighterFatterMaskGrowSize = pexConfig.Field(
570 dtype=int,
571 doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
572 "when brighter-fatter correction is applied.",
573 default=2,
574 )
575 brighterFatterFwhmForInterpolation = pexConfig.Field(
576 dtype=float,
577 doc="FWHM of PSF in arcseconds used for interpolation in brighter-fatter correction "
578 "(currently unused).",
579 default=1.0,
580 )
581 growSaturationFootprintSize = pexConfig.Field(
582 dtype=int,
583 doc="Number of pixels by which to grow the saturation footprints.",
584 default=1,
585 )
587 # Dark subtraction.
588 doDark = pexConfig.Field(
589 dtype=bool,
590 doc="Apply dark frame correction.",
591 default=True,
592 )
594 # Flat correction.
595 doFlat = pexConfig.Field(
596 dtype=bool,
597 doc="Apply flat field correction.",
598 default=True,
599 )
600 flatScalingType = pexConfig.ChoiceField(
601 dtype=str,
602 doc="The method for scaling the flat on the fly.",
603 default='USER',
604 allowed={
605 "USER": "Scale by flatUserScale",
606 "MEAN": "Scale by the inverse of the mean",
607 "MEDIAN": "Scale by the inverse of the median",
608 },
609 )
610 flatUserScale = pexConfig.Field(
611 dtype=float,
612 doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise.",
613 default=1.0,
614 )
616 # Calculate image quality statistics?
617 doStandardStatistics = pexConfig.Field(
618 dtype=bool,
619 doc="Should standard image quality statistics be calculated?",
620 default=True,
621 )
622 # Calculate additional statistics?
623 doCalculateStatistics = pexConfig.Field(
624 dtype=bool,
625 doc="Should additional ISR statistics be calculated?",
626 default=True,
627 )
628 isrStats = pexConfig.ConfigurableField(
629 target=IsrStatisticsTask,
630 doc="Task to calculate additional statistics.",
631 )
633 # Make binned images?
634 doBinnedExposures = pexConfig.Field(
635 dtype=bool,
636 doc="Should binned exposures be calculated?",
637 default=False,
638 )
639 binning = pexConfig.ConfigurableField(
640 target=BinImageDataTask,
641 doc="Task to bin the exposure.",
642 )
643 binFactor1 = pexConfig.Field(
644 dtype=int,
645 doc="Binning factor for first binned exposure. This is intended for a finely binned output.",
646 default=8,
647 check=lambda x: x > 1,
648 )
649 binFactor2 = pexConfig.Field(
650 dtype=int,
651 doc="Binning factor for second binned exposure. This is intended for a coarsely binned output.",
652 default=64,
653 check=lambda x: x > 1,
654 )
656 def validate(self):
657 super().validate()
659 if self.doBootstrap:
660 # Additional checks in bootstrap (no PTC/gains) mode.
661 if self.doApplyGains:
662 raise ValueError("Cannot run task with doBootstrap=True and doApplyGains=True.")
663 if self.doCorrectGains:
664 raise ValueError("Cannot run task with doBootstrap=True and doCorrectGains=True.")
665 if self.doCrosstalk and self.crosstalk.doQuadraticCrosstalkCorrection:
666 raise ValueError("Cannot apply quadratic crosstalk correction with doBootstrap=True.")
667 if numpy.isfinite(self.serialOverscanMedianShiftSigmaThreshold):
668 raise ValueError("Cannot do amp overscan level checks with doBootstrap=True.")
669 if numpy.isfinite(self.ampNoiseThreshold):
670 raise ValueError("Cannot do amp noise thresholds with doBootstrap=True.")
672 if self.doITLEdgeBleedMask and not self.doSaturation:
673 raise ValueError("Cannot do ITL edge bleed masking when doSaturation=False.")
674 if self.doE2VEdgeBleedMask and not self.doSaturation:
675 raise ValueError("Cannot do e2v edge bleed masking when doSaturation=False.")
677 def setDefaults(self):
678 super().setDefaults()
681class IsrTaskLSST(pipeBase.PipelineTask):
682 ConfigClass = IsrTaskLSSTConfig
683 _DefaultName = "isrLSST"
685 def __init__(self, **kwargs):
686 super().__init__(**kwargs)
687 self.makeSubtask("assembleCcd")
688 self.makeSubtask("deferredChargeCorrection")
689 self.makeSubtask("crosstalk")
690 self.makeSubtask("masking")
691 self.makeSubtask("isrStats")
692 self.makeSubtask("ampOffset")
693 self.makeSubtask("binning")
695 def runQuantum(self, butlerQC, inputRefs, outputRefs):
697 inputs = butlerQC.get(inputRefs)
698 self.validateInput(inputs)
700 if self.config.doHeaderProvenance:
701 # Add calibration provenanace info to header.
702 exposureMetadata = inputs['ccdExposure'].metadata
703 for inputName in sorted(list(inputs.keys())):
704 reference = getattr(inputRefs, inputName, None)
705 if reference is not None and hasattr(reference, "run"):
706 runKey = f"LSST CALIB RUN {inputName.upper()}"
707 runValue = reference.run
708 idKey = f"LSST CALIB UUID {inputName.upper()}"
709 idValue = str(reference.id)
710 dateKey = f"LSST CALIB DATE {inputName.upper()}"
711 dateValue = self.extractCalibDate(inputs[inputName])
712 if dateValue != "Unknown Unknown":
713 butlerQC.add_additional_provenance(reference, {"calib date": dateValue})
715 exposureMetadata[runKey] = runValue
716 exposureMetadata[idKey] = idValue
717 exposureMetadata[dateKey] = dateValue
719 outputs = self.run(**inputs)
720 butlerQC.put(outputs, outputRefs)
722 def validateInput(self, inputs):
723 """
724 This is a check that all the inputs required by the config
725 are available.
726 """
728 inputMap = {'dnlLUT': self.config.doDiffNonLinearCorrection,
729 'bias': self.config.doBias,
730 'deferredChargeCalib': self.config.doDeferredCharge,
731 'linearizer': self.config.doLinearize,
732 'ptc': self.config.doApplyGains,
733 'crosstalk': self.config.doCrosstalk,
734 'defects': self.config.doDefect,
735 'bfKernel': self.config.doBrighterFatter,
736 'dark': self.config.doDark,
737 }
739 for calibrationFile, configValue in inputMap.items():
740 if configValue and inputs[calibrationFile] is None:
741 raise RuntimeError("Must supply ", calibrationFile)
743 def diffNonLinearCorrection(self, ccdExposure, dnlLUT, **kwargs):
744 # TODO DM 36636
745 # isrFunctions.diffNonLinearCorrection
746 pass
748 def maskFullAmplifiers(self, ccdExposure, detector, defects, gains=None):
749 """
750 Check for fully masked bad amplifiers and mask them.
752 This includes defects which cover full amplifiers, as well
753 as amplifiers with nan gain values which should be used
754 if self.config.doApplyGains=True.
756 Full defect masking happens later to allow for defects which
757 cross amplifier boundaries.
759 Parameters
760 ----------
761 ccdExposure : `lsst.afw.image.Exposure`
762 Input exposure to be masked.
763 detector : `lsst.afw.cameraGeom.Detector`
764 Detector object.
765 defects : `lsst.ip.isr.Defects`
766 List of defects. Used to determine if an entire
767 amplifier is bad.
768 gains : `dict` [`str`, `float`], optional
769 Dictionary of gains to check if
770 self.config.doApplyGains=True.
772 Returns
773 -------
774 badAmpDict : `str`[`bool`]
775 Dictionary of amplifiers, keyed by name, value is True if
776 amplifier is fully masked.
777 """
778 badAmpDict = {}
780 maskedImage = ccdExposure.getMaskedImage()
782 for amp in detector:
783 ampName = amp.getName()
784 badAmpDict[ampName] = False
786 # Check if entire amp region is defined as a defect
787 # NB: need to use amp.getBBox() for correct comparison with current
788 # defects definition.
789 if defects is not None:
790 badAmpDict[ampName] = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
792 if badAmpDict[ampName]:
793 self.log.warning("Amplifier %s is bad (completely covered with defects)", ampName)
795 if gains is not None and self.config.doApplyGains:
796 if not math.isfinite(gains[ampName]):
797 badAmpDict[ampName] = True
799 self.log.warning("Amplifier %s is bad (non-finite gain)", ampName)
801 # In the case of a bad amp, we will set mask to "BAD"
802 # (here use amp.getRawBBox() for correct association with pixels in
803 # current ccdExposure).
804 if badAmpDict[ampName]:
805 dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
806 afwImage.PARENT)
807 maskView = dataView.getMask()
808 maskView |= maskView.getPlaneBitMask("BAD")
809 del maskView
811 return badAmpDict
813 def checkAllBadAmps(self, badAmpDict, detector):
814 """Check if all amps are marked as bad.
816 Parameters
817 ----------
818 badAmpDict : `str`[`bool`]
819 Dictionary of amplifiers, keyed by name, value is True if
820 amplifier is fully masked.
821 detector : `lsst.afw.cameraGeom.Detector`
822 Detector object.
824 Raises
825 ------
826 UnprocessableDataError if all amps are bad and doCheckUnprocessableData
827 configuration is True.
828 """
829 if not self.config.doCheckUnprocessableData:
830 return
832 for amp in detector:
833 if not badAmpDict.get(amp.getName(), False):
834 return
836 raise UnprocessableDataError(f"All amps in the exposure {detector.getName()} are bad; skipping ISR.")
838 def checkAmpOverscanLevel(self, badAmpDict, exposure, ptc):
839 """Check if the amplifier overscan levels have changed.
841 Any amplifier that has an overscan median level that has changed
842 significantly will be masked as BAD and added to toe badAmpDict.
844 Parameters
845 ----------
846 badAmpDict : `str` [`bool`]
847 Dictionary of amplifiers, keyed by name, value is True if
848 amplifier is fully masked.
849 exposure : `lsst.afw.image.Exposure`
850 Input exposure to be masked (untrimmed).
851 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
852 PTC dataset with gains/read noises.
854 Returns
855 -------
856 badAmpDict : `str`[`bool`]
857 Dictionary of amplifiers, keyed by name.
859 """
860 if isTrimmedExposure(exposure):
861 raise RuntimeError("checkAmpOverscanLevel must be run on an untrimmed exposure.")
863 # We want to consolidate all the amps into one warning if necessary.
864 # This config should not be set to a finite threshold if the necessary
865 # data is not in the PTC.
866 missingWarnString = "No PTC overscan information for amplifier "
867 missingWarnFlag = False
868 for amp in exposure.getDetector():
869 ampName = amp.getName()
871 if not numpy.isfinite(ptc.overscanMedian[ampName]) or \
872 not numpy.isfinite(ptc.overscanMedianSigma[ampName]):
873 missingWarnString += f"{ampName},"
874 missingWarnFlag = True
875 else:
876 key = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
877 # If it is missing, just return the PTC value and it
878 # will be skipped.
879 overscanLevel = exposure.metadata.get(key, ptc.overscanMedian[ampName])
880 pull = (overscanLevel - ptc.overscanMedian[ampName])/ptc.overscanMedianSigma[ampName]
881 if numpy.abs(pull) > self.config.serialOverscanMedianShiftSigmaThreshold:
882 self.log.warning(
883 "Amplifier %s has an overscan level that is %.2f sigma from the expected level; "
884 "masking it as BAD.",
885 ampName,
886 pull,
887 )
889 badAmpDict[ampName] = True
890 exposure.mask[amp.getRawBBox()] |= exposure.mask.getPlaneBitMask("BAD")
892 if missingWarnFlag:
893 self.log.warning(missingWarnString)
895 return badAmpDict
897 def checkAmpNoise(self, badAmpDict, exposure, ptc):
898 """Check if amplifier noise levels are above threshold.
900 Any amplifier that is above the noise level will be masked as BAD
901 and added to the badAmpDict.
903 Parameters
904 ----------
905 badAmpDict : `str` [`bool`]
906 Dictionary of amplifiers, keyed by name, value is True if
907 amplifier is fully masked.
908 exposure : `lsst.afw.image.Exposure`
909 Input exposure to be masked (untrimmed).
910 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
911 PTC dataset with gains/read noises.
913 Returns
914 -------
915 badAmpDict : `str`[`bool`]
916 Dictionary of amplifiers, keyed by name.
917 """
919 if isTrimmedExposure(exposure):
920 raise RuntimeError("checkAmpNoise must be run on an untrimmed exposure.")
922 for amp in exposure.getDetector():
923 ampName = amp.getName()
925 doMask = False
926 if ptc.noise[ampName] > self.config.ampNoiseThreshold:
927 self.log.info(
928 "Amplifier %s has a PTC noise level of %.2f e-, above threshold.",
929 ampName,
930 ptc.noise[ampName],
931 )
932 doMask = True
933 else:
934 key = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"
935 overscanNoise = exposure.metadata.get(key, numpy.nan)
936 if overscanNoise * ptc.gain[ampName] > self.config.ampNoiseThreshold:
937 self.log.warning(
938 "Amplifier %s has an overscan read noise level of %.2f e-, above threshold.",
939 ampName,
940 overscanNoise * ptc.gain[ampName],
941 )
942 doMask = True
944 if doMask:
945 badAmpDict[ampName] = True
947 exposure.mask[amp.getRawBBox()] |= exposure.mask.getPlaneBitMask("BAD")
949 return badAmpDict
951 def maskSaturatedPixels(self, badAmpDict, ccdExposure, detector, detectorConfig, ptc=None):
952 """
953 Mask SATURATED and SUSPECT pixels and check if any amplifiers
954 are fully masked.
956 Parameters
957 ----------
958 badAmpDict : `str` [`bool`]
959 Dictionary of amplifiers, keyed by name, value is True if
960 amplifier is fully masked.
961 ccdExposure : `lsst.afw.image.Exposure`
962 Input exposure to be masked.
963 detector : `lsst.afw.cameraGeom.Detector`
964 Detector object.
965 defects : `lsst.ip.isr.Defects`
966 List of defects. Used to determine if an entire
967 amplifier is bad.
968 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
969 Per-amplifier configurations.
970 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
971 PTC dataset (used if configured to use PTCTURNOFF).
973 Returns
974 -------
975 badAmpDict : `str`[`bool`]
976 Dictionary of amplifiers, keyed by name.
977 """
978 maskedImage = ccdExposure.getMaskedImage()
980 metadata = ccdExposure.metadata
982 if self.config.doSaturation and self.config.defaultSaturationSource == "PTCTURNOFF" and ptc is None:
983 raise RuntimeError("Must provide ptc if using PTCTURNOFF as saturation source.")
984 if self.config.doSuspect and self.config.defaultSuspectSource == "PTCTURNOFF" and ptc is None:
985 raise RuntimeError("Must provide ptc if using PTCTURNOFF as suspect source.")
987 for amp in detector:
988 ampName = amp.getName()
990 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
992 if badAmpDict[ampName]:
993 # No need to check fully bad amplifiers.
994 continue
996 # Mask saturated and suspect pixels.
997 limits = {}
998 if self.config.doSaturation:
999 if self.config.defaultSaturationSource == "PTCTURNOFF":
1000 limits.update({self.config.saturatedMaskName: ptc.ptcTurnoff[amp.getName()]})
1001 elif self.config.defaultSaturationSource == "CAMERAMODEL":
1002 # Set to the default from the camera model.
1003 limits.update({self.config.saturatedMaskName: amp.getSaturation()})
1004 elif self.config.defaultSaturationSource == "NONE":
1005 limits.update({self.config.saturatedMaskName: numpy.inf})
1007 # And update if it is set in the config.
1008 if math.isfinite(ampConfig.saturation):
1009 limits.update({self.config.saturatedMaskName: ampConfig.saturation})
1010 metadata[f"LSST ISR SATURATION LEVEL {ampName}"] = limits[self.config.saturatedMaskName]
1012 if self.config.doSuspect:
1013 if self.config.defaultSuspectSource == "PTCTURNOFF":
1014 limits.update({self.config.suspectMaskName: ptc.ptcTurnoff[amp.getName()]})
1015 elif self.config.defaultSuspectSource == "CAMERAMODEL":
1016 # Set to the default from the camera model.
1017 limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
1018 elif self.config.defaultSuspectSource == "NONE":
1019 limits.update({self.config.suspectMaskName: numpy.inf})
1021 # And update if it set in the config.
1022 if math.isfinite(ampConfig.suspectLevel):
1023 limits.update({self.config.suspectMaskName: ampConfig.suspectLevel})
1024 metadata[f"LSST ISR SUSPECT LEVEL {ampName}"] = limits[self.config.suspectMaskName]
1026 for maskName, maskThreshold in limits.items():
1027 if not math.isnan(maskThreshold):
1028 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
1029 toMask = (dataView.image.array >= maskThreshold)
1030 dataView.mask.array[toMask] |= dataView.mask.getPlaneBitMask(maskName)
1032 # Determine if we've fully masked this amplifier with SUSPECT and
1033 # SAT pixels.
1034 maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
1035 afwImage.PARENT)
1036 maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
1037 self.config.suspectMaskName])
1038 if numpy.all(maskView.getArray() & maskVal > 0):
1039 self.log.warning("Amplifier %s is bad (completely SATURATED or SUSPECT)", ampName)
1040 badAmpDict[ampName] = True
1041 maskView |= maskView.getPlaneBitMask("BAD")
1043 return badAmpDict
1045 def maskITLSatEdgesAndColumns(self, exposure, badAmpDict):
1047 # The following steps will rely on the footprint of saturated
1048 # cores with large areas.
1049 thresh = afwDetection.Threshold(exposure.mask.getPlaneBitMask("SAT"),
1050 afwDetection.Threshold.BITMASK
1051 )
1052 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
1054 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
1055 largeAreas, = numpy.where((satAreas >= self.config.itlEdgeBleedSatMinArea)
1056 & (satAreas < self.config.itlEdgeBleedSatMaxArea))
1058 for largeAreasIndex in largeAreas:
1060 fpCore = fpList[largeAreasIndex]
1062 # Edge bleed masking
1063 if self.config.doITLEdgeBleedMask:
1064 isrFunctions.maskITLEdgeBleed(ccdExposure=exposure,
1065 badAmpDict=badAmpDict,
1066 fpCore=fpCore,
1067 itlEdgeBleedSatMinArea=self.config.itlEdgeBleedSatMinArea,
1068 itlEdgeBleedSatMaxArea=self.config.itlEdgeBleedSatMaxArea,
1069 itlEdgeBleedThreshold=self.config.itlEdgeBleedThreshold,
1070 itlEdgeBleedModelConstant=self.config.itlEdgeBleedModelConstant,
1071 saturatedMaskName=self.config.saturatedMaskName,
1072 log=self.log
1073 )
1074 if self.config.doITLSatSagMask:
1075 isrFunctions.maskITLSatSag(ccdExposure=exposure, fpCore=fpCore,
1076 saturatedMaskName=self.config.saturatedMaskName)
1078 def overscanCorrection(self, mode, detectorConfig, detector, badAmpDict, ccdExposure):
1079 """Apply serial overscan correction in place to all amps.
1081 The actual overscan subtraction is performed by the
1082 `lsst.ip.isr.overscan.OverscanTask`, which is called here.
1084 Parameters
1085 ----------
1086 mode : `str`
1087 Must be `SERIAL` or `PARALLEL`.
1088 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1089 Per-amplifier configurations.
1090 detector : `lsst.afw.cameraGeom.Detector`
1091 Detector object.
1092 badAmpDict : `dict`
1093 Dictionary of amp name to whether it is a bad amp.
1094 ccdExposure : `lsst.afw.image.Exposure`
1095 Exposure to have overscan correction performed.
1097 Returns
1098 -------
1099 overscans : `list` [`lsst.pipe.base.Struct` or None]
1100 Overscan measurements (always in adu).
1101 Each result struct has components:
1103 ``imageFit``
1104 Value or fit subtracted from the amplifier image data.
1105 (scalar or `lsst.afw.image.Image`)
1106 ``overscanFit``
1107 Value or fit subtracted from the overscan image data.
1108 (scalar or `lsst.afw.image.Image`)
1109 ``overscanImage``
1110 Image of the overscan region with the overscan
1111 correction applied. This quantity is used to estimate
1112 the amplifier read noise empirically.
1113 (`lsst.afw.image.Image`)
1114 ``overscanMean``
1115 Mean overscan fit value. (`float`)
1116 ``overscanMedian``
1117 Median overscan fit value. (`float`)
1118 ``overscanSigma``
1119 Clipped standard deviation of the overscan fit. (`float`)
1120 ``residualMean``
1121 Mean of the overscan after fit subtraction. (`float`)
1122 ``residualMedian``
1123 Median of the overscan after fit subtraction. (`float`)
1124 ``residualSigma``
1125 Clipped standard deviation of the overscan after fit
1126 subtraction. (`float`)
1128 See Also
1129 --------
1130 lsst.ip.isr.overscan.OverscanTask
1131 """
1132 if mode not in ["SERIAL", "PARALLEL"]:
1133 raise ValueError("Mode must be SERIAL or PARALLEL")
1135 # This returns a list in amp order, with None for uncorrected amps.
1136 overscans = []
1138 for i, amp in enumerate(detector):
1139 ampName = amp.getName()
1141 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
1143 if mode == "SERIAL" and not ampConfig.doSerialOverscan:
1144 self.log.debug(
1145 "ISR_OSCAN: Amplifier %s/%s configured to skip serial overscan.",
1146 detector.getName(),
1147 ampName,
1148 )
1149 results = None
1150 elif mode == "PARALLEL" and not ampConfig.doParallelOverscan:
1151 self.log.debug(
1152 "ISR_OSCAN: Amplifier %s configured to skip parallel overscan.",
1153 detector.getName(),
1154 ampName,
1155 )
1156 results = None
1157 elif badAmpDict[ampName] or not ccdExposure.getBBox().contains(amp.getBBox()):
1158 results = None
1159 else:
1160 # This check is to confirm that we are not trying to run
1161 # overscan on an already trimmed image.
1162 if isTrimmedExposure(ccdExposure):
1163 self.log.warning(
1164 "ISR_OSCAN: No overscan region for amp %s. Not performing overscan correction.",
1165 ampName,
1166 )
1167 results = None
1168 else:
1169 if mode == "SERIAL":
1170 # We need to set up the subtask here with a custom
1171 # configuration.
1172 serialOverscan = SerialOverscanCorrectionTask(config=ampConfig.serialOverscanConfig)
1173 results = serialOverscan.run(ccdExposure, amp)
1174 else:
1175 config = ampConfig.parallelOverscanConfig
1176 parallelOverscan = ParallelOverscanCorrectionTask(
1177 config=config,
1178 )
1180 metadata = ccdExposure.metadata
1182 # We need to know the saturation level that was used
1183 # for the parallel overscan masking. If it isn't set
1184 # then the configured parallelOverscanSaturationLevel
1185 # will be used instead (assuming
1186 # doParallelOverscanSaturation is True). Note that
1187 # this will have the correct units (adu or electron)
1188 # depending on whether the gain has been applied.
1189 if self.config.doSaturation:
1190 saturationLevel = metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
1191 saturationLevel *= config.parallelOverscanSaturationLevelAdjustmentFactor
1192 else:
1193 saturationLevel = config.parallelOverscanSaturationLevel
1194 if ccdExposure.metadata["LSST ISR UNITS"] == "electron":
1195 # Need to convert to electron from adu.
1196 saturationLevel *= metadata[f"LSST ISR GAIN {amp.getName()}"]
1198 self.log.debug(
1199 "Using saturation level of %.2f for parallel overscan amp %s",
1200 saturationLevel,
1201 amp.getName(),
1202 )
1204 parallelOverscan.maskParallelOverscanAmp(
1205 ccdExposure,
1206 amp,
1207 saturationLevel=saturationLevel,
1208 )
1210 results = parallelOverscan.run(ccdExposure, amp)
1212 metadata = ccdExposure.metadata
1213 keyBase = "LSST ISR OVERSCAN"
1215 # The overscan is always in adu for the serial mode,
1216 # but, it may be electron in the parallel mode if
1217 # doApplyGains==True. If doApplyGains==True, then the
1218 # gains are applied to the untrimmed image, so the
1219 # overscan statistics units here will always match the
1220 # units of the image at this point.
1221 metadata[f"{keyBase} {mode} UNITS"] = ccdExposure.metadata["LSST ISR UNITS"]
1222 metadata[f"{keyBase} {mode} MEAN {ampName}"] = results.overscanMean
1223 metadata[f"{keyBase} {mode} MEDIAN {ampName}"] = results.overscanMedian
1224 metadata[f"{keyBase} {mode} STDEV {ampName}"] = results.overscanSigma
1226 metadata[f"{keyBase} RESIDUAL {mode} MEAN {ampName}"] = results.residualMean
1227 metadata[f"{keyBase} RESIDUAL {mode} MEDIAN {ampName}"] = results.residualMedian
1228 metadata[f"{keyBase} RESIDUAL {mode} STDEV {ampName}"] = results.residualSigma
1230 overscans.append(results)
1232 # Question: should this be finer grained?
1233 ccdExposure.metadata.set("OVERSCAN", "Overscan corrected")
1235 return overscans
1237 def correctGains(self, exposure, ptc, gains):
1238 # TODO DM 36639
1239 gains = []
1240 readNoise = []
1242 return gains, readNoise
1244 def maskNegativeVariance(self, exposure):
1245 """Identify and mask pixels with negative variance values.
1247 Parameters
1248 ----------
1249 exposure : `lsst.afw.image.Exposure`
1250 Exposure to process.
1252 See Also
1253 --------
1254 lsst.ip.isr.isrFunctions.updateVariance
1255 """
1256 maskPlane = exposure.getMask().getPlaneBitMask(self.config.negativeVarianceMaskName)
1257 bad = numpy.where(exposure.getVariance().getArray() <= 0.0)
1258 exposure.mask.array[bad] |= maskPlane
1260 def addVariancePlane(self, exposure, detector):
1261 """Add the variance plane to the image.
1263 The gain and read noise per amp must have been set in the
1264 exposure metadata as ``LSST ISR GAIN ampName`` and
1265 ``LSST ISR READNOISE ampName`` with the units of the image.
1266 Unit conversions for the variance plane will be done as
1267 necessary based on the exposure units.
1269 The units of the variance plane will always be of the same
1270 type as the units of the input image itself
1271 (``LSST ISR UNITS``^2).
1273 Parameters
1274 ----------
1275 exposure : `lsst.afw.image.Exposure`
1276 The exposure to add the variance plane.
1277 detector : `lsst.afw.cameraGeom.Detector`
1278 Detector with geometry info.
1279 """
1280 # NOTE: this will fail if the exposure is not trimmed.
1281 if not isTrimmedExposure(exposure):
1282 raise RuntimeError("Exposure must be trimmed to add variance plane.")
1284 isElectrons = (exposure.metadata["LSST ISR UNITS"] == "electron")
1286 for amp in detector:
1287 if exposure.getBBox().contains(amp.getBBox()):
1288 self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1289 ampExposure = exposure.Factory(exposure, amp.getBBox())
1291 # The effective gain is 1.0 if we are in electron units.
1292 # The metadata read noise is in the same units as the image.
1293 gain = exposure.metadata[f"LSST ISR GAIN {amp.getName()}"] if not isElectrons else 1.0
1294 readNoise = exposure.metadata[f"LSST ISR READNOISE {amp.getName()}"]
1296 isrFunctions.updateVariance(
1297 maskedImage=ampExposure.maskedImage,
1298 gain=gain,
1299 readNoise=readNoise,
1300 replace=False,
1301 )
1303 if self.config.qa is not None and self.config.qa.saveStats is True:
1304 qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1305 afwMath.MEDIAN | afwMath.STDEVCLIP)
1306 self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1307 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1308 qaStats.getValue(afwMath.STDEVCLIP))
1310 if self.config.maskNegativeVariance:
1311 self.maskNegativeVariance(exposure)
1313 def maskDefects(self, exposure, defectBaseList):
1314 """Mask defects using mask plane "BAD", in place.
1316 Parameters
1317 ----------
1318 exposure : `lsst.afw.image.Exposure`
1319 Exposure to process.
1321 defectBaseList : defect-type
1322 List of defects to mask. Can be of type `lsst.ip.isr.Defects`
1323 or `list` of `lsst.afw.image.DefectBase`.
1324 """
1325 maskedImage = exposure.getMaskedImage()
1326 if not isinstance(defectBaseList, Defects):
1327 # Promotes DefectBase to Defect
1328 defectList = Defects(defectBaseList)
1329 else:
1330 defectList = defectBaseList
1331 defectList.maskPixels(maskedImage, maskName="BAD")
1333 if len(self.config.badAmps) == 0:
1334 return
1336 detector = exposure.getDetector()
1337 mask = maskedImage.mask
1338 for badAmp in self.config.badAmps:
1339 if badAmp.startswith(detector.getName()):
1340 # Split on the full detector name plus _, which
1341 # gives us an empty string and the amp name.
1342 ampName = badAmp.split(detector.getName() + "_")[-1]
1343 self.log.info("Masking amplifier %s as bad via config.", ampName)
1344 mask[detector[ampName].getBBox()].array[:, :] |= mask.getPlaneBitMask("BAD")
1346 def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
1347 """Mask edge pixels with applicable mask plane.
1349 Parameters
1350 ----------
1351 exposure : `lsst.afw.image.Exposure`
1352 Exposure to process.
1353 numEdgePixels : `int`, optional
1354 Number of edge pixels to mask.
1355 maskPlane : `str`, optional
1356 Mask plane name to use.
1357 level : `str`, optional
1358 Level at which to mask edges.
1359 """
1360 maskedImage = exposure.getMaskedImage()
1361 maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
1363 if numEdgePixels > 0:
1364 if level == 'DETECTOR':
1365 boxes = [maskedImage.getBBox()]
1366 elif level == 'AMP':
1367 boxes = [amp.getBBox() for amp in exposure.getDetector()]
1369 for box in boxes:
1370 # This makes a bbox numEdgeSuspect pixels smaller than the
1371 # image on each side
1372 subImage = maskedImage[box]
1373 box.grow(-numEdgePixels)
1374 # Mask pixels outside box
1375 SourceDetectionTask.setEdgeBits(
1376 subImage,
1377 box,
1378 maskBitMask)
1380 def maskNan(self, exposure):
1381 """Mask NaNs using mask plane "UNMASKEDNAN", in place.
1383 Parameters
1384 ----------
1385 exposure : `lsst.afw.image.Exposure`
1386 Exposure to process.
1388 Notes
1389 -----
1390 We mask over all non-finite values (NaN, inf), including those
1391 that are masked with other bits (because those may or may not be
1392 interpolated over later, and we want to remove all NaN/infs).
1393 Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
1394 preserve the historical name.
1395 """
1396 maskedImage = exposure.getMaskedImage()
1398 # Find and mask NaNs
1399 maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
1400 numNans = maskNans(maskedImage, maskVal)
1401 self.metadata["NUMNANS"] = numNans
1402 if numNans > 0:
1403 self.log.warning("There were %d unmasked NaNs.", numNans)
1405 def setBadRegions(self, exposure):
1406 """Set bad regions from large contiguous regions.
1408 Parameters
1409 ----------
1410 exposure : `lsst.afw.Exposure`
1411 Exposure to set bad regions.
1413 Notes
1414 -----
1415 Reset and interpolate bad pixels.
1417 Large contiguous bad regions (which should have the BAD mask
1418 bit set) should have their values set to the image median.
1419 This group should include defects and bad amplifiers. As the
1420 area covered by these defects are large, there's little
1421 reason to expect that interpolation would provide a more
1422 useful value.
1424 Smaller defects can be safely interpolated after the larger
1425 regions have had their pixel values reset. This ensures
1426 that the remaining defects adjacent to bad amplifiers (as an
1427 example) do not attempt to interpolate extreme values.
1428 """
1429 badPixelCount, badPixelValue = isrFunctions.setBadRegions(exposure)
1430 if badPixelCount > 0:
1431 self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1433 @contextmanager
1434 def flatContext(self, exp, flat, dark=None):
1435 """Context manager that applies and removes flats and darks,
1436 if the task is configured to apply them.
1438 Parameters
1439 ----------
1440 exp : `lsst.afw.image.Exposure`
1441 Exposure to process.
1442 flat : `lsst.afw.image.Exposure`
1443 Flat exposure the same size as ``exp``.
1444 dark : `lsst.afw.image.Exposure`, optional
1445 Dark exposure the same size as ``exp``.
1447 Yields
1448 ------
1449 exp : `lsst.afw.image.Exposure`
1450 The flat and dark corrected exposure.
1451 """
1452 if self.config.doDark and dark is not None:
1453 self.darkCorrection(exp, dark)
1454 if self.config.doFlat and flat is not None:
1455 self.flatCorrection(exp, flat)
1456 try:
1457 yield exp
1458 finally:
1459 if self.config.doFlat and flat is not None:
1460 self.flatCorrection(exp, flat, invert=True)
1461 if self.config.doDark and dark is not None:
1462 self.darkCorrection(exp, dark, invert=True)
1464 def getBrighterFatterKernel(self, detector, bfKernel):
1465 detName = detector.getName()
1467 # This is expected to be a dictionary of amp-wise gains.
1468 bfGains = bfKernel.gain
1469 if bfKernel.level == 'DETECTOR':
1470 if detName in bfKernel.detKernels:
1471 bfKernelOut = bfKernel.detKernels[detName]
1472 return bfKernelOut, bfGains
1473 else:
1474 raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1475 elif bfKernel.level == 'AMP':
1476 self.log.info("Making DETECTOR level kernel from AMP based brighter "
1477 "fatter kernels.")
1478 bfKernel.makeDetectorKernelFromAmpwiseKernels(detName)
1479 bfKernelOut = bfKernel.detKernels[detName]
1480 return bfKernelOut, bfGains
1482 def applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, brighterFatterApplyGain,
1483 bfGains):
1484 """Apply a brighter fatter correction to the image using the
1485 method defined in Coulton et al. 2019.
1487 Note that this correction requires that the image is in units
1488 electrons.
1490 Parameters
1491 ----------
1492 ccdExposure : `lsst.afw.image.Exposure`
1493 Exposure to process.
1494 flat : `lsst.afw.image.Exposure`
1495 Flat exposure the same size as ``exp``.
1496 dark : `lsst.afw.image.Exposure`, optional
1497 Dark exposure the same size as ``exp``.
1498 bfKernel : `lsst.ip.isr.BrighterFatterKernel`
1499 The brighter-fatter kernel.
1500 brighterFatterApplyGain : `bool`
1501 Apply the gain to convert the image to electrons?
1502 bfGains : `dict`
1503 The gains to use if brighterFatterApplyGain = True.
1505 Yields
1506 ------
1507 exp : `lsst.afw.image.Exposure`
1508 The flat and dark corrected exposure.
1509 """
1510 interpExp = ccdExposure.clone()
1512 # We need to interpolate before we do B-F. Note that
1513 # brighterFatterFwhmForInterpolation is currently unused.
1514 isrFunctions.interpolateFromMask(
1515 maskedImage=interpExp.getMaskedImage(),
1516 fwhm=self.config.brighterFatterFwhmForInterpolation,
1517 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1518 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1519 useLegacyInterp=self.config.useLegacyInterp,
1520 )
1521 bfExp = interpExp.clone()
1522 bfResults = isrFunctions.brighterFatterCorrection(bfExp, bfKernel,
1523 self.config.brighterFatterMaxIter,
1524 self.config.brighterFatterThreshold,
1525 brighterFatterApplyGain,
1526 bfGains)
1527 bfCorrIters = bfResults[1]
1528 if bfCorrIters == self.config.brighterFatterMaxIter:
1529 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1530 bfResults[0])
1531 else:
1532 self.log.info("Finished brighter-fatter correction in %d iterations.",
1533 bfResults[1])
1535 image = ccdExposure.getMaskedImage().getImage()
1536 bfCorr = bfExp.getMaskedImage().getImage()
1537 bfCorr -= interpExp.getMaskedImage().getImage()
1538 image += bfCorr
1540 # Applying the brighter-fatter correction applies a
1541 # convolution to the science image. At the edges this
1542 # convolution may not have sufficient valid pixels to
1543 # produce a valid correction. Mark pixels within the size
1544 # of the brighter-fatter kernel as EDGE to warn of this
1545 # fact.
1546 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1547 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1548 maskPlane="EDGE")
1550 if self.config.brighterFatterMaskGrowSize > 0:
1551 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1552 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1553 isrFunctions.growMasks(ccdExposure.getMask(),
1554 radius=self.config.brighterFatterMaskGrowSize,
1555 maskNameList=maskPlane,
1556 maskValue=maskPlane)
1558 return ccdExposure, bfCorrIters
1560 def darkCorrection(self, exposure, darkExposure, invert=False):
1561 """Apply dark correction in place.
1563 Parameters
1564 ----------
1565 exposure : `lsst.afw.image.Exposure`
1566 Exposure to process.
1567 darkExposure : `lsst.afw.image.Exposure`
1568 Dark exposure of the same size as ``exposure``.
1569 invert : `Bool`, optional
1570 If True, re-add the dark to an already corrected image.
1572 Raises
1573 ------
1574 RuntimeError
1575 Raised if either ``exposure`` or ``darkExposure`` do not
1576 have their dark time defined.
1578 See Also
1579 --------
1580 lsst.ip.isr.isrFunctions.darkCorrection
1581 """
1582 expScale = exposure.visitInfo.darkTime
1583 if math.isnan(expScale):
1584 raise RuntimeError("Exposure darktime is NAN.")
1585 if darkExposure.visitInfo is not None \
1586 and not math.isnan(darkExposure.visitInfo.darkTime):
1587 darkScale = darkExposure.visitInfo.darkTime
1588 else:
1589 # DM-17444: darkExposure.visitInfo is None
1590 # so darkTime does not exist.
1591 self.log.warning("darkExposure.visitInfo does not exist. Using darkScale = 1.0.")
1592 darkScale = 1.0
1594 isrFunctions.darkCorrection(
1595 maskedImage=exposure.maskedImage,
1596 darkMaskedImage=darkExposure.maskedImage,
1597 expScale=expScale,
1598 darkScale=darkScale,
1599 invert=invert,
1600 )
1602 @staticmethod
1603 def extractCalibDate(calib):
1604 """Extract common calibration metadata values that will be written to
1605 output header.
1607 Parameters
1608 ----------
1609 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1610 Calibration to pull date information from.
1612 Returns
1613 -------
1614 dateString : `str`
1615 Calibration creation date string to add to header.
1616 """
1617 if hasattr(calib, "getMetadata"):
1618 if 'CALIB_CREATION_DATE' in calib.metadata:
1619 return " ".join((calib.metadata.get("CALIB_CREATION_DATE", "Unknown"),
1620 calib.metadata.get("CALIB_CREATION_TIME", "Unknown")))
1621 else:
1622 return " ".join((calib.metadata.get("CALIB_CREATE_DATE", "Unknown"),
1623 calib.metadata.get("CALIB_CREATE_TIME", "Unknown")))
1624 else:
1625 return "Unknown Unknown"
1627 def compareUnits(self, calibMetadata, calibName):
1628 """Compare units from calibration to ISR units.
1630 This compares calibration units (adu or electron) to whether
1631 doApplyGain is set.
1633 Parameters
1634 ----------
1635 calibMetadata : `lsst.daf.base.PropertyList`
1636 Calibration metadata from header.
1637 calibName : `str`
1638 Calibration name for log message.
1639 """
1640 calibUnits = calibMetadata.get("LSST ISR UNITS", "adu")
1641 isrUnits = "electron" if self.config.doApplyGains else "adu"
1642 if calibUnits != isrUnits:
1643 if self.config.doRaiseOnCalibMismatch:
1644 raise RuntimeError(
1645 "Unit mismatch: isr has %s units but %s has %s units",
1646 isrUnits,
1647 calibName,
1648 calibUnits,
1649 )
1650 else:
1651 self.log.warning(
1652 "Unit mismatch: isr has %s units but %s has %s units",
1653 isrUnits,
1654 calibName,
1655 calibUnits,
1656 )
1658 def convertIntToFloat(self, exposure):
1659 """Convert exposure image from uint16 to float.
1661 If the exposure does not need to be converted, the input is
1662 immediately returned. For exposures that are converted to use
1663 floating point pixels, the variance is set to unity and the
1664 mask to zero.
1666 Parameters
1667 ----------
1668 exposure : `lsst.afw.image.Exposure`
1669 The raw exposure to be converted.
1671 Returns
1672 -------
1673 newexposure : `lsst.afw.image.Exposure`
1674 The input ``exposure``, converted to floating point pixels.
1676 Raises
1677 ------
1678 RuntimeError
1679 Raised if the exposure type cannot be converted to float.
1681 """
1682 if isinstance(exposure, afwImage.ExposureF):
1683 # Nothing to be done
1684 self.log.debug("Exposure already of type float.")
1685 return exposure
1686 if not hasattr(exposure, "convertF"):
1687 raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
1689 newexposure = exposure.convertF()
1690 newexposure.variance[:] = 1
1691 newexposure.mask[:] = 0x0
1693 return newexposure
1695 def ditherCounts(self, exposure, detectorConfig, fallbackSeed=12345):
1696 """Dither the counts in the exposure.
1698 Parameters
1699 ----------
1700 exposure : `lsst.afw.image.Exposure`
1701 The raw exposure to be dithered.
1702 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1703 Configuration for overscan/etc for this detector.
1704 fallbackSeed : `int`, optional
1705 Random seed to fall back to if exposure.getInfo().getId() is
1706 not set.
1707 """
1708 if detectorConfig.integerDitherMode == "NONE":
1709 # Nothing to do here.
1710 return
1712 # This ID is a unique combination of {exposure, detector} for a raw
1713 # image as we have here. We additionally need to take the lower
1714 # 32 bits to be used as a random seed.
1715 if exposure.info.id is not None:
1716 seed = exposure.info.id & 0xFFFFFFFF
1717 else:
1718 seed = fallbackSeed
1719 self.log.warning("No exposure ID found; using fallback random seed.")
1721 self.log.info("Seeding dithering random number generator with %d.", seed)
1722 rng = numpy.random.RandomState(seed=seed)
1724 if detectorConfig.integerDitherMode == "POSITIVE":
1725 low = 0.0
1726 high = 1.0
1727 elif detectorConfig.integerDitherMode == "NEGATIVE":
1728 low = -1.0
1729 high = 0.0
1730 elif detectorConfig.integerDitherMode == "SYMMETRIC":
1731 low = -0.5
1732 high = 0.5
1733 else:
1734 raise RuntimeError("Invalid config")
1736 exposure.image.array[:, :] += rng.uniform(low=low, high=high, size=exposure.image.array.shape)
1738 def checkBssVoltage(self, exposure):
1739 """Check the back-side bias voltage to see if the detector is on.
1741 Parameters
1742 ----------
1743 exposure : `lsst.afw.image.ExposureF`
1744 Input exposure.
1746 Raises
1747 ------
1748 `UnprocessableDataError` if voltage is off.
1749 """
1750 voltage = exposure.metadata.get(self.config.bssVoltageKeyword, None)
1751 if voltage is None or not numpy.isfinite(voltage):
1752 self.log.warning(
1753 "Back-side bias voltage %s not found in metadata.",
1754 self.config.bssVoltageKeyword,
1755 )
1756 return
1758 hv = exposure.metadata.get(self.config.hvBiasKeyword, None)
1759 if hv is None:
1760 self.log.warning(
1761 "HV bias on %s not found in metadata.",
1762 self.config.hvBiasKeyword,
1763 )
1764 return
1766 if voltage < self.config.bssVoltageMinimum or hv == "OFF":
1767 detector = exposure.getDetector()
1768 raise UnprocessableDataError(
1769 f"Back-side bias voltage is turned off for {detector.getName()}; skipping ISR.",
1770 )
1772 @deprecated(
1773 reason=(
1774 "makeBinnedImages is no longer used. "
1775 "Please subtask lsst.ip.isr.BinImageDataTask instead."
1776 ),
1777 version="v28", category=FutureWarning
1778 )
1779 def makeBinnedImages(self, exposure):
1780 """Make visualizeVisit style binned exposures.
1782 Parameters
1783 ----------
1784 exposure : `lsst.afw.image.Exposure`
1785 Exposure to bin.
1787 Returns
1788 -------
1789 bin1 : `lsst.afw.image.Exposure`
1790 Binned exposure using binFactor1.
1791 bin2 : `lsst.afw.image.Exposure`
1792 Binned exposure using binFactor2.
1793 """
1794 mi = exposure.getMaskedImage()
1796 bin1 = afwMath.binImage(mi, self.config.binFactor1)
1797 bin2 = afwMath.binImage(mi, self.config.binFactor2)
1799 bin1 = afwImage.makeExposure(bin1)
1800 bin2 = afwImage.makeExposure(bin2)
1802 bin1.setInfo(exposure.getInfo())
1803 bin2.setInfo(exposure.getInfo())
1805 return bin1, bin2
1807 def run(self, ccdExposure, *, dnlLUT=None, bias=None, deferredChargeCalib=None, linearizer=None,
1808 ptc=None, crosstalk=None, defects=None, bfKernel=None, bfGains=None, dark=None,
1809 flat=None, camera=None, **kwargs
1810 ):
1812 detector = ccdExposure.getDetector()
1814 overscanDetectorConfig = self.config.overscanCamera.getOverscanDetectorConfig(detector)
1816 if self.config.doBootstrap and ptc is not None:
1817 self.log.warning("Task configured with doBootstrap=True. Ignoring provided PTC.")
1818 ptc = None
1820 if not self.config.doBootstrap:
1821 if ptc is None:
1822 raise RuntimeError("A PTC must be supplied if config.doBootstrap is False.")
1824 # Validation step: check inputs match exposure configuration.
1825 exposureMetadata = ccdExposure.metadata
1826 doRaise = self.config.doRaiseOnCalibMismatch
1827 keywords = self.config.cameraKeywordsToCompare
1828 if not self.config.doBootstrap:
1829 compareCameraKeywords(doRaise, keywords, exposureMetadata, ptc, "PTC", log=self.log)
1830 else:
1831 if self.config.doCorrectGains:
1832 raise RuntimeError("doCorrectGains is True but no ptc provided.")
1833 if self.config.doDiffNonLinearCorrection:
1834 if dnlLUT is None:
1835 raise RuntimeError("doDiffNonLinearCorrection is True but no dnlLUT provided.")
1836 compareCameraKeywords(doRaise, keywords, exposureMetadata, dnlLUT, "dnlLUT", log=self.log)
1837 if self.config.doLinearize:
1838 if linearizer is None:
1839 raise RuntimeError("doLinearize is True but no linearizer provided.")
1840 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer, "linearizer", log=self.log)
1841 if self.config.doBias:
1842 if bias is None:
1843 raise RuntimeError("doBias is True but no bias provided.")
1844 compareCameraKeywords(doRaise, keywords, exposureMetadata, bias, "bias", log=self.log)
1845 self.compareUnits(bias.metadata, "bias")
1846 if self.config.doCrosstalk:
1847 if crosstalk is None:
1848 raise RuntimeError("doCrosstalk is True but no crosstalk provided.")
1849 compareCameraKeywords(doRaise, keywords, exposureMetadata, crosstalk, "crosstalk", log=self.log)
1850 if self.config.doDeferredCharge:
1851 if deferredChargeCalib is None:
1852 raise RuntimeError("doDeferredCharge is True but no deferredChargeCalib provided.")
1853 compareCameraKeywords(
1854 doRaise,
1855 keywords,
1856 exposureMetadata,
1857 deferredChargeCalib,
1858 "CTI",
1859 log=self.log,
1860 )
1861 if self.config.doDefect:
1862 if defects is None:
1863 raise RuntimeError("doDefect is True but no defects provided.")
1864 compareCameraKeywords(doRaise, keywords, exposureMetadata, defects, "defects", log=self.log)
1865 if self.config.doDark:
1866 if dark is None:
1867 raise RuntimeError("doDark is True but no dark frame provided.")
1868 compareCameraKeywords(doRaise, keywords, exposureMetadata, dark, "dark", log=self.log)
1869 self.compareUnits(bias.metadata, "dark")
1870 if self.config.doBrighterFatter:
1871 if bfKernel is None:
1872 raise RuntimeError("doBrighterFatter is True not no bfKernel provided.")
1873 compareCameraKeywords(doRaise, keywords, exposureMetadata, bfKernel, "bf", log=self.log)
1874 if self.config.doFlat:
1875 if flat is None:
1876 raise RuntimeError("doFlat is True but no flat provided.")
1877 compareCameraKeywords(doRaise, keywords, exposureMetadata, flat, "flat", log=self.log)
1879 if self.config.doSaturation:
1880 if self.config.defaultSaturationSource in ["PTCTURNOFF",]:
1881 if ptc is None:
1882 raise RuntimeError(
1883 "doSaturation is True and defaultSaturationSource is "
1884 f"{self.config.defaultSaturationSource}, but no ptc provided."
1885 )
1886 if self.config.doSuspect:
1887 if self.config.defaultSuspectSource in ["PTCTURNOFF",]:
1888 if ptc is None:
1889 raise RuntimeError(
1890 "doSuspect is True and defaultSuspectSource is "
1891 f"{self.config.defaultSuspectSource}, but no ptc provided."
1892 )
1894 if self.config.doCheckUnprocessableData and self.config.bssVoltageMinimum > 0.0:
1895 self.checkBssVoltage(ccdExposure)
1897 # FIXME: Make sure that if linearity is done then it is matched
1898 # with the right PTC.
1900 # We keep track of units: start in adu.
1901 exposureMetadata["LSST ISR UNITS"] = "adu"
1902 exposureMetadata["LSST ISR CROSSTALK APPLIED"] = False
1903 exposureMetadata["LSST ISR OVERSCANLEVEL CHECKED"] = False
1904 exposureMetadata["LSST ISR NOISE CHECKED"] = False
1905 exposureMetadata["LSST ISR LINEARIZER APPLIED"] = False
1906 exposureMetadata["LSST ISR CTI APPLIED"] = False
1907 exposureMetadata["LSST ISR BIAS APPLIED"] = False
1908 exposureMetadata["LSST ISR DARK APPLIED"] = False
1909 exposureMetadata["LSST ISR BF APPLIED"] = False
1910 exposureMetadata["LSST ISR FLAT APPLIED"] = False
1911 exposureMetadata["LSST ISR DEFECTS APPLIED"] = False
1913 if self.config.doBootstrap:
1914 self.log.info("Configured using doBootstrap=True; using gain of 1.0 (adu units)")
1915 ptc = PhotonTransferCurveDataset([amp.getName() for amp in detector], "NOMINAL_PTC", 1)
1916 for amp in detector:
1917 ptc.gain[amp.getName()] = 1.0
1918 ptc.noise[amp.getName()] = 0.0
1920 exposureMetadata["LSST ISR BOOTSTRAP"] = self.config.doBootstrap
1922 # Set which gains to use
1923 gains = ptc.gain
1925 # And check if we have configured gains to override. This is
1926 # also a warning, since it should not be typical usage.
1927 for amp in detector:
1928 if not math.isnan(gain := overscanDetectorConfig.getOverscanAmpConfig(amp).gain):
1929 gains[amp.getName()] = gain
1930 self.log.warning(
1931 "Overriding gain for amp %s with configured value of %.3f.",
1932 amp.getName(),
1933 gain,
1934 )
1936 # First we convert the exposure to floating point values
1937 # (if necessary).
1938 self.log.debug("Converting exposure to floating point values.")
1939 ccdExposure = self.convertIntToFloat(ccdExposure)
1941 # Then we mark which amplifiers are completely bad from defects.
1942 badAmpDict = self.maskFullAmplifiers(ccdExposure, detector, defects, gains=gains)
1944 self.checkAllBadAmps(badAmpDict, detector)
1946 # Now we go through ISR steps.
1948 # Differential non-linearity correction.
1949 # Units: adu
1950 if self.config.doDiffNonLinearCorrection:
1951 self.diffNonLinearCorrection(ccdExposure, dnlLUT)
1953 # Dither the integer counts.
1954 # Input units: integerized adu
1955 # Output units: floating-point adu
1956 self.ditherCounts(ccdExposure, overscanDetectorConfig)
1958 # Serial overscan correction.
1959 # Input units: adu
1960 # Output units: adu
1961 if overscanDetectorConfig.doAnySerialOverscan:
1962 serialOverscans = self.overscanCorrection(
1963 "SERIAL",
1964 overscanDetectorConfig,
1965 detector,
1966 badAmpDict,
1967 ccdExposure,
1968 )
1970 if self.config.doBootstrap:
1971 # Get the empirical read noise
1972 for amp, serialOverscan in zip(detector, serialOverscans):
1973 if serialOverscan is None:
1974 ptc.noise[amp.getName()] = 0.0
1975 else:
1976 # All PhotonTransferCurveDataset objects should contain
1977 # noise attributes in units of electrons. The read
1978 # noise measured from overscans is always in adu, so we
1979 # scale it by the gain.
1980 # Note that in bootstrap mode, these gains will always
1981 # be 1.0, but we put this conversion here for clarity.
1982 ptc.noise[amp.getName()] = serialOverscan.residualSigma * gains[amp.getName()]
1983 else:
1984 serialOverscans = [None]*len(detector)
1986 # After serial overscan correction, we can mask SATURATED and
1987 # SUSPECT pixels. This updates badAmpDict if any amplifier
1988 # is fully saturated after serial overscan correction.
1990 # The saturation is currently assumed to be recorded in
1991 # overscan-corrected adu.
1992 badAmpDict = self.maskSaturatedPixels(
1993 badAmpDict,
1994 ccdExposure,
1995 detector,
1996 overscanDetectorConfig,
1997 ptc=ptc,
1998 )
2000 self.checkAllBadAmps(badAmpDict, detector)
2002 if self.config.doCorrectGains:
2003 # TODO: DM-36639
2004 # This requires the PTC (tbd) with the temperature dependence.
2005 self.log.info("Apply temperature dependence to the gains.")
2006 gains, readNoise = self.correctGains(ccdExposure, ptc, gains)
2008 # Do gain normalization.
2009 # Input units: adu
2010 # Output units: electron
2011 if self.config.doApplyGains:
2012 self.log.info("Using gain values to convert from adu to electron units.")
2013 isrFunctions.applyGains(ccdExposure, normalizeGains=False, ptcGains=gains, isTrimmed=False)
2014 # The units are now electron.
2015 exposureMetadata["LSST ISR UNITS"] = "electron"
2017 # Update the saturation units in the metadata if there.
2018 # These will always have the same units as the image.
2019 for amp in detector:
2020 ampName = amp.getName()
2021 if (key := f"LSST ISR SATURATION LEVEL {ampName}") in exposureMetadata:
2022 exposureMetadata[key] *= gains[ampName]
2023 if (key := f"LSST ISR SUSPECT LEVEL {ampName}") in exposureMetadata:
2024 exposureMetadata[key] *= gains[ampName]
2026 # Record gain and read noise in header.
2027 metadata = ccdExposure.metadata
2028 metadata["LSST ISR READNOISE UNITS"] = "electron"
2029 for amp in detector:
2030 # This includes any gain correction (if applied).
2031 metadata[f"LSST ISR GAIN {amp.getName()}"] = gains[amp.getName()]
2033 # At this stage, the read noise is always in electrons.
2034 noise = ptc.noise[amp.getName()]
2035 metadata[f"LSST ISR READNOISE {amp.getName()}"] = noise
2037 # Do crosstalk correction in the full region.
2038 # Output units: electron (adu if doBootstrap=True)
2039 if self.config.doCrosstalk:
2040 self.log.info("Applying crosstalk corrections to full amplifier region.")
2041 if self.config.doBootstrap and numpy.any(crosstalk.fitGains != 0):
2042 crosstalkGains = None
2043 else:
2044 crosstalkGains = gains
2045 self.crosstalk.run(
2046 ccdExposure,
2047 crosstalk=crosstalk,
2048 gains=crosstalkGains,
2049 fullAmplifier=True,
2050 badAmpDict=badAmpDict,
2051 ignoreVariance=True,
2052 )
2053 ccdExposure.metadata["LSST ISR CROSSTALK APPLIED"] = True
2055 # After crosstalk, we check for amplifier noise and state changes.
2056 if numpy.isfinite(self.config.serialOverscanMedianShiftSigmaThreshold):
2057 badAmpDict = self.checkAmpOverscanLevel(badAmpDict, ccdExposure, ptc)
2058 ccdExposure.metadata["LSST ISR OVERSCANLEVEL CHECKED"] = True
2060 if numpy.isfinite(self.config.ampNoiseThreshold):
2061 badAmpDict = self.checkAmpNoise(badAmpDict, ccdExposure, ptc)
2062 ccdExposure.metadata["LSST ISR NOISE CHECKED"] = True
2064 if numpy.isfinite(self.config.serialOverscanMedianShiftSigmaThreshold) or \
2065 numpy.isfinite(self.config.ampNoiseThreshold):
2066 self.checkAllBadAmps(badAmpDict, detector)
2068 # Parallel overscan correction.
2069 # Output units: electron (adu if doBootstrap=True)
2070 parallelOverscans = None
2071 if overscanDetectorConfig.doAnyParallelOverscan:
2072 # At the moment we do not use the return values from this task.
2073 parallelOverscans = self.overscanCorrection(
2074 "PARALLEL",
2075 overscanDetectorConfig,
2076 detector,
2077 badAmpDict,
2078 ccdExposure,
2079 )
2081 # Linearity correction
2082 # Output units: electron (adu if doBootstrap=True)
2083 if self.config.doLinearize:
2084 self.log.info("Applying linearizer.")
2085 # The linearizer is in units of adu.
2086 # If our units are electron, then pass in the gains
2087 # for conversion.
2088 if exposureMetadata["LSST ISR UNITS"] == "electron":
2089 linearityGains = gains
2090 else:
2091 linearityGains = None
2092 linearizer.applyLinearity(
2093 image=ccdExposure.image,
2094 detector=detector,
2095 log=self.log,
2096 gains=linearityGains,
2097 )
2098 ccdExposure.metadata["LSST ISR LINEARIZER APPLIED"] = True
2100 # Serial CTI (deferred charge) correction
2101 # This will be performed in electron units
2102 # Output units: same as input units
2103 if self.config.doDeferredCharge:
2104 self.deferredChargeCorrection.run(
2105 ccdExposure,
2106 deferredChargeCalib,
2107 gains=gains,
2108 )
2109 ccdExposure.metadata["LSST ISR CTI APPLIED"] = True
2111 # Save the untrimmed version for later statistics,
2112 # which still contains the overscan information
2113 untrimmedCcdExposure = None
2114 if self.config.isrStats.doCtiStatistics:
2115 untrimmedCcdExposure = ccdExposure.clone()
2117 # Assemble/trim
2118 # Output units: electron (adu if doBootstrap=True)
2119 if self.config.doAssembleCcd:
2120 self.log.info("Assembling CCD from amplifiers.")
2121 ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
2123 if self.config.expectWcs and not ccdExposure.getWcs():
2124 self.log.warning("No WCS found in input exposure.")
2126 # E2V edge bleed
2127 if self.config.doE2VEdgeBleedMask and detector.getPhysicalType() == "E2V":
2128 isrFunctions.maskE2VEdgeBleed(
2129 exposure=ccdExposure,
2130 e2vEdgeBleedSatMinArea=self.config.e2vEdgeBleedSatMinArea,
2131 e2vEdgeBleedSatMaxArea=self.config.e2vEdgeBleedSatMaxArea,
2132 e2vEdgeBleedYMax=self.config.e2vEdgeBleedYMax,
2133 saturatedMaskName=self.config.saturatedMaskName,
2134 log=self.log,
2135 )
2137 # ITL Dip Masking
2138 for maskPlane in self.config.itlDipMaskPlanes:
2139 if maskPlane not in ccdExposure.mask.getMaskPlaneDict():
2140 self.log.info("Adding %s mask plane to image.", maskPlane)
2141 ccdExposure.mask.addMaskPlane(maskPlane)
2143 if self.config.doITLDipMask:
2144 isrFunctions.maskITLDip(
2145 exposure=ccdExposure,
2146 detectorConfig=overscanDetectorConfig,
2147 log=self.log,
2148 maskPlaneNames=self.config.itlDipMaskPlanes,
2149 )
2151 if (self.config.doITLSatSagMask or self.config.doITLEdgeBleedMask) \
2152 and detector.getPhysicalType() == 'ITL':
2153 self.maskITLSatEdgesAndColumns(exposure=ccdExposure,
2154 badAmpDict=badAmpDict)
2156 # Bias subtraction
2157 # Output units: electron (adu if doBootstrap=True)
2158 if self.config.doBias:
2159 self.log.info("Applying bias correction.")
2160 # Bias frame and ISR unit consistency is checked at the top of
2161 # the run method.
2162 isrFunctions.biasCorrection(ccdExposure.maskedImage, bias.maskedImage)
2163 ccdExposure.metadata["LSST ISR BIAS APPLIED"] = True
2165 # Dark subtraction
2166 # Output units: electron (adu if doBootstrap=True)
2167 if self.config.doDark:
2168 self.log.info("Applying dark subtraction.")
2169 # Dark frame and ISR unit consistency is checked at the top of
2170 # the run method.
2171 self.darkCorrection(ccdExposure, dark)
2172 ccdExposure.metadata["LSST ISR DARK APPLIED"] = True
2174 # Defect masking
2175 # Masking block (defects, NAN pixels and trails).
2176 # Saturated and suspect pixels have already been masked.
2177 # Output units: electron (adu if doBootstrap=True)
2178 if self.config.doDefect:
2179 self.log.info("Applying defect masking.")
2180 self.maskDefects(ccdExposure, defects)
2181 ccdExposure.metadata["LSST ISR DEFECTS APPLIED"] = True
2183 self.log.info("Adding UNMASKEDNAN mask plane to image.")
2184 ccdExposure.mask.addMaskPlane("UNMASKEDNAN")
2185 if self.config.doNanMasking:
2186 self.log.info("Masking non-finite (NAN, inf) value pixels.")
2187 self.maskNan(ccdExposure)
2189 if self.config.doWidenSaturationTrails:
2190 self.log.info("Widening saturation trails.")
2191 isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
2193 # Brighter-Fatter
2194 # Output units: electron (adu if doBootstrap=True)
2195 if self.config.doBrighterFatter:
2196 self.log.info("Applying brighter-fatter correction.")
2198 bfKernelOut, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
2200 # Needs to be done in electrons; applyBrighterFatterCorrection
2201 # will convert the image if necessary.
2202 if exposureMetadata["LSST ISR UNITS"] == "electron":
2203 brighterFatterApplyGain = False
2204 else:
2205 brighterFatterApplyGain = True
2207 if brighterFatterApplyGain and (ptc is not None) and (bfGains != gains):
2208 # The supplied ptc should be the same as the ptc used to
2209 # generate the bfKernel, in which case they will have the
2210 # same stored amp-keyed dictionary of gains. If not, there
2211 # is a mismatch in the calibrations being used. This should
2212 # not be always be a fatal error, but ideally, everything
2213 # should to be consistent.
2214 self.log.warning("Need to apply gain for brighter-fatter, but the stored"
2215 "gains in the kernel are not the same as the gains stored"
2216 "in the PTC. Using the kernel gains.")
2218 ccdExposure, bfCorrIters = self.applyBrighterFatterCorrection(
2219 ccdExposure,
2220 flat,
2221 dark,
2222 bfKernelOut,
2223 brighterFatterApplyGain,
2224 bfGains,
2225 )
2227 ccdExposure.metadata["LSST ISR BF APPLIED"] = True
2228 metadata["LSST ISR BF ITERS"] = bfCorrIters
2230 # Variance plane creation
2231 # Output units: electron (adu if doBootstrap=True)
2232 if self.config.doVariance:
2233 self.addVariancePlane(ccdExposure, detector)
2235 # Flat-fielding
2236 # This may move elsewhere, but this is the most convenient
2237 # location for simple flat-fielding for attractive backgrounds.
2238 # Output units: electron (adu if doBootstrap=True)
2239 if self.config.doFlat:
2240 self.log.info("Applying flat correction.")
2241 isrFunctions.flatCorrection(
2242 maskedImage=ccdExposure.maskedImage,
2243 flatMaskedImage=flat.maskedImage,
2244 scalingType=self.config.flatScalingType,
2245 userScale=self.config.flatUserScale,
2246 )
2248 # Copy over valid polygon from flat if it is
2249 # available, and set NO_DATA to 0.0 (which may
2250 # be inherited from the flat in the fully
2251 # vignetted region).
2252 if (validPolygon := flat.info.getValidPolygon()) is not None:
2253 ccdExposure.info.setValidPolygon(validPolygon)
2255 noData = (ccdExposure.mask.array & ccdExposure.mask.getPlaneBitMask("NO_DATA")) > 0
2256 ccdExposure.image.array[noData] = 0.0
2257 ccdExposure.variance.array[noData] = 0.0
2259 ccdExposure.metadata["LSST ISR FLAT APPLIED"] = True
2260 ccdExposure.metadata["LSST ISR FLAT SOURCE"] = flat.metadata.get("FLATSRC", "UNKNOWN")
2262 # Pixel values for masked regions are set here
2263 preInterpExp = None
2264 if self.config.doSaveInterpPixels:
2265 preInterpExp = ccdExposure.clone()
2267 if self.config.doSetBadRegions:
2268 self.log.info('Setting values in large contiguous bad regions.')
2269 self.setBadRegions(ccdExposure)
2271 if self.config.doInterpolate:
2272 self.log.info("Interpolating masked pixels.")
2273 isrFunctions.interpolateFromMask(
2274 maskedImage=ccdExposure.getMaskedImage(),
2275 fwhm=self.config.brighterFatterFwhmForInterpolation,
2276 growSaturatedFootprints=self.config.growSaturationFootprintSize,
2277 maskNameList=list(self.config.maskListToInterpolate),
2278 useLegacyInterp=self.config.useLegacyInterp,
2279 )
2281 # Calculate amp offset corrections within the CCD.
2282 if self.config.doAmpOffset:
2283 if self.config.ampOffset.doApplyAmpOffset:
2284 self.log.info("Measuring and applying amp offset corrections.")
2285 else:
2286 self.log.info("Measuring amp offset corrections only, without applying them.")
2287 self.ampOffset.run(ccdExposure)
2289 # Calculate standard image quality statistics
2290 if self.config.doStandardStatistics:
2291 metadata = ccdExposure.metadata
2292 for amp in detector:
2293 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
2294 ampName = amp.getName()
2295 metadata[f"LSST ISR MASK SAT {ampName}"] = isrFunctions.countMaskedPixels(
2296 ampExposure.getMaskedImage(),
2297 [self.config.saturatedMaskName]
2298 )
2299 metadata[f"LSST ISR MASK BAD {ampName}"] = isrFunctions.countMaskedPixels(
2300 ampExposure.getMaskedImage(),
2301 ["BAD"]
2302 )
2303 metadata[f"LSST ISR MASK SUSPECT {ampName}"] = isrFunctions.countMaskedPixels(
2304 ampExposure.getMaskedImage(),
2305 ["SUSPECT"],
2306 )
2307 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
2308 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP)
2310 metadata[f"LSST ISR FINAL MEAN {ampName}"] = qaStats.getValue(afwMath.MEAN)
2311 metadata[f"LSST ISR FINAL MEDIAN {ampName}"] = qaStats.getValue(afwMath.MEDIAN)
2312 metadata[f"LSST ISR FINAL STDEV {ampName}"] = qaStats.getValue(afwMath.STDEVCLIP)
2314 k1 = f"LSST ISR FINAL MEDIAN {ampName}"
2315 k2 = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
2316 if overscanDetectorConfig.doAnySerialOverscan and k1 in metadata and k2 in metadata:
2317 metadata[f"LSST ISR LEVEL {ampName}"] = metadata[k1] - metadata[k2]
2318 else:
2319 metadata[f"LSST ISR LEVEL {ampName}"] = numpy.nan
2321 # calculate additional statistics.
2322 outputStatistics = None
2323 if self.config.doCalculateStatistics:
2324 outputStatistics = self.isrStats.run(ccdExposure,
2325 untrimmedInputExposure=untrimmedCcdExposure,
2326 serialOverscanResults=serialOverscans,
2327 parallelOverscanResults=parallelOverscans,
2328 bias=bias, dark=dark, flat=flat,
2329 ptc=ptc, defects=defects).results
2331 # do image binning.
2332 outputBin1Exposure = None
2333 outputBin2Exposure = None
2334 if self.config.doBinnedExposures:
2335 self.log.info("Creating binned exposures.")
2336 outputBin1Exposure = self.binning.run(
2337 ccdExposure,
2338 binFactor=self.config.binFactor1,
2339 ).outputData
2340 outputBin2Exposure = self.binning.run(
2341 ccdExposure,
2342 binFactor=self.config.binFactor2,
2343 ).outputData
2345 return pipeBase.Struct(
2346 exposure=ccdExposure,
2348 outputBin1Exposure=outputBin1Exposure,
2349 outputBin2Exposure=outputBin2Exposure,
2351 preInterpExposure=preInterpExp,
2352 outputExposure=ccdExposure,
2353 outputStatistics=outputStatistics,
2354 )