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