Coverage for python/lsst/analysis/tools/tasks/wholeTractImageAnalysis.py: 48%

131 statements  

« prev     ^ index     » next       coverage.py v7.14.3, created at 2026-06-29 23:06 +0000

1# This file is part of analysis_tools. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22__all__ = ( 

23 "WholeTractImageAnalysisConfig", 

24 "WholeTractImageAnalysisTask", 

25 "WholeTractMaskFractionAnalysisTask", 

26 "MakeBinnedCoaddConfig", 

27 "MakeBinnedCoaddTask", 

28) 

29 

30from collections import defaultdict 

31from collections.abc import Mapping 

32from typing import Any 

33 

34import numpy as np 

35 

36import lsst.pipe.base as pipeBase 

37from lsst.daf.butler import DataCoordinate 

38from lsst.ip.isr.binImageDataTask import binImageData 

39from lsst.pex.config import Field, ListField 

40from lsst.pex.exceptions import InvalidParameterError 

41from lsst.pipe.base import ( 

42 InputQuantizedConnection, 

43 OutputQuantizedConnection, 

44 PipelineTask, 

45 PipelineTaskConfig, 

46 PipelineTaskConnections, 

47 QuantumContext, 

48) 

49from lsst.pipe.base import connectionTypes as ct 

50from lsst.skymap import BaseSkyMap 

51 

52from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask 

53 

54 

55class WholeTractImageAnalysisConnections( 

56 AnalysisBaseConnections, 

57 dimensions=("skymap", "tract", "band"), 

58 defaultTemplates={ 

59 "coaddName": "deep", 

60 }, 

61): 

62 data = ct.Input( 

63 doc="Binned coadd image data to read from the butler.", 

64 name="{coaddName}Coadd_calexp_bin", 

65 storageClass="ExposureF", 

66 deferLoad=True, 

67 dimensions=( 

68 "skymap", 

69 "tract", 

70 "patch", 

71 "band", 

72 ), 

73 multiple=True, 

74 ) 

75 

76 skymap = ct.Input( 

77 doc="The skymap that covers the tract that the data is from.", 

78 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

79 storageClass="SkyMap", 

80 dimensions=("skymap",), 

81 ) 

82 

83 def __init__(self, *, config=None): 

84 """Customize the storageClass for a specific instance. This enables it 

85 to be dynamically set at runtime, allowing the task to work with 

86 different types of image-like data. 

87 

88 Parameters 

89 ---------- 

90 config : `WholeTractImageAnalysisConfig` 

91 A config for `WholeTractImageAnalysisTask`. 

92 """ 

93 super().__init__(config=config) 

94 if config and config.dataStorageClass != self.data.storageClass: 

95 self.data = ct.Input( 

96 name=self.data.name, 

97 doc=self.data.doc, 

98 storageClass=config.dataStorageClass, 

99 dimensions=self.data.dimensions, 

100 deferLoad=self.data.deferLoad, 

101 multiple=self.data.multiple, 

102 ) 

103 

104 if config.outputDimensions: 

105 self.dimensions.clear() 

106 self.dimensions.update(frozenset(sorted(config.outputDimensions))) 

107 

108 

109class WholeTractImageAnalysisConfig( 

110 AnalysisBaseConfig, pipelineConnections=WholeTractImageAnalysisConnections 

111): 

112 dataStorageClass = Field[str]( 

113 default="ExposureF", 

114 doc=( 

115 "Override the storageClass of the input data. " 

116 "Must be of type `Image`, `MaskedImage` or `Exposure`, or one of their subtypes." 

117 ), 

118 ) 

119 outputDimensions = ListField[str]( 

120 default=(), doc=("Override the dimensions of the output data." "Also overrides the task dimensions.") 

121 ) 

122 

123 

124class WholeTractImageAnalysisTask(AnalysisPipelineTask): 

125 

126 ConfigClass = WholeTractImageAnalysisConfig 

127 _DefaultName = "wholeTractImageAnalysis" 

128 

129 def runQuantum( 

130 self, 

131 butlerQC: QuantumContext, 

132 inputRefs: InputQuantizedConnection, 

133 outputRefs: OutputQuantizedConnection, 

134 ) -> None: 

135 inputs = butlerQC.get(inputRefs) 

136 dataId = butlerQC.quantum.dataId 

137 plotInfo = self.parsePlotInfo(inputs, dataId) 

138 

139 try: 

140 inputData = inputs.pop("data") 

141 except KeyError: 

142 raise RuntimeError("'data' is a required input connection, but is not defined.") 

143 

144 keyedData = dict() 

145 if "Exposure" in self.config.dataStorageClass: 

146 inputNames = {"mask"} 

147 inputNames.update(self.collectInputNames()) 

148 for inputName in inputNames: 

149 keyedData[inputName] = dict() 

150 for handle in inputData: 

151 keyedData[inputName][handle.dataId["patch"]] = handle.get(component=inputName) 

152 elif "Image" in self.config.dataStorageClass: 

153 keyedData["image"] = dict() 

154 for handle in inputData: 

155 image = handle.get() 

156 keyedData["image"][handle.dataId["patch"]] = image 

157 else: 

158 raise TypeError("'data' must be of type Image, MaskedImage, Exposure, or one of their subtypes") 

159 

160 outputs = self.run( 

161 data=keyedData, 

162 plotInfo=plotInfo, 

163 tractId=dataId["tract"], 

164 skymap=inputs["skymap"], 

165 bands=dataId["band"], 

166 ) 

167 

168 self.putByBand(butlerQC, outputs, outputRefs) 

169 

170 def parsePlotInfo( 

171 self, inputs: Mapping[str, Any] | None, dataId: DataCoordinate | None, connectionName: str = "data" 

172 ) -> Mapping[str, str]: 

173 """Parse the inputs and dataId to get the information needed to 

174 to add to the figure. The parent class parsePlotInfo cannot be 

175 used becuase it assumes a single input dataset, as opposed to the 

176 multiple datasets used by this analysis task. 

177 

178 Parameters 

179 ---------- 

180 inputs: `dict` 

181 The inputs to the task 

182 dataCoordinate: `lsst.daf.butler.DataCoordinate` 

183 The dataId that the task is being run on. 

184 connectionName: `str`, optional 

185 Name of the input connection to use for determining table name. 

186 

187 Returns 

188 ------- 

189 plotInfo : `dict` 

190 """ 

191 

192 if inputs is None: 

193 tableName = "" 

194 run = "" 

195 else: 

196 tableName = inputs[connectionName][0].ref.datasetType.name 

197 run = inputs[connectionName][0].ref.run 

198 

199 # Initialize the plot info dictionary 

200 plotInfo = {"tableName": tableName, "run": run} 

201 

202 self._populatePlotInfoWithDataId(plotInfo, dataId) 

203 return plotInfo 

204 

205 

206class WholeTractMaskFractionAnalysisConfig( 

207 WholeTractImageAnalysisConfig, pipelineConnections=WholeTractImageAnalysisConnections 

208): 

209 maskPlanes = ListField[str]( 

210 doc="Mask plane names to aggregate fractions of.", 

211 default=[ 

212 "BAD", 

213 "CLIPPED", 

214 "CR", 

215 "DETECTED", 

216 "DETECTED_NEGATIVE", 

217 "EDGE", 

218 "INEXACT_PSF", 

219 "INTRP", 

220 "NO_DATA", 

221 "NOT_DEBLENDED", 

222 "REJECTED", 

223 "SAT", 

224 "STREAK", 

225 "SUSPECT", 

226 "UNMASKEDNAN", 

227 "VIGNETTED", 

228 ], 

229 ) 

230 

231 

232class WholeTractMaskFractionAnalysisTask(AnalysisPipelineTask): 

233 """Computes per-patch mask plane pixel fractions.""" 

234 

235 ConfigClass = WholeTractMaskFractionAnalysisConfig 

236 _DefaultName = "wholeTractMaskFractionAnalysis" 

237 

238 def runQuantum( 

239 self, 

240 butlerQC: QuantumContext, 

241 inputRefs: InputQuantizedConnection, 

242 outputRefs: OutputQuantizedConnection, 

243 ) -> None: 

244 inputs = butlerQC.get(inputRefs) 

245 dataId = butlerQC.quantum.dataId 

246 plotInfo = self.parsePlotInfo(None, dataId) 

247 

248 handles = inputs["data"] 

249 

250 all_planes = set(self.config.maskPlanes) | {"NO_DATA"} 

251 

252 keyedData = self._computeMaskFractions(handles, all_planes) 

253 

254 outputs = self.run(data=keyedData, plotInfo=plotInfo, bands=dataId["band"]) 

255 self.putByBand(butlerQC, outputs, outputRefs) 

256 

257 def _computeMaskFractions(self, handles, all_planes): 

258 """Compute per-patch mask plane fractions across a set of patch 

259 handles. 

260 

261 Parameters 

262 ---------- 

263 handles : `list` 

264 Deferred dataset handles, one per patch. Each must return an 

265 exposure-like object with a ``getMask()`` method. 

266 all_planes : `set` [`str`] 

267 Mask plane names to compute fractions for. Must include 

268 ``"NO_DATA"``. 

269 

270 Returns 

271 ------- 

272 keyedData : `dict` [`str`, `numpy.ndarray`] 

273 Arrays of per-patch fractions, keyed by ``"{plane}_fraction"`` 

274 and ``"{plane}_valid_data_fraction"``, plus 

275 ``"valid_data_pixel_count"``. 

276 """ 

277 fractions: defaultdict[str, list[float]] = defaultdict(lambda: [float("nan")] * len(handles)) 

278 valid_data_pixel_counts: list[int] = [float("nan")] * len(handles) 

279 

280 for i, handle in enumerate(handles): 

281 exp = handle.get() 

282 mask = exp.getMask() 

283 arr = mask.array 

284 

285 no_data_bit = mask.getPlaneBitMask("NO_DATA") 

286 valid_data = (arr & no_data_bit) == 0 

287 

288 total = arr.size 

289 n_valid_data = int(np.sum(valid_data)) 

290 valid_data_pixel_counts[i] = n_valid_data 

291 for plane in all_planes: 

292 try: 

293 bit = mask.getPlaneBitMask(plane) 

294 except InvalidParameterError: 

295 continue 

296 flagged = (arr & bit) != 0 

297 fractions[f"{plane}_fraction"][i] = float(np.sum(flagged)) / total 

298 if plane != "NO_DATA": 

299 value = ( 

300 float(np.sum(flagged & valid_data)) / n_valid_data 

301 if n_valid_data > 0 

302 else float("nan") 

303 ) 

304 fractions[f"{plane}_valid_data_fraction"][i] = value 

305 

306 keyedData = {k: np.array(v) for k, v in fractions.items() if v} 

307 keyedData["valid_data_pixel_count"] = np.array(valid_data_pixel_counts) 

308 return keyedData 

309 

310 

311class MakeBinnedCoaddConnections( 

312 PipelineTaskConnections, 

313 dimensions=("skymap", "tract", "patch", "band"), 

314 defaultTemplates={"coaddName": "deep"}, 

315): 

316 

317 coadd = ct.Input( 

318 doc="Input coadd image data to bin.", 

319 name="{coaddName}Coadd_calexp", 

320 storageClass="ExposureF", 

321 dimensions=("skymap", "tract", "patch", "band"), 

322 deferLoad=True, 

323 ) 

324 skymap = ct.Input( 

325 doc="The skymap that covers the tract that the data is from.", 

326 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

327 storageClass="SkyMap", 

328 dimensions=("skymap",), 

329 ) 

330 binnedCoadd = ct.Output( 

331 doc="Binned coadd image data.", 

332 name="{coaddName}Coadd_calexp_bin", 

333 storageClass="ExposureF", 

334 dimensions=("skymap", "tract", "patch", "band"), 

335 ) 

336 

337 def __init__(self, *, config=None): 

338 """Customize the storageClass for a specific instance. 

339 This enables it to be dynamically set at runtime, allowing 

340 the task to work with different types of image-like data. 

341 

342 Parameters 

343 ---------- 

344 config : `MakeBinnedCoaddConfig` 

345 A config for `MakeBinnedCoaddTask`. 

346 """ 

347 super().__init__(config=config) 

348 if config and config.coaddStorageClass != self.coadd.storageClass: 

349 self.coadd = ct.Input( 

350 name=self.coadd.name, 

351 doc=self.coadd.doc, 

352 storageClass=config.coaddStorageClass, 

353 dimensions=self.coadd.dimensions, 

354 deferLoad=self.coadd.deferLoad, 

355 ) 

356 self.binnedCoadd = ct.Output( 

357 name=self.binnedCoadd.name, 

358 doc=self.binnedCoadd.doc, 

359 storageClass=config.coaddStorageClass, 

360 dimensions=self.binnedCoadd.dimensions, 

361 ) 

362 

363 

364class MakeBinnedCoaddConfig(PipelineTaskConfig, pipelineConnections=MakeBinnedCoaddConnections): 

365 """Config for MakeBinnedCoaddTask""" 

366 

367 doBinInnerBBox = Field[bool]( 

368 doc=( 

369 "Retrieve and bin the coadd image data within the patch Inner Bounding Box, ", 

370 "thereby excluding the regions that overlap neighboring patches.", 

371 ), 

372 default=False, 

373 ) 

374 binFactor = Field[int]( 

375 doc="Binning factor applied to both spatial dimensions.", 

376 default=8, 

377 check=lambda x: x > 1, 

378 ) 

379 coaddStorageClass = Field( 

380 default="ExposureF", 

381 dtype=str, 

382 doc=( 

383 "Override the storageClass of the input and binned coadd image data. " 

384 "Must be of type `Image`, `MaskedImage`, or `Exposure`, or one of their subtypes." 

385 ), 

386 ) 

387 

388 

389class MakeBinnedCoaddTask(PipelineTask): 

390 

391 ConfigClass = MakeBinnedCoaddConfig 

392 _DefaultName = "makeBinnedCoadd" 

393 

394 def runQuantum( 

395 self, 

396 butlerQC: QuantumContext, 

397 inputRefs: InputQuantizedConnection, 

398 outputRefs: OutputQuantizedConnection, 

399 ) -> None: 

400 """Takes coadd image data and bins it by the factor specified in 

401 self.config.binFactor. This task uses the binImageData function 

402 defined in ip_isr, but adds the option to only retrieve and bin the 

403 data contained within the patch's inner bounding box. 

404 

405 Parameters 

406 ---------- 

407 butlerQC : `lsst.pipe.base.QuantumContext` 

408 A butler which is specialized to operate in the context of a 

409 `lsst.daf.butler.Quantum`. 

410 inputRefs : `lsst.pipe.base.InputQuantizedConnection` 

411 Data structure containing named attributes 'coadd' and 'skymap'. 

412 The values of these attributes are the corresponding 

413 `lsst.daf.butler.DatasetRef` objects defined in the corresponding 

414 `PipelineTaskConnections` class. 

415 outputRefs : `lsst.pipe.base.OutputQuantizedConnection` 

416 Datastructure containing named attribute 'binnedCoadd'. 

417 The value of this attribute is the corresponding 

418 `lsst.daf.butler.DatasetRef` object defined in the corresponding 

419 `PipelineTaskConnections` class. 

420 """ 

421 

422 inputs = butlerQC.get(inputRefs) 

423 coaddRef = inputs["coadd"] 

424 

425 if self.config.doBinInnerBBox: 

426 skymap = inputs["skymap"] 

427 tractId = butlerQC.quantum.dataId["tract"] 

428 patchId = butlerQC.quantum.dataId["patch"] 

429 tractInfo = skymap.generateTract(tractId) 

430 bbox = tractInfo.getPatchInfo(patchId).getInnerBBox() 

431 

432 coadd = coaddRef.get(parameters={"bbox": bbox}) 

433 else: 

434 coadd = coaddRef.get() 

435 

436 binnedCoadd = binImageData(coadd, self.config.binFactor) 

437 

438 butlerQC.put(pipeBase.Struct(binnedCoadd=binnedCoadd), outputRefs)