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 badMaskPlanes = [mp for mp in self.config.badMaskPlanes 

110 if mp in afwImage.Mask().getMaskPlaneDict().keys()] 

111 return afwImage.Mask.getPlaneBitMask(badMaskPlanes) 

112 

113 

114def makeSkyInfo(skyMap, tractId, patchId): 

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

116 patchId formats. 

117 

118 Parameters 

119 ---------- 

120 skyMap : `lsst.skyMap.SkyMap` 

121 Sky map. 

122 tractId : `int` 

123 The ID of the tract. 

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

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

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

127 

128 Returns 

129 ------- 

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

131 pipe_base Struct with attributes: 

132 

133 ``skyMap`` 

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

135 ``tractInfo`` 

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

137 ``patchInfo`` 

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

139 ``wcs`` 

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

141 ``bbox`` 

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

143 """ 

144 tractInfo = skyMap[tractId] 

145 

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

147 # patch format is "xIndex,yIndex" 

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

149 else: 

150 patchIndex = patchId 

151 

152 patchInfo = tractInfo.getPatchInfo(patchIndex) 

153 

154 return pipeBase.Struct( 

155 skyMap=skyMap, 

156 tractInfo=tractInfo, 

157 patchInfo=patchInfo, 

158 wcs=tractInfo.getWcs(), 

159 bbox=patchInfo.getOuterBBox(), 

160 ) 

161 

162 

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

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

165 

166 Parameters 

167 ---------- 

168 inputList : `list` 

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

170 inputKeys : `iterable` 

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

172 outputKeys : `iterable` 

173 Iterable of values to be compared with inputKeys. 

174 padWith : `Unknown` 

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

176 

177 Returns 

178 ------- 

179 outputList : `list` 

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

181 so that the length matches length of outputKeys. 

182 """ 

183 outputList = [] 

184 for d in outputKeys: 

185 if d in inputKeys: 

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

187 else: 

188 outputList.append(padWith) 

189 return outputList 

190 

191 

192def reorderRefs(inputRefs, outputSortKeyOrder, dataIdKey): 

193 """Reorder inputRefs per outputSortKeyOrder. 

194 

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

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

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

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

199 

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

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

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

203 outputSortKeyOrder it will be removed. 

204 

205 Parameters 

206 ---------- 

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

208 Input references to be reordered and padded. 

209 outputSortKeyOrder : `iterable` 

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

211 dataIdKey : `str` 

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

213 

214 Returns 

215 ------- 

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

217 Quantized Connection with sorted DatasetRef values sorted if iterable. 

218 """ 

219 for connectionName, refs in inputRefs: 

220 if isinstance(refs, Iterable): 

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

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

223 else: 

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

225 if inputSortKeyOrder != outputSortKeyOrder: 

226 setattr(inputRefs, connectionName, 

227 reorderAndPadList(refs, inputSortKeyOrder, outputSortKeyOrder)) 

228 return inputRefs 

229 

230 

231def subBBoxIter(bbox, subregionSize): 

232 """Iterate over subregions of a bbox. 

233 

234 Parameters 

235 ---------- 

236 bbox : `lsst.geom.Box2I` 

237 Bounding box over which to iterate. 

238 subregionSize : `lsst.geom.Extent2I` 

239 Size of sub-bboxes. 

240 

241 Yields 

242 ------ 

243 subBBox : `lsst.geom.Box2I` 

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

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

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

247 

248 Raises 

249 ------ 

250 RuntimeError 

251 Raised if any of the following occur: 

252 - The given bbox is empty. 

253 - The subregionSize is 0. 

254 """ 

255 if bbox.isEmpty(): 

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

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

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

259 

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

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

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

263 subBBox.clip(bbox) 

264 if subBBox.isEmpty(): 

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

266 "colShift=%s, rowShift=%s" % 

267 (bbox, subregionSize, colShift, rowShift)) 

268 yield subBBox 

269 

270 

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

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

273# without creating any relationships between the two classes. 

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

275# DM-47916. 

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

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

278 

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

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

281 

282 Parameters 

283 ---------- 

284 coaddInputs : `lsst.afw.image.coaddInputs` 

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

286 growBy : `int` 

287 The value to grow the valid polygons by. 

288 

289 Notes 

290 ----- 

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

292 """ 

293 for ccd in coaddInputs.ccds: 

294 polyOrig = ccd.getValidPolygon() 

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

296 validPolyBBox.grow(growBy) 

297 if polyOrig: 

298 validPolygon = polyOrig.intersectionSingle(validPolyBBox) 

299 else: 

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

301 

302 ccd.validPolygon = validPolygon 

303 

304 

305def removeMaskPlanes( 

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

307): 

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

309 

310 Parameters 

311 ---------- 

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

313 The mask to be modified. 

314 mask_planes : `list` 

315 The list of mask planes to be removed. 

316 logger : `logging.Logger`, optional 

317 Logger to log messages. 

318 """ 

319 for maskPlane in mask_planes: 

320 try: 

321 mask &= ~mask.getPlaneBitMask(maskPlane) 

322 except InvalidParameterError: 

323 if logger: 

324 logger.warning( 

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

326 maskPlane, 

327 ) 

328 

329 

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

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

332 

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

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

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

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

337 

338 Parameters 

339 ---------- 

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

341 Statistics control object for coadd. 

342 

343 Returns 

344 ------- 

345 maskMap : `list` of `tuple` of `int` 

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

347 mask planes of the coadd. 

348 """ 

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

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

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

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

353 maskMap = [ 

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

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

356 (clipped, clipped), 

357 ] 

358 return maskMap