Coverage for python/lsst/pipe/tasks/multiBand.py: 21%

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__all__ = ["DetectCoaddSourcesConfig", "DetectCoaddSourcesTask", 

23 "MeasureMergedCoaddSourcesConfig", "MeasureMergedCoaddSourcesTask", 

24 "DEEP_COADD_BACKGROUND_DOCSTRING", 

25 ] 

26 

27import dataclasses 

28 

29import astropy.units 

30import numpy as np 

31 

32from lsst.geom import Extent2I 

33from lsst.pipe.base import ( 

34 AnnotatedPartialOutputsError, 

35 Struct, 

36 PipelineTask, 

37 PipelineTaskConfig, 

38 PipelineTaskConnections 

39) 

40import lsst.pipe.base.connectionTypes as cT 

41from lsst.pex.config import Field, ChoiceField, ConfigurableField 

42from lsst.cell_coadds import MultipleCellCoadd 

43from lsst.images import Mask, get_legacy_deep_coadd_mask_planes 

44from lsst.images.cells import CellCoadd 

45from lsst.images.fields import field_from_legacy_background 

46from lsst.meas.algorithms import ( 

47 DynamicDetectionTask, 

48 ExceedsMaxVarianceScaleError, 

49 InsufficientSourcesError, 

50 PsfGenerationError, 

51 ScaleVarianceTask, 

52 SetPrimaryFlagsTask, 

53 TooManyMaskedPixelsError, 

54 ZeroFootprintError, 

55) 

56from lsst.meas.base import ( 

57 SingleFrameMeasurementTask, 

58 ApplyApCorrTask, 

59 CatalogCalculationTask, 

60 SkyMapIdGeneratorConfig, 

61) 

62from lsst.meas.extensions.scarlet.io import updateCatalogFootprints 

63from lsst.pipe.tasks.propagateSourceFlags import PropagateSourceFlagsTask 

64import lsst.afw.image as afwImage 

65import lsst.afw.math as afwMath 

66import lsst.afw.table as afwTable 

67from lsst.daf.base import PropertyList 

68from lsst.skymap import BaseSkyMap 

69 

70# NOTE: these imports are a convenience so multiband users only have to import this file. 

71from .mergeDetections import MergeDetectionsConfig, MergeDetectionsTask # noqa: F401 

72from .mergeMeasurements import MergeMeasurementsConfig, MergeMeasurementsTask # noqa: F401 

73from .multiBandUtils import CullPeaksConfig # noqa: F401 

74from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiConfig # noqa: F401 

75from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiTask # noqa: F401 

76 

77 

78""" 

79New set types: 

80* deepCoadd_det: detections from what used to be processCoadd (tract, patch, filter) 

81* deepCoadd_mergeDet: merged detections (tract, patch) 

82* deepCoadd_meas: measurements of merged detections (tract, patch, filter) 

83* deepCoadd_ref: reference sources (tract, patch) 

84All of these have associated *_schema catalogs that require no data ID and hold no records. 

85 

86In addition, we have a schema-only dataset, which saves the schema for the PeakRecords in 

87the mergeDet, meas, and ref dataset Footprints: 

88* deepCoadd_peak_schema 

89""" 

90 

91 

92############################################################################################################## 

93 

94# The default for DetectCoaddSourcesConfig.backgroundDescription: 

95DEEP_COADD_BACKGROUND_DOCSTRING = ( 

96 "Background subtracted from the image when generating the Object catalog. " 

97 "This intentionally oversubtracts the background to reduce blending and ensure " 

98 "scattered light is subtracted. " 

99 "Restoring this background does not restore all original backgrounds, " 

100 "as the coadd was built from background-subtracted visit images; in most " 

101 "cases this background term is actually quite small " 

102) 

103 

104 

105class DetectCoaddSourcesConnections(PipelineTaskConnections, 

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

107 defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deep"}): 

108 detectionSchema = cT.InitOutput( 

109 doc="Schema of the detection catalog", 

110 name="{outputCoaddName}Coadd_det_schema", 

111 storageClass="SourceCatalog", 

112 ) 

113 exposure = cT.Input( 

114 doc="Exposure on which detections are to be performed (if useCellCoadd=False). ", 

115 name="{inputCoaddName}Coadd", 

116 storageClass="ExposureF", 

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

118 ) 

119 exposure_cells = cT.Input( 

120 doc="Exposure on which detections are to be performed (if useCellCoadd=True). ", 

121 name="{inputCoaddName}CoaddCell", 

122 storageClass="MultipleCellCoadd", 

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

124 ) 

125 skyMap = cT.Input( 

126 doc="Description of the skymap's tracts and patches.", 

127 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

128 storageClass="SkyMap", 

129 dimensions=("skymap",), 

130 ) 

131 outputBackgrounds = cT.Output( 

132 doc="Output Backgrounds used in detection", 

133 name="{outputCoaddName}Coadd_calexp_background", 

134 storageClass="Background", 

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

136 ) 

137 outputSources = cT.Output( 

138 doc="Detected sources catalog", 

139 name="{outputCoaddName}Coadd_det", 

140 storageClass="SourceCatalog", 

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

142 ) 

143 outputExposure = cT.Output( 

144 doc="Exposure post detection", 

145 name="{outputCoaddName}Coadd_calexp", 

146 storageClass="ExposureF", 

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

148 ) 

149 

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

151 super().__init__(config=config) 

152 assert isinstance(config, DetectCoaddSourcesConfig) 

153 

154 if config.imageType == "future": 

155 if config.useCellCoadds: 

156 self.exposure_cells = dataclasses.replace(self.exposure_cells, storageClass="CellCoadd") 

157 else: 

158 self.exposure = dataclasses.replace(self.exposure, storageClass="CellCoadd") 

159 self.outputExposure = dataclasses.replace(self.outputExposure, storageClass="CellCoadd") 

160 del self.outputBackgrounds 

161 if config.useCellCoadds: 

162 del self.exposure 

163 else: 

164 del self.exposure_cells 

165 if not self.config.forceExactBinning: 

166 del self.skyMap 

167 if self.config.writeOnlyBackgrounds: 

168 del self.outputExposure 

169 del self.outputSources 

170 del self.detectionSchema 

171 

172 

173class DetectCoaddSourcesConfig(PipelineTaskConfig, pipelineConnections=DetectCoaddSourcesConnections): 

174 """Configuration parameters for the DetectCoaddSourcesTask 

175 """ 

176 

177 doScaleVariance = Field(dtype=bool, default=True, doc="Scale variance plane using empirical noise?") 

178 scaleVariance = ConfigurableField(target=ScaleVarianceTask, doc="Variance rescaling") 

179 detection = ConfigurableField(target=DynamicDetectionTask, doc="Source detection") 

180 coaddName = Field(dtype=str, default="deep", doc="Name of coadd") 

181 useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?") 

182 hasFakes = Field( 

183 dtype=bool, 

184 default=False, 

185 doc="Should be set to True if fake sources have been inserted into the input data.", 

186 ) 

187 idGenerator = SkyMapIdGeneratorConfig.make_field() 

188 forceExactBinning = Field( 

189 dtype=bool, 

190 default=False, 

191 doc=( 

192 "Check that the background bin size evenly divides the patch inner region, and " 

193 "crop the outer region to an integer number of bins." 

194 ) 

195 ) 

196 writeOnlyBackgrounds = Field(dtype=bool, default=False, doc="If true, only save the background models.") 

197 writeEmptyBackgrounds = Field( 

198 dtype=bool, 

199 default=True, 

200 doc=( 

201 "If true, save a placeholder background with NaNs in all bins (but the right geometry) when " 

202 "there are no pixels to compute a background from. This can be useful if a later task combines " 

203 "backgrounds from multiple patches as input." 

204 ) 

205 ) 

206 imageType = ChoiceField( 

207 "Which image type to use for the input and output coadd. " 

208 "This option only directly affects connection storage classes and hence 'runQuantum'; the 'run' " 

209 "method behavior is determined by which type is actually passed in.", 

210 allowed={ 

211 "legacy": ( 

212 "Read a lsst.cell_coadds.MultipleCellCoadd (if useCellCoadd) or " 

213 "lsst.afw.image.Exposure (if not useCellCoadd) and write an lsst.afw.image.Exposure." 

214 ), 

215 "future": ( 

216 "Read and write lsst.images.cells.CellCoadd (useCellCoadd just " 

217 "sets which of 'exposure' or 'exposure_cells' will be used for inputs). " 

218 "The 'outputBackground' is deleted, writeEmptyBackgrounds is ignored, and " 

219 "writeOnlyBackgrounds=True is invalid." 

220 ), 

221 }, 

222 dtype=str, 

223 optional=False, 

224 default="legacy", 

225 ) 

226 backgroundName = Field( 

227 "Name of the subtracted background, to be stored with the image when the input and " 

228 "output images are lsst.images.cells.CellCoadd.", 

229 dtype=str, 

230 default="object", 

231 ) 

232 backgroundDescription = Field( 

233 "Description of the subtracted background, to be stored with the image when the input and " 

234 "output images are lsst.images.cells.CellCoadd.", 

235 dtype=str, 

236 default=DEEP_COADD_BACKGROUND_DOCSTRING, 

237 ) 

238 

239 def setDefaults(self): 

240 super().setDefaults() 

241 self.detection.thresholdType = "pixel_stdev" 

242 self.detection.isotropicGrow = True 

243 # Coadds are made from background-subtracted CCDs, so any background subtraction should be very basic 

244 self.detection.reEstimateBackground = False 

245 self.detection.background.useApprox = False 

246 self.detection.background.binSize = 4096 

247 self.detection.background.undersampleStyle = 'REDUCE_INTERP_ORDER' 

248 self.detection.doTempWideBackground = True # Suppress large footprints that overwhelm the deblender 

249 # Include band in packed data IDs that go into object IDs (None -> "as 

250 # many bands as are defined", rather than the default of zero). 

251 self.idGenerator.packer.n_bands = None 

252 

253 def validate(self): 

254 super().validate() 

255 if self.imageType == "future": 

256 if self.doScaleVariance: 

257 raise ValueError("doScaleVariance=True is not compatible with imageType='future'") 

258 if self.forceExactBinning: 

259 raise ValueError("forceExactBinning=True is not compatible with imageType='future'") 

260 if self.writeOnlyBackgrounds: 

261 raise ValueError("writeOnlyBackgrounds=True is not compatible with imageType='future'") 

262 

263 

264class DetectCoaddSourcesTask(PipelineTask): 

265 """Detect sources on a single filter coadd. 

266 

267 Coadding individual visits requires each exposure to be warped. This 

268 introduces covariance in the noise properties across pixels. Before 

269 detection, we correct the coadd variance by scaling the variance plane in 

270 the coadd to match the observed variance. This is an approximate 

271 approach -- strictly, we should propagate the full covariance matrix -- 

272 but it is simple and works well in practice. 

273 

274 After scaling the variance plane, we detect sources and generate footprints 

275 by delegating to the @ref SourceDetectionTask_ "detection" subtask. 

276 

277 DetectCoaddSourcesTask is meant to be run after assembling a coadded image 

278 in a given band. The purpose of the task is to update the background, 

279 detect all sources in a single band and generate a set of parent 

280 footprints. Subsequent tasks in the multi-band processing procedure will 

281 merge sources across bands and, eventually, perform forced photometry. 

282 

283 Parameters 

284 ---------- 

285 schema : `lsst.afw.table.Schema`, optional 

286 Initial schema for the output catalog, modified-in place to include all 

287 fields set by this task. If None, the source minimal schema will be used. 

288 **kwargs 

289 Additional keyword arguments. 

290 """ 

291 

292 _DefaultName = "detectCoaddSources" 

293 ConfigClass = DetectCoaddSourcesConfig 

294 

295 def __init__(self, schema=None, **kwargs): 

296 # N.B. Super is used here to handle the multiple inheritance of PipelineTasks, the init tree 

297 # call structure has been reviewed carefully to be sure super will work as intended. 

298 super().__init__(**kwargs) 

299 if schema is None: 

300 schema = afwTable.SourceTable.makeMinimalSchema() 

301 self.schema = schema 

302 self.makeSubtask("detection", schema=self.schema) 

303 if self.config.doScaleVariance: 

304 self.makeSubtask("scaleVariance") 

305 

306 self.detectionSchema = afwTable.SourceCatalog(self.schema) 

307 

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

309 inputs = butlerQC.get(inputRefs) 

310 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) 

311 

312 if self.config.useCellCoadds: 

313 multiple_cell_coadd = inputs.pop("exposure_cells") 

314 match self.config.imageType: 

315 case "legacy": 

316 exposure = multiple_cell_coadd.stitch().asExposure() 

317 case "future": 

318 exposure = multiple_cell_coadd # conversion deferred to run(). 

319 case _: 

320 raise AssertionError(f"Invalid choice {self.config.imageType!r} for imageType.") 

321 else: 

322 exposure = inputs.pop("exposure") 

323 

324 skyMap = inputs.pop("skyMap", None) 

325 if skyMap is not None: 

326 patchInfo = skyMap[butlerQC.quantum.dataId["tract"]][butlerQC.quantum.dataId["patch"]] 

327 else: 

328 patchInfo = None 

329 

330 assert not inputs, "runQuantum got more inputs than expected." 

331 try: 

332 outputs = self.run( 

333 exposure=exposure, 

334 idFactory=idGenerator.make_table_id_factory(), 

335 expId=idGenerator.catalog_id, 

336 patchInfo=patchInfo, 

337 ) 

338 except ( 

339 TooManyMaskedPixelsError, 

340 ExceedsMaxVarianceScaleError, 

341 InsufficientSourcesError, 

342 PsfGenerationError, 

343 ZeroFootprintError, 

344 ) as e: 

345 if self.config.writeEmptyBackgrounds: 

346 butlerQC.put(self._makeEmptyBackground(exposure, patchInfo), outputRefs.outputBackgrounds) 

347 # Detection failed, so clear any leftover detected mask planes. 

348 if isinstance(exposure, CellCoadd): 

349 # If we passed a CellCoadd in, it won't be modified until 'run' 

350 # is about to exit, and hence it can't have picked up a 

351 # DETECTED_NEGATIVE plane, which can only be temporary here. 

352 # But it might have a preexisting DETECTED plane (i.e. a union 

353 # of the DETECTED plane from the warps) that we'd still want to 

354 # clear. 

355 exposure.mask.clear("DETECTED") 

356 else: 

357 for maskName in ["DETECTED", "DETECTED_NEGATIVE"]: 

358 if maskName in exposure.mask.getMaskPlaneDict().keys(): 

359 detectedMask = exposure.mask.getMaskPlane(maskName) 

360 exposure.mask.clearMaskPlane(detectedMask) 

361 error = AnnotatedPartialOutputsError.annotate( 

362 e, 

363 self, 

364 exposure, 

365 log=self.log, 

366 ) 

367 butlerQC.put(exposure, outputRefs.outputExposure) 

368 raise error from e 

369 

370 butlerQC.put(outputs, outputRefs) 

371 

372 def run(self, exposure, idFactory, expId, patchInfo=None): 

373 """Run detection on an exposure. 

374 

375 First scale the variance plane to match the observed variance 

376 using ``ScaleVarianceTask``. Then invoke the ``SourceDetectionTask_`` "detection" subtask to 

377 detect sources. 

378 

379 Parameters 

380 ---------- 

381 exposure : `lsst.afw.image.Exposure` or `lsst.images.cells.CellCoadd`. 

382 Exposure on which to detect (may be background-subtracted and scaled, 

383 depending on configuration). 

384 idFactory : `lsst.afw.table.IdFactory` 

385 IdFactory to set source identifiers. 

386 expId : `int` 

387 Exposure identifier (integer) for RNG seed. 

388 patchInfo : `lsst.skymap.PatchInfo`, optional 

389 Description of the patch geometry. Only needed if 

390 `~DetectCoaddSourceConfig.forceExactBinning` is `True`. 

391 

392 Returns 

393 ------- 

394 result : `lsst.pipe.base.Struct` 

395 Results as a struct with attributes: 

396 

397 ``outputSources`` 

398 Catalog of detections (`lsst.afw.table.SourceCatalog`). 

399 ``outputBackgrounds`` 

400 List of backgrounds (`list`). 

401 ``outputExposure`` 

402 The background-subtracted coadd image, with its mask plane 

403 updated to include detections. This will have the same type 

404 as ``exposure``. 

405 """ 

406 cell_coadd = None 

407 if isinstance(exposure, CellCoadd): 

408 cell_coadd = exposure 

409 exposure = cell_coadd.to_legacy() 

410 if self.config.forceExactBinning: 

411 if cell_coadd is not None: 

412 raise ValueError("forceExactBinning=True is not compatible with CellCoadd inputs") 

413 exposure = self._cropToExactBinning(exposure, patchInfo) 

414 if self.config.doScaleVariance: 

415 if cell_coadd is not None: 

416 raise ValueError("doScaleVariance=True is not compatible with CellCoadd inputs") 

417 varScale = self.scaleVariance.run(exposure.maskedImage) 

418 exposure.getMetadata().add("VARIANCE_SCALE", varScale) 

419 backgrounds = afwMath.BackgroundList() 

420 table = afwTable.SourceTable.make(self.schema, idFactory) 

421 detections = self.detection.run(table, exposure, expId=expId) 

422 sources = detections.sources 

423 if hasattr(detections, "background") and detections.background: 

424 for bg in detections.background: 

425 backgrounds.append(bg) 

426 if len(backgrounds) == 0: 

427 # Persist a constant background with value of NaN to get around 

428 # inability to persist empty BackgroundList. 

429 emptyBg = self._makeEmptyBackground(exposure, patchInfo) 

430 backgrounds.append(emptyBg) 

431 if cell_coadd is not None: 

432 cell_coadd.image.array[...] = exposure.image.array 

433 cell_coadd.mask = Mask.from_legacy( 

434 exposure.mask, plane_map=get_legacy_deep_coadd_mask_planes() 

435 ).view(sky_projection=cell_coadd.sky_projection) 

436 if backgrounds: 

437 cell_coadd.backgrounds.add( 

438 self.config.backgroundName, 

439 field_from_legacy_background(backgrounds, unit=astropy.units.nJy), 

440 self.config.backgroundDescription, 

441 is_subtracted=True, 

442 ) 

443 exposure = cell_coadd 

444 return Struct(outputSources=sources, outputBackgrounds=backgrounds, outputExposure=exposure) 

445 

446 def _cropToExactBinning(self, exposure, patchInfo): 

447 """Crop a coadd `~lsst.afw.image.Exposure` instance to ensure exact 

448 background binning. 

449 

450 Parameters 

451 ---------- 

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

453 Exposure to crop, assumed to cover the patch outer bounding box. 

454 patchInfo : `lsst.skymap.PatchInfo` 

455 Description of the patch geometry. 

456 

457 Returns 

458 ------- 

459 cropped : `lsst.afw.image.Exposure` 

460 View of ``exposure`` with background bins that evenly divide both 

461 the full cropped image and the patch inner region. The bounding 

462 box is guaranteed to contain the patch inner bounding box and be 

463 contained by the patch outer bounding box. 

464 

465 Raises 

466 ------ 

467 ValueError 

468 Raised if the patch inner region width or height is not a multiple 

469 of the background bin size. 

470 """ 

471 bbox = patchInfo.getInnerBBox() 

472 if bbox.width % self.detection.background.binSizeX: 

473 raise ValueError( 

474 f"Patch inner width {bbox.width} does not evenly " 

475 f"divide bin width {self.detection.background.binSizeX}." 

476 ) 

477 if bbox.height % self.detection.background.binSizeY: 

478 raise ValueError( 

479 f"Patch inner height {bbox.height} does not evenly " 

480 f"divide bin height {self.detection.background.binSizeY}." 

481 ) 

482 outer_bbox = patchInfo.getOuterBBox() 

483 n_bins_grow_x = (bbox.x.begin - outer_bbox.x.begin) // self.detection.background.binSizeX 

484 n_bins_grow_y = (bbox.y.begin - outer_bbox.y.begin) // self.detection.background.binSizeY 

485 bbox.grow( 

486 Extent2I( 

487 n_bins_grow_x*self.detection.background.binSizeX, 

488 n_bins_grow_y*self.detection.background.binSizeY, 

489 ) 

490 ) 

491 assert outer_bbox.contains(bbox) 

492 assert bbox.contains(patchInfo.getInnerBBox()) 

493 assert bbox.width % self.detection.background.binSizeX == 0 

494 assert bbox.height % self.detection.background.binSizeY == 0 

495 return exposure[bbox] 

496 

497 def _makeEmptyBackground(self, exposure, patchInfo=None): 

498 """Construct an empty `lsst.afw.math.BackgroundList` with NaN values. 

499 

500 Parameters 

501 ---------- 

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

503 Exposure that the background should correspond to. 

504 patchInfo : `lsst.skymap.PatchInfo`, optional 

505 Description of the patch geometry. Only needed if 

506 `~DetectCoaddSourceConfig.forceExactBinning` is `True`. 

507 

508 Returns 

509 ------- 

510 background : `lsst.afw.math.BackgroundList` 

511 A background object with a single layer and the same bin geometry 

512 that a background for that exposure would have had if it had enough 

513 usable pixels. This object cannot actually be used for background 

514 subtraction. 

515 """ 

516 # Create a backgroundList with one entry whose "stats image" is NaNs 

517 # and has all pixels set as NO_DATA. 

518 if self.config.forceExactBinning: 

519 exposure = self._cropToExactBinning(exposure, patchInfo).clone() 

520 

521 bgLevel = np.nan 

522 bgStats = afwImage.MaskedImageF(1, 1) 

523 bgStats.set(bgLevel, 0, bgLevel) 

524 bg = afwMath.BackgroundMI(exposure.getBBox(), bgStats) 

525 bgData = (bg, afwMath.Interpolate.LINEAR, afwMath.REDUCE_INTERP_ORDER, 

526 afwMath.ApproximateControl.UNKNOWN, 0, 0, False) 

527 background = afwMath.BackgroundList() 

528 background.append(bgData) 

529 for bg, *_ in background: 

530 stats = bg.getStatsImage() 

531 stats.mask.array[:, :] = stats.mask.getPlaneBitMask("NO_DATA") 

532 stats.variance.array[:, :] = 0.0 

533 return background 

534 

535 

536class MeasureMergedCoaddSourcesConnections( 

537 PipelineTaskConnections, 

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

539 defaultTemplates={ 

540 "inputCoaddName": "deep", 

541 "outputCoaddName": "deep", 

542 }, 

543): 

544 inputSchema = cT.InitInput( 

545 doc="Input schema for measure merged task produced by a deblender or detection task", 

546 name="{inputCoaddName}Coadd_deblendedFlux_schema", 

547 storageClass="SourceCatalog" 

548 ) 

549 outputSchema = cT.InitOutput( 

550 doc="Output schema after all new fields are added by task", 

551 name="{inputCoaddName}Coadd_meas_schema", 

552 storageClass="SourceCatalog" 

553 ) 

554 exposure = cT.Input( 

555 doc="Input non-cell-based coadd image", 

556 name="{inputCoaddName}Coadd_calexp", 

557 storageClass="ExposureF", 

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

559 ) 

560 exposure_cells = cT.Input( 

561 doc="Input cell-based coadd image", 

562 name="{inputCoaddName}CoaddCell", 

563 storageClass="MultipleCellCoadd", 

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

565 ) 

566 background = cT.Input( 

567 doc="Background to subtract from cell-based coadd image", 

568 name="{inputCoaddName}Coadd_calexp_background", 

569 storageClass="Background", 

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

571 ) 

572 skyMap = cT.Input( 

573 doc="SkyMap to use in processing", 

574 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

575 storageClass="SkyMap", 

576 dimensions=("skymap",), 

577 ) 

578 sourceTableHandles = cT.Input( 

579 doc=("Source tables that are derived from the ``CalibrateTask`` sources. " 

580 "These tables contain astrometry and photometry flags, and optionally " 

581 "PSF flags."), 

582 name="sourceTable_visit", 

583 storageClass="ArrowAstropy", 

584 dimensions=("instrument", "visit"), 

585 multiple=True, 

586 deferLoad=True, 

587 ) 

588 finalizedSourceTableHandles = cT.Input( 

589 doc=("Finalized source tables from ``FinalizeCalibrationTask``. These " 

590 "tables contain PSF flags from the finalized PSF estimation."), 

591 name="finalized_src_table", 

592 storageClass="ArrowAstropy", 

593 dimensions=("instrument", "visit"), 

594 multiple=True, 

595 deferLoad=True, 

596 ) 

597 finalVisitSummaryHandles = cT.Input( 

598 doc="Final visit summary table", 

599 name="finalVisitSummary", 

600 storageClass="ExposureCatalog", 

601 dimensions=("instrument", "visit"), 

602 multiple=True, 

603 deferLoad=True, 

604 ) 

605 scarletCatalog = cT.Input( 

606 doc="Catalogs produced by multiband deblending", 

607 name="{inputCoaddName}Coadd_deblendedCatalog", 

608 storageClass="SourceCatalog", 

609 dimensions=("tract", "patch", "skymap"), 

610 ) 

611 scarletModels = cT.Input( 

612 doc="Multiband scarlet models produced by the deblender", 

613 name="{inputCoaddName}Coadd_scarletModelData", 

614 storageClass="LsstScarletModelData", 

615 dimensions=("tract", "patch", "skymap"), 

616 ) 

617 outputSources = cT.Output( 

618 doc="Source catalog containing all the measurement information generated in this task", 

619 name="{outputCoaddName}Coadd_meas", 

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

621 storageClass="SourceCatalog", 

622 ) 

623 

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

625 super().__init__(config=config) 

626 if not config.doPropagateFlags: 

627 del self.sourceTableHandles 

628 del self.finalizedSourceTableHandles 

629 del self.finalVisitSummaryHandles 

630 else: 

631 # Check for types of flags required. 

632 if not config.propagateFlags.source_flags: 

633 del self.sourceTableHandles 

634 if not config.propagateFlags.finalized_source_flags: 

635 del self.finalizedSourceTableHandles 

636 if not config.doAddFootprints: 

637 del self.scarletModels 

638 if self.config.imageType == "future": 

639 self.exposure = dataclasses.replace(self.exposure, storageClass="CellCoadd") 

640 del self.exposure_cells 

641 del self.background 

642 elif self.config.useCellCoadds: 

643 del self.exposure 

644 else: 

645 del self.exposure_cells 

646 del self.background 

647 

648 

649class MeasureMergedCoaddSourcesConfig(PipelineTaskConfig, 

650 pipelineConnections=MeasureMergedCoaddSourcesConnections): 

651 """Configuration parameters for the MeasureMergedCoaddSourcesTask 

652 """ 

653 doAddFootprints = Field(dtype=bool, 

654 default=True, 

655 doc="Whether or not to add footprints to the input catalog from scarlet models. " 

656 "This should be true whenever using the multi-band deblender, " 

657 "otherwise this should be False.") 

658 doConserveFlux = Field(dtype=bool, default=True, 

659 doc="Whether to use the deblender models as templates to re-distribute the flux " 

660 "from the 'exposure' (True), or to perform measurements on the deblender " 

661 "model footprints.") 

662 doStripFootprints = Field(dtype=bool, default=True, 

663 doc="Whether to strip footprints from the output catalog before " 

664 "saving to disk. " 

665 "This is usually done when using scarlet models to save disk space.") 

666 useCellCoadds = Field(dtype=bool, default=False, doc="Whether to use cell coadds?") 

667 measurement = ConfigurableField(target=SingleFrameMeasurementTask, doc="Source measurement") 

668 setPrimaryFlags = ConfigurableField(target=SetPrimaryFlagsTask, doc="Set flags for primary tract/patch") 

669 doPropagateFlags = Field( 

670 dtype=bool, default=True, 

671 doc="Whether to match sources to CCD catalogs to propagate flags (to e.g. identify PSF stars)" 

672 ) 

673 propagateFlags = ConfigurableField(target=PropagateSourceFlagsTask, doc="Propagate source flags to coadd") 

674 coaddName = Field(dtype=str, default="deep", doc="Name of coadd") 

675 psfCache = Field(dtype=int, default=100, doc="Size of psfCache") 

676 checkUnitsParseStrict = Field( 

677 doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'", 

678 dtype=str, 

679 default="raise", 

680 ) 

681 doApCorr = Field( 

682 dtype=bool, 

683 default=True, 

684 doc="Apply aperture corrections" 

685 ) 

686 applyApCorr = ConfigurableField( 

687 target=ApplyApCorrTask, 

688 doc="Subtask to apply aperture corrections" 

689 ) 

690 doRunCatalogCalculation = Field( 

691 dtype=bool, 

692 default=True, 

693 doc='Run catalogCalculation task' 

694 ) 

695 catalogCalculation = ConfigurableField( 

696 target=CatalogCalculationTask, 

697 doc="Subtask to run catalogCalculation plugins on catalog" 

698 ) 

699 

700 hasFakes = Field( 

701 dtype=bool, 

702 default=False, 

703 doc="Should be set to True if fake sources have been inserted into the input data." 

704 ) 

705 idGenerator = SkyMapIdGeneratorConfig.make_field() 

706 imageType = ChoiceField( 

707 "Which image type to expect for the input coadd. " 

708 "This option only directly affects connection storage classes and hence 'runQuantum'; the 'run' " 

709 "method behavior is determined by which type is actually passed in.", 

710 allowed={ 

711 "legacy": ( 

712 "Read a lsst.cell_coadds.MultipleCellCoadd via 'exposure_cells` and restore 'background' " 

713 "(if useCellCoadd) or lsst.afw.image.Exposure via `exposure` (if not useCellCoadd)." 

714 ), 

715 "future": ( 

716 "Read lsst.images.cells.CellCoadd via the 'exposure' connection. useCellCoadd is ignored." 

717 ), 

718 }, 

719 dtype=str, 

720 optional=False, 

721 default="legacy", 

722 ) 

723 

724 def setDefaults(self): 

725 super().setDefaults() 

726 self.measurement.plugins.names |= ['base_InputCount', 

727 'base_Variance', 

728 'base_LocalPhotoCalib', 

729 'base_LocalWcs'] 

730 

731 # TODO: Remove STREAK in DM-44658, streak masking to happen only in 

732 # ip_diffim; if we can propagate the streak mask from diffim, we can 

733 # still set flags with it here. 

734 self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['CLIPPED', 'SENSOR_EDGE', 

735 'INEXACT_PSF'] 

736 self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['CLIPPED', 'SENSOR_EDGE', 

737 'INEXACT_PSF'] 

738 

739 

740class MeasureMergedCoaddSourcesTask(PipelineTask): 

741 """Deblend sources from main catalog in each coadd seperately and measure. 

742 

743 Use peaks and footprints from a master catalog to perform deblending and 

744 measurement in each coadd. 

745 

746 Given a master input catalog of sources (peaks and footprints) or deblender 

747 outputs(including a HeavyFootprint in each band), measure each source on 

748 the coadd. Repeating this procedure with the same master catalog across 

749 multiple coadds will generate a consistent set of child sources. 

750 

751 The deblender retains all peaks and deblends any missing peaks (dropouts in 

752 that band) as PSFs. Source properties are measured and the @c is-primary 

753 flag (indicating sources with no children) is set. Visit flags are 

754 propagated to the coadd sources. 

755 

756 After MeasureMergedCoaddSourcesTask has been run on multiple coadds, we 

757 have a set of per-band catalogs. The next stage in the multi-band 

758 processing procedure will merge these measurements into a suitable catalog 

759 for driving forced photometry. 

760 

761 Parameters 

762 ---------- 

763 schema : ``lsst.afw.table.Schema`, optional 

764 The schema of the merged detection catalog used as input to this one. 

765 peakSchema : ``lsst.afw.table.Schema`, optional 

766 The schema of the PeakRecords in the Footprints in the merged detection catalog. 

767 initInputs : `dict`, optional 

768 Dictionary that can contain a key ``inputSchema`` containing the 

769 input schema. If present will override the value of ``schema``. 

770 **kwargs 

771 Additional keyword arguments. 

772 """ 

773 

774 _DefaultName = "measureCoaddSources" 

775 ConfigClass = MeasureMergedCoaddSourcesConfig 

776 

777 def __init__(self, schema=None, peakSchema=None, initInputs=None, **kwargs): 

778 super().__init__(**kwargs) 

779 if initInputs is not None: 

780 schema = initInputs['inputSchema'].schema 

781 if schema is None: 

782 raise ValueError("Schema must be defined.") 

783 self.schemaMapper = afwTable.SchemaMapper(schema) 

784 self.schemaMapper.addMinimalSchema(schema) 

785 self.schema = self.schemaMapper.getOutputSchema() 

786 self.algMetadata = PropertyList() 

787 self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata) 

788 self.makeSubtask("setPrimaryFlags", schema=self.schema) 

789 if self.config.doPropagateFlags: 

790 self.makeSubtask("propagateFlags", schema=self.schema) 

791 self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict) 

792 if self.config.doApCorr: 

793 self.makeSubtask("applyApCorr", schema=self.schema) 

794 if self.config.doRunCatalogCalculation: 

795 self.makeSubtask("catalogCalculation", schema=self.schema) 

796 

797 self.outputSchema = afwTable.SourceCatalog(self.schema) 

798 

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

800 inputs = butlerQC.get(inputRefs) 

801 

802 if self.config.imageType == "future": 

803 coadd = inputs.pop("exposure") 

804 band = inputRefs.exposure.dataId["band"] 

805 # Instead of going directly from lsst.images.cells.CellCoadd to 

806 # Exposure, it's cleaner for now to go through MultipleCellCoadd 

807 # because the apCorrMap and ccdInputs need special handling - the 

808 # cell-based versions can't be attached to Exposure. Eventually 

809 # we'll rewrite the lower-level code to use the lsst.images 

810 # equivalents natively. 

811 coadd = coadd.to_legacy_cell_coadd() 

812 elif self.config.useCellCoadds: 

813 coadd = inputs.pop("exposure_cells") 

814 band = inputRefs.exposure_cells.dataId["band"] 

815 else: 

816 coadd = inputs.pop("exposure") 

817 band = inputRefs.exposure.dataId["band"] 

818 if isinstance(coadd, MultipleCellCoadd): 

819 stitched_coadd = coadd.stitch() 

820 exposure = stitched_coadd.asExposure() 

821 if self.config.imageType == "legacy": 

822 background = inputs.pop("background") 

823 exposure.image -= background.getImage() 

824 ccdInputs = stitched_coadd.ccds 

825 apCorrMap = stitched_coadd.ap_corr_map 

826 else: 

827 exposure = coadd 

828 # Set psfcache only when we don't have a cell-based coadd. 

829 exposure.getPsf().setCacheCapacity(self.config.psfCache) 

830 

831 ccdInputs = exposure.getInfo().getCoaddInputs().ccds 

832 apCorrMap = exposure.getInfo().getApCorrMap() 

833 

834 # Get unique integer ID for IdFactory and RNG seeds; only the latter 

835 # should really be used as the IDs all come from the input catalog. 

836 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) 

837 

838 # Transform inputCatalog 

839 table = afwTable.SourceTable.make(self.schema, idGenerator.make_table_id_factory()) 

840 sources = afwTable.SourceCatalog(table) 

841 # Load the correct input catalog 

842 inputCatalog = inputs.pop("scarletCatalog") 

843 catalogRef = inputRefs.scarletCatalog 

844 sources.extend(inputCatalog, self.schemaMapper) 

845 del inputCatalog 

846 # Add the HeavyFootprints to the deblended sources 

847 if self.config.doAddFootprints: 

848 modelData = inputs.pop('scarletModels') 

849 if self.config.doConserveFlux: 

850 imageForRedistribution = exposure 

851 else: 

852 imageForRedistribution = None 

853 updateCatalogFootprints( 

854 modelData=modelData, 

855 catalog=sources, 

856 band=band, 

857 imageForRedistribution=imageForRedistribution, 

858 removeScarletData=True, 

859 updateFluxColumns=True, 

860 ) 

861 table = sources.getTable()