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

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 ] 

25 

26import numpy as np 

27 

28from lsst.geom import Extent2I 

29from lsst.pipe.base import ( 

30 AnnotatedPartialOutputsError, 

31 Struct, 

32 PipelineTask, 

33 PipelineTaskConfig, 

34 PipelineTaskConnections 

35) 

36import lsst.pipe.base.connectionTypes as cT 

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

38from lsst.meas.algorithms import DynamicDetectionTask, ReferenceObjectLoader, ScaleVarianceTask, \ 

39 SetPrimaryFlagsTask 

40from lsst.meas.algorithms.subtractBackground import TooManyMaskedPixelsError 

41from lsst.meas.base import ( 

42 SingleFrameMeasurementTask, 

43 ApplyApCorrTask, 

44 CatalogCalculationTask, 

45 SkyMapIdGeneratorConfig, 

46) 

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

48from lsst.meas.astrom import DirectMatchTask, denormalizeMatches 

49from lsst.pipe.tasks.propagateSourceFlags import PropagateSourceFlagsTask 

50import lsst.afw.table as afwTable 

51import lsst.afw.math as afwMath 

52from lsst.daf.base import PropertyList 

53from lsst.skymap import BaseSkyMap 

54 

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

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

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

58from .multiBandUtils import CullPeaksConfig # noqa: F401 

59from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiConfig # noqa: F401 

60from .deblendCoaddSourcesPipeline import DeblendCoaddSourcesMultiTask # noqa: F401 

61 

62 

63""" 

64New set types: 

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

66* deepCoadd_mergeDet: merged detections (tract, patch) 

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

68* deepCoadd_ref: reference sources (tract, patch) 

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

70 

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

72the mergeDet, meas, and ref dataset Footprints: 

73* deepCoadd_peak_schema 

74""" 

75 

76 

77############################################################################################################## 

78class DetectCoaddSourcesConnections(PipelineTaskConnections, 

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

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

81 detectionSchema = cT.InitOutput( 

82 doc="Schema of the detection catalog", 

83 name="{outputCoaddName}Coadd_det_schema", 

84 storageClass="SourceCatalog", 

85 ) 

86 exposure = cT.Input( 

87 doc="Exposure on which detections are to be performed. ", 

88 name="{inputCoaddName}Coadd", 

89 storageClass="ExposureF", 

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

91 ) 

92 exposure_cells = cT.Input( 

93 doc="Exposure on which detections are to be performed. ", 

94 name="{inputCoaddName}CoaddCell", 

95 storageClass="MultipleCellCoadd", 

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

97 ) 

98 skyMap = cT.Input( 

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

100 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

101 storageClass="SkyMap", 

102 dimensions=("skymap",), 

103 ) 

104 outputBackgrounds = cT.Output( 

105 doc="Output Backgrounds used in detection", 

106 name="{outputCoaddName}Coadd_calexp_background", 

107 storageClass="Background", 

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

109 ) 

110 outputSources = cT.Output( 

111 doc="Detected sources catalog", 

112 name="{outputCoaddName}Coadd_det", 

113 storageClass="SourceCatalog", 

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

115 ) 

116 outputExposure = cT.Output( 

117 doc="Exposure post detection", 

118 name="{outputCoaddName}Coadd_calexp", 

119 storageClass="ExposureF", 

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

121 ) 

122 

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

124 super().__init__(config=config) 

125 assert isinstance(config, DetectCoaddSourcesConfig) 

126 

127 if config.useCellCoadds: 

128 del self.exposure 

129 else: 

130 del self.exposure_cells 

131 

132 if not self.config.forceExactBinning: 

133 del self.skyMap 

134 if self.config.writeOnlyBackgrounds: 

135 del self.outputExposure 

136 del self.outputSources 

137 del self.detectionSchema 

138 

139 

140class DetectCoaddSourcesConfig(PipelineTaskConfig, pipelineConnections=DetectCoaddSourcesConnections): 

141 """Configuration parameters for the DetectCoaddSourcesTask 

142 """ 

143 

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

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

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

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

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

149 hasFakes = Field( 

150 dtype=bool, 

151 default=False, 

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

153 ) 

154 idGenerator = SkyMapIdGeneratorConfig.make_field() 

155 forceExactBinning = Field( 

156 dtype=bool, 

157 default=False, 

158 doc=( 

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

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

161 ) 

162 ) 

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

164 writeEmptyBackgrounds = Field( 

165 dtype=bool, 

166 default=False, 

167 doc=( 

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

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

170 "backgrounds from multiple patches as input." 

171 ) 

172 ) 

173 

174 def setDefaults(self): 

175 super().setDefaults() 

176 self.detection.thresholdType = "pixel_stdev" 

177 self.detection.isotropicGrow = True 

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

179 self.detection.reEstimateBackground = False 

180 self.detection.background.useApprox = False 

181 self.detection.background.binSize = 4096 

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

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

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

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

186 self.idGenerator.packer.n_bands = None 

187 

188 

189class DetectCoaddSourcesTask(PipelineTask): 

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

191 

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

193 introduces covariance in the noise properties across pixels. Before 

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

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

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

197 but it is simple and works well in practice. 

198 

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

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

201 

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

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

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

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

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

207 

208 Parameters 

209 ---------- 

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

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

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

213 **kwargs 

214 Additional keyword arguments. 

215 """ 

216 

217 _DefaultName = "detectCoaddSources" 

218 ConfigClass = DetectCoaddSourcesConfig 

219 

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

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

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

223 super().__init__(**kwargs) 

224 if schema is None: 

225 schema = afwTable.SourceTable.makeMinimalSchema() 

226 self.schema = schema 

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

228 if self.config.doScaleVariance: 

229 self.makeSubtask("scaleVariance") 

230 

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

232 

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

234 inputs = butlerQC.get(inputRefs) 

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

236 

237 if self.config.useCellCoadds: 

238 multiple_cell_coadd = inputs.pop("exposure_cells") 

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

240 else: 

241 exposure = inputs.pop("exposure") 

242 

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

244 if skyMap is not None: 

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

246 else: 

247 patchInfo = None 

248 

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

250 try: 

251 outputs = self.run( 

252 exposure=exposure, 

253 idFactory=idGenerator.make_table_id_factory(), 

254 expId=idGenerator.catalog_id, 

255 patchInfo=patchInfo, 

256 ) 

257 except TooManyMaskedPixelsError as e: 

258 if self.config.writeEmptyBackgrounds: 

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

260 error = AnnotatedPartialOutputsError.annotate( 

261 e, 

262 self, 

263 exposure, 

264 log=self.log, 

265 ) 

266 raise error from e 

267 

268 butlerQC.put(outputs, outputRefs) 

269 

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

271 """Run detection on an exposure. 

272 

273 First scale the variance plane to match the observed variance 

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

275 detect sources. 

276 

277 Parameters 

278 ---------- 

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

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

281 depending on configuration). 

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

283 IdFactory to set source identifiers. 

284 expId : `int` 

285 Exposure identifier (integer) for RNG seed. 

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

287 Description of the patch geometry. Only needed if 

288 `~DetectCoaddSourceConfig.forceExactBinning` is `True`. 

289 

290 Returns 

291 ------- 

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

293 Results as a struct with attributes: 

294 

295 ``sources`` 

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

297 ``backgrounds`` 

298 List of backgrounds (`list`). 

299 """ 

300 if self.config.forceExactBinning: 

301 exposure = self._cropToExactBinning(exposure, patchInfo) 

302 if self.config.doScaleVariance: 

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

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

305 backgrounds = afwMath.BackgroundList() 

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

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

308 sources = detections.sources 

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

310 for bg in detections.background: 

311 backgrounds.append(bg) 

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

313 

314 def _cropToExactBinning(self, exposure, patchInfo): 

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

316 background binning. 

317 

318 Parameters 

319 ---------- 

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

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

322 patchInfo : `lsst.skymap.PatchInfo` 

323 Description of the patch geometry. 

324 

325 Returns 

326 ------- 

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

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

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

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

331 contained by the patch outer bounding box. 

332 

333 Raises 

334 ------ 

335 ValueError 

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

337 of the background bin size. 

338 """ 

339 bbox = patchInfo.getInnerBBox() 

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

341 raise ValueError( 

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

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

344 ) 

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

346 raise ValueError( 

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

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

349 ) 

350 outer_bbox = patchInfo.getOuterBBox() 

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

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

353 bbox.grow( 

354 Extent2I( 

355 n_bins_grow_x*self.detection.background.binSizeX, 

356 n_bins_grow_y*self.detection.background.binSizeY, 

357 ) 

358 ) 

359 assert outer_bbox.contains(bbox) 

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

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

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

363 return exposure[bbox] 

364 

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

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

367 

368 Parameters 

369 ---------- 

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

371 Exposure that the background should correspond to. 

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

373 Description of the patch geometry. Only needed if 

374 `~DetectCoaddSourceConfig.forceExactBinning` is `True`. 

375 

376 Returns 

377 ------- 

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

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

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

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

382 subtraction. 

383 """ 

384 # BackgroundList objects are a huge pain to construct without actually 

385 # measuring the background, so that's what we do - on an image with 

386 # 0s everywhere and no masks. And then we replace the background 

387 # "stats image" with NaNs and NO_DATA. 

388 if self.config.forceExactBinning: 

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

390 exposure.image.array[:, :] = 0.0 

391 exposure.mask[:, :] = 0 

392 background = self.detection.background.run(exposure).background 

393 for bg, *_ in background: 

394 stats = bg.getStatsImage() 

395 stats.image.array[:, :] = np.nan 

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

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

398 return background 

399 

400 

401class MeasureMergedCoaddSourcesConnections( 

402 PipelineTaskConnections, 

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

404 defaultTemplates={ 

405 "inputCoaddName": "deep", 

406 "outputCoaddName": "deep", 

407 "deblendedCatalog": "deblendedFlux", 

408 }, 

409 deprecatedTemplates={ 

410 # TODO[DM-47797]: remove this deprecated connection template. 

411 "deblendedCatalog": "Support for old deblender outputs will be removed after v29." 

412 }, 

413): 

414 inputSchema = cT.InitInput( 

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

416 name="{inputCoaddName}Coadd_deblendedFlux_schema", 

417 storageClass="SourceCatalog" 

418 ) 

419 outputSchema = cT.InitOutput( 

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

421 name="{inputCoaddName}Coadd_meas_schema", 

422 storageClass="SourceCatalog" 

423 ) 

424 # TODO[DM-47797]: remove this deprecated connection. 

425 refCat = cT.PrerequisiteInput( 

426 doc="Reference catalog used to match measured sources against known sources", 

427 name="ref_cat", 

428 storageClass="SimpleCatalog", 

429 dimensions=("skypix",), 

430 deferLoad=True, 

431 multiple=True, 

432 deprecated="Reference matching in measureCoaddSources will be removed after v29.", 

433 ) 

434 exposure = cT.Input( 

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

436 name="{inputCoaddName}Coadd_calexp", 

437 storageClass="ExposureF", 

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

439 ) 

440 exposure_cells = cT.Input( 

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

442 name="{inputCoaddName}CoaddCell", 

443 storageClass="MultipleCellCoadd", 

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

445 ) 

446 background = cT.Input( 

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

448 name="{inputCoaddName}Coadd_calexp_background", 

449 storageClass="Background", 

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

451 ) 

452 skyMap = cT.Input( 

453 doc="SkyMap to use in processing", 

454 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

455 storageClass="SkyMap", 

456 dimensions=("skymap",), 

457 ) 

458 # TODO[DM-47424]: remove this deprecated connection. 

459 visitCatalogs = cT.Input( 

460 doc="Deprecated and unused.", 

461 name="src", 

462 dimensions=("instrument", "visit", "detector"), 

463 storageClass="SourceCatalog", 

464 multiple=True, 

465 deprecated="Deprecated and unused. Will be removed after v29.", 

466 ) 

467 sourceTableHandles = cT.Input( 

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

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

470 "PSF flags."), 

471 name="sourceTable_visit", 

472 storageClass="ArrowAstropy", 

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

474 multiple=True, 

475 deferLoad=True, 

476 ) 

477 finalizedSourceTableHandles = cT.Input( 

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

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

480 name="finalized_src_table", 

481 storageClass="ArrowAstropy", 

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

483 multiple=True, 

484 deferLoad=True, 

485 ) 

486 finalVisitSummaryHandles = cT.Input( 

487 doc="Final visit summary table", 

488 name="finalVisitSummary", 

489 storageClass="ExposureCatalog", 

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

491 multiple=True, 

492 deferLoad=True, 

493 ) 

494 # TODO[DM-47797]: remove this deprecated connection. 

495 inputCatalog = cT.Input( 

496 doc=("Name of the input catalog to use." 

497 "If the single band deblender was used this should be 'deblendedFlux." 

498 "If the multi-band deblender was used this should be 'deblendedModel, " 

499 "or deblendedFlux if the multiband deblender was configured to output " 

500 "deblended flux catalogs. If no deblending was performed this should " 

501 "be 'mergeDet'"), 

502 name="{inputCoaddName}Coadd_{deblendedCatalog}", 

503 storageClass="SourceCatalog", 

504 deprecated="Support for old deblender outputs will be removed after v29.", 

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

506 ) 

507 scarletCatalog = cT.Input( 

508 doc="Catalogs produced by multiband deblending", 

509 name="{inputCoaddName}Coadd_deblendedCatalog", 

510 storageClass="SourceCatalog", 

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

512 ) 

513 scarletModels = cT.Input( 

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

515 name="{inputCoaddName}Coadd_scarletModelData", 

516 storageClass="ScarletModelData", 

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

518 ) 

519 outputSources = cT.Output( 

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

521 name="{outputCoaddName}Coadd_meas", 

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

523 storageClass="SourceCatalog", 

524 ) 

525 # TODO[DM-47797]: remove this deprecated connection. 

526 matchResult = cT.Output( 

527 doc="Match catalog produced by configured matcher, optional on doMatchSources", 

528 name="{outputCoaddName}Coadd_measMatch", 

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

530 storageClass="Catalog", 

531 deprecated="Reference matching in measureCoaddSources will be removed after v29.", 

532 ) 

533 # TODO[DM-47797]: remove this deprecated connection. 

534 denormMatches = cT.Output( 

535 doc="Denormalized Match catalog produced by configured matcher, optional on " 

536 "doWriteMatchesDenormalized", 

537 name="{outputCoaddName}Coadd_measMatchFull", 

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

539 storageClass="Catalog", 

540 deprecated="Reference matching in measureCoaddSources will be removed after v29.", 

541 ) 

542 

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

544 super().__init__(config=config) 

545 del self.visitCatalogs 

546 if not config.doPropagateFlags: 

547 del self.sourceTableHandles 

548 del self.finalizedSourceTableHandles 

549 else: 

550 # Check for types of flags required. 

551 if not config.propagateFlags.source_flags: 

552 del self.sourceTableHandles 

553 if not config.propagateFlags.finalized_source_flags: 

554 del self.finalizedSourceTableHandles 

555 # TODO[DM-47797]: only the 'if' block contents here should survive. 

556 if config.inputCatalog == "deblendedCatalog": 

557 del self.inputCatalog 

558 if not config.doAddFootprints: 

559 del self.scarletModels 

560 else: 

561 del self.deblendedCatalog 

562 del self.scarletModels 

563 

564 # TODO[DM-47797]: delete the conditionals below. 

565 if not config.doMatchSources: 

566 del self.refCat 

567 del self.matchResult 

568 

569 if not config.doWriteMatchesDenormalized: 

570 del self.denormMatches 

571 

572 if config.useCellCoadds: 

573 del self.exposure 

574 else: 

575 del self.exposure_cells 

576 del self.background 

577 

578 

579class MeasureMergedCoaddSourcesConfig(PipelineTaskConfig, 

580 pipelineConnections=MeasureMergedCoaddSourcesConnections): 

581 """Configuration parameters for the MeasureMergedCoaddSourcesTask 

582 """ 

583 inputCatalog = ChoiceField( 

584 dtype=str, 

585 default="deblendedCatalog", 

586 allowed={ 

587 "deblendedCatalog": "Output catalog from ScarletDeblendTask", 

588 "deblendedFlux": "Output catalog from SourceDeblendTask", 

589 "mergeDet": "The merged detections before deblending." 

590 }, 

591 doc="The name of the input catalog.", 

592 # TODO[DM-47797]: remove this config option and anything using it. 

593 deprecated="Support for old deblender outputs will be removed after v29.", 

594 ) 

595 doAddFootprints = Field(dtype=bool, 

596 default=True, 

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

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

599 "otherwise this should be False.") 

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

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

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

603 "model footprints.") 

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

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

606 "saving to disk. " 

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

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

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

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

611 doPropagateFlags = Field( 

612 dtype=bool, default=True, 

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

614 ) 

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

616 doMatchSources = Field( 

617 dtype=bool, 

618 default=False, 

619 doc="Match sources to reference catalog?", 

620 deprecated="Reference matching in measureCoaddSources will be removed after v29.", 

621 ) 

622 match = ConfigurableField( 

623 target=DirectMatchTask, 

624 doc="Matching to reference catalog", 

625 deprecated="Reference matching in measureCoaddSources will be removed after v29.", 

626 ) 

627 doWriteMatchesDenormalized = Field( 

628 dtype=bool, 

629 default=False, 

630 doc=("Write reference matches in denormalized format? " 

631 "This format uses more disk space, but is more convenient to read."), 

632 deprecated="Reference matching in measureCoaddSources will be removed after v29.", 

633 ) 

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

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

636 checkUnitsParseStrict = Field( 

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

638 dtype=str, 

639 default="raise", 

640 ) 

641 doApCorr = Field( 

642 dtype=bool, 

643 default=True, 

644 doc="Apply aperture corrections" 

645 ) 

646 applyApCorr = ConfigurableField( 

647 target=ApplyApCorrTask, 

648 doc="Subtask to apply aperture corrections" 

649 ) 

650 doRunCatalogCalculation = Field( 

651 dtype=bool, 

652 default=True, 

653 doc='Run catalogCalculation task' 

654 ) 

655 catalogCalculation = ConfigurableField( 

656 target=CatalogCalculationTask, 

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

658 ) 

659 

660 hasFakes = Field( 

661 dtype=bool, 

662 default=False, 

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

664 ) 

665 idGenerator = SkyMapIdGeneratorConfig.make_field() 

666 

667 @property 

668 def refObjLoader(self): 

669 return self.match.refObjLoader 

670 

671 def setDefaults(self): 

672 super().setDefaults() 

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

674 'base_Variance', 

675 'base_LocalPhotoCalib', 

676 'base_LocalWcs'] 

677 

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

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

680 # still set flags with it here. 

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

682 'INEXACT_PSF'] 

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

684 'INEXACT_PSF'] 

685 

686 def validate(self): 

687 super().validate() 

688 

689 if not self.doMatchSources and self.doWriteMatchesDenormalized: 

690 raise ValueError("Cannot set doWriteMatchesDenormalized if doMatchSources is False.") 

691 

692 

693class MeasureMergedCoaddSourcesTask(PipelineTask): 

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

695 

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

697 measurement in each coadd. 

698 

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

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

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

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

703 

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

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

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

707 propagated to the coadd sources. 

708 

709 Optionally, we can match the coadd sources to an external reference 

710 catalog. 

711 

712 After MeasureMergedCoaddSourcesTask has been run on multiple coadds, we 

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

714 processing procedure will merge these measurements into a suitable catalog 

715 for driving forced photometry. 

716 

717 Parameters 

718 ---------- 

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

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

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

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

723 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`, optional 

724 An instance of ReferenceObjectLoader that supplies an external reference 

725 catalog. May be None if the loader can be constructed from the butler argument or all steps 

726 requiring a reference catalog are disabled. 

727 initInputs : `dict`, optional 

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

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

730 **kwargs 

731 Additional keyword arguments. 

732 """ 

733 

734 _DefaultName = "measureCoaddSources" 

735 ConfigClass = MeasureMergedCoaddSourcesConfig 

736 

737 def __init__(self, schema=None, peakSchema=None, refObjLoader=None, initInputs=None, 

738 **kwargs): 

739 super().__init__(**kwargs) 

740 self.deblended = self.config.inputCatalog.startswith("deblended") 

741 self.inputCatalog = "Coadd_" + self.config.inputCatalog 

742 if initInputs is not None: 

743 schema = initInputs['inputSchema'].schema 

744 if schema is None: 

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

746 self.schemaMapper = afwTable.SchemaMapper(schema) 

747 self.schemaMapper.addMinimalSchema(schema) 

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

749 self.algMetadata = PropertyList() 

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

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

752 # TODO[DM-47797]: remove match subtask 

753 if self.config.doMatchSources: 

754 self.makeSubtask("match", refObjLoader=refObjLoader) 

755 if self.config.doPropagateFlags: 

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

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

758 if self.config.doApCorr: 

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

760 if self.config.doRunCatalogCalculation: 

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

762 

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

764 

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

766 inputs = butlerQC.get(inputRefs) 

767 

768 # TODO[DM-47797]: remove this block 

769 if self.config.doMatchSources: 

770 refObjLoader = ReferenceObjectLoader([ref.datasetRef.dataId for ref in inputRefs.refCat], 

771 inputs.pop('refCat'), 

772 name=self.config.connections.refCat, 

773 config=self.config.refObjLoader, 

774 log=self.log) 

775 self.match.setRefObjLoader(refObjLoader) 

776 

777 if self.config.useCellCoadds: 

778 multiple_cell_coadd = inputs.pop("exposure_cells") 

779 stitched_coadd = multiple_cell_coadd.stitch() 

780 exposure = stitched_coadd.asExposure() 

781 background = inputs.pop("background") 

782 exposure.image -= background.getImage() 

783 

784 ccdInputs = stitched_coadd.ccds 

785 apCorrMap = stitched_coadd.ap_corr_map 

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

787 else: 

788 exposure = inputs.pop("exposure") 

789 # Set psfcache 

790 # move this to run after gen2 deprecation 

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

792 

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

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

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

796 

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

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

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

800 

801 # Transform inputCatalog 

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

803 sources = afwTable.SourceCatalog(table) 

804 # Load the correct input catalog 

805 if "scarletCatalog" in inputs: 

806 inputCatalog = inputs.pop("scarletCatalog") 

807 catalogRef = inputRefs.scarletCatalog 

808 else: 

809 inputCatalog = inputs.pop("inputCatalog") 

810 catalogRef = inputRefs.inputCatalog 

811 sources.extend(inputCatalog, self.schemaMapper) 

812 del inputCatalog 

813 # Add the HeavyFootprints to the deblended sources 

814 if self.config.doAddFootprints: 

815 modelData = inputs.pop('scarletModels') 

816 if self.config.doConserveFlux: 

817 imageForRedistribution = exposure 

818 else: 

819 imageForRedistribution = None 

820 updateCatalogFootprints( 

821 modelData=modelData, 

822 catalog=sources, 

823 band=band, 

824 imageForRedistribution=imageForRedistribution, 

825 removeScarletData=True, 

826 updateFluxColumns=True, 

827 ) 

828 table = sources.getTable() 

829 table.setMetadata(self.algMetadata) # Capture algorithm metadata to write out to the source catalog. 

830 

831 skyMap = inputs.pop('skyMap') 

832 tractNumber = catalogRef.dataId['tract'] 

833 tractInfo = skyMap[tractNumber] 

834 patchInfo = tractInfo.getPatchInfo(catalogRef.dataId['patch']) 

835 skyInfo = Struct( 

836 skyMap=skyMap, 

837 tractInfo=tractInfo, 

838 patchInfo=patchInfo, 

839 wcs=tractInfo.getWcs(), 

840 bbox=patchInfo.getOuterBBox() 

841 ) 

842 

843 if self.config.doPropagateFlags: 

844 if "sourceTableHandles" in inputs: 

845 sourceTableHandles = inputs.pop("sourceTableHandles") 

846 sourceTableHandleDict = {handle.dataId["visit"]: handle for handle in sourceTableHandles} 

847 else: 

848 sourceTableHandleDict = None 

849 if "finalizedSourceTableHandles" in inputs: 

850 finalizedSourceTableHandles = inputs.pop("finalizedSourceTableHandles") 

851 finalizedSourceTableHandleDict = {handle.dataId["visit"]: handle 

852 for handle in finalizedSourceTableHandles} 

853 else: 

854 finalizedSourceTableHandleDict = None 

855 if "finalVisitSummaryHandles" in inputs: 

856 finalVisitSummaryHandles = inputs.pop("finalVisitSummaryHandles") 

857 finalVisitSummaryHandleDict = {handle.dataId["visit"]: handle