Coverage for python/lsst/ap/association/diaPipe.py: 85%

1# 

2# LSST Data Management System 

3# Copyright 2008-2016 AURA/LSST. 

4# 

5# This product includes software developed by the 

6# LSST Project (http://www.lsst.org/). 

7# 

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

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

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

11# (at your option) any later version. 

12# 

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

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

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

16# GNU General Public License for more details. 

17# 

18# You should have received a copy of the LSST License Statement and 

19# the GNU General Public License along with this program. If not, 

20# see <https://www.lsstcorp.org/LegalNotices/>. 

21# 

22 

23"""PipelineTask for associating DiaSources with previous DiaObjects. 

24 

25Additionally performs forced photometry on the calibrated and difference 

26images at the updated locations of DiaObjects. 

27""" 

28 

29__all__ = ("DiaPipelineConfig", 

30 "DiaPipelineTask", 

31 "DiaPipelineConnections") 

32 

33 

34import dataclasses 

35 

36import lsst.dax.apdb as daxApdb 

37import lsst.pex.config as pexConfig 

38import lsst.pipe.base as pipeBase 

39import lsst.pipe.base.connectionTypes as connTypes 

40import lsst.sphgeom 

41 

42from astropy.table import Table 

43import numpy as np 

44import pandas as pd 

45from lsst.ap.association import ( 

46 AssociationTask, 

47 DiaForcedSourceTask, 

48 PackageAlertsTask) 

49 

50from lsst.ap.association.loadDiaCatalogs import loadDiaObjectsFromApdb, loadDiaSourcesFromApdb, \ 

51 loadDiaForcedSourcesFromApdb 

52from lsst.ap.association.utils import makeEmptyForcedSourceTable, getRegion, paddedRegion, readSchemaFromApdb 

53from lsst.daf.base import DateTime 

54from lsst.meas.base import DetectorVisitIdGeneratorConfig, \ 

55 DiaObjectCalculationTask 

56from lsst.pipe.tasks.schemaUtils import convertDataFrameToSdmSchema, checkSdmSchemaColumns, \ 

57 dropEmptyColumns, make_empty_catalog 

58from lsst.pipe.tasks.ssoAssociation import SolarSystemAssociationTask 

59from lsst.utils.timer import timeMethod, duration_from_timeMethod 

60 

61 

62class TooManyDiaObjectsError(pipeBase.AlgorithmError): 

63 """Raised if there are an unusually large number of unassociated DiaSources. 

64 This is usually indicative of an image subtraction error, and needs to be 

65 caught before updating the APDB with a large number of spurious entries. 

66 """ 

67 def __init__(self, *, nNewDiaObjects, threshold): 

68 msg = ("Aborting processing before writing to the APDB." 

69 f" {nNewDiaObjects} new DiaObjects would be created, which exceeds the" 

70 f" configured maximum of {threshold}") 

71 super().__init__(msg) 

72 self.nNewDiaObjects = nNewDiaObjects 

73 self.threshold = threshold 

74 

75 @property 

76 def metadata(self): 

77 return {"nNewDiaObjects": self.nNewDiaObjects, 

78 "threshold": self.threshold 

79 } 

80 

81 

82class PostApdbUpdateError(pipeBase.AlgorithmError): 

83 """Raised for any error that occurs after the APDB has been updated. 

84 This allows partial outputs to be written, and signals that processing 

85 can't be retried. 

86 """ 

87 def __init__(self, *, errorMsg): 

88 msg = ("Aborting processing after writing to the APDB. " 

89 "This image cannot be retried! " 

90 "The original error was:\n" 

91 f"{errorMsg}") 

92 super().__init__(msg) 

93 

94 @property 

95 def metadata(self): 

96 return {} 

97 

98 

99class DiaPipelineConnections( 

100 pipeBase.PipelineTaskConnections, 

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

102 defaultTemplates={"coaddName": "deep", "fakesType": ""}): 

103 """Butler connections for DiaPipelineTask. 

104 """ 

105 diaSourceTable = connTypes.Input( 

106 doc="Catalog of calibrated DiaSources.", 

107 name="{fakesType}{coaddName}Diff_diaSrcTable", 

108 storageClass="DataFrame", 

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

110 ) 

111 solarSystemObjectTable = connTypes.Input( 

112 doc="Catalog of SolarSolarSystem objects expected to be observable in " 

113 "this detectorVisit.", 

114 name="preloaded_SsObjects", 

115 storageClass="ArrowAstropy", 

116 dimensions=("instrument", "group", "detector"), 

117 minimum=0, 

118 ) 

119 diffIm = connTypes.Input( 

120 doc="Difference image on which the DiaSources were detected.", 

121 name="{fakesType}{coaddName}Diff_differenceExp", 

122 storageClass="ExposureF", 

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

124 ) 

125 exposure = connTypes.Input( 

126 doc="Calibrated exposure differenced with a template image during " 

127 "image differencing.", 

128 name="{fakesType}calexp", 

129 storageClass="ExposureF", 

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

131 ) 

132 template = connTypes.Input( 

133 doc="Warped template used to create `subtractedExposure`. Not PSF " 

134 "matched.", 

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

136 storageClass="ExposureF", 

137 name="{fakesType}{coaddName}Diff_templateExp", 

138 ) 

139 preloadedDiaObjects = connTypes.Input( 

140 doc="DiaObjects preloaded from the APDB.", 

141 name="preloaded_diaObjects", 

142 storageClass="DataFrame", 

143 dimensions=("instrument", "group", "detector"), 

144 ) 

145 preloadedDiaSources = connTypes.Input( 

146 doc="DiaSources preloaded from the APDB.", 

147 name="preloaded_diaSources", 

148 storageClass="DataFrame", 

149 dimensions=("instrument", "group", "detector"), 

150 ) 

151 preloadedDiaForcedSources = connTypes.Input( 

152 doc="DiaForcedSources preloaded from the APDB.", 

153 name="preloaded_diaForcedSources", 

154 storageClass="DataFrame", 

155 dimensions=("instrument", "group", "detector"), 

156 ) 

157 apdbMarker = connTypes.Output( 

158 doc="Marker dataset storing the configuration of the Apdb for each " 

159 "visit/detector. Used to signal the completion of the pipeline.", 

160 name="apdb_marker", 

161 storageClass="Config", 

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

163 ) 

164 associatedDiaSources = connTypes.Output( 

165 doc="Optional output storing the DiaSource catalog after matching, " 

166 "calibration, and standardization for insertion into the Apdb.", 

167 name="{fakesType}{coaddName}Diff_assocDiaSrc", 

168 storageClass="ArrowAstropy", 

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

170 ) 

171 associatedSsSources = connTypes.Output( 

172 doc="Optional output storing ssSource data computed during association.", 

173 name="{fakesType}{coaddName}Diff_associatedSsSources", 

174 storageClass="ArrowAstropy", 

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

176 ) 

177 unassociatedSsObjects = connTypes.Output( 

178 doc="Expected locations of an ssObject with no source", 

179 name="ssUnassociatedObjects", 

180 storageClass="ArrowAstropy", 

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

182 ) 

183 

184 diaForcedSources = connTypes.Output( 

185 doc="Optional output storing the forced sources computed at the diaObject positions.", 

186 name="{fakesType}{coaddName}Diff_diaForcedSrc", 

187 storageClass="ArrowAstropy", 

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

189 ) 

190 diaObjects = connTypes.Output( 

191 doc="Optional output storing the updated diaObjects associated to these sources.", 

192 name="{fakesType}{coaddName}Diff_diaObject", 

193 storageClass="ArrowAstropy", 

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

195 ) 

196 newDiaSources = connTypes.Output( 

197 doc="New diaSources not associated with an existing diaObject that" 

198 " were used to create a new diaObject", 

199 name="{fakesType}new_dia_source", 

200 storageClass="ArrowAstropy", 

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

202 ) 

203 marginalDiaSources = connTypes.Output( 

204 doc="Low SNR diaSources not associated with an existing diaObject that" 

205 " were rejected instead of creating a new diaObject", 

206 name="{fakesType}marginal_new_dia_source", 

207 storageClass="ArrowAstropy", 

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

209 ) 

210 

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

212 super().__init__(config=config) 

213 

214 if config.input_image_type == "future": 

215 self.diffIm = dataclasses.replace(self.diffIm, storageClass="DifferenceImage") 

216 self.exposure = dataclasses.replace(self.exposure, storageClass="VisitImage") 

217 self.template = dataclasses.replace(self.template, storageClass="DifferenceImage") 

218 

219 if not config.doWriteAssociatedSources: 219 ↛ 220line 219 didn't jump to line 220 because the condition on line 219 was never true

220 self.outputs.remove("associatedDiaSources") 

221 self.outputs.remove("diaForcedSources") 

222 self.outputs.remove("diaObjects") 

223 self.outputs.remove("newDiaSources") 

224 self.outputs.remove("marginalDiaSources") 

225 else: 

226 if not config.doRunForcedMeasurement: 

227 self.outputs.remove("diaForcedSources") 

228 if not config.filterUnAssociatedSources: 228 ↛ 229line 228 didn't jump to line 229 because the condition on line 228 was never true

229 self.outputs.remove("newDiaSources") 

230 self.outputs.remove("marginalDiaSources") 

231 if not config.doSolarSystemAssociation: 

232 self.inputs.remove("solarSystemObjectTable") 

233 if (not config.doWriteAssociatedSources) or (not config.doSolarSystemAssociation): 

234 self.outputs.remove("associatedSsSources") 

235 self.outputs.remove("unassociatedSsObjects") 

236 if config.doReloadAllApdbCatalogs: 

237 # If this is set, the complete history will be read from the APDB 

238 # during association and the preloaded catalogs will not be used. 

239 self.inputs.remove("preloadedDiaObjects") 

240 self.inputs.remove("preloadedDiaSources") 

241 self.inputs.remove("preloadedDiaForcedSources") 

242 

243 def adjustQuantum(self, inputs, outputs, label, dataId): 

244 """Override to make adjustments to `lsst.daf.butler.DatasetRef` objects 

245 in the `lsst.daf.butler.core.Quantum` during the graph generation stage 

246 of the activator. 

247 

248 This implementation checks to make sure that the filters in the dataset 

249 are compatible with AP processing as set by the Apdb/DPDD schema. 

250 

251 Parameters 

252 ---------- 

253 inputs : `dict` 

254 Dictionary whose keys are an input (regular or prerequisite) 

255 connection name and whose values are a tuple of the connection 

256 instance and a collection of associated `DatasetRef` objects. 

257 The exact type of the nested collections is unspecified; it can be 

258 assumed to be multi-pass iterable and support `len` and ``in``, but 

259 it should not be mutated in place. In contrast, the outer 

260 dictionaries are guaranteed to be temporary copies that are true 

261 `dict` instances, and hence may be modified and even returned; this 

262 is especially useful for delegating to `super` (see notes below). 

263 outputs : `dict` 

264 Dict of output datasets, with the same structure as ``inputs``. 

265 label : `str` 

266 Label for this task in the pipeline (should be used in all 

267 diagnostic messages). 

268 data_id : `lsst.daf.butler.DataCoordinate` 

269 Data ID for this quantum in the pipeline (should be used in all 

270 diagnostic messages). 

271 

272 Returns 

273 ------- 

274 adjusted_inputs : `dict` 

275 Dict of the same form as ``inputs`` with updated containers of 

276 input `DatasetRef` objects. Connections that are not changed 

277 should not be returned at all. Datasets may only be removed, not 

278 added. Nested collections may be of any multi-pass iterable type, 

279 and the order of iteration will set the order of iteration within 

280 `PipelineTask.runQuantum`. 

281 adjusted_outputs : `dict` 

282 Dict of updated output datasets, with the same structure and 

283 interpretation as ``adjusted_inputs``. 

284 

285 Raises 

286 ------ 

287 ScalarError 

288 Raised if any `Input` or `PrerequisiteInput` connection has 

289 ``multiple`` set to `False`, but multiple datasets. 

290 NoWorkFound 

291 Raised to indicate that this quantum should not be run; not enough 

292 datasets were found for a regular `Input` connection, and the 

293 quantum should be pruned or skipped. 

294 FileNotFoundError 

295 Raised to cause QuantumGraph generation to fail (with the message 

296 included in this exception); not enough datasets were found for a 

297 `PrerequisiteInput` connection. 

298 """ 

299 _, refs = inputs["diffIm"] 

300 for ref in refs: 

301 if ref.dataId["band"] not in self.config.validBands: 

302 raise ValueError( 

303 f"Requested '{ref.dataId['band']}' not in " 

304 "DiaPipelineConfig.validBands. To process bands not in " 

305 "the standard Rubin set (ugrizy) you must add the band to " 

306 "the validBands list in DiaPipelineConfig and add the " 

307 "appropriate columns to the Apdb schema.") 

308 return super().adjustQuantum(inputs, outputs, label, dataId) 

309 

310 

311class DiaPipelineConfig(pipeBase.PipelineTaskConfig, 

312 pipelineConnections=DiaPipelineConnections): 

313 """Config for DiaPipelineTask. 

314 """ 

315 input_image_type = pexConfig.ChoiceField[str]( 

316 doc="Type of the input images. run() takes lsst.afw.image.ExposureF" 

317 " either way.", 

318 allowed={ 

319 "legacy": "Read `lsst.afw.image.ExposureF`.", 

320 "future": "Read `lsst.images.VisitImage` for the science image and" 

321 " `lsst.images.DifferenceImage` for the difference image and" 

322 " template.", 

323 }, 

324 optional=False, 

325 default="legacy", 

326 ) 

327 coaddName = pexConfig.Field( 

328 doc="coadd name: typically one of deep, goodSeeing, or dcr", 

329 dtype=str, 

330 default="deep", 

331 ) 

332 apdb_config_url = pexConfig.Field( 

333 dtype=str, 

334 default=None, 

335 optional=False, 

336 doc="A config file specifying the APDB and its connection parameters, " 

337 "typically written by the apdb-cli command-line utility. " 

338 "The database must already be initialized.", 

339 ) 

340 validBands = pexConfig.ListField( 

341 dtype=str, 

342 default=["u", "g", "r", "i", "z", "y"], 

343 doc="List of bands that are valid for AP processing. To process a " 

344 "band not on this list, the appropriate band specific columns " 

345 "must be added to the Apdb schema in dax_apdb.", 

346 ) 

347 associator = pexConfig.ConfigurableField( 

348 target=AssociationTask, 

349 doc="Task used to associate DiaSources with DiaObjects.", 

350 ) 

351 doSolarSystemAssociation = pexConfig.Field( 

352 dtype=bool, 

353 default=True, 

354 doc="Process SolarSystem objects through the pipeline.", 

355 ) 

356 solarSystemAssociator = pexConfig.ConfigurableField( 

357 target=SolarSystemAssociationTask, 

358 doc="Task used to associate DiaSources with SolarSystemObjects.", 

359 ) 

360 diaCalculation = pexConfig.ConfigurableField( 

361 target=DiaObjectCalculationTask, 

362 doc="Task to compute summary statistics for DiaObjects.", 

363 ) 

364 doReloadDiaObjects = pexConfig.Field( 

365 dtype=bool, 

366 default=True, 

367 doc="Drop preloaded DiaObjects and reload them from the APDB?" 

368 "Used in production when the very latest objects from the APDB " 

369 "are needed. Ignored and superceded if `doReloadAllApdbCatalogs` " 

370 "is set.", 

371 ) 

372 doReloadAllApdbCatalogs = pexConfig.Field( 

373 dtype=bool, 

374 default=False, 

375 doc="Read the complete DiaObject, DiaSource, and DiaForcedSource " 

376 "history from the APDB during association? Use only during " 

377 "reprocessing, and never in Prompt Processing.", 

378 ) 

379 angleMargin = pexConfig.RangeField( 

380 doc="Padding to add when reloading catalogs from the APDB. " 

381 "Smaller than the padding used for ``LoadDiaCatalogsTask`` because " 

382 "this task is after astrometric calibration.", 

383 dtype=float, 

384 default=2, 

385 min=0, 

386 ) 

387 doRunForcedMeasurement = pexConfig.Field( 

388 dtype=bool, 

389 default=True, 

390 deprecated="Added to allow disabling forced sources for performance " 

391 "reasons during the ops rehearsal. " 

392 "It is expected to be removed.", 

393 doc="Run forced measurement on all of the diaObjects? " 

394 "This should only be turned off for debugging purposes.", 

395 ) 

396 diaForcedSource = pexConfig.ConfigurableField( 

397 target=DiaForcedSourceTask, 

398 doc="Task used for force photometer DiaObject locations in direct and " 

399 "difference images.", 

400 ) 

401 forcedReliabilityThreshold = pexConfig.Field( 

402 dtype=float, 

403 default=0.5, 

404 doc="Minimum reliability score of constituent diaSources to use when " 

405 "determining whether to calculate forced photometry for a " 

406 "diaObject.", 

407 ) 

408 forcedTrailLengthThreshold = pexConfig.Field( 

409 dtype=float, 

410 default=1.4, 

411 doc="Maximum trail length (arseconds) of constituent diaSources to use " 

412 "when determining whether to calculate forced photometry for a " 

413 "diaObject.", 

414 ) 

415 forcedBadFlags = pexConfig.ListField( 

416 dtype=str, 

417 default=["glint_trail", "isDipole"], 

418 doc="Flags of constituent diaSources to exclude when determining " 

419 "whether to calculate forced photometry for a diaObject.", 

420 ) 

421 alertPackager = pexConfig.ConfigurableField( 

422 target=PackageAlertsTask, 

423 doc="Subtask for packaging Ap data into alerts.", 

424 ) 

425 doPackageAlerts = pexConfig.Field( 

426 dtype=bool, 

427 default=False, 

428 doc="Package Dia-data into serialized alerts for distribution and " 

429 "write them to disk.", 

430 ) 

431 doWriteAssociatedSources = pexConfig.Field( 

432 dtype=bool, 

433 default=True, 

434 doc="Write out associated DiaSources, DiaForcedSources, and DiaObjects, " 

435 "formatted following the Science Data Model.", 

436 ) 

437 imagePixelMargin = pexConfig.RangeField( 

438 dtype=int, 

439 default=10, 

440 min=0, 

441 doc="Pad the image by this many pixels before removing off-image " 

442 "diaObjects for association.", 

443 ) 

444 filterUnAssociatedSources = pexConfig.Field( 

445 dtype=bool, 

446 default=True, 

447 doc="Check unassociated diaSources for quality before creating new ." 

448 "diaObjects. DiaSources that are associated with existing diaObjects " 

449 "will not be affected." 

450 ) 

451 newObjectBadFlags = pexConfig.ListField( 

452 dtype=str, 

453 default=("centroid_flag", 

454 "pixelFlags_crCenter", 

455 "pixelFlags_nodataCenter", 

456 "pixelFlags_interpolatedCenter", 

457 "pixelFlags_saturatedCenter", 

458 "pixelFlags_suspectCenter", 

459 "pixelFlags_streakCenter", 

460 "glint_trail"), 

461 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

462 "these flags set before creating new diaObjects.", 

463 ) 

464 maxNewDiaObjects = pexConfig.RangeField( 

465 dtype=float, 

466 default=0, 

467 min=0, 

468 doc="Maximum number of new DiaObjects to create before raising an error." 

469 "Set to zero to disable.", 

470 ) 

471 newObjectSnrThreshold = pexConfig.RangeField( 

472 dtype=float, 

473 default=0, 

474 min=0, 

475 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

476 "Abs(flux/fluxErr) less than this threshold before creating new " 

477 "diaObjects." 

478 "Set to zero to disable.", 

479 ) 

480 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

481 dtype=float, 

482 default=0, 

483 min=0, 

484 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

485 "signal-to-noise ratios less than this threshold if they have" 

486 " low reliability scores before creating new diaObjects." 

487 "Set to zero to disable.", 

488 ) 

489 newObjectReliabilityThreshold = pexConfig.RangeField( 

490 dtype=float, 

491 default=0.1, 

492 min=0, 

493 max=1, 

494 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

495 "reliability scores less than this threshold before creating new " 

496 "diaObjects." 

497 "Set to zero to disable.", 

498 ) 

499 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

500 dtype=float, 

501 default=0.1, 

502 min=0, 

503 max=1, 

504 doc="If `filterUnAssociatedSources` is set, exclude diaSources with " 

505 "low signal-to-noise and reliability scores below this threshold " 

506 "before creating new diaObjects." 

507 "Set to zero to disable.", 

508 ) 

509 newObjectFluxField = pexConfig.Field( 

510 dtype=str, 

511 default="apFlux", 

512 doc="Name of the flux field in the standardized diaSource catalog to " 

513 "use for checking the signal-to-noise before creating new diaObjects.", 

514 ) 

515 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

516 

517 def setDefaults(self): 

518 self.diaCalculation.plugins = ["ap_meanPosition", 

519 "ap_nDiaSources", 

520 "ap_meanFlux", 

521 "ap_sigmaFlux", 

522 "ap_minMaxFlux", 

523 "ap_maxSlopeFlux", 

524 "ap_meanErrFlux", 

525 "ap_meanTotFlux"] 

526 

527 

528class DiaPipelineTask(pipeBase.PipelineTask): 

529 """Task for loading, associating and storing Difference Image Analysis 

530 (DIA) Objects and Sources. 

531 """ 

532 ConfigClass = DiaPipelineConfig 

533 _DefaultName = "diaPipe" 

534 

535 def __init__(self, initInputs=None, **kwargs): 

536 super().__init__(**kwargs) 

537 self.apdb = daxApdb.Apdb.from_uri(self.config.apdb_config_url) 

538 self.schema = readSchemaFromApdb(self.apdb) 

539 self.makeSubtask("associator") 

540 self.makeSubtask("diaCalculation") 

541 if self.config.doRunForcedMeasurement: 

542 self.makeSubtask("diaForcedSource") 

543 if self.config.doPackageAlerts: 

544 self.makeSubtask("alertPackager") 

545 if self.config.doSolarSystemAssociation: 

546 self.makeSubtask("solarSystemAssociator") 

547 if self.config.filterUnAssociatedSources: 

548 columns = [self.config.newObjectFluxField, 

549 self.config.newObjectFluxField + "Err", 

550 "reliability"] 

551 columns += self.config.newObjectBadFlags 

552 

553 missing = checkSdmSchemaColumns(self.schema, columns, "DiaSource") 

554 if missing: 

555 raise pipeBase.InvalidQuantumError("Field %s not in the DiaSource schema" % missing) 

556 

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

558 inputs = butlerQC.get(inputRefs) 

559 if self.config.input_image_type == "future": 

560 for name in ("diffIm", "exposure", "template"): 

561 inputs[name] = inputs[name].to_legacy() 

562 inputs["idGenerator"] = self.config.idGenerator.apply(butlerQC.quantum.dataId) 

563 inputs["band"] = butlerQC.quantum.dataId["band"] 

564 inputs["legacySolarSystemTable"] = None 

565 if not self.config.doSolarSystemAssociation: 565 ↛ 566line 565 didn't jump to line 566 because the condition on line 565 was never true

566 inputs["solarSystemObjectTable"] = None 

567 

568 associationResults = pipeBase.Struct( 

569 apdbMarker=None, 

570 associatedDiaSources=None, 

571 diaForcedSources=None, 

572 diaObjects=None, 

573 associatedSsSources=None, 

574 unassociatedSsObjects=None, 

575 newDiaSources=None, 

576 marginalDiaSources=None, 

577 ) 

578 try: 

579 self.run(**inputs, associationResults=associationResults) 

580 except pipeBase.AlgorithmError as e: 

581 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

582 e, 

583 self, 

584 log=self.log 

585 ) 

586 butlerQC.put(associationResults, outputRefs) 

587 raise error from e 

588 butlerQC.put(associationResults, outputRefs) 

589 

590 @timeMethod 

591 def run(self, 

592 diaSourceTable, 

593 legacySolarSystemTable, 

594 diffIm, 

595 exposure, 

596 template, 

597 preloadedDiaObjects=None, 

598 preloadedDiaSources=None, 

599 preloadedDiaForcedSources=None, 

600 band=None, 

601 idGenerator=None, 

602 solarSystemObjectTable=None, 

603 associationResults=None): 

604 """Process DiaSources and DiaObjects. 

605 

606 Load previous DiaObjects and their DiaSource history. Calibrate the 

607 values in the diaSourceCat. Associate new DiaSources with previous 

608 DiaObjects. Run forced photometry at the updated DiaObject locations. 

609 Store the results in the Alert Production Database (Apdb). 

610 

611 Parameters 

612 ---------- 

613 diaSourceTable : `pandas.DataFrame` 

614 Newly detected DiaSources. 

615 legacySolarSystemTable : `pandas.DataFrame` 

616 Not used 

617 diffIm : `lsst.afw.image.ExposureF` 

618 Difference image exposure in which the sources in ``diaSourceCat`` 

619 were detected. 

620 exposure : `lsst.afw.image.ExposureF` 

621 Calibrated exposure differenced with a template to create 

622 ``diffIm``. 

623 template : `lsst.afw.image.ExposureF` 

624 Template exposure used to create diffIm. 

625 preloadedDiaObjects : `pandas.DataFrame`, optional 

626 Previously detected DiaObjects, loaded from the APDB. 

627 `None` if ``doReloadAllApdbCatalogs`` is set. 

628 preloadedDiaSources : `pandas.DataFrame`, optional 

629 Previously detected DiaSources, loaded from the APDB. 

630 `None` if ``doReloadAllApdbCatalogs`` is set. 

631 preloadedDiaForcedSources : `pandas.DataFrame`, optional 

632 Catalog of previously detected forced DiaSources, from the APDB. 

633 `None` if ``doReloadAllApdbCatalogs`` is set. 

634 band : `str` 

635 The band in which the new DiaSources were detected. Required, 

636 despite the `None` default. 

637 idGenerator : `lsst.meas.base.IdGenerator` 

638 Object that generates source IDs and random number generator seeds. 

639 Required, despite the `None` default. 

640 solarSystemObjectTable : `astropy.table.Table`, optional 

641 Preloaded Solar System objects expected to be visible in the image. 

642 associationResults : `lsst.pipe.base.Struct`, optional 

643 Result struct that is modified to allow saving of partial outputs 

644 for some failure conditions. If the task completes successfully, 

645 this is also returned. 

646 

647 Returns 

648 ------- 

649 associationResults : `lsst.pipe.base.Struct` 

650 Results struct with components. 

651 

652 - ``apdbMarker`` : Marker dataset to store in the Butler indicating 

653 that this ccdVisit has completed successfully. 

654 (`lsst.dax.apdb.ApdbConfig`) 

655 - ``associatedDiaSources`` : Catalog of newly associated 

656 DiaSources. (`pandas.DataFrame`) 

657 - ``diaForcedSources`` : Catalog of new and previously detected 

658 forced DiaSources. (`pandas.DataFrame`) 

659 - ``diaObjects`` : Updated table of DiaObjects. (`pandas.DataFrame`) 

660 - ``associatedSsSources`` : Catalog of ssSource records. 

661 (`pandas.DataFrame`) 

662 

663 Raises 

664 ------ 

665 RuntimeError 

666 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

667 ValueError 

668 Raised if ``band`` or ``idGenerator`` is `None`. 

669 """ 

670 # These default to `None` only so that the preloaded catalogs, which 

671 # precede them, can be omitted. 

672 if band is None: 

673 raise ValueError("band is required.") 

674 if idGenerator is None: 

675 raise ValueError("idGenerator is required.") 

676 

677 if associationResults is None: 677 ↛ 679line 677 didn't jump to line 679 because the condition on line 677 was always true

678 associationResults = pipeBase.Struct() 

679 self._validateExposure(exposure, "Science") 

680 self._validateExposure(diffIm, "Difference") 

681 self._validateExposure(template, "Template") 

682 

683 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

684 if legacySolarSystemTable is not None and solarSystemObjectTable is None: 684 ↛ 685line 684 didn't jump to line 685 because the condition on line 684 was never true

685 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

686 region = getRegion(exposure) 

687 if self.config.doReloadDiaObjects or self.config.doReloadAllApdbCatalogs: 

688 try: 

689 preloadedDiaObjects = self.loadRefreshedDiaObjects(region, preloadedDiaObjects) 

690 except Exception as e: 

691 if self.config.doReloadAllApdbCatalogs: 

692 # We can't continue if the diaObjects were not loaded and 

693 # ``doReloadAllApdbCatalogs`` is set, because there will be 

694 # no preloaded diaObjects to fall back on. 

695 raise 

696 self.log.warning("Error encountered while attempting to load " 

697 "the latest diaObjects from the APDB. Processing " 

698 "will continue with only the diaObjects from " 

699 "preload.", exc_info=e) 

700 finally: 

701 self.metadata["loadDiaObjectsDuration"] = duration_from_timeMethod( 

702 self.metadata, "loadRefreshedDiaObjects", clock="Utc" 

703 ) 

704 self.log.verbose("Re-loading DiaObjects: Took %.4f seconds", 

705 self.metadata["loadDiaObjectsDuration"]) 

706 

707 else: 

708 self.metadata["loadDiaObjectsDuration"] = -1 

709 

710 if self.config.doReloadAllApdbCatalogs: 

711 # The metadata time records must use the same names as 

712 # LoadDiaCatalogsTask to keep the metrics upload consistent. 

713 visitTime = exposure.visitInfo.date.toAstropy() 

714 try: 

715 preloadedDiaSources = self.loadRefreshedDiaSources(region, preloadedDiaObjects, visitTime) 

716 finally: 

717 self.metadata["loadDiaSourcesDuration"] = duration_from_timeMethod( 

718 self.metadata, "loadRefreshedDiaSources", clock="Utc" 

719 ) 

720 self.log.verbose("Re-loading DiaSources: Took %.4f seconds", 

721 self.metadata["loadDiaSourcesDuration"]) 

722 try: 

723 preloadedDiaForcedSources = self.loadRefreshedDiaForcedSources( 

724 region, preloadedDiaObjects, visitTime) 

725 finally: 

726 self.metadata["loadDiaForcedSourcesDuration"] = duration_from_timeMethod( 

727 self.metadata, "loadRefreshedDiaForcedSources", clock="Utc" 

728 ) 

729 self.log.verbose("Re-loading DiaForcedSources: Took %.4f seconds", 

730 self.metadata["loadDiaForcedSourcesDuration"]) 

731 else: 

732 self.metadata["loadDiaSourcesDuration"] = -1 

733 self.metadata["loadDiaForcedSourcesDuration"] = -1 

734 

735 self.checkTableIndex(preloadedDiaSources, index=["diaObjectId", "band", "diaSourceId"]) 

736 self.checkTableIndex(preloadedDiaObjects, index="diaObjectId") 

737 self.checkTableIndex(preloadedDiaForcedSources, index=["diaObjectId", "diaForcedSourceId"]) 

738 

739 if not preloadedDiaObjects.empty: 739 ↛ 744line 739 didn't jump to line 744 because the condition on line 739 was always true

740 # Include a small buffer outside the image so that we can associate sources near the edge 

741 diaObjects, _ = self.purgeDiaObjects(diffIm.getBBox(), diffIm.getWcs(), preloadedDiaObjects, 

742 buffer=self.config.imagePixelMargin) 

743 else: 

744 self.log.info("Preloaded DiaObject table is empty.") 

745 diaObjects = preloadedDiaObjects 

746 

747 # Associate DiaSources with DiaObjects 

748 assocResults = self.associateDiaSources(diaSourceTable, solarSystemObjectTable, diffIm, diaObjects) 

749 

750 # Set unassociated diaObjectIds and ssObjectIds to NULL, and convert to SDM schema format 

751 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

752 assocResults.associatedDiaSources, 

753 "DiaSource", 

754 nullColumns=["diaObjectId", "ssObjectId"] 

755 ) 

756 

757 # Merge associated diaSources 

758 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

759 preloadedDiaSources, 

760 assocResults.associatedDiaSources, 

761 diaObjects, 

762 assocResults.newDiaObjects, 

763 diffIm 

764 ) 

765 

766 # Compute DiaObject Summary statistics from their full DiaSource 

767 # history. 

768 diaCalResult = self.diaCalculation.run( 

769 mergedDiaObjects, 

770 mergedDiaSourceHistory, 

771 updatedDiaObjectIds, 

772 self.config.validBands) 

773 updatedDiaObjects = convertDataFrameToSdmSchema(self.schema, diaCalResult.updatedDiaObjects, 

774 tableName="DiaObject", skipIndex=True) 

775 

776 # Test for duplication in the updated DiaObjects. 

777 if self.testDataFrameIndex(diaCalResult.diaObjectCat): 777 ↛ 778line 777 didn't jump to line 778 because the condition on line 777 was never true

778 raise RuntimeError( 

779 "Duplicate DiaObjects (loaded + updated) created after " 

780 "DiaCalculation. This is unexpected behavior and should be " 

781 "reported. Exiting.") 

782 if self.testDataFrameIndex(updatedDiaObjects): 782 ↛ 783line 782 didn't jump to line 783 because the condition on line 782 was never true

783 raise RuntimeError( 

784 "Duplicate DiaObjects (updated) created after " 

785 "DiaCalculation. This is unexpected behavior and should be " 

786 "reported. Exiting.") 

787 

788 # Forced source measurement 

789 if self.config.doRunForcedMeasurement: 

790 diaObjectsForced = self._selectGoodDiaObjects(diaCalResult.diaObjectCat, mergedDiaSourceHistory) 

791 diaForcedSources = self.runForcedMeasurement( 

792 diaObjectsForced, updatedDiaObjects, exposure, diffIm, idGenerator 

793 ) 

794 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

795 else: 

796 # alertPackager needs correct columns 

797 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

798 forcedSourceHistoryThreshold = 0 

799 

800 # Write results to Alert Production Database (APDB) 

801 try: 

802 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

803 diaForcedSources) 

804 finally: 

805 self.metadata["writeToApdbDuration"] = duration_from_timeMethod(self.metadata, "writeToApdb", 

806 clock="Utc") 

807 # A single log message is easier for Loki to parse than timeMethod's start+end pairs. 

808 self.log.verbose("writeToApdb: Took %.4f seconds", self.metadata["writeToApdbDuration"]) 

809 

810 # For historical reasons, apdbMarker is a Config even if it's not meant to be read. 

811 # A default Config is the cheapest way to satisfy the storage class. 

812 associationResults.apdbMarker = pexConfig.Config() 

813 associationResults.associatedDiaSources = assocResults.associatedDiaSources 

814 associationResults.diaForcedSources = diaForcedSources 

815 associationResults.diaObjects = diaCalResult.diaObjectCat 

816 associationResults.unassociatedSsObjects = assocResults.unassociatedSsObjects 

817 associationResults.newDiaSources = assocResults.newDiaSources 

818 associationResults.marginalDiaSources = assocResults.marginalDiaSources 

819 

820 # Catch *any* error after we have updated the APDB 

821 try: 

822 associatedSsSources = assocResults.associatedSsSources 

823 associatedSsSourcesPlusMpcorb = None 

824 if associatedSsSources is not None: 

825 mpcorbColumns = [col for col in associatedSsSources.columns if col[:7] == 'MPCORB_'] 

826 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

827 associatedSsSources = self.standardizeTable(associatedSsSources, "SSSource", nullColumns=[]) 

828 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

829 for mpcorbColumn in mpcorbColumns: 829 ↛ 830line 829 didn't jump to line 830 because the loop on line 829 never started

830 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

831 associationResults.associatedSsSources = associatedSsSources 

832 

833 # patch the otherwise-empty validityStart field for the alerts 

834 updatedDiaObjects['validityStartMjdTai'] = validityStart.get(system=DateTime.MJD, 

835 scale=DateTime.TAI) 

836 

837 # Package alerts 

838 if self.config.doPackageAlerts: 

839 # Append new forced sources to the full history 

840 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

841 tableName="DiaForcedSource") 

842 if self.testDataFrameIndex(diaForcedSourcesFull): 842 ↛ 843line 842 didn't jump to line 843 because the condition on line 842 was never true

843 self.log.warning( 

844 "Duplicate DiaForcedSources created after merge with " 

845 "history and new sources. This may cause downstream " 

846 "problems. Dropping duplicates.") 

847 # Drop duplicates via index and keep the first appearance. 

848 diaForcedSourcesFull = diaForcedSourcesFull[ 

849 ~diaForcedSourcesFull.index.duplicated(keep="first")] 

850 

851 self.alertPackager.run(assocResults.associatedDiaSources, 

852 updatedDiaObjects, 

853 preloadedDiaSources, 

854 diaForcedSourcesFull, 

855 diffIm, 

856 exposure, 

857 template, 

858 ssSrc=associatedSsSourcesPlusMpcorb, 

859 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

860 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

861 ) 

862 except Exception as e: 

863 # Catch *any* error after we have updated the APDB 

864 raise PostApdbUpdateError(errorMsg=repr(e)) from e 

865 

866 return associationResults 

867 

868 def _validateExposure(self, exposure, label): 

869 """Validate exposure metadata early to avoid failures after updating 

870 the APDB. 

871 

872 Parameters 

873 ---------- 

874 exposure : `lsst.afw.image.ExposureF` 

875 Calibrated science image that was subtracted. Corresponds to the 

876 ``exposure`` Butler connection (``{fakesType}calexp``). 

877 label : `str` 

878 Specification of the image being validated, for logging messages. 

879 

880 Raises 

881 ------ 

882 ValueError 

883 Raised when any required metadata field is missing or invalid. 

884 """ 

885 errors = [] 

886 

887 photo_calib = exposure.getPhotoCalib() 

888 if photo_calib is None: 888 ↛ 889line 888 didn't jump to line 889 because the condition on line 888 was never true

889 errors.append( 

890 f"{label}: PhotoCalib is None; flux-to-nJy calibration " 

891 "in alert packaging will fail." 

892 ) 

893 else: 

894 mean = photo_calib.getCalibrationMean() 

895 if not np.isfinite(mean) or mean <= 0.0: 895 ↛ 896line 895 didn't jump to line 896 because the condition on line 895 was never true

896 errors.append( 

897 f"{label}: PhotoCalib calibration mean is {mean!r}; " 

898 "must be a finite, positive value." 

899 ) 

900 

901 if exposure.getWcs() is None: 901 ↛ 902line 901 didn't jump to line 902 because the condition on line 901 was never true

902 errors.append( 

903 f"{label}: WCS is None; sky-to-pixel projection and " 

904 "DIA Object association will fail." 

905 ) 

906 

907 # Skip visitInfo checks for the template 

908 if label != "Template":