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._backgrounds.clear() 

438 cell_coadd.backgrounds.add( 

439 self.config.backgroundName, 

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

441 self.config.backgroundDescription, 

442 is_subtracted=True, 

443 ) 

444 exposure = cell_coadd 

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

446 

447 def _cropToExactBinning(self, exposure, patchInfo): 

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

449 background binning. 

450 

451 Parameters 

452 ---------- 

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

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

455 patchInfo : `lsst.skymap.PatchInfo` 

456 Description of the patch geometry. 

457 

458 Returns 

459 ------- 

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

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

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

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

464 contained by the patch outer bounding box. 

465 

466 Raises 

467 ------ 

468 ValueError 

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

470 of the background bin size. 

471 """ 

472 bbox = patchInfo.getInnerBBox() 

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

474 raise ValueError( 

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

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

477 ) 

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

479 raise ValueError( 

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

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

482 ) 

483 outer_bbox = patchInfo.getOuterBBox() 

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

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

486 bbox.grow( 

487 Extent2I( 

488 n_bins_grow_x*self.detection.background.binSizeX, 

489 n_bins_grow_y*self.detection.background.binSizeY, 

490 ) 

491 ) 

492 assert outer_bbox.contains(bbox) 

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

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

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

496 return exposure[bbox] 

497 

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

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

500 

501 Parameters 

502 ---------- 

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

504 Exposure that the background should correspond to. 

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

506 Description of the patch geometry. Only needed if 

507 `~DetectCoaddSourceConfig.forceExactBinning` is `True`. 

508 

509 Returns 

510 ------- 

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

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

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

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

515 subtraction. 

516 """ 

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

518 # and has all pixels set as NO_DATA. 

519 if self.config.forceExactBinning: 

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

521 

522 bgLevel = np.nan 

523 bgStats = afwImage.MaskedImageF(1, 1) 

524 bgStats.set(bgLevel, 0, bgLevel) 

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

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

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

528 background = afwMath.BackgroundList() 

529 background.append(bgData) 

530 for bg, *_ in background: 

531 stats = bg.getStatsImage() 

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

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

534 return background 

535 

536 

537class MeasureMergedCoaddSourcesConnections( 

538 PipelineTaskConnections, 

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

540 defaultTemplates={ 

541 "inputCoaddName": "deep", 

542 "outputCoaddName": "deep", 

543 }, 

544): 

545 inputSchema = cT.InitInput( 

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

547 name="{inputCoaddName}Coadd_deblendedFlux_schema", 

548 storageClass="SourceCatalog" 

549 ) 

550 outputSchema = cT.InitOutput( 

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

552 name="{inputCoaddName}Coadd_meas_schema", 

553 storageClass="SourceCatalog" 

554 ) 

555 exposure = cT.Input( 

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

557 name="{inputCoaddName}Coadd_calexp", 

558 storageClass="ExposureF", 

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

560 ) 

561 exposure_cells = cT.Input( 

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

563 name="{inputCoaddName}CoaddCell", 

564 storageClass="MultipleCellCoadd", 

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

566 ) 

567 background = cT.Input( 

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

569 name="{inputCoaddName}Coadd_calexp_background", 

570 storageClass="Background", 

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

572 ) 

573 skyMap = cT.Input( 

574 doc="SkyMap to use in processing", 

575 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

576 storageClass="SkyMap", 

577 dimensions=("skymap",), 

578 ) 

579 sourceTableHandles = cT.Input( 

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

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

582 "PSF flags."), 

583 name="sourceTable_visit", 

584 storageClass="ArrowAstropy", 

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

586 multiple=True, 

587 deferLoad=True, 

588 ) 

589 finalizedSourceTableHandles = cT.Input( 

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

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

592 name="finalized_src_table", 

593 storageClass="ArrowAstropy", 

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

595 multiple=True, 

596 deferLoad=True, 

597 ) 

598 finalVisitSummaryHandles = cT.Input( 

599 doc="Final visit summary table", 

600 name="finalVisitSummary", 

601 storageClass="ExposureCatalog", 

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

603 multiple=True, 

604 deferLoad=True, 

605 ) 

606 scarletCatalog = cT.Input( 

607 doc="Catalogs produced by multiband deblending", 

608 name="{inputCoaddName}Coadd_deblendedCatalog", 

609 storageClass="SourceCatalog", 

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

611 ) 

612 scarletModels = cT.Input( 

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

614 name="{inputCoaddName}Coadd_scarletModelData", 

615 storageClass="LsstScarletModelData", 

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

617 ) 

618 outputSources = cT.Output( 

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

620 name="{outputCoaddName}Coadd_meas", 

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

622 storageClass="SourceCatalog", 

623 ) 

624 

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

626 super().__init__(config=config) 

627 if not config.doPropagateFlags: 

628 del self.sourceTableHandles 

629 del self.finalizedSourceTableHandles 

630 del self.finalVisitSummaryHandles 

631 else: 

632 # Check for types of flags required. 

633 if not config.propagateFlags.source_flags: 

634 del self.sourceTableHandles 

635 if not config.propagateFlags.finalized_source_flags: 

636 del self.finalizedSourceTableHandles 

637 if not config.doAddFootprints: 

638 del self.scarletModels 

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

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

641 del self.exposure_cells 

642 del self.background 

643 elif self.config.useCellCoadds: 

644 del self.exposure 

645 else: 

646 del self.exposure_cells 

647 del self.background 

648 

649 

650class MeasureMergedCoaddSourcesConfig(PipelineTaskConfig, 

651 pipelineConnections=MeasureMergedCoaddSourcesConnections): 

652 """Configuration parameters for the MeasureMergedCoaddSourcesTask 

653 """ 

654 doAddFootprints = Field(dtype=bool, 

655 default=True, 

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

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

658 "otherwise this should be False.") 

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

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

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

662 "model footprints.") 

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

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

665 "saving to disk. " 

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

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

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

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

670 doPropagateFlags = Field( 

671 dtype=bool, default=True, 

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

673 ) 

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

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

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

677 checkUnitsParseStrict = Field( 

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

679 dtype=str, 

680 default="raise", 

681 ) 

682 doApCorr = Field( 

683 dtype=bool, 

684 default=True, 

685 doc="Apply aperture corrections" 

686 ) 

687 applyApCorr = ConfigurableField( 

688 target=ApplyApCorrTask, 

689 doc="Subtask to apply aperture corrections" 

690 ) 

691 doRunCatalogCalculation = Field( 

692 dtype=bool, 

693 default=True, 

694 doc='Run catalogCalculation task' 

695 ) 

696 catalogCalculation = ConfigurableField( 

697 target=CatalogCalculationTask, 

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

699 ) 

700 

701 hasFakes = Field( 

702 dtype=bool, 

703 default=False, 

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

705 ) 

706 idGenerator = SkyMapIdGeneratorConfig.make_field() 

707 imageType = ChoiceField( 

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

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

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

711 allowed={ 

712 "legacy": ( 

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

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

715 ), 

716 "future": ( 

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

718 ), 

719 }, 

720 dtype=str, 

721 optional=False, 

722 default="legacy", 

723 ) 

724 

725 def setDefaults(self): 

726 super().setDefaults() 

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

728 'base_Variance', 

729 'base_LocalPhotoCalib', 

730 'base_LocalWcs'] 

731 

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

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

734 # still set flags with it here. 

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

736 'INEXACT_PSF'] 

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

738 'INEXACT_PSF'] 

739 

740 

741class MeasureMergedCoaddSourcesTask(PipelineTask): 

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

743 

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

745 measurement in each coadd. 

746 

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

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

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

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

751 

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

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

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

755 propagated to the coadd sources. 

756 

757 After MeasureMergedCoaddSourcesTask has been run on multiple coadds, we 

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

759 processing procedure will merge these measurements into a suitable catalog 

760 for driving forced photometry. 

761 

762 Parameters 

763 ---------- 

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

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

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

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

768 initInputs : `dict`, optional 

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

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

771 **kwargs 

772 Additional keyword arguments. 

773 """ 

774 

775 _DefaultName = "measureCoaddSources" 

776 ConfigClass = MeasureMergedCoaddSourcesConfig 

777 

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

779 super().__init__(**kwargs) 

780 if initInputs is not None: 

781 schema = initInputs['inputSchema'].schema 

782 if schema is None: 

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

784 self.schemaMapper = afwTable.SchemaMapper(schema) 

785 self.schemaMapper.addMinimalSchema(schema) 

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

787 self.algMetadata = PropertyList() 

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

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

790 if self.config.doPropagateFlags: 

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

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

793 if self.config.doApCorr: 

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

795 if self.config.doRunCatalogCalculation: 

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

797 

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

799 

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

801 inputs = butlerQC.get(inputRefs) 

802 

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

804 coadd = inputs.pop("exposure") 

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

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

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

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

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

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

811 # equivalents natively. 

812 coadd = coadd.to_legacy_cell_coadd() 

813 elif self.config.useCellCoadds: 

814 coadd = inputs.pop("exposure_cells") 

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

816 else: 

817 coadd = inputs.pop("exposure") 

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

819 if isinstance(coadd, MultipleCellCoadd): 

820 stitched_coadd = coadd.stitch() 

821 exposure = stitched_coadd.asExposure() 

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

823 background = inputs.pop("background") 

824 exposure.image -= background.getImage() 

825 ccdInputs = stitched_coadd.ccds 

826 apCorrMap = stitched_coadd.ap_corr_map 

827 else: 

828 exposure = coadd 

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

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

831 

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

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

834 

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

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

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

838 

839 # Transform inputCatalog 

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

841 sources = afwTable.SourceCatalog(table) 

842 # Load the correct input catalog 

843 inputCatalog = inputs.pop("scarletCatalog") 

844 catalogRef = inputRefs.scarletCatalog 

845 sources.extend(inputCatalog, self.schemaMapper) 

846 del inputCatalog 

847 # Add the HeavyFootprints to the deblended sources 

848 if self.config.doAddFootprints: 

849 modelData = inputs.pop('scarletModels') 

850 if self.config.doConserveFlux: 

851 imageForRedistribution = exposure 

852 else: 

853 imageForRedistribution = None 

854 updateCatalogFootprints( 

855 modelData=modelData, 

856 catalog=sources, 

857 band=band, 

858 imageForRedistribution=imageForRedistribution, 

859 removeScarletData=True, 

860 updateFluxColumns=True,