Coverage for python / lsst / cp / pipe / cpSky.py: 41%

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 numpy as np 

22 

23import lsst.pex.config as pexConfig 

24import lsst.pipe.base as pipeBase 

25import lsst.pipe.base.connectionTypes as cT 

26from lsst.daf.base import PropertyList 

27from lsst.pipe.tasks.background import ( 

28 FocalPlaneBackground, 

29 FocalPlaneBackgroundConfig, 

30 MaskObjectsTask, 

31 SkyMeasurementTask, 

32) 

33 

34from .cpCombine import CalibCombineTask 

35 

36__all__ = [ 

37 "CpSkyImageTask", 

38 "CpSkyImageConfig", 

39 "CpSkyScaleMeasureTask", 

40 "CpSkyScaleMeasureConfig", 

41 "CpSkySubtractBackgroundTask", 

42 "CpSkySubtractBackgroundConfig", 

43 "CpSkyCombineTask", 

44 "CpSkyCombineConfig", 

45] 

46 

47 

48class CpSkyImageConnections( 

49 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter", "exposure", "detector") 

50): 

51 inputExp = cT.Input( 

52 name="cpSkyIsrExp", 

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

54 storageClass="Exposure", 

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

56 ) 

57 camera = cT.PrerequisiteInput( 

58 name="camera", 

59 doc="Input camera to use for geometry.", 

60 storageClass="Camera", 

61 dimensions=("instrument",), 

62 isCalibration=True, 

63 ) 

64 

65 maskedExp = cT.Output( 

66 name="cpSkyMaskedIsrExp", 

67 doc="Output masked post-ISR exposure.", 

68 storageClass="Exposure", 

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

70 ) 

71 maskedBkg = cT.Output( 

72 name="cpSkyDetectorBackground", 

73 doc="Initial background model from one image.", 

74 storageClass="FocalPlaneBackground", 

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

76 ) 

77 

78 

79class CpSkyImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CpSkyImageConnections): 

80 maskTask = pexConfig.ConfigurableField( 

81 target=MaskObjectsTask, 

82 doc="Object masker to use.", 

83 ) 

84 

85 maskThresh = pexConfig.Field( 

86 dtype=float, 

87 default=3.0, 

88 doc="k-sigma threshold for masking pixels.", 

89 ) 

90 maskList = pexConfig.ListField( 

91 dtype=str, 

92 default=["DETECTED", "BAD", "NO_DATA", "SAT"], 

93 doc="Mask planes to reject.", 

94 ) 

95 

96 largeScaleBackground = pexConfig.ConfigField( 

97 dtype=FocalPlaneBackgroundConfig, 

98 doc="Large-scale background configuration.", 

99 ) 

100 

101 

102class CpSkyImageTask(pipeBase.PipelineTask): 

103 """Mask the detections on the postISRCCD. 

104 

105 This task maps the MaskObjectsTask across all of the initial ISR 

106 processed cpSkyIsr images to create cpSkyMaskedIsr products for 

107 all (exposure, detector) values. 

108 """ 

109 

110 ConfigClass = CpSkyImageConfig 

111 _DefaultName = "CpSkyImage" 

112 

113 def __init__(self, **kwargs): 

114 super().__init__(**kwargs) 

115 self.makeSubtask("maskTask") 

116 

117 def run(self, inputExp, camera): 

118 """Mask the detections on the postISRCCD. 

119 

120 Parameters 

121 ---------- 

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

123 An ISR processed exposure that will have detections 

124 masked. 

125 camera : `lsst.afw.cameraGeom.Camera` 

126 The camera geometry for this exposure. This is needed to 

127 create the background model. 

128 

129 Returns 

130 ------- 

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

132 The results struct containing: 

133 

134 ``maskedExp`` : `lsst.afw.image.Exposure` 

135 The detection-masked version of the ``inputExp``. 

136 ``maskedBkg`` : `lsst.pipe.tasks.background.FocalPlaneBackground` 

137 The partial focal plane background containing only 

138 this exposure/detector's worth of data. 

139 """ 

140 # As constructCalibs.py SkyTask.processSingleBackground() 

141 # Except: check if a detector is fully masked to avoid 

142 # self.maskTask raising. 

143 currentMask = inputExp.getMask() 

144 badMask = currentMask.getPlaneBitMask(self.config.maskList) 

145 if (currentMask.getArray() & badMask).all(): 

146 self.log.warning("All pixels are masked!") 

147 else: 

148 self.maskTask.run(inputExp, self.config.maskList) 

149 

150 # As constructCalibs.py SkyTask.measureBackground() 

151 bgModel = FocalPlaneBackground.fromCamera(self.config.largeScaleBackground, camera) 

152 bgModel.addCcd(inputExp) 

153 

154 return pipeBase.Struct( 

155 maskedExp=inputExp, 

156 maskedBkg=bgModel, 

157 ) 

158 

159 

160class CpSkyScaleMeasureConnections( 

161 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter", "exposure") 

162): 

163 camera = cT.PrerequisiteInput( 

164 name="camera", 

165 doc="Input camera to use for geometry.", 

166 storageClass="Camera", 

167 dimensions=("instrument",), 

168 isCalibration=True, 

169 ) 

170 inputBkgs = cT.Input( 

171 name="cpSkyDetectorBackground", 

172 doc="Initial background model from one exposure/detector", 

173 storageClass="FocalPlaneBackground", 

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

175 multiple=True, 

176 ) 

177 

178 outputBkg = cT.Output( 

179 name="cpSkyExpBackground", 

180 doc="Background model for a full exposure.", 

181 storageClass="FocalPlaneBackground", 

182 dimensions=("instrument", "exposure"), 

183 ) 

184 outputBkgAlternate = cT.Output( 

185 name="cpSkyExpBackgroundAlternate", 

186 doc="Background model for a full exposure from an alternate physical type.", 

187 storageClass="FocalPlaneBackground", 

188 dimensions=("instrument", "exposure"), 

189 ) 

190 outputScale = cT.Output( 

191 name="cpSkyExpScale", 

192 doc="Scale for the full exposure.", 

193 storageClass="PropertyList", 

194 dimensions=("instrument", "exposure"), 

195 ) 

196 

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

198 super().__init__(config=config) 

199 if not config.includeAlternateBackground: 

200 del self.camera 

201 del self.outputBkgAlternate 

202 

203 

204class CpSkyScaleMeasureConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CpSkyScaleMeasureConnections): 

205 includeAlternateBackground = pexConfig.Field( 

206 dtype=bool, 

207 default=False, 

208 doc="Include an alternate focal plane background for detectors of a different physical type?", 

209 ) 

210 

211 

212class CpSkyScaleMeasureTask(pipeBase.PipelineTask): 

213 """Measure per-exposure scale factors and merge focal plane backgrounds. 

214 

215 Merge all the per-detector partial backgrounds to a full focal 

216 plane background for each exposure, and measure the scale factor 

217 from that full background. 

218 """ 

219 

220 ConfigClass = CpSkyScaleMeasureConfig 

221 _DefaultName = "cpSkyScaleMeasure" 

222 

223 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

224 """Ensure that the input and output dimensions are passed along. 

225 

226 Parameters 

227 ---------- 

228 butlerQC : `lsst.daf.butler.QuantumContext` 

229 Butler to operate on. 

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

231 Input data refs to load. 

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

233 Output data refs to persist. 

234 """ 

235 inputs = butlerQC.get(inputRefs) 

236 

237 inputs["inputDims"] = [dict(bkg.dataId.required) for bkg in inputRefs.inputBkgs] 

238 

239 outputs = self.run(**inputs) 

240 butlerQC.put(outputs, outputRefs) 

241 

242 def run(self, inputBkgs, camera=None, inputDims=None): 

243 """Merge focal plane backgrounds and measure the scale factor. 

244 

245 Parameters 

246 ---------- 

247 inputBkgs : `list` [`lsst.pipe.tasks.background.FocalPlaneBackground`] 

248 A list of all of the partial focal plane backgrounds, one 

249 from each detector in this exposure. 

250 camera : `lsst.afw.cameraGeom.Camera`, optional 

251 The camera geometry for this exposure. This is needed to 

252 create the background model. 

253 inputDims : `list` [`dict`], optional 

254 The data IDs for each of the input backgrounds. This is 

255 used to set provenance information on the output background. 

256 

257 Returns 

258 ------- 

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

260 The results struct containing: 

261 

262 ``outputBkg`` : `lsst.pipe.tasks.background.FocalPlaneBackground` 

263 The full merged background for the entire focal plane. 

264 ``outputScale`` : `lsst.daf.base.PropertyList` 

265 A metadata containing the median level of the 

266 background, stored in the key 'scale'. 

267 """ 

268 if self.config.includeAlternateBackground and camera is not None and inputDims is not None: 

269 detectorTypesAll = [camera[ref["detector"]].getPhysicalType() for ref in inputDims] 

270 else: 

271 detectorTypesAll = ["homogeneous" for _ in inputBkgs] 

272 detectorTypes = sorted(set(detectorTypesAll)) 

273 

274 backgrounds = {} 

275 scales = [] 

276 for detectorType in detectorTypes: 

277 inputBkgsSingleType = [bg for bg, dt in zip(inputBkgs, detectorTypesAll) if dt == detectorType] 

278 

279 # Merge per-detector backgrounds into a full focal plane background 

280 background = inputBkgsSingleType[0] 

281 for bg in inputBkgsSingleType[1:]: 

282 background.merge(bg) 

283 backgrounds[detectorType] = background 

284 

285 backgroundPixels = background.getStatsImage().getArray() 

286 self.log.info( 

287 "Background model%s min/max: %f %f. Scale: %f", 

288 "" if detectorType == "homogeneous" else f" ({detectorType})", 

289 np.min(backgroundPixels), 

290 np.max(backgroundPixels), 

291 np.median(backgroundPixels), 

292 ) 

293 

294 # TODO: Ultimately, we should modify FocalPlaneBackground to 

295 # store metadata directly and set up a storage class which allows 

296 # us to persist multiple FocalPlaneBackground types per exposure. 

297 # For now, we store the scale and type data in a PropertyList. 

298 scales.append(np.median(background.getStatsImage().getArray())) 

299 

300 scaleMD = PropertyList() 

301 scaleMD.set("scale", float(scales[0])) 

302 if len(detectorTypes) > 1: 

303 scaleMD.set("detectorType", detectorTypes[0]) # TODO: this is not technically needed 

304 scaleMD.set("scaleAlternate", float(scales[1])) 

305 scaleMD.set("detectorTypeAlternate", detectorTypes[1]) 

306 

307 return pipeBase.Struct( 

308 outputBkg=backgrounds[detectorTypes[0]], 

309 outputBkgAlternate=backgrounds.get(scaleMD.get("detectorTypeAlternate")), 

310 outputScale=scaleMD, 

311 ) 

312 

313 

314class CpSkySubtractBackgroundConnections( 

315 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter", "exposure", "detector") 

316): 

317 inputExp = cT.Input( 

318 name="cpSkyMaskedIsrExp", 

319 doc="Masked post-ISR image.", 

320 storageClass="Exposure", 

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

322 ) 

323 inputBkg = cT.Input( 

324 name="cpSkyExpBackground", 

325 doc="Background model for the full exposure.", 

326 storageClass="FocalPlaneBackground", 

327 dimensions=("instrument", "exposure"), 

328 ) 

329 inputBkgAlternate = cT.Input( 

330 name="cpSkyExpBackgroundAlternate", 

331 doc="Background model for a full exposure from an alternate physical type.", 

332 storageClass="FocalPlaneBackground", 

333 dimensions=("instrument", "exposure"), 

334 ) 

335 inputScale = cT.Input( 

336 name="cpSkyExpScale", 

337 doc="Scale for the full exposure.", 

338 storageClass="PropertyList", 

339 dimensions=("instrument", "exposure"), 

340 ) 

341 

342 outputBkg = cT.Output( 

343 name="cpSkyExpResidualBackground", 

344 doc="Normalized, static background.", 

345 storageClass="Background", 

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

347 ) 

348 

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

350 super().__init__(config=config) 

351 if not config.includeAlternateBackground: 

352 del self.inputBkgAlternate 

353 

354 

355class CpSkySubtractBackgroundConfig( 

356 pipeBase.PipelineTaskConfig, pipelineConnections=CpSkySubtractBackgroundConnections 

357): 

358 sky = pexConfig.ConfigurableField( 

359 target=SkyMeasurementTask, 

360 doc="Sky measurement", 

361 ) 

362 includeAlternateBackground = pexConfig.Field( 

363 dtype=bool, 

364 default=False, 

365 doc="Use an alternate focal plane background for detectors of a different physical type?", 

366 ) 

367 

368 

369class CpSkySubtractBackgroundTask(pipeBase.PipelineTask): 

370 """Subtract per-exposure background from individual detector masked images. 

371 

372 The cpSkyMaskedIsr images constructed by CpSkyImageTask have the 

373 scaled background constructed by CpSkyScaleMeasureTask subtracted, 

374 and new background models are constructed for the remaining 

375 signal. 

376 

377 The output was called `icExpBackground` in gen2, but the product 

378 created here has definition clashes that prevent that from being 

379 reused. 

380 """ 

381 

382 ConfigClass = CpSkySubtractBackgroundConfig 

383 _DefaultName = "cpSkySubtractBkg" 

384 

385 def __init__(self, **kwargs): 

386 super().__init__(**kwargs) 

387 self.makeSubtask("sky") 

388 

389 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

390 """Ensure that the input and output dimensions are passed along. 

391 

392 Parameters 

393 ---------- 

394 butlerQC : `lsst.daf.butler.QuantumContext` 

395 Butler to operate on. 

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

397 Input data refs to load. 

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

399 Output data refs to persist. 

400 """ 

401 inputs = butlerQC.get(inputRefs) 

402 

403 if self.config.includeAlternateBackground: 

404 detectorType = inputs["inputExp"].getDetector().getPhysicalType() 

405 scaleMD = inputs["inputScale"] 

406 # Swap in the alternate background and scale if appropriate 

407 if detectorType == scaleMD.get("detectorTypeAlternate"): 

408 inputs["inputBkg"] = inputs["inputBkgAlternate"] 

409 scaleMDAlternate = PropertyList() 

410 scaleMDAlternate.set("scale", scaleMD.get("scaleAlternate")) 

411 inputs["inputScale"] = scaleMDAlternate 

412 _ = inputs.pop("inputBkgAlternate", None) 

413 

414 outputs = self.run(**inputs) 

415 butlerQC.put(outputs, outputRefs) 

416 

417 def run(self, inputExp, inputBkg, inputScale): 

418 """Subtract per-exposure background from individual detector masked 

419 images. 

420 

421 Parameters 

422 ---------- 

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

424 The ISR processed, detection masked image. 

425 inputBkg : `lsst.pipe.tasks.background.FocalPlaneBackground. 

426 Full focal plane background for this exposure. 

427 inputScale : `lsst.daf.base.PropertyList` 

428 Metadata containing the scale factor. 

429 

430 Returns 

431 ------- 

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

433 The results struct containing: 

434 

435 ``outputBkg`` 

436 Remnant sky background with the full-exposure 

437 component removed. (`lsst.afw.math.BackgroundList`) 

438 """ 

439 # As constructCalibs.py SkyTask.processSingle() 

440 image = inputExp.getMaskedImage() 

441 detector = inputExp.getDetector() 

442 bbox = image.getBBox() 

443 

444 scale = inputScale.get("scale") 

445 background = inputBkg.toCcdBackground(detector, bbox) 

446 image -= background.getImage() 

447 image /= scale 

448 

449 newBackground = self.sky.measureBackground(image) 

450 return pipeBase.Struct(outputBkg=newBackground) 

451 

452 

453class CpSkyCombineConnections( 

454 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter", "detector") 

455): 

456 inputBkgs = cT.Input( 

457 name="cpSkyExpResidualBackground", 

458 doc="Normalized, static background.", 

459 storageClass="Background", 

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

461 multiple=True, 

462 ) 

463 inputExpHandles = cT.Input( 

464 name="cpSkyMaskedIsrExp", 

465 doc="Masked post-ISR image.", 

466 storageClass="Exposure", 

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

468 multiple=True, 

469 deferLoad=True, 

470 ) 

471 

472 outputCalib = cT.Output( 

473 name="sky", 

474 doc="Averaged static background.", 

475 storageClass="ExposureF", 

476 dimensions=("instrument", "detector", "physical_filter"), 

477 isCalibration=True, 

478 ) 

479 

480 

481class CpSkyCombineConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CpSkyCombineConnections): 

482 sky = pexConfig.ConfigurableField( 

483 target=SkyMeasurementTask, 

484 doc="Sky measurement", 

485 ) 

486 

487 

488class CpSkyCombineTask(pipeBase.PipelineTask): 

489 """Merge per-exposure measurements into a detector level calibration. 

490 

491 Each of the per-detector results from all input exposures are 

492 averaged to produce the final SKY calibration. 

493 

494 As before, this is written to a skyCalib instead of a SKY to avoid 

495 definition classes in gen3. 

496 """ 

497 

498 ConfigClass = CpSkyCombineConfig 

499 _DefaultName = "cpSkyCombine" 

500 

501 def __init__(self, **kwargs): 

502 super().__init__(**kwargs) 

503 self.makeSubtask("sky") 

504 

505 def run(self, inputBkgs, inputExpHandles): 

506 """Merge per-exposure measurements into a detector level calibration. 

507 

508 Parameters 

509 ---------- 

510 inputBkgs : `list` [`lsst.afw.math.BackgroundList`] 

511 Remnant backgrounds from each exposure. 

512 inputHandles : `list` [`lsst.daf.butler.DeferredDatasetHandles`] 

513 The Butler handles to the ISR processed, detection masked images. 

514 

515 Returns 

516 ------- 

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

518 The results struct containing: 

519 

520 `outputCalib` : `lsst.afw.image.Exposure` 

521 The final sky calibration product. 

522 """ 

523 skyCalib = self.sky.averageBackgrounds(inputBkgs) 

524 skyCalib.setDetector(inputExpHandles[0].get(component="detector")) 

525 skyCalib.setFilter(inputExpHandles[0].get(component="filter")) 

526 

527 CalibCombineTask().combineHeaders(inputExpHandles, skyCalib, calibType="SKY") 

528 

529 return pipeBase.Struct( 

530 outputCalib=skyCalib, 

531 )