LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
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 import math
23 import numpy
24 
25 import lsst.geom
26 import lsst.afw.image as afwImage
27 import lsst.afw.math as afwMath
28 import lsst.pex.config as pexConfig
29 import lsst.pipe.base as pipeBase
31 
32 from contextlib import contextmanager
33 from lsstDebug import getDebugFrame
34 
35 from lsst.afw.cameraGeom import (PIXELS, FOCAL_PLANE, NullLinearityType,
36  ReadoutCorner)
37 from lsst.afw.display import getDisplay
38 from lsst.afw.geom import Polygon
39 from lsst.daf.persistence import ButlerDataRef
40 from lsst.daf.persistence.butler import NoResults
41 from lsst.meas.algorithms.detection import SourceDetectionTask
42 
43 from . import isrFunctions
44 from . import isrQa
45 from . import linearize
46 from .defects import Defects
47 
48 from .assembleCcdTask import AssembleCcdTask
49 from .crosstalk import CrosstalkTask, CrosstalkCalib
50 from .fringe import FringeTask
51 from .isr import maskNans
52 from .masking import MaskingTask
53 from .overscan import OverscanCorrectionTask
54 from .straylight import StrayLightTask
55 from .vignette import VignetteTask
56 from lsst.daf.butler import DimensionGraph
57 
58 
59 __all__ = ["IsrTask", "IsrTaskConfig", "RunIsrTask", "RunIsrConfig"]
60 
61 
62 def 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  This will be unused until DM-25348 resolves the quantum graph
70  generation issue.
71 
72  Parameters
73  ----------
74  datasetType : `str`
75  Dataset to lookup.
76  registry : `lsst.daf.butler.Registry`
77  Butler registry to query.
78  quantumDataId : `lsst.daf.butler.ExpandedDataCoordinate`
79  Data id to transform to identify crosstalkSources. The
80  ``detector`` entry will be stripped.
81  collections : `lsst.daf.butler.CollectionSearch`
82  Collections to search through.
83 
84  Returns
85  -------
86  results : `list` [`lsst.daf.butler.DatasetRef`]
87  List of datasets that match the query that will be used as
88  crosstalkSources.
89  """
90  newDataId = quantumDataId.subset(DimensionGraph(registry.dimensions, names=["instrument", "exposure"]))
91  results = list(registry.queryDatasets(datasetType,
92  collections=collections,
93  dataId=newDataId,
94  findFirst=True,
95  ).expanded())
96  return results
97 
98 
99 class 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  )
123  # TODO: DM-25348. This does not work yet to correctly load
124  # possible crosstalk sources.
125  crosstalkSources = cT.PrerequisiteInput(
126  name="isrOverscanCorrected",
127  doc="Overscan corrected input images.",
128  storageClass="Exposure",
129  dimensions=["instrument", "exposure", "detector"],
130  deferLoad=True,
131  multiple=True,
132  lookupFunction=crosstalkSourceLookup,
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  fringes = cT.PrerequisiteInput(
156  name="fringe",
157  doc="Input fringe calibration.",
158  storageClass="ExposureF",
159  dimensions=["instrument", "physical_filter", "detector"],
160  isCalibration=True,
161  )
162  strayLightData = cT.PrerequisiteInput(
163  name='yBackground',
164  doc="Input stray light calibration.",
165  storageClass="StrayLightData",
166  dimensions=["instrument", "physical_filter", "detector"],
167  isCalibration=True,
168  )
169  bfKernel = cT.PrerequisiteInput(
170  name='bfKernel',
171  doc="Input brighter-fatter kernel.",
172  storageClass="NumpyArray",
173  dimensions=["instrument"],
174  isCalibration=True,
175  )
176  newBFKernel = cT.PrerequisiteInput(
177  name='brighterFatterKernel',
178  doc="Newer complete kernel + gain solutions.",
179  storageClass="BrighterFatterKernel",
180  dimensions=["instrument", "detector"],
181  isCalibration=True,
182  )
183  defects = cT.PrerequisiteInput(
184  name='defects',
185  doc="Input defect tables.",
186  storageClass="Defects",
187  dimensions=["instrument", "detector"],
188  isCalibration=True,
189  )
190  linearizer = cT.PrerequisiteInput(
191  name='linearizer',
192  storageClass="Linearizer",
193  doc="Linearity correction calibration.",
194  dimensions=["instrument", "detector"],
195  isCalibration=True,
196  )
197  opticsTransmission = cT.PrerequisiteInput(
198  name="transmission_optics",
199  storageClass="TransmissionCurve",
200  doc="Transmission curve due to the optics.",
201  dimensions=["instrument"],
202  isCalibration=True,
203  )
204  filterTransmission = cT.PrerequisiteInput(
205  name="transmission_filter",
206  storageClass="TransmissionCurve",
207  doc="Transmission curve due to the filter.",
208  dimensions=["instrument", "physical_filter"],
209  isCalibration=True,
210  )
211  sensorTransmission = cT.PrerequisiteInput(
212  name="transmission_sensor",
213  storageClass="TransmissionCurve",
214  doc="Transmission curve due to the sensor.",
215  dimensions=["instrument", "detector"],
216  isCalibration=True,
217  )
218  atmosphereTransmission = cT.PrerequisiteInput(
219  name="transmission_atmosphere",
220  storageClass="TransmissionCurve",
221  doc="Transmission curve due to the atmosphere.",
222  dimensions=["instrument"],
223  isCalibration=True,
224  )
225  illumMaskedImage = cT.PrerequisiteInput(
226  name="illum",
227  doc="Input illumination correction.",
228  storageClass="MaskedImageF",
229  dimensions=["instrument", "physical_filter", "detector"],
230  isCalibration=True,
231  )
232 
233  outputExposure = cT.Output(
234  name='postISRCCD',
235  doc="Output ISR processed exposure.",
236  storageClass="Exposure",
237  dimensions=["instrument", "exposure", "detector"],
238  )
239  preInterpExposure = cT.Output(
240  name='preInterpISRCCD',
241  doc="Output ISR processed exposure, with pixels left uninterpolated.",
242  storageClass="ExposureF",
243  dimensions=["instrument", "exposure", "detector"],
244  )
245  outputOssThumbnail = cT.Output(
246  name="OssThumb",
247  doc="Output Overscan-subtracted thumbnail image.",
248  storageClass="Thumbnail",
249  dimensions=["instrument", "exposure", "detector"],
250  )
251  outputFlattenedThumbnail = cT.Output(
252  name="FlattenedThumb",
253  doc="Output flat-corrected thumbnail image.",
254  storageClass="Thumbnail",
255  dimensions=["instrument", "exposure", "detector"],
256  )
257 
258  def __init__(self, *, config=None):
259  super().__init__(config=config)
260 
261  if config.doBias is not True:
262  self.prerequisiteInputs.discard("bias")
263  if config.doLinearize is not True:
264  self.prerequisiteInputs.discard("linearizer")
265  if config.doCrosstalk is not True:
266  self.inputs.discard("crosstalkSources")
267  self.prerequisiteInputs.discard("crosstalk")
268  if config.doBrighterFatter is not True:
269  self.prerequisiteInputs.discard("bfKernel")
270  self.prerequisiteInputs.discard("newBFKernel")
271  if config.doDefect is not True:
272  self.prerequisiteInputs.discard("defects")
273  if config.doDark is not True:
274  self.prerequisiteInputs.discard("dark")
275  if config.doFlat is not True:
276  self.prerequisiteInputs.discard("flat")
277  if config.doAttachTransmissionCurve is not True:
278  self.prerequisiteInputs.discard("opticsTransmission")
279  self.prerequisiteInputs.discard("filterTransmission")
280  self.prerequisiteInputs.discard("sensorTransmission")
281  self.prerequisiteInputs.discard("atmosphereTransmission")
282  if config.doUseOpticsTransmission is not True:
283  self.prerequisiteInputs.discard("opticsTransmission")
284  if config.doUseFilterTransmission is not True:
285  self.prerequisiteInputs.discard("filterTransmission")
286  if config.doUseSensorTransmission is not True:
287  self.prerequisiteInputs.discard("sensorTransmission")
288  if config.doUseAtmosphereTransmission is not True:
289  self.prerequisiteInputs.discard("atmosphereTransmission")
290  if config.doIlluminationCorrection is not True:
291  self.prerequisiteInputs.discard("illumMaskedImage")
292 
293  if config.doWrite is not True:
294  self.outputs.discard("outputExposure")
295  self.outputs.discard("preInterpExposure")
296  self.outputs.discard("outputFlattenedThumbnail")
297  self.outputs.discard("outputOssThumbnail")
298  if config.doSaveInterpPixels is not True:
299  self.outputs.discard("preInterpExposure")
300  if config.qa.doThumbnailOss is not True:
301  self.outputs.discard("outputOssThumbnail")
302  if config.qa.doThumbnailFlattened is not True:
303  self.outputs.discard("outputFlattenedThumbnail")
304 
305 
306 class IsrTaskConfig(pipeBase.PipelineTaskConfig,
307  pipelineConnections=IsrTaskConnections):
308  """Configuration parameters for IsrTask.
309 
310  Items are grouped in the order in which they are executed by the task.
311  """
312  datasetType = pexConfig.Field(
313  dtype=str,
314  doc="Dataset type for input data; users will typically leave this alone, "
315  "but camera-specific ISR tasks will override it",
316  default="raw",
317  )
318 
319  fallbackFilterName = pexConfig.Field(
320  dtype=str,
321  doc="Fallback default filter name for calibrations.",
322  optional=True
323  )
324  useFallbackDate = pexConfig.Field(
325  dtype=bool,
326  doc="Pass observation date when using fallback filter.",
327  default=False,
328  )
329  expectWcs = pexConfig.Field(
330  dtype=bool,
331  default=True,
332  doc="Expect input science images to have a WCS (set False for e.g. spectrographs)."
333  )
334  fwhm = pexConfig.Field(
335  dtype=float,
336  doc="FWHM of PSF in arcseconds.",
337  default=1.0,
338  )
339  qa = pexConfig.ConfigField(
340  dtype=isrQa.IsrQaConfig,
341  doc="QA related configuration options.",
342  )
343 
344  # Image conversion configuration
345  doConvertIntToFloat = pexConfig.Field(
346  dtype=bool,
347  doc="Convert integer raw images to floating point values?",
348  default=True,
349  )
350 
351  # Saturated pixel handling.
352  doSaturation = pexConfig.Field(
353  dtype=bool,
354  doc="Mask saturated pixels? NB: this is totally independent of the"
355  " interpolation option - this is ONLY setting the bits in the mask."
356  " To have them interpolated make sure doSaturationInterpolation=True",
357  default=True,
358  )
359  saturatedMaskName = pexConfig.Field(
360  dtype=str,
361  doc="Name of mask plane to use in saturation detection and interpolation",
362  default="SAT",
363  )
364  saturation = pexConfig.Field(
365  dtype=float,
366  doc="The saturation level to use if no Detector is present in the Exposure (ignored if NaN)",
367  default=float("NaN"),
368  )
369  growSaturationFootprintSize = pexConfig.Field(
370  dtype=int,
371  doc="Number of pixels by which to grow the saturation footprints",
372  default=1,
373  )
374 
375  # Suspect pixel handling.
376  doSuspect = pexConfig.Field(
377  dtype=bool,
378  doc="Mask suspect pixels?",
379  default=False,
380  )
381  suspectMaskName = pexConfig.Field(
382  dtype=str,
383  doc="Name of mask plane to use for suspect pixels",
384  default="SUSPECT",
385  )
386  numEdgeSuspect = pexConfig.Field(
387  dtype=int,
388  doc="Number of edge pixels to be flagged as untrustworthy.",
389  default=0,
390  )
391  edgeMaskLevel = pexConfig.ChoiceField(
392  dtype=str,
393  doc="Mask edge pixels in which coordinate frame: DETECTOR or AMP?",
394  default="DETECTOR",
395  allowed={
396  'DETECTOR': 'Mask only the edges of the full detector.',
397  'AMP': 'Mask edges of each amplifier.',
398  },
399  )
400 
401  # Initial masking options.
402  doSetBadRegions = pexConfig.Field(
403  dtype=bool,
404  doc="Should we set the level of all BAD patches of the chip to the chip's average value?",
405  default=True,
406  )
407  badStatistic = pexConfig.ChoiceField(
408  dtype=str,
409  doc="How to estimate the average value for BAD regions.",
410  default='MEANCLIP',
411  allowed={
412  "MEANCLIP": "Correct using the (clipped) mean of good data",
413  "MEDIAN": "Correct using the median of the good data",
414  },
415  )
416 
417  # Overscan subtraction configuration.
418  doOverscan = pexConfig.Field(
419  dtype=bool,
420  doc="Do overscan subtraction?",
421  default=True,
422  )
423  overscan = pexConfig.ConfigurableField(
424  target=OverscanCorrectionTask,
425  doc="Overscan subtraction task for image segments.",
426  )
427 
428  overscanFitType = pexConfig.ChoiceField(
429  dtype=str,
430  doc="The method for fitting the overscan bias level.",
431  default='MEDIAN',
432  allowed={
433  "POLY": "Fit ordinary polynomial to the longest axis of the overscan region",
434  "CHEB": "Fit Chebyshev polynomial to the longest axis of the overscan region",
435  "LEG": "Fit Legendre polynomial to the longest axis of the overscan region",
436  "NATURAL_SPLINE": "Fit natural spline to the longest axis of the overscan region",
437  "CUBIC_SPLINE": "Fit cubic spline to the longest axis of the overscan region",
438  "AKIMA_SPLINE": "Fit Akima spline to the longest axis of the overscan region",
439  "MEAN": "Correct using the mean of the overscan region",
440  "MEANCLIP": "Correct using a clipped mean of the overscan region",
441  "MEDIAN": "Correct using the median of the overscan region",
442  "MEDIAN_PER_ROW": "Correct using the median per row of the overscan region",
443  },
444  deprecated=("Please configure overscan via the OverscanCorrectionConfig interface."
445  " This option will no longer be used, and will be removed after v20.")
446  )
447  overscanOrder = pexConfig.Field(
448  dtype=int,
449  doc=("Order of polynomial or to fit if overscan fit type is a polynomial, "
450  "or number of spline knots if overscan fit type is a spline."),
451  default=1,
452  deprecated=("Please configure overscan via the OverscanCorrectionConfig interface."
453  " This option will no longer be used, and will be removed after v20.")
454  )
455  overscanNumSigmaClip = pexConfig.Field(
456  dtype=float,
457  doc="Rejection threshold (sigma) for collapsing overscan before fit",
458  default=3.0,
459  deprecated=("Please configure overscan via the OverscanCorrectionConfig interface."
460  " This option will no longer be used, and will be removed after v20.")
461  )
462  overscanIsInt = pexConfig.Field(
463  dtype=bool,
464  doc="Treat overscan as an integer image for purposes of overscan.FitType=MEDIAN"
465  " and overscan.FitType=MEDIAN_PER_ROW.",
466  default=True,
467  deprecated=("Please configure overscan via the OverscanCorrectionConfig interface."
468  " This option will no longer be used, and will be removed after v20.")
469  )
470  # These options do not get deprecated, as they define how we slice up the image data.
471  overscanNumLeadingColumnsToSkip = pexConfig.Field(
472  dtype=int,
473  doc="Number of columns to skip in overscan, i.e. those closest to amplifier",
474  default=0,
475  )
476  overscanNumTrailingColumnsToSkip = pexConfig.Field(
477  dtype=int,
478  doc="Number of columns to skip in overscan, i.e. those farthest from amplifier",
479  default=0,
480  )
481  overscanMaxDev = pexConfig.Field(
482  dtype=float,
483  doc="Maximum deviation from the median for overscan",
484  default=1000.0, check=lambda x: x > 0
485  )
486  overscanBiasJump = pexConfig.Field(
487  dtype=bool,
488  doc="Fit the overscan in a piecewise-fashion to correct for bias jumps?",
489  default=False,
490  )
491  overscanBiasJumpKeyword = pexConfig.Field(
492  dtype=str,
493  doc="Header keyword containing information about devices.",
494  default="NO_SUCH_KEY",
495  )
496  overscanBiasJumpDevices = pexConfig.ListField(
497  dtype=str,
498  doc="List of devices that need piecewise overscan correction.",
499  default=(),
500  )
501  overscanBiasJumpLocation = pexConfig.Field(
502  dtype=int,
503  doc="Location of bias jump along y-axis.",
504  default=0,
505  )
506 
507  # Amplifier to CCD assembly configuration
508  doAssembleCcd = pexConfig.Field(
509  dtype=bool,
510  default=True,
511  doc="Assemble amp-level exposures into a ccd-level exposure?"
512  )
513  assembleCcd = pexConfig.ConfigurableField(
514  target=AssembleCcdTask,
515  doc="CCD assembly task",
516  )
517 
518  # General calibration configuration.
519  doAssembleIsrExposures = pexConfig.Field(
520  dtype=bool,
521  default=False,
522  doc="Assemble amp-level calibration exposures into ccd-level exposure?"
523  )
524  doTrimToMatchCalib = pexConfig.Field(
525  dtype=bool,
526  default=False,
527  doc="Trim raw data to match calibration bounding boxes?"
528  )
529 
530  # Bias subtraction.
531  doBias = pexConfig.Field(
532  dtype=bool,
533  doc="Apply bias frame correction?",
534  default=True,
535  )
536  biasDataProductName = pexConfig.Field(
537  dtype=str,
538  doc="Name of the bias data product",
539  default="bias",
540  )
541  doBiasBeforeOverscan = pexConfig.Field(
542  dtype=bool,
543  doc="Reverse order of overscan and bias correction.",
544  default=False
545  )
546 
547  # Variance construction
548  doVariance = pexConfig.Field(
549  dtype=bool,
550  doc="Calculate variance?",
551  default=True
552  )
553  gain = pexConfig.Field(
554  dtype=float,
555  doc="The gain to use if no Detector is present in the Exposure (ignored if NaN)",
556  default=float("NaN"),
557  )
558  readNoise = pexConfig.Field(
559  dtype=float,
560  doc="The read noise to use if no Detector is present in the Exposure",
561  default=0.0,
562  )
563  doEmpiricalReadNoise = pexConfig.Field(
564  dtype=bool,
565  default=False,
566  doc="Calculate empirical read noise instead of value from AmpInfo data?"
567  )
568 
569  # Linearization.
570  doLinearize = pexConfig.Field(
571  dtype=bool,
572  doc="Correct for nonlinearity of the detector's response?",
573  default=True,
574  )
575 
576  # Crosstalk.
577  doCrosstalk = pexConfig.Field(
578  dtype=bool,
579  doc="Apply intra-CCD crosstalk correction?",
580  default=False,
581  )
582  doCrosstalkBeforeAssemble = pexConfig.Field(
583  dtype=bool,
584  doc="Apply crosstalk correction before CCD assembly, and before trimming?",
585  default=False,
586  )
587  crosstalk = pexConfig.ConfigurableField(
588  target=CrosstalkTask,
589  doc="Intra-CCD crosstalk correction",
590  )
591 
592  # Masking options.
593  doDefect = pexConfig.Field(
594  dtype=bool,
595  doc="Apply correction for CCD defects, e.g. hot pixels?",
596  default=True,
597  )
598  doNanMasking = pexConfig.Field(
599  dtype=bool,
600  doc="Mask non-finite (NAN, inf) pixels?",
601  default=True,
602  )
603  doWidenSaturationTrails = pexConfig.Field(
604  dtype=bool,
605  doc="Widen bleed trails based on their width?",
606  default=True
607  )
608 
609  # Brighter-Fatter correction.
610  doBrighterFatter = pexConfig.Field(
611  dtype=bool,
612  default=False,
613  doc="Apply the brighter-fatter correction?"
614  )
615  brighterFatterLevel = pexConfig.ChoiceField(
616  dtype=str,
617  default="DETECTOR",
618  doc="The level at which to correct for brighter-fatter.",
619  allowed={
620  "AMP": "Every amplifier treated separately.",
621  "DETECTOR": "One kernel per detector",
622  }
623  )
624  brighterFatterMaxIter = pexConfig.Field(
625  dtype=int,
626  default=10,
627  doc="Maximum number of iterations for the brighter-fatter correction"
628  )
629  brighterFatterThreshold = pexConfig.Field(
630  dtype=float,
631  default=1000,
632  doc="Threshold used to stop iterating the brighter-fatter correction. It is the "
633  "absolute value of the difference between the current corrected image and the one "
634  "from the previous iteration summed over all the pixels."
635  )
636  brighterFatterApplyGain = pexConfig.Field(
637  dtype=bool,
638  default=True,
639  doc="Should the gain be applied when applying the brighter-fatter correction?"
640  )
641  brighterFatterMaskListToInterpolate = pexConfig.ListField(
642  dtype=str,
643  doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
644  "correction.",
645  default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
646  )
647  brighterFatterMaskGrowSize = pexConfig.Field(
648  dtype=int,
649  default=0,
650  doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
651  "when brighter-fatter correction is applied."
652  )
653 
654  # Dark subtraction.
655  doDark = pexConfig.Field(
656  dtype=bool,
657  doc="Apply dark frame correction?",
658  default=True,
659  )
660  darkDataProductName = pexConfig.Field(
661  dtype=str,
662  doc="Name of the dark data product",
663  default="dark",
664  )
665 
666  # Camera-specific stray light removal.
667  doStrayLight = pexConfig.Field(
668  dtype=bool,
669  doc="Subtract stray light in the y-band (due to encoder LEDs)?",
670  default=False,
671  )
672  strayLight = pexConfig.ConfigurableField(
673  target=StrayLightTask,
674  doc="y-band stray light correction"
675  )
676 
677  # Flat correction.
678  doFlat = pexConfig.Field(
679  dtype=bool,
680  doc="Apply flat field correction?",
681  default=True,
682  )
683  flatDataProductName = pexConfig.Field(
684  dtype=str,
685  doc="Name of the flat data product",
686  default="flat",
687  )
688  flatScalingType = pexConfig.ChoiceField(
689  dtype=str,
690  doc="The method for scaling the flat on the fly.",
691  default='USER',
692  allowed={
693  "USER": "Scale by flatUserScale",
694  "MEAN": "Scale by the inverse of the mean",
695  "MEDIAN": "Scale by the inverse of the median",
696  },
697  )
698  flatUserScale = pexConfig.Field(
699  dtype=float,
700  doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise",
701  default=1.0,
702  )
703  doTweakFlat = pexConfig.Field(
704  dtype=bool,
705  doc="Tweak flats to match observed amplifier ratios?",
706  default=False
707  )
708 
709  # Amplifier normalization based on gains instead of using flats configuration.
710  doApplyGains = pexConfig.Field(
711  dtype=bool,
712  doc="Correct the amplifiers for their gains instead of applying flat correction",
713  default=False,
714  )
715  normalizeGains = pexConfig.Field(
716  dtype=bool,
717  doc="Normalize all the amplifiers in each CCD to have the same median value.",
718  default=False,
719  )
720 
721  # Fringe correction.
722  doFringe = pexConfig.Field(
723  dtype=bool,
724  doc="Apply fringe correction?",
725  default=True,
726  )
727  fringe = pexConfig.ConfigurableField(
728  target=FringeTask,
729  doc="Fringe subtraction task",
730  )
731  fringeAfterFlat = pexConfig.Field(
732  dtype=bool,
733  doc="Do fringe subtraction after flat-fielding?",
734  default=True,
735  )
736 
737  # Initial CCD-level background statistics options.
738  doMeasureBackground = pexConfig.Field(
739  dtype=bool,
740  doc="Measure the background level on the reduced image?",
741  default=False,
742  )
743 
744  # Camera-specific masking configuration.
745  doCameraSpecificMasking = pexConfig.Field(
746  dtype=bool,
747  doc="Mask camera-specific bad regions?",
748  default=False,
749  )
750  masking = pexConfig.ConfigurableField(
751  target=MaskingTask,
752  doc="Masking task."
753  )
754 
755  # Interpolation options.
756 
757  doInterpolate = pexConfig.Field(
758  dtype=bool,
759  doc="Interpolate masked pixels?",
760  default=True,
761  )
762  doSaturationInterpolation = pexConfig.Field(
763  dtype=bool,
764  doc="Perform interpolation over pixels masked as saturated?"
765  " NB: This is independent of doSaturation; if that is False this plane"
766  " will likely be blank, resulting in a no-op here.",
767  default=True,
768  )
769  doNanInterpolation = pexConfig.Field(
770  dtype=bool,
771  doc="Perform interpolation over pixels masked as NaN?"
772  " NB: This is independent of doNanMasking; if that is False this plane"
773  " will likely be blank, resulting in a no-op here.",
774  default=True,
775  )
776  doNanInterpAfterFlat = pexConfig.Field(
777  dtype=bool,
778  doc=("If True, ensure we interpolate NaNs after flat-fielding, even if we "
779  "also have to interpolate them before flat-fielding."),
780  default=False,
781  )
782  maskListToInterpolate = pexConfig.ListField(
783  dtype=str,
784  doc="List of mask planes that should be interpolated.",
785  default=['SAT', 'BAD'],
786  )
787  doSaveInterpPixels = pexConfig.Field(
788  dtype=bool,
789  doc="Save a copy of the pre-interpolated pixel values?",
790  default=False,
791  )
792 
793  # Default photometric calibration options.
794  fluxMag0T1 = pexConfig.DictField(
795  keytype=str,
796  itemtype=float,
797  doc="The approximate flux of a zero-magnitude object in a one-second exposure, per filter.",
798  default=dict((f, pow(10.0, 0.4*m)) for f, m in (("Unknown", 28.0),
799  ))
800  )
801  defaultFluxMag0T1 = pexConfig.Field(
802  dtype=float,
803  doc="Default value for fluxMag0T1 (for an unrecognized filter).",
804  default=pow(10.0, 0.4*28.0)
805  )
806 
807  # Vignette correction configuration.
808  doVignette = pexConfig.Field(
809  dtype=bool,
810  doc="Apply vignetting parameters?",
811  default=False,
812  )
813  vignette = pexConfig.ConfigurableField(
814  target=VignetteTask,
815  doc="Vignetting task.",
816  )
817 
818  # Transmission curve configuration.
819  doAttachTransmissionCurve = pexConfig.Field(
820  dtype=bool,
821  default=False,
822  doc="Construct and attach a wavelength-dependent throughput curve for this CCD image?"
823  )
824  doUseOpticsTransmission = pexConfig.Field(
825  dtype=bool,
826  default=True,
827  doc="Load and use transmission_optics (if doAttachTransmissionCurve is True)?"
828  )
829  doUseFilterTransmission = pexConfig.Field(
830  dtype=bool,
831  default=True,
832  doc="Load and use transmission_filter (if doAttachTransmissionCurve is True)?"
833  )
834  doUseSensorTransmission = pexConfig.Field(
835  dtype=bool,
836  default=True,
837  doc="Load and use transmission_sensor (if doAttachTransmissionCurve is True)?"
838  )
839  doUseAtmosphereTransmission = pexConfig.Field(
840  dtype=bool,
841  default=True,
842  doc="Load and use transmission_atmosphere (if doAttachTransmissionCurve is True)?"
843  )
844 
845  # Illumination correction.
846  doIlluminationCorrection = pexConfig.Field(
847  dtype=bool,
848  default=False,
849  doc="Perform illumination correction?"
850  )
851  illuminationCorrectionDataProductName = pexConfig.Field(
852  dtype=str,
853  doc="Name of the illumination correction data product.",
854  default="illumcor",
855  )
856  illumScale = pexConfig.Field(
857  dtype=float,
858  doc="Scale factor for the illumination correction.",
859  default=1.0,
860  )
861  illumFilters = pexConfig.ListField(
862  dtype=str,
863  default=[],
864  doc="Only perform illumination correction for these filters."
865  )
866 
867  # Write the outputs to disk. If ISR is run as a subtask, this may not be needed.
868  doWrite = pexConfig.Field(
869  dtype=bool,
870  doc="Persist postISRCCD?",
871  default=True,
872  )
873 
874  def validate(self):
875  super().validate()
876  if self.doFlatdoFlat and self.doApplyGainsdoApplyGains:
877  raise ValueError("You may not specify both doFlat and doApplyGains")
878  if self.doBiasBeforeOverscandoBiasBeforeOverscan and self.doTrimToMatchCalibdoTrimToMatchCalib:
879  raise ValueError("You may not specify both doBiasBeforeOverscan and doTrimToMatchCalib")
880  if self.doSaturationInterpolationdoSaturationInterpolation and self.saturatedMaskNamesaturatedMaskName not in self.maskListToInterpolatemaskListToInterpolate:
881  self.maskListToInterpolatemaskListToInterpolate.append(self.saturatedMaskNamesaturatedMaskName)
882  if not self.doSaturationInterpolationdoSaturationInterpolation and self.saturatedMaskNamesaturatedMaskName in self.maskListToInterpolatemaskListToInterpolate:
883  self.maskListToInterpolatemaskListToInterpolate.remove(self.saturatedMaskNamesaturatedMaskName)
884  if self.doNanInterpolationdoNanInterpolation and "UNMASKEDNAN" not in self.maskListToInterpolatemaskListToInterpolate:
885  self.maskListToInterpolatemaskListToInterpolate.append("UNMASKEDNAN")
886 
887 
888 class IsrTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
889  """Apply common instrument signature correction algorithms to a raw frame.
890 
891  The process for correcting imaging data is very similar from
892  camera to camera. This task provides a vanilla implementation of
893  doing these corrections, including the ability to turn certain
894  corrections off if they are not needed. The inputs to the primary
895  method, `run()`, are a raw exposure to be corrected and the
896  calibration data products. The raw input is a single chip sized
897  mosaic of all amps including overscans and other non-science
898  pixels. The method `runDataRef()` identifies and defines the
899  calibration data products, and is intended for use by a
900  `lsst.pipe.base.cmdLineTask.CmdLineTask` and takes as input only a
901  `daf.persistence.butlerSubset.ButlerDataRef`. This task may be
902  subclassed for different camera, although the most camera specific
903  methods have been split into subtasks that can be redirected
904  appropriately.
905 
906  The __init__ method sets up the subtasks for ISR processing, using
907  the defaults from `lsst.ip.isr`.
908 
909  Parameters
910  ----------
911  args : `list`
912  Positional arguments passed to the Task constructor. None used at this time.
913  kwargs : `dict`, optional
914  Keyword arguments passed on to the Task constructor. None used at this time.
915  """
916  ConfigClass = IsrTaskConfig
917  _DefaultName = "isr"
918 
919  def __init__(self, **kwargs):
920  super().__init__(**kwargs)
921  self.makeSubtask("assembleCcd")
922  self.makeSubtask("crosstalk")
923  self.makeSubtask("strayLight")
924  self.makeSubtask("fringe")
925  self.makeSubtask("masking")
926  self.makeSubtask("overscan")
927  self.makeSubtask("vignette")
928 
929  def runQuantum(self, butlerQC, inputRefs, outputRefs):
930  inputs = butlerQC.get(inputRefs)
931 
932  try:
933  inputs['detectorNum'] = inputRefs.ccdExposure.dataId['detector']
934  except Exception as e:
935  raise ValueError("Failure to find valid detectorNum value for Dataset %s: %s." %
936  (inputRefs, e))
937 
938  inputs['isGen3'] = True
939 
940  detector = inputs['ccdExposure'].getDetector()
941 
942  if self.config.doCrosstalk is True:
943  # Crosstalk sources need to be defined by the pipeline
944  # yaml if they exist.
945  if 'crosstalk' in inputs and inputs['crosstalk'] is not None:
946  if not isinstance(inputs['crosstalk'], CrosstalkCalib):
947  inputs['crosstalk'] = CrosstalkCalib.fromTable(inputs['crosstalk'])
948  else:
949  coeffVector = (self.config.crosstalk.crosstalkValues
950  if self.config.crosstalk.useConfigCoefficients else None)
951  crosstalkCalib = CrosstalkCalib().fromDetector(detector, coeffVector=coeffVector)
952  inputs['crosstalk'] = crosstalkCalib
953  if inputs['crosstalk'].interChip and len(inputs['crosstalk'].interChip) > 0:
954  if 'crosstalkSources' not in inputs:
955  self.log.warn("No crosstalkSources found for chip with interChip terms!")
956 
957  if self.doLinearizedoLinearize(detector) is True:
958  if 'linearizer' in inputs:
959  if isinstance(inputs['linearizer'], dict):
960  linearizer = linearize.Linearizer(detector=detector, log=self.log)
961  linearizer.fromYaml(inputs['linearizer'])
962  self.log.warn("Dictionary linearizers will be deprecated in DM-28741.")
963  elif isinstance(inputs['linearizer'], numpy.ndarray):
964  linearizer = linearize.Linearizer(table=inputs.get('linearizer', None),
965  detector=detector,
966  log=self.log)
967  self.log.warn("Bare lookup table linearizers will be deprecated in DM-28741.")
968  else:
969  linearizer = inputs['linearizer']
970  linearizer.log = self.log
971  inputs['linearizer'] = linearizer
972  else:
973  inputs['linearizer'] = linearize.Linearizer(detector=detector, log=self.log)
974  self.log.warn("Constructing linearizer from cameraGeom information.")
975 
976  if self.config.doDefect is True:
977  if "defects" in inputs and inputs['defects'] is not None:
978  # defects is loaded as a BaseCatalog with columns x0, y0, width, height.
979  # masking expects a list of defects defined by their bounding box
980  if not isinstance(inputs["defects"], Defects):
981  inputs["defects"] = Defects.fromTable(inputs["defects"])
982 
983  # Load the correct style of brighter-fatter kernel, and repack
984  # the information as a numpy array.
985  if self.config.doBrighterFatter:
986  brighterFatterKernel = inputs.pop('newBFKernel', None)
987  if brighterFatterKernel is None:
988  brighterFatterKernel = inputs.get('bfKernel', None)
989 
990  if brighterFatterKernel is not None and not isinstance(brighterFatterKernel, numpy.ndarray):
991  detId = detector.getId()
992  inputs['bfGains'] = brighterFatterKernel.gain
993  # If the kernel is not an ndarray, it's the cp_pipe version
994  # so extract the kernel for this detector, or raise an error
995  if self.config.brighterFatterLevel == 'DETECTOR':
996  if brighterFatterKernel.detectorKernel:
997  inputs['bfKernel'] = brighterFatterKernel.detectorKernel[detId]
998  elif brighterFatterKernel.detectorKernelFromAmpKernels:
999  inputs['bfKernel'] = brighterFatterKernel.detectorKernelFromAmpKernels[detId]
1000  else:
1001  raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1002  else:
1003  # TODO DM-15631 for implementing this
1004  raise NotImplementedError("Per-amplifier brighter-fatter correction not implemented")
1005 
1006  if self.config.doFringe is True and self.fringe.checkFilter(inputs['ccdExposure']):
1007  expId = inputs['ccdExposure'].getInfo().getVisitInfo().getExposureId()
1008  inputs['fringes'] = self.fringe.loadFringes(inputs['fringes'],
1009  expId=expId,
1010  assembler=self.assembleCcd
1011  if self.config.doAssembleIsrExposures else None)
1012  else:
1013  inputs['fringes'] = pipeBase.Struct(fringes=None)
1014 
1015  if self.config.doStrayLight is True and self.strayLight.checkFilter(inputs['ccdExposure']):
1016  if 'strayLightData' not in inputs:
1017  inputs['strayLightData'] = None
1018 
1019  outputs = self.runrun(**inputs)
1020  butlerQC.put(outputs, outputRefs)
1021 
1022  def readIsrData(self, dataRef, rawExposure):
1023  """Retrieve necessary frames for instrument signature removal.
1024 
1025  Pre-fetching all required ISR data products limits the IO
1026  required by the ISR. Any conflict between the calibration data
1027  available and that needed for ISR is also detected prior to
1028  doing processing, allowing it to fail quickly.
1029 
1030  Parameters
1031  ----------
1032  dataRef : `daf.persistence.butlerSubset.ButlerDataRef`
1033  Butler reference of the detector data to be processed
1034  rawExposure : `afw.image.Exposure`
1035  The raw exposure that will later be corrected with the
1036  retrieved calibration data; should not be modified in this
1037  method.
1038 
1039  Returns
1040  -------
1041  result : `lsst.pipe.base.Struct`
1042  Result struct with components (which may be `None`):
1043  - ``bias``: bias calibration frame (`afw.image.Exposure`)
1044  - ``linearizer``: functor for linearization (`ip.isr.linearize.LinearizeBase`)
1045  - ``crosstalkSources``: list of possible crosstalk sources (`list`)
1046  - ``dark``: dark calibration frame (`afw.image.Exposure`)
1047  - ``flat``: flat calibration frame (`afw.image.Exposure`)
1048  - ``bfKernel``: Brighter-Fatter kernel (`numpy.ndarray`)
1049  - ``defects``: list of defects (`lsst.ip.isr.Defects`)
1050  - ``fringes``: `lsst.pipe.base.Struct` with components:
1051  - ``fringes``: fringe calibration frame (`afw.image.Exposure`)
1052  - ``seed``: random seed derived from the ccdExposureId for random
1053  number generator (`uint32`).
1054  - ``opticsTransmission``: `lsst.afw.image.TransmissionCurve`
1055  A ``TransmissionCurve`` that represents the throughput of the optics,
1056  to be evaluated in focal-plane coordinates.
1057  - ``filterTransmission`` : `lsst.afw.image.TransmissionCurve`
1058  A ``TransmissionCurve`` that represents the throughput of the filter
1059  itself, to be evaluated in focal-plane coordinates.
1060  - ``sensorTransmission`` : `lsst.afw.image.TransmissionCurve`
1061  A ``TransmissionCurve`` that represents the throughput of the sensor
1062  itself, to be evaluated in post-assembly trimmed detector coordinates.
1063  - ``atmosphereTransmission`` : `lsst.afw.image.TransmissionCurve`
1064  A ``TransmissionCurve`` that represents the throughput of the
1065  atmosphere, assumed to be spatially constant.
1066  - ``strayLightData`` : `object`
1067  An opaque object containing calibration information for
1068  stray-light correction. If `None`, no correction will be
1069  performed.
1070  - ``illumMaskedImage`` : illumination correction image (`lsst.afw.image.MaskedImage`)
1071 
1072  Raises
1073  ------
1074  NotImplementedError :
1075  Raised if a per-amplifier brighter-fatter kernel is requested by the configuration.
1076  """
1077  try:
1078  dateObs = rawExposure.getInfo().getVisitInfo().getDate()
1079  dateObs = dateObs.toPython().isoformat()
1080  except RuntimeError:
1081  self.log.warn("Unable to identify dateObs for rawExposure.")
1082  dateObs = None
1083 
1084  ccd = rawExposure.getDetector()
1085  filterLabel = rawExposure.getFilterLabel()
1086  rawExposure.mask.addMaskPlane("UNMASKEDNAN") # needed to match pre DM-15862 processing.
1087  biasExposure = (self.getIsrExposuregetIsrExposure(dataRef, self.config.biasDataProductName)
1088  if self.config.doBias else None)
1089  # immediate=True required for functors and linearizers are functors; see ticket DM-6515
1090  linearizer = (dataRef.get("linearizer", immediate=True)
1091  if self.doLinearizedoLinearize(ccd) else None)
1092  if linearizer is not None and not isinstance(linearizer, numpy.ndarray):
1093  linearizer.log = self.log
1094  if isinstance(linearizer, numpy.ndarray):
1095  linearizer = linearize.Linearizer(table=linearizer, detector=ccd)
1096 
1097  crosstalkCalib = None
1098  if self.config.doCrosstalk:
1099  try:
1100  crosstalkCalib = dataRef.get("crosstalk", immediate=True)
1101  except NoResults:
1102  coeffVector = (self.config.crosstalk.crosstalkValues
1103  if self.config.crosstalk.useConfigCoefficients else None)
1104  crosstalkCalib = CrosstalkCalib().fromDetector(ccd, coeffVector=coeffVector)
1105  crosstalkSources = (self.crosstalk.prepCrosstalk(dataRef, crosstalkCalib)
1106  if self.config.doCrosstalk else None)
1107 
1108  darkExposure = (self.getIsrExposuregetIsrExposure(dataRef, self.config.darkDataProductName)
1109  if self.config.doDark else None)
1110  flatExposure = (self.getIsrExposuregetIsrExposure(dataRef, self.config.flatDataProductName,
1111  dateObs=dateObs)
1112  if self.config.doFlat else None)
1113 
1114  brighterFatterKernel = None
1115  brighterFatterGains = None
1116  if self.config.doBrighterFatter is True:
1117  try:
1118  # Use the new-style cp_pipe version of the kernel if it exists
1119  # If using a new-style kernel, always use the self-consistent
1120  # gains, i.e. the ones inside the kernel object itself
1121  brighterFatterKernel = dataRef.get("brighterFatterKernel")
1122  brighterFatterGains = brighterFatterKernel.gain
1123  self.log.info("New style brighter-fatter kernel (brighterFatterKernel) loaded")
1124  except NoResults:
1125  try: # Fall back to the old-style numpy-ndarray style kernel if necessary.
1126  brighterFatterKernel = dataRef.get("bfKernel")
1127  self.log.info("Old style brighter-fatter kernel (np.array) loaded")
1128  except NoResults:
1129  brighterFatterKernel = None
1130  if brighterFatterKernel is not None and not isinstance(brighterFatterKernel, numpy.ndarray):
1131  # If the kernel is not an ndarray, it's the cp_pipe version
1132  # so extract the kernel for this detector, or raise an error
1133  if self.config.brighterFatterLevel == 'DETECTOR':
1134  if brighterFatterKernel.detectorKernel:
1135  brighterFatterKernel = brighterFatterKernel.detectorKernel[ccd.getId()]
1136  elif brighterFatterKernel.detectorKernelFromAmpKernels:
1137  brighterFatterKernel = brighterFatterKernel.detectorKernelFromAmpKernels[ccd.getId()]
1138  else:
1139  raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1140  else:
1141  # TODO DM-15631 for implementing this
1142  raise NotImplementedError("Per-amplifier brighter-fatter correction not implemented")
1143 
1144  defectList = (dataRef.get("defects")
1145  if self.config.doDefect else None)
1146  fringeStruct = (self.fringe.readFringes(dataRef, assembler=self.assembleCcd
1147  if self.config.doAssembleIsrExposures else None)
1148  if self.config.doFringe and self.fringe.checkFilter(rawExposure)
1149  else pipeBase.Struct(fringes=None))
1150 
1151  if self.config.doAttachTransmissionCurve:
1152  opticsTransmission = (dataRef.get("transmission_optics")
1153  if self.config.doUseOpticsTransmission else None)
1154  filterTransmission = (dataRef.get("transmission_filter")
1155  if self.config.doUseFilterTransmission else None)
1156  sensorTransmission = (dataRef.get("transmission_sensor")
1157  if self.config.doUseSensorTransmission else None)
1158  atmosphereTransmission = (dataRef.get("transmission_atmosphere")
1159  if self.config.doUseAtmosphereTransmission else None)
1160  else:
1161  opticsTransmission = None
1162  filterTransmission = None
1163  sensorTransmission = None
1164  atmosphereTransmission = None
1165 
1166  if self.config.doStrayLight:
1167  strayLightData = self.strayLight.readIsrData(dataRef, rawExposure)
1168  else:
1169  strayLightData = None
1170 
1171  illumMaskedImage = (self.getIsrExposuregetIsrExposure(dataRef,
1172  self.config.illuminationCorrectionDataProductName).getMaskedImage()
1173  if (self.config.doIlluminationCorrection
1174  and filterLabel in self.config.illumFilters)
1175  else None)
1176 
1177  # Struct should include only kwargs to run()
1178  return pipeBase.Struct(bias=biasExposure,
1179  linearizer=linearizer,
1180  crosstalk=crosstalkCalib,
1181  crosstalkSources=crosstalkSources,
1182  dark=darkExposure,
1183  flat=flatExposure,
1184  bfKernel=brighterFatterKernel,
1185  bfGains=brighterFatterGains,
1186  defects=defectList,
1187  fringes=fringeStruct,
1188  opticsTransmission=opticsTransmission,
1189  filterTransmission=filterTransmission,
1190  sensorTransmission=sensorTransmission,
1191  atmosphereTransmission=atmosphereTransmission,
1192  strayLightData=strayLightData,
1193  illumMaskedImage=illumMaskedImage
1194  )
1195 
1196  @pipeBase.timeMethod
1197  def run(self, ccdExposure, camera=None, bias=None, linearizer=None,
1198  crosstalk=None, crosstalkSources=None,
1199  dark=None, flat=None, bfKernel=None, bfGains=None, defects=None,
1200  fringes=pipeBase.Struct(fringes=None), opticsTransmission=None, filterTransmission=None,
1201  sensorTransmission=None, atmosphereTransmission=None,
1202  detectorNum=None, strayLightData=None, illumMaskedImage=None,
1203  isGen3=False,
1204  ):
1205  """Perform instrument signature removal on an exposure.
1206 
1207  Steps included in the ISR processing, in order performed, are:
1208  - saturation and suspect pixel masking
1209  - overscan subtraction
1210  - CCD assembly of individual amplifiers
1211  - bias subtraction
1212  - variance image construction
1213  - linearization of non-linear response
1214  - crosstalk masking
1215  - brighter-fatter correction
1216  - dark subtraction
1217  - fringe correction
1218  - stray light subtraction
1219  - flat correction
1220  - masking of known defects and camera specific features
1221  - vignette calculation
1222  - appending transmission curve and distortion model
1223 
1224  Parameters
1225  ----------
1226  ccdExposure : `lsst.afw.image.Exposure`
1227  The raw exposure that is to be run through ISR. The
1228  exposure is modified by this method.
1229  camera : `lsst.afw.cameraGeom.Camera`, optional
1230  The camera geometry for this exposure. Required if ``isGen3`` is
1231  `True` and one or more of ``ccdExposure``, ``bias``, ``dark``, or
1232  ``flat`` does not have an associated detector.
1233  bias : `lsst.afw.image.Exposure`, optional
1234  Bias calibration frame.
1235  linearizer : `lsst.ip.isr.linearize.LinearizeBase`, optional
1236  Functor for linearization.
1237  crosstalk : `lsst.ip.isr.crosstalk.CrosstalkCalib`, optional
1238  Calibration for crosstalk.
1239  crosstalkSources : `list`, optional
1240  List of possible crosstalk sources.
1241  dark : `lsst.afw.image.Exposure`, optional
1242  Dark calibration frame.
1243  flat : `lsst.afw.image.Exposure`, optional
1244  Flat calibration frame.
1245  bfKernel : `numpy.ndarray`, optional
1246  Brighter-fatter kernel.
1247  bfGains : `dict` of `float`, optional
1248  Gains used to override the detector's nominal gains for the
1249  brighter-fatter correction. A dict keyed by amplifier name for
1250  the detector in question.
1251  defects : `lsst.ip.isr.Defects`, optional
1252  List of defects.
1253  fringes : `lsst.pipe.base.Struct`, optional
1254  Struct containing the fringe correction data, with
1255  elements:
1256  - ``fringes``: fringe calibration frame (`afw.image.Exposure`)
1257  - ``seed``: random seed derived from the ccdExposureId for random
1258  number generator (`uint32`)
1259  opticsTransmission: `lsst.afw.image.TransmissionCurve`, optional
1260  A ``TransmissionCurve`` that represents the throughput of the optics,
1261  to be evaluated in focal-plane coordinates.
1262  filterTransmission : `lsst.afw.image.TransmissionCurve`
1263  A ``TransmissionCurve`` that represents the throughput of the filter
1264  itself, to be evaluated in focal-plane coordinates.
1265  sensorTransmission : `lsst.afw.image.TransmissionCurve`
1266  A ``TransmissionCurve`` that represents the throughput of the sensor
1267  itself, to be evaluated in post-assembly trimmed detector coordinates.
1268  atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1269  A ``TransmissionCurve`` that represents the throughput of the
1270  atmosphere, assumed to be spatially constant.
1271  detectorNum : `int`, optional
1272  The integer number for the detector to process.
1273  isGen3 : bool, optional
1274  Flag this call to run() as using the Gen3 butler environment.
1275  strayLightData : `object`, optional
1276  Opaque object containing calibration information for stray-light
1277  correction. If `None`, no correction will be performed.
1278  illumMaskedImage : `lsst.afw.image.MaskedImage`, optional
1279  Illumination correction image.
1280 
1281  Returns
1282  -------
1283  result : `lsst.pipe.base.Struct`
1284  Result struct with component:
1285  - ``exposure`` : `afw.image.Exposure`
1286  The fully ISR corrected exposure.
1287  - ``outputExposure`` : `afw.image.Exposure`
1288  An alias for `exposure`
1289  - ``ossThumb`` : `numpy.ndarray`
1290  Thumbnail image of the exposure after overscan subtraction.
1291  - ``flattenedThumb`` : `numpy.ndarray`
1292  Thumbnail image of the exposure after flat-field correction.
1293 
1294  Raises
1295  ------
1296  RuntimeError
1297  Raised if a configuration option is set to True, but the
1298  required calibration data has not been specified.
1299 
1300  Notes
1301  -----
1302  The current processed exposure can be viewed by setting the
1303  appropriate lsstDebug entries in the `debug.display`
1304  dictionary. The names of these entries correspond to some of
1305  the IsrTaskConfig Boolean options, with the value denoting the
1306  frame to use. The exposure is shown inside the matching
1307  option check and after the processing of that step has
1308  finished. The steps with debug points are:
1309 
1310  doAssembleCcd
1311  doBias
1312  doCrosstalk
1313  doBrighterFatter
1314  doDark
1315  doFringe
1316  doStrayLight
1317  doFlat
1318 
1319  In addition, setting the "postISRCCD" entry displays the
1320  exposure after all ISR processing has finished.
1321 
1322  """
1323 
1324  if isGen3 is True:
1325  # Gen3 currently cannot automatically do configuration overrides.
1326  # DM-15257 looks to discuss this issue.
1327  # Configure input exposures;
1328  if detectorNum is None:
1329  raise RuntimeError("Must supply the detectorNum if running as Gen3.")
1330 
1331  ccdExposure = self.ensureExposureensureExposure(ccdExposure, camera, detectorNum)
1332  bias = self.ensureExposureensureExposure(bias, camera, detectorNum)
1333  dark = self.ensureExposureensureExposure(dark, camera, detectorNum)
1334  flat = self.ensureExposureensureExposure(flat, camera, detectorNum)
1335  else:
1336  if isinstance(ccdExposure, ButlerDataRef):
1337  return self.runDataRefrunDataRef(ccdExposure)
1338 
1339  ccd = ccdExposure.getDetector()
1340  filterLabel = ccdExposure.getFilterLabel()
1341 
1342  if not ccd:
1343  assert not self.config.doAssembleCcd, "You need a Detector to run assembleCcd."
1344  ccd = [FakeAmp(ccdExposure, self.config)]
1345 
1346  # Validate Input
1347  if self.config.doBias and bias is None:
1348  raise RuntimeError("Must supply a bias exposure if config.doBias=True.")
1349  if self.doLinearizedoLinearize(ccd) and linearizer is None:
1350  raise RuntimeError("Must supply a linearizer if config.doLinearize=True for this detector.")
1351  if self.config.doBrighterFatter and bfKernel is None:
1352  raise RuntimeError("Must supply a kernel if config.doBrighterFatter=True.")
1353  if self.config.doDark and dark is None:
1354  raise RuntimeError("Must supply a dark exposure if config.doDark=True.")
1355  if self.config.doFlat and flat is None:
1356  raise RuntimeError("Must supply a flat exposure if config.doFlat=True.")
1357  if self.config.doDefect and defects is None:
1358  raise RuntimeError("Must supply defects if config.doDefect=True.")
1359  if (self.config.doFringe and filterLabel in self.fringe.config.filters
1360  and fringes.fringes is None):
1361  # The `fringes` object needs to be a pipeBase.Struct, as
1362  # we use it as a `dict` for the parameters of
1363  # `FringeTask.run()`. The `fringes.fringes` `list` may
1364  # not be `None` if `doFringe=True`. Otherwise, raise.
1365  raise RuntimeError("Must supply fringe exposure as a pipeBase.Struct.")
1366  if (self.config.doIlluminationCorrection and filterLabel in self.config.illumFilters
1367  and illumMaskedImage is None):
1368  raise RuntimeError("Must supply an illumcor if config.doIlluminationCorrection=True.")
1369 
1370  # Begin ISR processing.
1371  if self.config.doConvertIntToFloat:
1372  self.log.info("Converting exposure to floating point values.")
1373  ccdExposure = self.convertIntToFloatconvertIntToFloat(ccdExposure)
1374 
1375  if self.config.doBias and self.config.doBiasBeforeOverscan:
1376  self.log.info("Applying bias correction.")
1377  isrFunctions.biasCorrection(ccdExposure.getMaskedImage(), bias.getMaskedImage(),
1378  trimToFit=self.config.doTrimToMatchCalib)
1379  self.debugViewdebugView(ccdExposure, "doBias")
1380 
1381  # Amplifier level processing.
1382  overscans = []
1383  for amp in ccd:
1384  # if ccdExposure is one amp, check for coverage to prevent performing ops multiple times
1385  if ccdExposure.getBBox().contains(amp.getBBox()):
1386  # Check for fully masked bad amplifiers, and generate masks for SUSPECT and SATURATED values.
1387  badAmp = self.maskAmplifiermaskAmplifier(ccdExposure, amp, defects)
1388 
1389  if self.config.doOverscan and not badAmp:
1390  # Overscan correction on amp-by-amp basis.
1391  overscanResults = self.overscanCorrectionoverscanCorrection(ccdExposure, amp)
1392  self.log.debug("Corrected overscan for amplifier %s.", amp.getName())
1393  if overscanResults is not None and \
1394  self.config.qa is not None and self.config.qa.saveStats is True:
1395  if isinstance(overscanResults.overscanFit, float):
1396  qaMedian = overscanResults.overscanFit
1397  qaStdev = float("NaN")
1398  else:
1399  qaStats = afwMath.makeStatistics(overscanResults.overscanFit,
1400  afwMath.MEDIAN | afwMath.STDEVCLIP)
1401  qaMedian = qaStats.getValue(afwMath.MEDIAN)
1402  qaStdev = qaStats.getValue(afwMath.STDEVCLIP)
1403 
1404  self.metadata.set(f"FIT MEDIAN {amp.getName()}", qaMedian)
1405  self.metadata.set(f"FIT STDEV {amp.getName()}", qaStdev)
1406  self.log.debug(" Overscan stats for amplifer %s: %f +/- %f",
1407  amp.getName(), qaMedian, qaStdev)
1408 
1409  # Residuals after overscan correction
1410  qaStatsAfter = afwMath.makeStatistics(overscanResults.overscanImage,
1411  afwMath.MEDIAN | afwMath.STDEVCLIP)
1412  qaMedianAfter = qaStatsAfter.getValue(afwMath.MEDIAN)
1413  qaStdevAfter = qaStatsAfter.getValue(afwMath.STDEVCLIP)
1414 
1415  self.metadata.set(f"RESIDUAL MEDIAN {amp.getName()}", qaMedianAfter)
1416  self.metadata.set(f"RESIDUAL STDEV {amp.getName()}", qaStdevAfter)
1417  self.log.debug(" Overscan stats for amplifer %s after correction: %f +/- %f",
1418  amp.getName(), qaMedianAfter, qaStdevAfter)
1419 
1420  ccdExposure.getMetadata().set('OVERSCAN', "Overscan corrected")
1421  else:
1422  if badAmp:
1423  self.log.warn("Amplifier %s is bad.", amp.getName())
1424  overscanResults = None
1425 
1426  overscans.append(overscanResults if overscanResults is not None else None)
1427  else:
1428  self.log.info("Skipped OSCAN for %s.", amp.getName())
1429 
1430  if self.config.doCrosstalk and self.config.doCrosstalkBeforeAssemble:
1431  self.log.info("Applying crosstalk correction.")
1432  self.crosstalk.run(ccdExposure, crosstalk=crosstalk,
1433  crosstalkSources=crosstalkSources, camera=camera)
1434  self.debugViewdebugView(ccdExposure, "doCrosstalk")
1435 
1436  if self.config.doAssembleCcd:
1437  self.log.info("Assembling CCD from amplifiers.")
1438  ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
1439 
1440  if self.config.expectWcs and not ccdExposure.getWcs():
1441  self.log.warn("No WCS found in input exposure.")
1442  self.debugViewdebugView(ccdExposure, "doAssembleCcd")
1443 
1444  ossThumb = None
1445  if self.config.qa.doThumbnailOss:
1446  ossThumb = isrQa.makeThumbnail(ccdExposure, isrQaConfig=self.config.qa)
1447 
1448  if self.config.doBias and not self.config.doBiasBeforeOverscan:
1449  self.log.info("Applying bias correction.")
1450  isrFunctions.biasCorrection(ccdExposure.getMaskedImage(), bias.getMaskedImage(),
1451  trimToFit=self.config.doTrimToMatchCalib)
1452  self.debugViewdebugView(ccdExposure, "doBias")
1453 
1454  if self.config.doVariance:
1455  for amp, overscanResults in zip(ccd, overscans):
1456  if ccdExposure.getBBox().contains(amp.getBBox()):
1457  self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1458  ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1459  if overscanResults is not None:
1460  self.updateVarianceupdateVariance(ampExposure, amp,
1461  overscanImage=overscanResults.overscanImage)
1462  else:
1463  self.updateVarianceupdateVariance(ampExposure, amp,
1464  overscanImage=None)
1465  if self.config.qa is not None and self.config.qa.saveStats is True:
1466  qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1467  afwMath.MEDIAN | afwMath.STDEVCLIP)
1468  self.metadata.set(f"ISR VARIANCE {amp.getName()} MEDIAN",
1469  qaStats.getValue(afwMath.MEDIAN))
1470  self.metadata.set(f"ISR VARIANCE {amp.getName()} STDEV",
1471  qaStats.getValue(afwMath.STDEVCLIP))
1472  self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1473  amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1474  qaStats.getValue(afwMath.STDEVCLIP))
1475 
1476  if self.doLinearizedoLinearize(ccd):
1477  self.log.info("Applying linearizer.")
1478  linearizer.applyLinearity(image=ccdExposure.getMaskedImage().getImage(),
1479  detector=ccd, log=self.log)
1480 
1481  if self.config.doCrosstalk and not self.config.doCrosstalkBeforeAssemble:
1482  self.log.info("Applying crosstalk correction.")
1483  self.crosstalk.run(ccdExposure, crosstalk=crosstalk,
1484  crosstalkSources=crosstalkSources, isTrimmed=True)
1485  self.debugViewdebugView(ccdExposure, "doCrosstalk")
1486 
1487  # Masking block. Optionally mask known defects, NAN/inf pixels, widen trails, and do
1488  # anything else the camera needs. Saturated and suspect pixels have already been masked.
1489  if self.config.doDefect:
1490  self.log.info("Masking defects.")
1491  self.maskDefectmaskDefect(ccdExposure, defects)
1492 
1493  if self.config.numEdgeSuspect > 0:
1494  self.log.info("Masking edges as SUSPECT.")
1495  self.maskEdgesmaskEdges(ccdExposure, numEdgePixels=self.config.numEdgeSuspect,
1496  maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
1497 
1498  if self.config.doNanMasking:
1499  self.log.info("Masking non-finite (NAN, inf) value pixels.")
1500  self.maskNanmaskNan(ccdExposure)
1501 
1502  if self.config.doWidenSaturationTrails:
1503  self.log.info("Widening saturation trails.")
1504  isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
1505 
1506  if self.config.doCameraSpecificMasking:
1507  self.log.info("Masking regions for camera specific reasons.")
1508  self.masking.run(ccdExposure)
1509 
1510  if self.config.doBrighterFatter:
1511  # We need to apply flats and darks before we can interpolate, and we
1512  # need to interpolate before we do B-F, but we do B-F without the
1513  # flats and darks applied so we can work in units of electrons or holes.
1514  # This context manager applies and then removes the darks and flats.
1515  #
1516  # We also do not want to interpolate values here, so operate on temporary
1517  # images so we can apply only the BF-correction and roll back the
1518  # interpolation.
1519  interpExp = ccdExposure.clone()
1520  with self.flatContextflatContext(interpExp, flat, dark):
1521  isrFunctions.interpolateFromMask(
1522  maskedImage=interpExp.getMaskedImage(),
1523  fwhm=self.config.fwhm,
1524  growSaturatedFootprints=self.config.growSaturationFootprintSize,
1525  maskNameList=list(self.config.brighterFatterMaskListToInterpolate)
1526  )
1527  bfExp = interpExp.clone()
1528 
1529  self.log.info("Applying brighter-fatter correction using kernel type %s / gains %s.",
1530  type(bfKernel), type(bfGains))
1531  bfResults = isrFunctions.brighterFatterCorrection(bfExp, bfKernel,
1532  self.config.brighterFatterMaxIter,
1533  self.config.brighterFatterThreshold,
1534  self.config.brighterFatterApplyGain,
1535  bfGains)
1536  if bfResults[1] == self.config.brighterFatterMaxIter:
1537  self.log.warn("Brighter-fatter correction did not converge, final difference %f.",
1538  bfResults[0])
1539  else:
1540  self.log.info("Finished brighter-fatter correction in %d iterations.",
1541  bfResults[1])
1542  image = ccdExposure.getMaskedImage().getImage()
1543  bfCorr = bfExp.getMaskedImage().getImage()
1544  bfCorr -= interpExp.getMaskedImage().getImage()
1545  image += bfCorr
1546 
1547  # Applying the brighter-fatter correction applies a
1548  # convolution to the science image. At the edges this
1549  # convolution may not have sufficient valid pixels to
1550  # produce a valid correction. Mark pixels within the size
1551  # of the brighter-fatter kernel as EDGE to warn of this
1552  # fact.
1553  self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1554  self.maskEdgesmaskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1555  maskPlane="EDGE")
1556 
1557  if self.config.brighterFatterMaskGrowSize > 0:
1558  self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1559  for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1560  isrFunctions.growMasks(ccdExposure.getMask(),
1561  radius=self.config.brighterFatterMaskGrowSize,
1562  maskNameList=maskPlane,
1563  maskValue=maskPlane)
1564 
1565  self.debugViewdebugView(ccdExposure, "doBrighterFatter")
1566 
1567  if self.config.doDark:
1568  self.log.info("Applying dark correction.")
1569  self.darkCorrectiondarkCorrection(ccdExposure, dark)
1570  self.debugViewdebugView(ccdExposure, "doDark")
1571 
1572  if self.config.doFringe and not self.config.fringeAfterFlat:
1573  self.log.info("Applying fringe correction before flat.")
1574  self.fringe.run(ccdExposure, **fringes.getDict())
1575  self.debugViewdebugView(ccdExposure, "doFringe")
1576 
1577  if self.config.doStrayLight and self.strayLight.check(ccdExposure):
1578  self.log.info("Checking strayLight correction.")
1579  self.strayLight.run(ccdExposure, strayLightData)
1580  self.debugViewdebugView(ccdExposure, "doStrayLight")
1581 
1582  if self.config.doFlat:
1583  self.log.info("Applying flat correction.")
1584  self.flatCorrectionflatCorrection(ccdExposure, flat)
1585  self.debugViewdebugView(ccdExposure, "doFlat")
1586 
1587  if self.config.doApplyGains:
1588  self.log.info("Applying gain correction instead of flat.")
1589  isrFunctions.applyGains(ccdExposure, self.config.normalizeGains)
1590 
1591  if self.config.doFringe and self.config.fringeAfterFlat:
1592  self.log.info("Applying fringe correction after flat.")
1593  self.fringe.run(ccdExposure, **fringes.getDict())
1594 
1595  if self.config.doVignette:
1596  self.log.info("Constructing Vignette polygon.")
1597  self.vignettePolygonvignettePolygon = self.vignette.run(ccdExposure)
1598 
1599  if self.config.vignette.doWriteVignettePolygon:
1600  self.setValidPolygonIntersectsetValidPolygonIntersect(ccdExposure, self.vignettePolygonvignettePolygon)
1601 
1602  if self.config.doAttachTransmissionCurve:
1603  self.log.info("Adding transmission curves.")
1604  isrFunctions.attachTransmissionCurve(ccdExposure, opticsTransmission=opticsTransmission,
1605  filterTransmission=filterTransmission,
1606  sensorTransmission=sensorTransmission,
1607  atmosphereTransmission=atmosphereTransmission)
1608 
1609  flattenedThumb = None
1610  if self.config.qa.doThumbnailFlattened:
1611  flattenedThumb = isrQa.makeThumbnail(ccdExposure, isrQaConfig=self.config.qa)
1612 
1613  if self.config.doIlluminationCorrection and filterLabel in self.config.illumFilters:
1614  self.log.info("Performing illumination correction.")
1615  isrFunctions.illuminationCorrection(ccdExposure.getMaskedImage(),
1616  illumMaskedImage, illumScale=self.config.illumScale,
1617  trimToFit=self.config.doTrimToMatchCalib)
1618 
1619  preInterpExp = None
1620  if self.config.doSaveInterpPixels:
1621  preInterpExp = ccdExposure.clone()
1622 
1623  # Reset and interpolate bad pixels.
1624  #
1625  # Large contiguous bad regions (which should have the BAD mask
1626  # bit set) should have their values set to the image median.
1627  # This group should include defects and bad amplifiers. As the
1628  # area covered by these defects are large, there's little
1629  # reason to expect that interpolation would provide a more
1630  # useful value.
1631  #
1632  # Smaller defects can be safely interpolated after the larger
1633  # regions have had their pixel values reset. This ensures
1634  # that the remaining defects adjacent to bad amplifiers (as an
1635  # example) do not attempt to interpolate extreme values.
1636  if self.config.doSetBadRegions:
1637  badPixelCount, badPixelValue = isrFunctions.setBadRegions(ccdExposure)
1638  if badPixelCount > 0:
1639  self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1640 
1641  if self.config.doInterpolate:
1642  self.log.info("Interpolating masked pixels.")
1643  isrFunctions.interpolateFromMask(
1644  maskedImage=ccdExposure.getMaskedImage(),
1645  fwhm=self.config.fwhm,
1646  growSaturatedFootprints=self.config.growSaturationFootprintSize,
1647  maskNameList=list(self.config.maskListToInterpolate)
1648  )
1649 
1650  self.roughZeroPointroughZeroPoint(ccdExposure)
1651 
1652  if self.config.doMeasureBackground:
1653  self.log.info("Measuring background level.")
1654  self.measureBackgroundmeasureBackground(ccdExposure, self.config.qa)
1655 
1656  if self.config.qa is not None and self.config.qa.saveStats is True:
1657  for amp in ccd:
1658  ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1659  qaStats = afwMath.makeStatistics(ampExposure.getImage(),
1660  afwMath.MEDIAN | afwMath.STDEVCLIP)
1661  self.metadata.set("ISR BACKGROUND {} MEDIAN".format(amp.getName()),
1662  qaStats.getValue(afwMath.MEDIAN))
1663  self.metadata.set("ISR BACKGROUND {} STDEV".format(amp.getName()),
1664  qaStats.getValue(afwMath.STDEVCLIP))
1665  self.log.debug(" Background stats for amplifer %s: %f +/- %f",
1666  amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1667  qaStats.getValue(afwMath.STDEVCLIP))
1668 
1669  self.debugViewdebugView(ccdExposure, "postISRCCD")
1670 
1671  return pipeBase.Struct(
1672  exposure=ccdExposure,
1673  ossThumb=ossThumb,
1674  flattenedThumb=flattenedThumb,
1675 
1676  preInterpolatedExposure=preInterpExp,
1677  outputExposure=ccdExposure,
1678  outputOssThumbnail=ossThumb,
1679  outputFlattenedThumbnail=flattenedThumb,
1680  )
1681 
1682  @pipeBase.timeMethod
1683  def runDataRef(self, sensorRef):
1684  """Perform instrument signature removal on a ButlerDataRef of a Sensor.
1685 
1686  This method contains the `CmdLineTask` interface to the ISR
1687  processing. All IO is handled here, freeing the `run()` method
1688  to manage only pixel-level calculations. The steps performed
1689  are:
1690  - Read in necessary detrending/isr/calibration data.
1691  - Process raw exposure in `run()`.
1692  - Persist the ISR-corrected exposure as "postISRCCD" if
1693  config.doWrite=True.
1694 
1695  Parameters
1696  ----------
1697  sensorRef : `daf.persistence.butlerSubset.ButlerDataRef`
1698  DataRef of the detector data to be processed
1699 
1700  Returns
1701  -------
1702  result : `lsst.pipe.base.Struct`
1703  Result struct with component:
1704  - ``exposure`` : `afw.image.Exposure`
1705  The fully ISR corrected exposure.
1706 
1707  Raises
1708  ------
1709  RuntimeError
1710  Raised if a configuration option is set to True, but the
1711  required calibration data does not exist.
1712 
1713  """
1714  self.log.info("Performing ISR on sensor %s.", sensorRef.dataId)
1715 
1716  ccdExposure = sensorRef.get(self.config.datasetType)
1717 
1718  camera = sensorRef.get("camera")
1719  isrData = self.readIsrDatareadIsrData(sensorRef, ccdExposure)
1720 
1721  result = self.runrun(ccdExposure, camera=camera, **isrData.getDict())
1722 
1723  if self.config.doWrite:
1724  sensorRef.put(result.exposure, "postISRCCD")
1725  if result.preInterpolatedExposure is not None:
1726  sensorRef.put(result.preInterpolatedExposure, "postISRCCD_uninterpolated")
1727  if result.ossThumb is not None:
1728  isrQa.writeThumbnail(sensorRef, result.ossThumb, "ossThumb")
1729  if result.flattenedThumb is not None:
1730  isrQa.writeThumbnail(sensorRef, result.flattenedThumb, "flattenedThumb")
1731 
1732  return result
1733 
1734  def getIsrExposure(self, dataRef, datasetType, dateObs=None, immediate=True):
1735  """Retrieve a calibration dataset for removing instrument signature.
1736 
1737  Parameters
1738  ----------
1739 
1740  dataRef : `daf.persistence.butlerSubset.ButlerDataRef`
1741  DataRef of the detector data to find calibration datasets
1742  for.
1743  datasetType : `str`
1744  Type of dataset to retrieve (e.g. 'bias', 'flat', etc).
1745  dateObs : `str`, optional
1746  Date of the observation. Used to correct butler failures
1747  when using fallback filters.
1748  immediate : `Bool`
1749  If True, disable butler proxies to enable error handling
1750  within this routine.
1751 
1752  Returns
1753  -------
1754  exposure : `lsst.afw.image.Exposure`
1755  Requested calibration frame.
1756 
1757  Raises
1758  ------
1759  RuntimeError
1760  Raised if no matching calibration frame can be found.
1761  """
1762  try:
1763  exp = dataRef.get(datasetType, immediate=immediate)
1764  except Exception as exc1:
1765  if not self.config.fallbackFilterName:
1766  raise RuntimeError("Unable to retrieve %s for %s: %s." % (datasetType, dataRef.dataId, exc1))
1767  try:
1768  if self.config.useFallbackDate and dateObs:
1769  exp = dataRef.get(datasetType, filter=self.config.fallbackFilterName,
1770  dateObs=dateObs, immediate=immediate)
1771  else:
1772  exp = dataRef.get(datasetType, filter=self.config.fallbackFilterName, immediate=immediate)
1773  except Exception as exc2:
1774  raise RuntimeError("Unable to retrieve %s for %s, even with fallback filter %s: %s AND %s." %
1775  (datasetType, dataRef.dataId, self.config.fallbackFilterName, exc1, exc2))
1776  self.log.warn("Using fallback calibration from filter %s.", self.config.fallbackFilterName)
1777 
1778  if self.config.doAssembleIsrExposures:
1779  exp = self.assembleCcd.assembleCcd(exp)
1780  return exp
1781 
1782  def ensureExposure(self, inputExp, camera, detectorNum):
1783  """Ensure that the data returned by Butler is a fully constructed exposure.
1784 
1785  ISR requires exposure-level image data for historical reasons, so if we did
1786  not recieve that from Butler, construct it from what we have, modifying the
1787  input in place.
1788 
1789  Parameters
1790  ----------
1791  inputExp : `lsst.afw.image.Exposure`, `lsst.afw.image.DecoratedImageU`, or
1792  `lsst.afw.image.ImageF`
1793  The input data structure obtained from Butler.
1794  camera : `lsst.afw.cameraGeom.camera`
1795  The camera associated with the image. Used to find the appropriate
1796  detector.
1797  detectorNum : `int`
1798  The detector this exposure should match.
1799 
1800  Returns
1801  -------
1802  inputExp : `lsst.afw.image.Exposure`
1803  The re-constructed exposure, with appropriate detector parameters.
1804 
1805  Raises
1806  ------
1807  TypeError
1808  Raised if the input data cannot be used to construct an exposure.
1809  """
1810  if isinstance(inputExp, afwImage.DecoratedImageU):
1811  inputExp = afwImage.makeExposure(afwImage.makeMaskedImage(inputExp))
1812  elif isinstance(inputExp, afwImage.ImageF):
1813  inputExp = afwImage.makeExposure(afwImage.makeMaskedImage(inputExp))
1814  elif isinstance(inputExp, afwImage.MaskedImageF):
1815  inputExp = afwImage.makeExposure(inputExp)
1816  elif isinstance(inputExp, afwImage.Exposure):
1817  pass
1818  elif inputExp is None:
1819  # Assume this will be caught by the setup if it is a problem.
1820  return inputExp
1821  else:
1822  raise TypeError("Input Exposure is not known type in isrTask.ensureExposure: %s." %
1823  (type(inputExp), ))
1824 
1825  if inputExp.getDetector() is None:
1826  inputExp.setDetector(camera[detectorNum])
1827 
1828  return inputExp
1829 
1830  def convertIntToFloat(self, exposure):
1831  """Convert exposure image from uint16 to float.
1832 
1833  If the exposure does not need to be converted, the input is
1834  immediately returned. For exposures that are converted to use
1835  floating point pixels, the variance is set to unity and the
1836  mask to zero.
1837 
1838  Parameters
1839  ----------
1840  exposure : `lsst.afw.image.Exposure`
1841  The raw exposure to be converted.
1842 
1843  Returns
1844  -------
1845  newexposure : `lsst.afw.image.Exposure`
1846  The input ``exposure``, converted to floating point pixels.
1847 
1848  Raises
1849  ------
1850  RuntimeError
1851  Raised if the exposure type cannot be converted to float.
1852 
1853  """
1854  if isinstance(exposure, afwImage.ExposureF):
1855  # Nothing to be done
1856  self.log.debug("Exposure already of type float.")
1857  return exposure
1858  if not hasattr(exposure, "convertF"):
1859  raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
1860 
1861  newexposure = exposure.convertF()
1862  newexposure.variance[:] = 1
1863  newexposure.mask[:] = 0x0
1864 
1865  return newexposure
1866 
1867  def maskAmplifier(self, ccdExposure, amp, defects):
1868  """Identify bad amplifiers, saturated and suspect pixels.
1869 
1870  Parameters
1871  ----------
1872  ccdExposure : `lsst.afw.image.Exposure`
1873  Input exposure to be masked.
1874  amp : `lsst.afw.table.AmpInfoCatalog`
1875  Catalog of parameters defining the amplifier on this
1876  exposure to mask.
1877  defects : `lsst.ip.isr.Defects`
1878  List of defects. Used to determine if the entire
1879  amplifier is bad.
1880 
1881  Returns
1882  -------
1883  badAmp : `Bool`
1884  If this is true, the entire amplifier area is covered by
1885  defects and unusable.
1886 
1887  """
1888  maskedImage = ccdExposure.getMaskedImage()
1889 
1890  badAmp = False
1891 
1892  # Check if entire amp region is defined as a defect (need to use amp.getBBox() for correct
1893  # comparison with current defects definition.
1894  if defects is not None:
1895  badAmp = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
1896 
1897  # In the case of a bad amp, we will set mask to "BAD" (here use amp.getRawBBox() for correct
1898  # association with pixels in current ccdExposure).
1899  if badAmp:
1900  dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
1901  afwImage.PARENT)
1902  maskView = dataView.getMask()
1903  maskView |= maskView.getPlaneBitMask("BAD")
1904  del maskView
1905  return badAmp
1906 
1907  # Mask remaining defects after assembleCcd() to allow for defects that cross amplifier boundaries.
1908  # Saturation and suspect pixels can be masked now, though.
1909  limits = dict()
1910  if self.config.doSaturation and not badAmp:
1911  limits.update({self.config.saturatedMaskName: amp.getSaturation()})
1912  if self.config.doSuspect and not badAmp:
1913  limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
1914  if math.isfinite(self.config.saturation):
1915  limits.update({self.config.saturatedMaskName: self.config.saturation})
1916 
1917  for maskName, maskThreshold in limits.items():
1918  if not math.isnan(maskThreshold):
1919  dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
1920  isrFunctions.makeThresholdMask(
1921  maskedImage=dataView,
1922  threshold=maskThreshold,
1923  growFootprints=0,
1924  maskName=maskName
1925  )
1926 
1927  # Determine if we've fully masked this amplifier with SUSPECT and SAT pixels.
1928  maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
1929  afwImage.PARENT)
1930  maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
1931  self.config.suspectMaskName])
1932  if numpy.all(maskView.getArray() & maskVal > 0):
1933  badAmp = True
1934  maskView |= maskView.getPlaneBitMask("BAD")
1935 
1936  return badAmp
1937 
1938  def overscanCorrection(self, ccdExposure, amp):
1939  """Apply overscan correction in place.
1940 
1941  This method does initial pixel rejection of the overscan
1942  region. The overscan can also be optionally segmented to
1943  allow for discontinuous overscan responses to be fit
1944  separately. The actual overscan subtraction is performed by
1945  the `lsst.ip.isr.isrFunctions.overscanCorrection` function,
1946  which is called here after the amplifier is preprocessed.
1947 
1948  Parameters
1949  ----------
1950  ccdExposure : `lsst.afw.image.Exposure`
1951  Exposure to have overscan correction performed.
1952  amp : `lsst.afw.cameraGeom.Amplifer`
1953  The amplifier to consider while correcting the overscan.
1954 
1955  Returns
1956  -------
1957  overscanResults : `lsst.pipe.base.Struct`
1958  Result struct with components:
1959  - ``imageFit`` : scalar or `lsst.afw.image.Image`
1960  Value or fit subtracted from the amplifier image data.
1961  - ``overscanFit`` : scalar or `lsst.afw.image.Image`
1962  Value or fit subtracted from the overscan image data.
1963  - ``overscanImage`` : `lsst.afw.image.Image`
1964  Image of the overscan region with the overscan
1965  correction applied. This quantity is used to estimate
1966  the amplifier read noise empirically.
1967 
1968  Raises
1969  ------
1970  RuntimeError
1971  Raised if the ``amp`` does not contain raw pixel information.
1972 
1973  See Also
1974  --------
1975  lsst.ip.isr.isrFunctions.overscanCorrection
1976  """
1977  if amp.getRawHorizontalOverscanBBox().isEmpty():
1978  self.log.info("ISR_OSCAN: No overscan region. Not performing overscan correction.")
1979  return None
1980 
1981  statControl = afwMath.StatisticsControl()
1982  statControl.setAndMask(ccdExposure.mask.getPlaneBitMask("SAT"))
1983 
1984  # Determine the bounding boxes
1985  dataBBox = amp.getRawDataBBox()
1986  oscanBBox = amp.getRawHorizontalOverscanBBox()
1987  dx0 = 0
1988  dx1 = 0
1989 
1990  prescanBBox = amp.getRawPrescanBBox()
1991  if (oscanBBox.getBeginX() > prescanBBox.getBeginX()): # amp is at the right
1992  dx0 += self.config.overscanNumLeadingColumnsToSkip
1993  dx1 -= self.config.overscanNumTrailingColumnsToSkip
1994  else:
1995  dx0 += self.config.overscanNumTrailingColumnsToSkip
1996  dx1 -= self.config.overscanNumLeadingColumnsToSkip
1997 
1998  # Determine if we need to work on subregions of the amplifier and overscan.
1999  imageBBoxes = []
2000  overscanBBoxes = []
2001 
2002  if ((self.config.overscanBiasJump
2003  and self.config.overscanBiasJumpLocation)
2004  and (ccdExposure.getMetadata().exists(self.config.overscanBiasJumpKeyword)
2005  and ccdExposure.getMetadata().getScalar(self.config.overscanBiasJumpKeyword) in
2006  self.config.overscanBiasJumpDevices)):
2007  if amp.getReadoutCorner() in (ReadoutCorner.LL, ReadoutCorner.LR):
2008  yLower = self.config.overscanBiasJumpLocation
2009  yUpper = dataBBox.getHeight() - yLower
2010  else:
2011  yUpper = self.config.overscanBiasJumpLocation
2012  yLower = dataBBox.getHeight() - yUpper
2013 
2014  imageBBoxes.append(lsst.geom.Box2I(dataBBox.getBegin(),
2015  lsst.geom.Extent2I(dataBBox.getWidth(), yLower)))
2016  overscanBBoxes.append(lsst.geom.Box2I(oscanBBox.getBegin() + lsst.geom.Extent2I(dx0, 0),
2017  lsst.geom.Extent2I(oscanBBox.getWidth() - dx0 + dx1,
2018  yLower)))
2019 
2020  imageBBoxes.append(lsst.geom.Box2I(dataBBox.getBegin() + lsst.geom.Extent2I(0, yLower),
2021  lsst.geom.Extent2I(dataBBox.getWidth(), yUpper)))
2022  overscanBBoxes.append(lsst.geom.Box2I(oscanBBox.getBegin() + lsst.geom.Extent2I(dx0, yLower),
2023  lsst.geom.Extent2I(oscanBBox.getWidth() - dx0 + dx1,
2024  yUpper)))
2025  else:
2026  imageBBoxes.append(lsst.geom.Box2I(dataBBox.getBegin(),
2027  lsst.geom.Extent2I(dataBBox.getWidth(), dataBBox.getHeight())))
2028  overscanBBoxes.append(lsst.geom.Box2I(oscanBBox.getBegin() + lsst.geom.Extent2I(dx0, 0),
2029  lsst.geom.Extent2I(oscanBBox.getWidth() - dx0 + dx1,
2030  oscanBBox.getHeight())))
2031 
2032  # Perform overscan correction on subregions, ensuring saturated pixels are masked.
2033  for imageBBox, overscanBBox in zip(imageBBoxes, overscanBBoxes):
2034  ampImage = ccdExposure.maskedImage[imageBBox]
2035  overscanImage = ccdExposure.maskedImage[overscanBBox]
2036 
2037  overscanArray = overscanImage.image.array
2038  median = numpy.ma.median(numpy.ma.masked_where(overscanImage.mask.array, overscanArray))
2039  bad = numpy.where(numpy.abs(overscanArray - median) > self.config.overscanMaxDev)
2040  overscanImage.mask.array[bad] = overscanImage.mask.getPlaneBitMask("SAT")
2041 
2042  statControl = afwMath.StatisticsControl()
2043  statControl.setAndMask(ccdExposure.mask.getPlaneBitMask("SAT"))
2044 
2045  overscanResults = self.overscan.run(ampImage.getImage(), overscanImage, amp)
2046 
2047  # Measure average overscan levels and record them in the metadata.
2048  levelStat = afwMath.MEDIAN
2049  sigmaStat = afwMath.STDEVCLIP
2050 
2051  sctrl = afwMath.StatisticsControl(self.config.qa.flatness.clipSigma,
2052  self.config.qa.flatness.nIter)
2053  metadata = ccdExposure.getMetadata()
2054  ampNum = amp.getName()
2055  # if self.config.overscanFitType in ("MEDIAN", "MEAN", "MEANCLIP"):
2056  if isinstance(overscanResults.overscanFit, float):
2057  metadata.set("ISR_OSCAN_LEVEL%s" % ampNum, overscanResults.overscanFit)
2058  metadata.set("ISR_OSCAN_SIGMA%s" % ampNum, 0.0)
2059  else:
2060  stats = afwMath.makeStatistics(overscanResults.overscanFit, levelStat | sigmaStat, sctrl)
2061  metadata.set("ISR_OSCAN_LEVEL%s" % ampNum, stats.getValue(levelStat))
2062  metadata.set("ISR_OSCAN_SIGMA%s" % ampNum, stats.getValue(sigmaStat))
2063 
2064  return overscanResults
2065 
2066  def updateVariance(self, ampExposure, amp, overscanImage=None):
2067  """Set the variance plane using the amplifier gain and read noise
2068 
2069  The read noise is calculated from the ``overscanImage`` if the
2070  ``doEmpiricalReadNoise`` option is set in the configuration; otherwise
2071  the value from the amplifier data is used.
2072 
2073  Parameters
2074  ----------
2075  ampExposure : `lsst.afw.image.Exposure`
2076  Exposure to process.
2077  amp : `lsst.afw.table.AmpInfoRecord` or `FakeAmp`
2078  Amplifier detector data.
2079  overscanImage : `lsst.afw.image.MaskedImage`, optional.
2080  Image of overscan, required only for empirical read noise.
2081 
2082  See also
2083  --------
2084  lsst.ip.isr.isrFunctions.updateVariance
2085  """
2086  maskPlanes = [self.config.saturatedMaskName, self.config.suspectMaskName]
2087  gain = amp.getGain()
2088 
2089  if math.isnan(gain):
2090  gain = 1.0
2091  self.log.warn("Gain set to NAN! Updating to 1.0 to generate Poisson variance.")
2092  elif gain <= 0:
2093  patchedGain = 1.0
2094  self.log.warn("Gain for amp %s == %g <= 0; setting to %f.",
2095  amp.getName(), gain, patchedGain)
2096  gain = patchedGain
2097 
2098  if self.config.doEmpiricalReadNoise and overscanImage is None:
2099  self.log.info("Overscan is none for EmpiricalReadNoise.")
2100 
2101  if self.config.doEmpiricalReadNoise and overscanImage is not None:
2102  stats = afwMath.StatisticsControl()
2103  stats.setAndMask(overscanImage.mask.getPlaneBitMask(maskPlanes))
2104  readNoise = afwMath.makeStatistics(overscanImage, afwMath.STDEVCLIP, stats).getValue()
2105  self.log.info("Calculated empirical read noise for amp %s: %f.",
2106  amp.getName(), readNoise)
2107  else:
2108  readNoise = amp.getReadNoise()
2109 
2110  isrFunctions.updateVariance(
2111  maskedImage=ampExposure.getMaskedImage(),
2112  gain=gain,
2113  readNoise=readNoise,
2114  )
2115 
2116  def darkCorrection(self, exposure, darkExposure, invert=False):
2117  """Apply dark correction in place.
2118 
2119  Parameters
2120  ----------
2121  exposure : `lsst.afw.image.Exposure`
2122  Exposure to process.
2123  darkExposure : `lsst.afw.image.Exposure`
2124  Dark exposure of the same size as ``exposure``.
2125  invert : `Bool`, optional
2126  If True, re-add the dark to an already corrected image.
2127 
2128  Raises
2129  ------
2130  RuntimeError
2131  Raised if either ``exposure`` or ``darkExposure`` do not
2132  have their dark time defined.
2133 
2134  See Also
2135  --------
2136  lsst.ip.isr.isrFunctions.darkCorrection
2137  """
2138  expScale = exposure.getInfo().getVisitInfo().getDarkTime()
2139  if math.isnan(expScale):
2140  raise RuntimeError("Exposure darktime is NAN.")
2141  if darkExposure.getInfo().getVisitInfo() is not None \
2142  and not math.isnan(darkExposure.getInfo().getVisitInfo().getDarkTime()):
2143  darkScale = darkExposure.getInfo().getVisitInfo().getDarkTime()
2144  else:
2145  # DM-17444: darkExposure.getInfo.getVisitInfo() is None
2146  # so getDarkTime() does not exist.
2147  self.log.warn("darkExposure.getInfo().getVisitInfo() does not exist. Using darkScale = 1.0.")
2148  darkScale = 1.0
2149 
2150  isrFunctions.darkCorrection(
2151  maskedImage=exposure.getMaskedImage(),
2152  darkMaskedImage=darkExposure.getMaskedImage(),
2153  expScale=expScale,
2154  darkScale=darkScale,
2155  invert=invert,
2156  trimToFit=self.config.doTrimToMatchCalib
2157  )
2158 
2159  def doLinearize(self, detector):
2160  """Check if linearization is needed for the detector cameraGeom.
2161 
2162  Checks config.doLinearize and the linearity type of the first
2163  amplifier.
2164 
2165  Parameters
2166  ----------
2167  detector : `lsst.afw.cameraGeom.Detector`
2168  Detector to get linearity type from.
2169 
2170  Returns
2171  -------
2172  doLinearize : `Bool`
2173  If True, linearization should be performed.
2174  """
2175  return self.config.doLinearize and \
2176  detector.getAmplifiers()[0].getLinearityType() != NullLinearityType
2177 
2178  def flatCorrection(self, exposure, flatExposure, invert=False):
2179  """Apply flat correction in place.
2180 
2181  Parameters
2182  ----------
2183  exposure : `lsst.afw.image.Exposure`
2184  Exposure to process.
2185  flatExposure : `lsst.afw.image.Exposure`
2186  Flat exposure of the same size as ``exposure``.
2187  invert : `Bool`, optional
2188  If True, unflatten an already flattened image.
2189 
2190  See Also
2191  --------
2192  lsst.ip.isr.isrFunctions.flatCorrection
2193  """
2194  isrFunctions.flatCorrection(
2195  maskedImage=exposure.getMaskedImage(),
2196  flatMaskedImage=flatExposure.getMaskedImage(),
2197  scalingType=self.config.flatScalingType,
2198  userScale=self.config.flatUserScale,
2199  invert=invert,
2200  trimToFit=self.config.doTrimToMatchCalib
2201  )
2202 
2203  def saturationDetection(self, exposure, amp):
2204  """Detect saturated pixels and mask them using mask plane config.saturatedMaskName, in place.
2205 
2206  Parameters
2207  ----------
2208  exposure : `lsst.afw.image.Exposure`
2209  Exposure to process. Only the amplifier DataSec is processed.
2210  amp : `lsst.afw.table.AmpInfoCatalog`
2211  Amplifier detector data.
2212 
2213  See Also
2214  --------
2215  lsst.ip.isr.isrFunctions.makeThresholdMask
2216  """
2217  if not math.isnan(amp.getSaturation()):
2218  maskedImage = exposure.getMaskedImage()
2219  dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2220  isrFunctions.makeThresholdMask(
2221  maskedImage=dataView,
2222  threshold=amp.getSaturation(),
2223  growFootprints=0,
2224  maskName=self.config.saturatedMaskName,
2225  )
2226 
2227  def saturationInterpolation(self, exposure):
2228  """Interpolate over saturated pixels, in place.
2229 
2230  This method should be called after `saturationDetection`, to
2231  ensure that the saturated pixels have been identified in the
2232  SAT mask. It should also be called after `assembleCcd`, since
2233  saturated regions may cross amplifier boundaries.
2234 
2235  Parameters
2236  ----------
2237  exposure : `lsst.afw.image.Exposure`
2238  Exposure to process.
2239 
2240  See Also
2241  --------
2242  lsst.ip.isr.isrTask.saturationDetection
2243  lsst.ip.isr.isrFunctions.interpolateFromMask
2244  """
2245  isrFunctions.interpolateFromMask(
2246  maskedImage=exposure.getMaskedImage(),
2247  fwhm=self.config.fwhm,
2248  growSaturatedFootprints=self.config.growSaturationFootprintSize,
2249  maskNameList=list(self.config.saturatedMaskName),
2250  )
2251 
2252  def suspectDetection(self, exposure, amp):
2253  """Detect suspect pixels and mask them using mask plane config.suspectMaskName, in place.
2254 
2255  Parameters
2256  ----------
2257  exposure : `lsst.afw.image.Exposure`
2258  Exposure to process. Only the amplifier DataSec is processed.
2259  amp : `lsst.afw.table.AmpInfoCatalog`
2260  Amplifier detector data.
2261 
2262  See Also
2263  --------
2264  lsst.ip.isr.isrFunctions.makeThresholdMask
2265 
2266  Notes
2267  -----
2268  Suspect pixels are pixels whose value is greater than amp.getSuspectLevel().
2269  This is intended to indicate pixels that may be affected by unknown systematics;
2270  for example if non-linearity corrections above a certain level are unstable
2271  then that would be a useful value for suspectLevel. A value of `nan` indicates
2272  that no such level exists and no pixels are to be masked as suspicious.
2273  """
2274  suspectLevel = amp.getSuspectLevel()
2275  if math.isnan(suspectLevel):
2276  return
2277 
2278  maskedImage = exposure.getMaskedImage()
2279  dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2280  isrFunctions.makeThresholdMask(
2281  maskedImage=dataView,
2282  threshold=suspectLevel,
2283  growFootprints=0,
2284  maskName=self.config.suspectMaskName,
2285  )
2286 
2287  def maskDefect(self, exposure, defectBaseList):
2288  """Mask defects using mask plane "BAD", in place.
2289 
2290  Parameters
2291  ----------
2292  exposure : `lsst.afw.image.Exposure`
2293  Exposure to process.
2294  defectBaseList : `lsst.ip.isr.Defects` or `list` of
2295  `lsst.afw.image.DefectBase`.
2296  List of defects to mask.
2297 
2298  Notes
2299  -----
2300  Call this after CCD assembly, since defects may cross amplifier boundaries.
2301  """
2302  maskedImage = exposure.getMaskedImage()
2303  if not isinstance(defectBaseList, Defects):
2304  # Promotes DefectBase to Defect
2305  defectList = Defects(defectBaseList)
2306  else:
2307  defectList = defectBaseList
2308  defectList.maskPixels(maskedImage, maskName="BAD")
2309 
2310  def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
2311  """Mask edge pixels with applicable mask plane.
2312 
2313  Parameters
2314  ----------
2315  exposure : `lsst.afw.image.Exposure`
2316  Exposure to process.
2317  numEdgePixels : `int`, optional
2318  Number of edge pixels to mask.
2319  maskPlane : `str`, optional
2320  Mask plane name to use.
2321  level : `str`, optional
2322  Level at which to mask edges.
2323  """
2324  maskedImage = exposure.getMaskedImage()
2325  maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
2326 
2327  if numEdgePixels > 0:
2328  if level == 'DETECTOR':
2329  boxes = [maskedImage.getBBox()]
2330  elif level == 'AMP':
2331  boxes = [amp.getBBox() for amp in exposure.getDetector()]
2332 
2333  for box in boxes:
2334  # This makes a bbox numEdgeSuspect pixels smaller than the image on each side
2335  subImage = maskedImage[box]
2336  box.grow(-numEdgePixels)
2337  # Mask pixels outside box
2338  SourceDetectionTask.setEdgeBits(
2339  subImage,
2340  box,
2341  maskBitMask)
2342 
2343  def maskAndInterpolateDefects(self, exposure, defectBaseList):
2344  """Mask and interpolate defects using mask plane "BAD", in place.
2345 
2346  Parameters
2347  ----------
2348  exposure : `lsst.afw.image.Exposure`
2349  Exposure to process.
2350  defectBaseList : `lsst.ip.isr.Defects` or `list` of
2351  `lsst.afw.image.DefectBase`.
2352  List of defects to mask and interpolate.
2353 
2354  See Also
2355  --------
2356  lsst.ip.isr.isrTask.maskDefect
2357  """
2358  self.maskDefectmaskDefect(exposure, defectBaseList)
2359  self.maskEdgesmaskEdges(exposure, numEdgePixels=self.config.numEdgeSuspect,
2360  maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
2361  isrFunctions.interpolateFromMask(
2362  maskedImage=exposure.getMaskedImage(),
2363  fwhm=self.config.fwhm,
2364  growSaturatedFootprints=0,
2365  maskNameList=["BAD"],
2366  )
2367 
2368  def maskNan(self, exposure):
2369  """Mask NaNs using mask plane "UNMASKEDNAN", in place.
2370 
2371  Parameters
2372  ----------
2373  exposure : `lsst.afw.image.Exposure`
2374  Exposure to process.
2375 
2376  Notes
2377  -----
2378  We mask over all non-finite values (NaN, inf), including those
2379  that are masked with other bits (because those may or may not be
2380  interpolated over later, and we want to remove all NaN/infs).
2381  Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
2382  preserve the historical name.
2383  """
2384  maskedImage = exposure.getMaskedImage()
2385 
2386  # Find and mask NaNs
2387  maskedImage.getMask().addMaskPlane("UNMASKEDNAN")
2388  maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
2389  numNans = maskNans(maskedImage, maskVal)
2390  self.metadata.set("NUMNANS", numNans)
2391  if numNans > 0:
2392  self.log.warn("There were %d unmasked NaNs.", numNans)
2393 
2394  def maskAndInterpolateNan(self, exposure):
2395  """"Mask and interpolate NaN/infs using mask plane "UNMASKEDNAN",
2396  in place.
2397 
2398  Parameters
2399  ----------
2400  exposure : `lsst.afw.image.Exposure`
2401  Exposure to process.
2402 
2403  See Also
2404  --------
2405  lsst.ip.isr.isrTask.maskNan
2406  """
2407  self.maskNanmaskNan(exposure)
2408  isrFunctions.interpolateFromMask(
2409  maskedImage=exposure.getMaskedImage(),
2410  fwhm=self.config.fwhm,
2411  growSaturatedFootprints=0,
2412  maskNameList=["UNMASKEDNAN"],
2413  )
2414 
2415  def measureBackground(self, exposure, IsrQaConfig=None):
2416  """Measure the image background in subgrids, for quality control purposes.
2417 
2418  Parameters
2419  ----------
2420  exposure : `lsst.afw.image.Exposure`
2421  Exposure to process.
2422  IsrQaConfig : `lsst.ip.isr.isrQa.IsrQaConfig`
2423  Configuration object containing parameters on which background
2424  statistics and subgrids to use.
2425  """
2426  if IsrQaConfig is not None:
2427  statsControl = afwMath.StatisticsControl(IsrQaConfig.flatness.clipSigma,
2428  IsrQaConfig.flatness.nIter)
2429  maskVal = exposure.getMaskedImage().getMask().getPlaneBitMask(["BAD", "SAT", "DETECTED"])
2430  statsControl.setAndMask(maskVal)
2431  maskedImage = exposure.getMaskedImage()
2432  stats = afwMath.makeStatistics(maskedImage, afwMath.MEDIAN | afwMath.STDEVCLIP, statsControl)
2433  skyLevel = stats.getValue(afwMath.MEDIAN)
2434  skySigma = stats.getValue(afwMath.STDEVCLIP)
2435  self.log.info("Flattened sky level: %f +/- %f.", skyLevel, skySigma)
2436  metadata = exposure.getMetadata()
2437  metadata.set('SKYLEVEL', skyLevel)
2438  metadata.set('SKYSIGMA', skySigma)
2439 
2440  # calcluating flatlevel over the subgrids
2441  stat = afwMath.MEANCLIP if IsrQaConfig.flatness.doClip else afwMath.MEAN
2442  meshXHalf = int(IsrQaConfig.flatness.meshX/2.)
2443  meshYHalf = int(IsrQaConfig.flatness.meshY/2.)
2444  nX = int((exposure.getWidth() + meshXHalf) / IsrQaConfig.flatness.meshX)
2445  nY = int((exposure.getHeight() + meshYHalf) / IsrQaConfig.flatness.meshY)
2446  skyLevels = numpy.zeros((nX, nY))
2447 
2448  for j in range(nY):
2449  yc = meshYHalf + j * IsrQaConfig.flatness.meshY
2450  for i in range(nX):
2451  xc = meshXHalf + i * IsrQaConfig.flatness.meshX
2452 
2453  xLLC = xc - meshXHalf
2454  yLLC = yc - meshYHalf
2455  xURC = xc + meshXHalf - 1
2456  yURC = yc + meshYHalf - 1
2457 
2458  bbox = lsst.geom.Box2I(lsst.geom.Point2I(xLLC, yLLC), lsst.geom.Point2I(xURC, yURC))
2459  miMesh = maskedImage.Factory(exposure.getMaskedImage(), bbox, afwImage.LOCAL)
2460 
2461  skyLevels[i, j] = afwMath.makeStatistics(miMesh, stat, statsControl).getValue()
2462 
2463  good = numpy.where(numpy.isfinite(skyLevels))
2464  skyMedian = numpy.median(skyLevels[good])
2465  flatness = (skyLevels[good] - skyMedian) / skyMedian
2466  flatness_rms = numpy.std(flatness)
2467  flatness_pp = flatness.max() - flatness.min() if len(flatness) > 0 else numpy.nan
2468 
2469  self.log.info("Measuring sky levels in %dx%d grids: %f.", nX, nY, skyMedian)
2470  self.log.info("Sky flatness in %dx%d grids - pp: %f rms: %f.",
2471  nX, nY, flatness_pp, flatness_rms)
2472 
2473  metadata.set('FLATNESS_PP', float(flatness_pp))
2474  metadata.set('FLATNESS_RMS', float(flatness_rms))
2475  metadata.set('FLATNESS_NGRIDS', '%dx%d' % (nX, nY))
2476  metadata.set('FLATNESS_MESHX', IsrQaConfig.flatness.meshX)
2477  metadata.set('FLATNESS_MESHY', IsrQaConfig.flatness.meshY)
2478 
2479  def roughZeroPoint(self, exposure):
2480  """Set an approximate magnitude zero point for the exposure.
2481 
2482  Parameters
2483  ----------
2484  exposure : `lsst.afw.image.Exposure`
2485  Exposure to process.
2486  """
2487  filterLabel = exposure.getFilterLabel()
2488  if filterLabel in self.config.fluxMag0T1:
2489  fluxMag0 = self.config.fluxMag0T1[filterLabel]
2490  else:
2491  self.log.warn("No rough magnitude zero point set for filter %s.", filterLabel)
2492  fluxMag0 = self.config.defaultFluxMag0T1
2493 
2494  expTime = exposure.getInfo().getVisitInfo().getExposureTime()
2495  if not expTime > 0: # handle NaN as well as <= 0
2496  self.log.warn("Non-positive exposure time; skipping rough zero point.")
2497  return
2498 
2499  self.log.info("Setting rough magnitude zero point: %f", 2.5*math.log10(fluxMag0*expTime))
2500  exposure.setPhotoCalib(afwImage.makePhotoCalibFromCalibZeroPoint(fluxMag0*expTime, 0.0))
2501 
2502  def setValidPolygonIntersect(self, ccdExposure, fpPolygon):
2503  """Set the valid polygon as the intersection of fpPolygon and the ccd corners.
2504 
2505  Parameters
2506  ----------
2507  ccdExposure : `lsst.afw.image.Exposure`
2508  Exposure to process.
2509  fpPolygon : `lsst.afw.geom.Polygon`
2510  Polygon in focal plane coordinates.
2511  """
2512  # Get ccd corners in focal plane coordinates
2513  ccd = ccdExposure.getDetector()
2514  fpCorners = ccd.getCorners(FOCAL_PLANE)
2515  ccdPolygon = Polygon(fpCorners)
2516 
2517  # Get intersection of ccd corners with fpPolygon
2518  intersect = ccdPolygon.intersectionSingle(fpPolygon)
2519 
2520  # Transform back to pixel positions and build new polygon
2521  ccdPoints = ccd.transform(intersect, FOCAL_PLANE, PIXELS)
2522  validPolygon = Polygon(ccdPoints)
2523  ccdExposure.getInfo().setValidPolygon(validPolygon)
2524 
2525  @contextmanager
2526  def flatContext(self, exp, flat, dark=None):
2527  """Context manager that applies and removes flats and darks,
2528  if the task is configured to apply them.
2529 
2530  Parameters
2531  ----------
2532  exp : `lsst.afw.image.Exposure`
2533  Exposure to process.
2534  flat : `lsst.afw.image.Exposure`
2535  Flat exposure the same size as ``exp``.
2536  dark : `lsst.afw.image.Exposure`, optional
2537  Dark exposure the same size as ``exp``.
2538 
2539  Yields
2540  ------
2541  exp : `lsst.afw.image.Exposure`
2542  The flat and dark corrected exposure.
2543  """
2544  if self.config.doDark and dark is not None:
2545  self.darkCorrectiondarkCorrection(exp, dark)
2546  if self.config.doFlat:
2547  self.flatCorrectionflatCorrection(exp, flat)
2548  try:
2549  yield exp
2550  finally:
2551  if self.config.doFlat:
2552  self.flatCorrectionflatCorrection(exp, flat, invert=True)
2553  if self.config.doDark and dark is not None:
2554  self.darkCorrectiondarkCorrection(exp, dark, invert=True)
2555 
2556  def debugView(self, exposure, stepname):
2557  """Utility function to examine ISR exposure at different stages.
2558 
2559  Parameters
2560  ----------
2561  exposure : `lsst.afw.image.Exposure`
2562  Exposure to view.
2563  stepname : `str`
2564  State of processing to view.
2565  """
2566  frame = getDebugFrame(self._display, stepname)
2567  if frame:
2568  display = getDisplay(frame)
2569  display.scale('asinh', 'zscale')
2570  display.mtv(exposure)
2571  prompt = "Press Enter to continue [c]... "
2572  while True:
2573  ans = input(prompt).lower()
2574  if ans in ("", "c",):
2575  break
2576 
2577 
2579  """A Detector-like object that supports returning gain and saturation level
2580 
2581  This is used when the input exposure does not have a detector.
2582 
2583  Parameters
2584  ----------
2585  exposure : `lsst.afw.image.Exposure`
2586  Exposure to generate a fake amplifier for.
2587  config : `lsst.ip.isr.isrTaskConfig`
2588  Configuration to apply to the fake amplifier.
2589  """
2590 
2591  def __init__(self, exposure, config):
2592  self._bbox_bbox = exposure.getBBox(afwImage.LOCAL)
2593  self._RawHorizontalOverscanBBox_RawHorizontalOverscanBBox = lsst.geom.Box2I()
2594  self._gain_gain = config.gain
2595  self._readNoise_readNoise = config.readNoise
2596  self._saturation_saturation = config.saturation
2597 
2598  def getBBox(self):
2599  return self._bbox_bbox
2600 
2601  def getRawBBox(self):
2602  return self._bbox_bbox
2603 
2605  return self._RawHorizontalOverscanBBox_RawHorizontalOverscanBBox
2606 
2607  def getGain(self):
2608  return self._gain_gain
2609 
2610  def getReadNoise(self):
2611  return self._readNoise_readNoise
2612 
2613  def getSaturation(self):
2614  return self._saturation_saturation
2615 
2616  def getSuspectLevel(self):
2617  return float("NaN")
2618 
2619 
2620 class RunIsrConfig(pexConfig.Config):
2621  isr = pexConfig.ConfigurableField(target=IsrTask, doc="Instrument signature removal")
2622 
2623 
2624 class RunIsrTask(pipeBase.CmdLineTask):
2625  """Task to wrap the default IsrTask to allow it to be retargeted.
2626 
2627  The standard IsrTask can be called directly from a command line
2628  program, but doing so removes the ability of the task to be
2629  retargeted. As most cameras override some set of the IsrTask
2630  methods, this would remove those data-specific methods in the
2631  output post-ISR images. This wrapping class fixes the issue,
2632  allowing identical post-ISR images to be generated by both the
2633  processCcd and isrTask code.
2634  """
2635  ConfigClass = RunIsrConfig
2636  _DefaultName = "runIsr"
2637 
2638  def __init__(self, *args, **kwargs):
2639  super().__init__(*args, **kwargs)
2640  self.makeSubtask("isr")
2641 
2642  def runDataRef(self, dataRef):
2643  """
2644  Parameters
2645  ----------
2646  dataRef : `lsst.daf.persistence.ButlerDataRef`
2647  data reference of the detector data to be processed
2648 
2649  Returns
2650  -------
2651  result : `pipeBase.Struct`
2652  Result struct with component:
2653 
2654  - exposure : `lsst.afw.image.Exposure`
2655  Post-ISR processed exposure.
2656  """
2657  return self.isr.runDataRef(dataRef)
table::Key< int > type
Definition: Detector.cc:163
Cartesian polygons.
Definition: Polygon.h:59
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:93
An integer coordinate rectangle.
Definition: Box.h:55
def getRawHorizontalOverscanBBox(self)
Definition: isrTask.py:2604
def __init__(self, exposure, config)
Definition: isrTask.py:2591
def __init__(self, *config=None)
Definition: isrTask.py:258
def flatCorrection(self, exposure, flatExposure, invert=False)
Definition: isrTask.py:2178
def maskAndInterpolateNan(self, exposure)
Definition: isrTask.py:2394
def saturationInterpolation(self, exposure)
Definition: isrTask.py:2227
def runDataRef(self, sensorRef)
Definition: isrTask.py:1683
def maskNan(self, exposure)
Definition: isrTask.py:2368
def maskAmplifier(self, ccdExposure, amp, defects)
Definition: isrTask.py:1867
def debugView(self, exposure, stepname)
Definition: isrTask.py:2556
def getIsrExposure(self, dataRef, datasetType, dateObs=None, immediate=True)
Definition: isrTask.py:1734
def saturationDetection(self, exposure, amp)
Definition: isrTask.py:2203
def maskDefect(self, exposure, defectBaseList)
Definition: isrTask.py:2287
def __init__(self, **kwargs)
Definition: isrTask.py:919
def runQuantum(self, butlerQC, inputRefs, outputRefs)
Definition: isrTask.py:929
def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR')
Definition: isrTask.py:2310
def overscanCorrection(self, ccdExposure, amp)
Definition: isrTask.py:1938
def measureBackground(self, exposure, IsrQaConfig=None)
Definition: isrTask.py:2415
def roughZeroPoint(self, exposure)
Definition: isrTask.py:2479
def maskAndInterpolateDefects(self, exposure, defectBaseList)
Definition: isrTask.py:2343
def setValidPolygonIntersect(self, ccdExposure, fpPolygon)
Definition: isrTask.py:2502
def readIsrData(self, dataRef, rawExposure)
Definition: isrTask.py:1022
def ensureExposure(self, inputExp, camera, detectorNum)
Definition: isrTask.py:1782
def updateVariance(self, ampExposure, amp, overscanImage=None)
Definition: isrTask.py:2066
def doLinearize(self, detector)
Definition: isrTask.py:2159
def flatContext(self, exp, flat, dark=None)
Definition: isrTask.py:2526
def convertIntToFloat(self, exposure)
Definition: isrTask.py:1830
def suspectDetection(self, exposure, amp)
Definition: isrTask.py:2252
def run(self, ccdExposure, camera=None, bias=None, linearizer=None, crosstalk=None, crosstalkSources=None, dark=None, flat=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, isGen3=False)
Definition: isrTask.py:1204
def darkCorrection(self, exposure, darkExposure, invert=False)
Definition: isrTask.py:2116
def __init__(self, *args, **kwargs)
Definition: isrTask.py:2638
def runDataRef(self, dataRef)
Definition: isrTask.py:2642
daf::base::PropertyList * list
Definition: fits.cc:913
daf::base::PropertySet * set
Definition: fits.cc:912
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
Definition: PhotoCalib.cc:614
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:462
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.
Definition: MaskedImage.h:1268
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:354
def checkFilter(exposure, filterList, log)
def crosstalkSourceLookup(datasetType, registry, quantumDataId, collections)
Definition: isrTask.py:62
size_t maskNans(afw::image::MaskedImage< PixelT > const &mi, afw::image::MaskPixel maskVal, afw::image::MaskPixel allow=0)
Mask NANs in an image.
Definition: Isr.cc:35
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
def getDebugFrame(debugDisplay, name)
Definition: lsstDebug.py:95