Coverage for python/lsst/pipe/tasks/coaddBase.py: 25%

1# This file is part of pipe_tasks. 

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 

23from __future__ import annotations 

24 

25__all__ = ["CoaddBaseTask", "makeSkyInfo"] 

26 

27from collections.abc import Iterable 

28from typing import TYPE_CHECKING 

29 

30import lsst.afw.image as afwImage 

31import lsst.geom as geom 

32import lsst.pex.config as pexConfig 

33import lsst.pipe.base as pipeBase 

34from lsst.afw.geom import Polygon 

35from lsst.pex.exceptions import InvalidParameterError 

36 

37from .coaddInputRecorder import CoaddInputRecorderTask 

38from .selectImages import PsfWcsSelectImagesTask 

39 

40if TYPE_CHECKING: 

41 from logging import Logger 

42 

43 from lsst.afw.math import StatisticsControl 

44 

45 

46class CoaddBaseConfig(pexConfig.Config): 

47 """Configuration parameters for CoaddBaseTask 

48 

49 Configuration parameters shared between MakeCoaddTempExp and AssembleCoadd 

50 """ 

51 

52 coaddName = pexConfig.Field( 

53 doc="Coadd name: typically one of deep or goodSeeing.", 

54 dtype=str, 

55 default="deep", 

56 ) 

57 select = pexConfig.ConfigurableField( 

58 doc="Image selection subtask.", 

59 target=PsfWcsSelectImagesTask, 

60 ) 

61 badMaskPlanes = pexConfig.ListField( 

62 dtype=str, 

63 doc="Mask planes that, if set, the associated pixel should not be included in the coaddTempExp.", 

64 default=("NO_DATA",), 

65 ) 

66 inputRecorder = pexConfig.ConfigurableField( 

67 doc="Subtask that helps fill CoaddInputs catalogs added to the final Exposure", 

68 target=CoaddInputRecorderTask, 

69 ) 

70 # TODO[DM-49400]: remove this field (it already does nothing). 

71 includeCalibVar = pexConfig.Field( 

72 dtype=bool, 

73 doc="Add photometric calibration variance to warp variance plane.", 

74 default=False, 

75 deprecated="Deprecated and ignored. Will be removed after v29.", 

76 ) 

77 

78 

79class CoaddBaseTask(pipeBase.PipelineTask): 

80 """Base class for coaddition. 

81 

82 Subclasses must specify _DefaultName 

83 """ 

84 

85 ConfigClass = CoaddBaseConfig 

86 

87 def __init__(self, **kwargs): 

88 super().__init__(**kwargs) 

89 self.makeSubtask("select") 

90 self.makeSubtask("inputRecorder") 

91 

92 def getTempExpDatasetName(self, warpType="direct"): 

93 """Return warp name for given warpType and task config 

94 

95 Parameters 

96 ---------- 

97 warpType : `str` 

98 Either 'direct' or 'psfMatched'. 

99 

100 Returns 

101 ------- 

102 WarpDatasetName : `str` 

103 """ 

104 return self.config.coaddName + "Coadd_" + warpType + "Warp" 

105 

106 def getBadPixelMask(self): 

107 """Convenience method to provide the bitmask from the mask plane names 

108 """ 

109 return afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes) 

110 

111 

112def makeSkyInfo(skyMap, tractId, patchId): 

113 """Constructs SkyInfo used by coaddition tasks for multiple 

114 patchId formats. 

115 

116 Parameters 

117 ---------- 

118 skyMap : `lsst.skyMap.SkyMap` 

119 Sky map. 

120 tractId : `int` 

121 The ID of the tract. 

122 patchId : `str` or `int` or `tuple` of `int` 

123 Either Gen2-style comma delimited string (e.g. '4,5'), 

124 tuple of integers (e.g (4, 5), Gen3-style integer. 

125 

126 Returns 

127 ------- 

128 makeSkyInfo : `lsst.pipe.base.Struct` 

129 pipe_base Struct with attributes: 

130 

131 ``skyMap`` 

132 Sky map (`lsst.skyMap.SkyMap`). 

133 ``tractInfo`` 

134 Information for chosen tract of sky map (`lsst.skyMap.TractInfo`). 

135 ``patchInfo`` 

136 Information about chosen patch of tract (`lsst.skyMap.PatchInfo`). 

137 ``wcs`` 

138 WCS of tract (`lsst.afw.image.SkyWcs`). 

139 ``bbox`` 

140 Outer bbox of patch, as an geom Box2I (`lsst.afw.geom.Box2I`). 

141 """ 

142 tractInfo = skyMap[tractId] 

143 

144 if isinstance(patchId, str) and ',' in patchId: 

145 # patch format is "xIndex,yIndex" 

146 patchIndex = tuple(int(i) for i in patchId.split(",")) 

147 else: 

148 patchIndex = patchId 

149 

150 patchInfo = tractInfo.getPatchInfo(patchIndex) 

151 

152 return pipeBase.Struct( 

153 skyMap=skyMap, 

154 tractInfo=tractInfo, 

155 patchInfo=patchInfo, 

156 wcs=tractInfo.getWcs(), 

157 bbox=patchInfo.getOuterBBox(), 

158 ) 

159 

160 

161def reorderAndPadList(inputList, inputKeys, outputKeys, padWith=None): 

162 """Match the order of one list to another, padding if necessary 

163 

164 Parameters 

165 ---------- 

166 inputList : `list` 

167 List to be reordered and padded. Elements can be any type. 

168 inputKeys : `iterable` 

169 Iterable of values to be compared with outputKeys. Length must match `inputList`. 

170 outputKeys : `iterable` 

171 Iterable of values to be compared with inputKeys. 

172 padWith : `Unknown` 

173 Any value to be inserted where inputKey not in outputKeys. 

174 

175 Returns 

176 ------- 

177 outputList : `list` 

178 Copy of inputList reordered per outputKeys and padded with `padWith` 

179 so that the length matches length of outputKeys. 

180 """ 

181 outputList = [] 

182 for d in outputKeys: 

183 if d in inputKeys: 

184 outputList.append(inputList[inputKeys.index(d)]) 

185 else: 

186 outputList.append(padWith) 

187 return outputList 

188 

189 

190def reorderRefs(inputRefs, outputSortKeyOrder, dataIdKey): 

191 """Reorder inputRefs per outputSortKeyOrder. 

192 

193 Any inputRefs which are lists will be resorted per specified key e.g., 

194 'detector.' Only iterables will be reordered, and values can be of type 

195 `lsst.pipe.base.connections.DeferredDatasetRef` or 

196 `lsst.daf.butler.core.datasets.ref.DatasetRef`. 

197 

198 Returned lists of refs have the same length as the outputSortKeyOrder. 

199 If an outputSortKey not in the inputRef, then it will be padded with None. 

200 If an inputRef contains an inputSortKey that is not in the 

201 outputSortKeyOrder it will be removed. 

202 

203 Parameters 

204 ---------- 

205 inputRefs : `lsst.pipe.base.connections.QuantizedConnection` 

206 Input references to be reordered and padded. 

207 outputSortKeyOrder : `iterable` 

208 Iterable of values to be compared with inputRef's dataId[dataIdKey]. 

209 dataIdKey : `str` 

210 The data ID key in the dataRefs to compare with the outputSortKeyOrder. 

211 

212 Returns 

213 ------- 

214 inputRefs : `lsst.pipe.base.connections.QuantizedConnection` 

215 Quantized Connection with sorted DatasetRef values sorted if iterable. 

216 """ 

217 for connectionName, refs in inputRefs: 

218 if isinstance(refs, Iterable): 

219 if hasattr(refs[0], "dataId"): 

220 inputSortKeyOrder = [ref.dataId[dataIdKey] for ref in refs] 

221 else: 

222 inputSortKeyOrder = [handle.datasetRef.dataId[dataIdKey] for handle in refs] 

223 if inputSortKeyOrder != outputSortKeyOrder: 

224 setattr(inputRefs, connectionName, 

225 reorderAndPadList(refs, inputSortKeyOrder, outputSortKeyOrder)) 

226 return inputRefs 

227 

228 

229def subBBoxIter(bbox, subregionSize): 

230 """Iterate over subregions of a bbox. 

231 

232 Parameters 

233 ---------- 

234 bbox : `lsst.geom.Box2I` 

235 Bounding box over which to iterate. 

236 subregionSize : `lsst.geom.Extent2I` 

237 Size of sub-bboxes. 

238 

239 Yields 

240 ------ 

241 subBBox : `lsst.geom.Box2I` 

242 Next sub-bounding box of size ``subregionSize`` or smaller; each ``subBBox`` 

243 is contained within ``bbox``, so it may be smaller than ``subregionSize`` at 

244 the edges of ``bbox``, but it will never be empty. 

245 

246 Raises 

247 ------ 

248 RuntimeError 

249 Raised if any of the following occur: 

250 - The given bbox is empty. 

251 - The subregionSize is 0. 

252 """ 

253 if bbox.isEmpty(): 

254 raise RuntimeError("bbox %s is empty" % (bbox,)) 

255 if subregionSize[0] < 1 or subregionSize[1] < 1: 

256 raise RuntimeError("subregionSize %s must be nonzero" % (subregionSize,)) 

257 

258 for rowShift in range(0, bbox.getHeight(), subregionSize[1]): 

259 for colShift in range(0, bbox.getWidth(), subregionSize[0]): 

260 subBBox = geom.Box2I(bbox.getMin() + geom.Extent2I(colShift, rowShift), subregionSize) 

261 subBBox.clip(bbox) 

262 if subBBox.isEmpty(): 

263 raise RuntimeError("Bug: empty bbox! bbox=%s, subregionSize=%s, " 

264 "colShift=%s, rowShift=%s" % 

265 (bbox, subregionSize, colShift, rowShift)) 

266 yield subBBox 

267 

268 

269# Note that this is implemented as a free-floating function to enable reuse in 

270# lsst.pipe.tasks.makeWarp and in lsst.drp.tasks.make_psf_matched_warp 

271# without creating any relationships between the two classes. 

272# This may be converted to a method after makeWarp.py is removed altogether in 

273# DM-47916. 

274def growValidPolygons(coaddInputs, growBy: int) -> None: 

275 """Grow coaddInputs' ccds' ValidPolygons in place. 

276 

277 Either modify each ccd's validPolygon in place, or if CoaddInputs 

278 does not have a validPolygon, create one from its bbox. 

279 

280 Parameters 

281 ---------- 

282 coaddInputs : `lsst.afw.image.coaddInputs` 

283 CoaddInputs object containing the ccds to grow the valid polygons of. 

284 growBy : `int` 

285 The value to grow the valid polygons by. 

286 

287 Notes 

288 ----- 

289 Negative values for ``growBy`` can shrink the polygons. 

290 """ 

291 for ccd in coaddInputs.ccds: 

292 polyOrig = ccd.getValidPolygon() 

293 validPolyBBox = polyOrig.getBBox() if polyOrig else ccd.getBBox() 

294 validPolyBBox.grow(growBy) 

295 if polyOrig: 

296 validPolygon = polyOrig.intersectionSingle(validPolyBBox) 

297 else: 

298 validPolygon = Polygon(geom.Box2D(validPolyBBox)) 

299 

300 ccd.validPolygon = validPolygon 

301 

302 

303def removeMaskPlanes( 

304 mask: afwImage.Mask, mask_planes: Iterable, logger: Logger | None = None 

305): 

306 """Unset the mask of an image for mask planes specified in the config. 

307 

308 Parameters 

309 ---------- 

310 mask : `lsst.afw.image.Mask` 

311 The mask to be modified. 

312 mask_planes : `list` 

313 The list of mask planes to be removed. 

314 logger : `logging.Logger`, optional 

315 Logger to log messages. 

316 """ 

317 for maskPlane in mask_planes: 

318 try: 

319 mask &= ~mask.getPlaneBitMask(maskPlane) 

320 except InvalidParameterError: 

321 if logger: 

322 logger.warning( 

323 "Unable to remove mask plane %s: no mask plane with that name was found.", 

324 maskPlane, 

325 ) 

326 

327 

328def setRejectedMaskMapping(statsCtrl: StatisticsControl) -> list[tuple[int, int]]: 

329 """Map certain mask planes of the warps to new planes for the coadd. 

330 

331 If a pixel is rejected due to a mask value other than EDGE, NO_DATA, 

332 or CLIPPED, set it to REJECTED on the coadd. 

333 If a pixel is rejected due to EDGE, set the coadd pixel to SENSOR_EDGE. 

334 If a pixel is rejected due to CLIPPED, set the coadd pixel to CLIPPED. 

335 

336 Parameters 

337 ---------- 

338 statsCtrl : `lsst.afw.math.StatisticsControl` 

339 Statistics control object for coadd. 

340 

341 Returns 

342 ------- 

343 maskMap : `list` of `tuple` of `int` 

344 A list of mappings of mask planes of the warped exposures to 

345 mask planes of the coadd. 

346 """ 

347 edge = 2 ** afwImage.Mask.addMaskPlane("EDGE") 

348 noData = 2 ** afwImage.Mask.addMaskPlane("NO_DATA") 

349 clipped = 2 ** afwImage.Mask.addMaskPlane("CLIPPED") 

350 toReject = statsCtrl.getAndMask() & (~noData) & (~edge) & (~clipped) 

351 maskMap = [ 

352 (toReject, 2 ** afwImage.Mask.addMaskPlane("REJECTED")), 

353 (edge, 2 ** afwImage.Mask.addMaskPlane("SENSOR_EDGE")), 

354 (clipped, clipped), 

355 ] 

356 return maskMap