Coverage for python / lsst / cp / pipe / cpDark.py: 37%

1# This file is part of cp_pipe. 

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 <http://www.gnu.org/licenses/>. 

21import math 

22import numpy as np 

23 

24import lsst.afw.geom as afwGeom 

25import lsst.meas.algorithms as measAlg 

26import lsst.pex.config as pexConfig 

27import lsst.pipe.base as pipeBase 

28import lsst.pipe.base.connectionTypes as cT 

29 

30from lsst.ip.isr import gainContext, interpolateFromMask 

31from lsst.pex.exceptions import LengthError 

32from lsst.pipe.tasks.repair import RepairTask 

33 

34 

35__all__ = ["CpDarkTask", "CpDarkTaskConfig"] 

36 

37 

38class CpDarkConnections(pipeBase.PipelineTaskConnections, 

39 dimensions=("instrument", "exposure", "detector")): 

40 inputExp = cT.Input( 

41 name="cpDarkISR", 

42 doc="Input pre-processed exposures to combine.", 

43 storageClass="Exposure", 

44 dimensions=("instrument", "exposure", "detector"), 

45 ) 

46 

47 outputExp = cT.Output( 

48 name="cpDarkProc", 

49 doc="Output combined proposed calibration.", 

50 storageClass="Exposure", 

51 dimensions=("instrument", "exposure", "detector"), 

52 ) 

53 

54 

55class CpDarkTaskConfig(pipeBase.PipelineTaskConfig, 

56 pipelineConnections=CpDarkConnections): 

57 psfFwhm = pexConfig.Field( 

58 dtype=float, 

59 default=3.0, 

60 doc="Repair PSF FWHM (pixels).", 

61 ) 

62 psfSize = pexConfig.Field( 

63 dtype=int, 

64 default=21, 

65 doc="Repair PSF size (pixels).", 

66 ) 

67 crGrow = pexConfig.Field( 

68 dtype=int, 

69 default=2, 

70 doc="Grow radius for CR (pixels).", 

71 ) 

72 repair = pexConfig.ConfigurableField( 

73 target=RepairTask, 

74 doc="Repair task to use.", 

75 ) 

76 maskListToInterpolate = pexConfig.ListField( 

77 dtype=str, 

78 doc="List of mask planes that should be interpolated.", 

79 default=["SAT", "BAD"], 

80 ) 

81 useLegacyInterp = pexConfig.Field( 

82 dtype=bool, 

83 doc="Use the legacy interpolation algorithm. If False use Gaussian Process.", 

84 default=True, 

85 ) 

86 

87 

88class CpDarkTask(pipeBase.PipelineTask): 

89 """Combine pre-processed dark frames into a proposed master calibration. 

90 """ 

91 

92 ConfigClass = CpDarkTaskConfig 

93 _DefaultName = "cpDark" 

94 

95 def __init__(self, **kwargs): 

96 super().__init__(**kwargs) 

97 self.makeSubtask("repair") 

98 

99 def run(self, inputExp): 

100 """Preprocess input exposures prior to DARK combination. 

101 

102 This task detects and repairs cosmic rays strikes. 

103 

104 Parameters 

105 ---------- 

106 inputExp : `lsst.afw.image.Exposure` 

107 Pre-processed dark frame data to combine. 

108 

109 Returns 

110 ------- 

111 results : `lsst.pipe.base.Struct` 

112 The results struct containing: 

113 

114 ``outputExp`` 

115 CR rejected, ISR processed Dark Frame 

116 (`lsst.afw.image.Exposure`). 

117 """ 

118 psf = measAlg.SingleGaussianPsf(self.config.psfSize, 

119 self.config.psfSize, 

120 self.config.psfFwhm/(2*math.sqrt(2*math.log(2)))) 

121 inputExp.setPsf(psf) 

122 

123 # Get the gain used to set the variance plane from the 

124 # exposure, if possible. Otherwise, use the cameraGeom value. 

125 gains = self._get_gains(inputExp) 

126 

127 # Is this gainContext still required? TODO DM-48754: 

128 # Investigate cosmic ray rejection during dark construction 

129 with gainContext(inputExp, inputExp.getVariance(), apply=True, gains=gains): 

130 # Scale the variance to match the image plane. A similar 

131 # scaling happens during flat-field correction for science 

132 # images. 

133 self.log.debug("Median image and variance: %f %f", 

134 np.median(inputExp.image.array), np.median(inputExp.variance.array)) 

135 crImage = inputExp.clone() 

136 # Interpolate the crImage, so the CR code can ignore 

137 # defects (which will now be interpolated). 

138 interpolateFromMask( 

139 maskedImage=crImage.getMaskedImage(), 

140 fwhm=self.config.psfFwhm, 

141 growSaturatedFootprints=0, 

142 maskNameList=list(self.config.maskListToInterpolate), 

143 useLegacyInterp=self.config.useLegacyInterp, 

144 ) 

145 

146 try: 

147 self.repair.run(crImage, keepCRs=False) 

148 except LengthError: 

149 self.log.warning("CR rejection failed!") 

150 

151 # Copy results to input frame. 

152 crBit = crImage.mask.getPlaneBitMask("CR") 

153 crPixels = (crImage.mask.array & crBit) > 0 

154 inputExp.mask.array[crPixels] |= crBit 

155 self.log.info("Number of CR pixels: %d", 

156 np.count_nonzero(crPixels)) 

157 

158 if self.config.crGrow > 0: 

159 crMask = inputExp.getMaskedImage().getMask().getPlaneBitMask("CR") 

160 spans = afwGeom.SpanSet.fromMask(inputExp.mask, crMask) 

161 spans = spans.dilated(self.config.crGrow) 

162 spans = spans.clippedTo(inputExp.getBBox()) 

163 spans.setMask(inputExp.mask, crMask) 

164 

165 # Clear the defect (BAD) mask plane, if it exists. 

166 planes = inputExp.mask.getMaskPlaneDict() 

167 if "BAD" in planes: 

168 inputExp.mask.clearMaskPlane(planes["BAD"]) 

169 

170 return pipeBase.Struct( 

171 outputExp=inputExp, 

172 ) 

173 

174 @staticmethod 

175 def _get_gains(exposure): 

176 """Get the per-amplifier gains used for this exposure. 

177 

178 Parameters 

179 ---------- 

180 exposure : `lsst.afw.image.Exposure` 

181 The exposure to find gains for. 

182 

183 Returns 

184 ------- 

185 gains : `dict` [`str` `float`] 

186 Dictionary of gain values, keyed by amplifier name. 

187 """ 

188 det = exposure.getDetector() 

189 metadata = exposure.getMetadata() 

190 gains = {} 

191 for amp in det: 

192 ampName = amp.getName() 

193 # The GAIN key may be the new LSST ISR GAIN or the old 

194 # LSST GAIN. 

195 if (key1 := f"LSST ISR GAIN {ampName}") in metadata: 

196 gains[ampName] = metadata[key1] 

197 elif (key2 := f"LSST GAIN {ampName}") in metadata: 

198 gains[ampName] = metadata[key2] 

199 else: 

200 gains[ampName] = amp.getGain() 

201 return gains