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

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 lsst.dax.apdb as daxApdb 

35import lsst.pex.config as pexConfig 

36import lsst.pipe.base as pipeBase 

37import lsst.pipe.base.connectionTypes as connTypes 

38import lsst.sphgeom 

39 

40from astropy.table import Table 

41import numpy as np 

42import pandas as pd 

43from lsst.ap.association import ( 

44 AssociationTask, 

45 DiaForcedSourceTask, 

46 PackageAlertsTask) 

47 

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

49from lsst.daf.base import DateTime 

50from lsst.meas.base import DetectorVisitIdGeneratorConfig, \ 

51 DiaObjectCalculationTask 

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

53 dropEmptyColumns, make_empty_catalog 

54from lsst.pipe.tasks.ssoAssociation import SolarSystemAssociationTask 

55from lsst.utils.timer import timeMethod, duration_from_timeMethod 

56 

57 

58class TooManyDiaObjectsError(pipeBase.AlgorithmError): 

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

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

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

62 """ 

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

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

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

66 f" configured maximum of {threshold}") 

67 super().__init__(msg) 

68 self.nNewDiaObjects = nNewDiaObjects 

69 self.threshold = threshold 

70 

71 @property 

72 def metadata(self): 

73 return {"nNewDiaObjects": self.nNewDiaObjects, 

74 "threshold": self.threshold 

75 } 

76 

77 

78class PostApdbUpdateError(pipeBase.AlgorithmError): 

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

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

81 can't be retried. 

82 """ 

83 def __init__(self, *, errorMsg): 

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

85 "This image cannot be retried! " 

86 "The original error was:\n" 

87 f"{errorMsg}") 

88 super().__init__(msg) 

89 

90 @property 

91 def metadata(self): 

92 return {} 

93 

94 

95class DiaPipelineConnections( 

96 pipeBase.PipelineTaskConnections, 

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

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

99 """Butler connections for DiaPipelineTask. 

100 """ 

101 diaSourceTable = connTypes.Input( 

102 doc="Catalog of calibrated DiaSources.", 

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

104 storageClass="DataFrame", 

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

106 ) 

107 solarSystemObjectTable = connTypes.Input( 

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

109 "this detectorVisit.", 

110 name="preloaded_SsObjects", 

111 storageClass="ArrowAstropy", 

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

113 minimum=0, 

114 ) 

115 diffIm = connTypes.Input( 

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

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

118 storageClass="ExposureF", 

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

120 ) 

121 exposure = connTypes.Input( 

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

123 "image differencing.", 

124 name="{fakesType}calexp", 

125 storageClass="ExposureF", 

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

127 ) 

128 template = connTypes.Input( 

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

130 "matched.", 

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

132 storageClass="ExposureF", 

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

134 ) 

135 preloadedDiaObjects = connTypes.Input( 

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

137 name="preloaded_diaObjects", 

138 storageClass="DataFrame", 

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

140 ) 

141 preloadedDiaSources = connTypes.Input( 

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

143 name="preloaded_diaSources", 

144 storageClass="DataFrame", 

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

146 ) 

147 preloadedDiaForcedSources = connTypes.Input( 

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

149 name="preloaded_diaForcedSources", 

150 storageClass="DataFrame", 

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

152 ) 

153 apdbMarker = connTypes.Output( 

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

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

156 name="apdb_marker", 

157 storageClass="Config", 

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

159 ) 

160 associatedDiaSources = connTypes.Output( 

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

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

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

164 storageClass="ArrowAstropy", 

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

166 ) 

167 associatedSsSources = connTypes.Output( 

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

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

170 storageClass="ArrowAstropy", 

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

172 ) 

173 unassociatedSsObjects = connTypes.Output( 

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

175 name="ssUnassociatedObjects", 

176 storageClass="ArrowAstropy", 

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

178 ) 

179 

180 diaForcedSources = connTypes.Output( 

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

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

183 storageClass="ArrowAstropy", 

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

185 ) 

186 diaObjects = connTypes.Output( 

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

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

189 storageClass="ArrowAstropy", 

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

191 ) 

192 newDiaSources = connTypes.Output( 

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

194 " were used to create a new diaObject", 

195 name="{fakesType}new_dia_source", 

196 storageClass="ArrowAstropy", 

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

198 ) 

199 marginalDiaSources = connTypes.Output( 

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

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

202 name="{fakesType}marginal_new_dia_source", 

203 storageClass="ArrowAstropy", 

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

205 ) 

206 

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

208 super().__init__(config=config) 

209 

210 if not config.doWriteAssociatedSources: 

211 self.outputs.remove("associatedDiaSources") 

212 self.outputs.remove("diaForcedSources") 

213 self.outputs.remove("diaObjects") 

214 self.outputs.remove("newDiaSources") 

215 self.outputs.remove("marginalDiaSources") 

216 else: 

217 if not config.doRunForcedMeasurement: 

218 self.outputs.remove("diaForcedSources") 

219 if not config.filterUnAssociatedSources: 

220 self.outputs.remove("newDiaSources") 

221 self.outputs.remove("marginalDiaSources") 

222 if not config.doSolarSystemAssociation: 

223 self.inputs.remove("solarSystemObjectTable") 

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

225 self.outputs.remove("associatedSsSources") 

226 self.outputs.remove("unassociatedSsObjects") 

227 

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

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

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

231 of the activator. 

232 

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

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

235 

236 Parameters 

237 ---------- 

238 inputs : `dict` 

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

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

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

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

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

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

245 dictionaries are guaranteed to be temporary copies that are true 

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

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

248 outputs : `dict` 

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

250 label : `str` 

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

252 diagnostic messages). 

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

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

255 diagnostic messages). 

256 

257 Returns 

258 ------- 

259 adjusted_inputs : `dict` 

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

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

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

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

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

265 `PipelineTask.runQuantum`. 

266 adjusted_outputs : `dict` 

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

268 interpretation as ``adjusted_inputs``. 

269 

270 Raises 

271 ------ 

272 ScalarError 

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

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

275 NoWorkFound 

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

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

278 quantum should be pruned or skipped. 

279 FileNotFoundError 

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

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

282 `PrerequisiteInput` connection. 

283 """ 

284 _, refs = inputs["diffIm"] 

285 for ref in refs: 

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

287 raise ValueError( 

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

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

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

291 "the validBands list in DiaPipelineConfig and add the " 

292 "appropriate columns to the Apdb schema.") 

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

294 

295 

296class DiaPipelineConfig(pipeBase.PipelineTaskConfig, 

297 pipelineConnections=DiaPipelineConnections): 

298 """Config for DiaPipelineTask. 

299 """ 

300 coaddName = pexConfig.Field( 

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

302 dtype=str, 

303 default="deep", 

304 ) 

305 apdb_config_url = pexConfig.Field( 

306 dtype=str, 

307 default=None, 

308 optional=False, 

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

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

311 "The database must already be initialized.", 

312 ) 

313 validBands = pexConfig.ListField( 

314 dtype=str, 

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

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

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

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

319 ) 

320 associator = pexConfig.ConfigurableField( 

321 target=AssociationTask, 

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

323 ) 

324 doSolarSystemAssociation = pexConfig.Field( 

325 dtype=bool, 

326 default=True, 

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

328 ) 

329 solarSystemAssociator = pexConfig.ConfigurableField( 

330 target=SolarSystemAssociationTask, 

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

332 ) 

333 diaCalculation = pexConfig.ConfigurableField( 

334 target=DiaObjectCalculationTask, 

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

336 ) 

337 doReloadDiaObjects = pexConfig.Field( 

338 dtype=bool, 

339 default=True, 

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

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

342 "are needed.", 

343 ) 

344 angleMargin = pexConfig.RangeField( 

345 doc="Padding to add when loading diaObjects if `doReloadDiaObjects=True`", 

346 dtype=float, 

347 default=2, 

348 min=0, 

349 ) 

350 doRunForcedMeasurement = pexConfig.Field( 

351 dtype=bool, 

352 default=True, 

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

354 "reasons during the ops rehearsal. " 

355 "It is expected to be removed.", 

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

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

358 ) 

359 diaForcedSource = pexConfig.ConfigurableField( 

360 target=DiaForcedSourceTask, 

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

362 "difference images.", 

363 ) 

364 forcedReliabilityThreshold = pexConfig.Field( 

365 dtype=float, 

366 default=0.5, 

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

368 "determining whether to calculate forced photometry for a " 

369 "diaObject.", 

370 ) 

371 forcedTrailLengthThreshold = pexConfig.Field( 

372 dtype=float, 

373 default=1.4, 

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

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

376 "diaObject.", 

377 ) 

378 forcedBadFlags = pexConfig.ListField( 

379 dtype=str, 

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

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

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

383 ) 

384 alertPackager = pexConfig.ConfigurableField( 

385 target=PackageAlertsTask, 

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

387 ) 

388 doPackageAlerts = pexConfig.Field( 

389 dtype=bool, 

390 default=False, 

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

392 "write them to disk.", 

393 ) 

394 doWriteAssociatedSources = pexConfig.Field( 

395 dtype=bool, 

396 default=True, 

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

398 "formatted following the Science Data Model.", 

399 ) 

400 imagePixelMargin = pexConfig.RangeField( 

401 dtype=int, 

402 default=10, 

403 min=0, 

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

405 "diaObjects for association.", 

406 ) 

407 filterUnAssociatedSources = pexConfig.Field( 

408 dtype=bool, 

409 default=True, 

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

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

412 "will not be affected." 

413 ) 

414 newObjectBadFlags = pexConfig.ListField( 

415 dtype=str, 

416 default=("centroid_flag", 

417 "pixelFlags_crCenter", 

418 "pixelFlags_nodataCenter", 

419 "pixelFlags_interpolatedCenter", 

420 "pixelFlags_saturatedCenter", 

421 "pixelFlags_suspectCenter", 

422 "pixelFlags_streakCenter", 

423 "glint_trail"), 

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

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

426 ) 

427 maxNewDiaObjects = pexConfig.RangeField( 

428 dtype=float, 

429 default=0, 

430 min=0, 

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

432 "Set to zero to disable.", 

433 ) 

434 newObjectSnrThreshold = pexConfig.RangeField( 

435 dtype=float, 

436 default=0, 

437 min=0, 

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

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

440 "diaObjects." 

441 "Set to zero to disable.", 

442 ) 

443 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

444 dtype=float, 

445 default=0, 

446 min=0, 

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

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

449 " low reliability scores before creating new diaObjects." 

450 "Set to zero to disable.", 

451 ) 

452 newObjectReliabilityThreshold = pexConfig.RangeField( 

453 dtype=float, 

454 default=0.1, 

455 min=0, 

456 max=1, 

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

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

459 "diaObjects." 

460 "Set to zero to disable.", 

461 ) 

462 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

463 dtype=float, 

464 default=0.1, 

465 min=0, 

466 max=1, 

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

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

469 "before creating new diaObjects." 

470 "Set to zero to disable.", 

471 ) 

472 newObjectFluxField = pexConfig.Field( 

473 dtype=str, 

474 default="apFlux", 

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

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

477 ) 

478 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

479 

480 def setDefaults(self): 

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

482 "ap_nDiaSources", 

483 "ap_meanFlux", 

484 "ap_sigmaFlux", 

485 "ap_minMaxFlux", 

486 "ap_maxSlopeFlux", 

487 "ap_meanErrFlux", 

488 "ap_meanTotFlux"] 

489 

490 

491class DiaPipelineTask(pipeBase.PipelineTask): 

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

493 (DIA) Objects and Sources. 

494 """ 

495 ConfigClass = DiaPipelineConfig 

496 _DefaultName = "diaPipe" 

497 

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

499 super().__init__(**kwargs) 

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

501 self.schema = readSchemaFromApdb(self.apdb) 

502 self.makeSubtask("associator") 

503 self.makeSubtask("diaCalculation") 

504 if self.config.doRunForcedMeasurement: 

505 self.makeSubtask("diaForcedSource") 

506 if self.config.doPackageAlerts: 

507 self.makeSubtask("alertPackager") 

508 if self.config.doSolarSystemAssociation: 

509 self.makeSubtask("solarSystemAssociator") 

510 if self.config.filterUnAssociatedSources: 

511 columns = [self.config.newObjectFluxField, 

512 self.config.newObjectFluxField + "Err", 

513 "reliability"] 

514 columns += self.config.newObjectBadFlags 

515 

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

517 if missing: 

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

519 

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

521 inputs = butlerQC.get(inputRefs) 

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

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

524 inputs["legacySolarSystemTable"] = None 

525 if not self.config.doSolarSystemAssociation: 

526 inputs["solarSystemObjectTable"] = None 

527 

528 associationResults = pipeBase.Struct( 

529 apdbMarker=None, 

530 associatedDiaSources=None, 

531 diaForcedSources=None, 

532 diaObjects=None, 

533 associatedSsSources=None, 

534 unassociatedSsObjects=None, 

535 newDiaSources=None, 

536 marginalDiaSources=None, 

537 ) 

538 try: 

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

540 except pipeBase.AlgorithmError as e: 

541 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

542 e, 

543 self, 

544 log=self.log 

545 ) 

546 butlerQC.put(associationResults, outputRefs) 

547 raise error from e 

548 butlerQC.put(associationResults, outputRefs) 

549 

550 @timeMethod 

551 def run(self, 

552 diaSourceTable, 

553 legacySolarSystemTable, 

554 diffIm, 

555 exposure, 

556 template, 

557 preloadedDiaObjects, 

558 preloadedDiaSources, 

559 preloadedDiaForcedSources, 

560 band, 

561 idGenerator, 

562 solarSystemObjectTable=None, 

563 associationResults=None): 

564 """Process DiaSources and DiaObjects. 

565 

566 Load previous DiaObjects and their DiaSource history. Calibrate the 

567 values in the diaSourceCat. Associate new DiaSources with previous 

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

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

570 

571 Parameters 

572 ---------- 

573 diaSourceTable : `pandas.DataFrame` 

574 Newly detected DiaSources. 

575 legacySolarSystemTable : `pandas.DataFrame` 

576 Not used 

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

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

579 were detected. 

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

581 Calibrated exposure differenced with a template to create 

582 ``diffIm``. 

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

584 Template exposure used to create diffIm. 

585 preloadedDiaObjects : `pandas.DataFrame` 

586 Previously detected DiaObjects, loaded from the APDB. 

587 preloadedDiaSources : `pandas.DataFrame` 

588 Previously detected DiaSources, loaded from the APDB. 

589 preloadedDiaForcedSources : `pandas.DataFrame` 

590 Catalog of previously detected forced DiaSources, from the APDB 

591 band : `str` 

592 The band in which the new DiaSources were detected. 

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

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

595 solarSystemObjectTable : `astropy.table.Table` 

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

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

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

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

600 this is also returned. 

601 

602 Returns 

603 ------- 

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

605 Results struct with components. 

606 

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

608 that this ccdVisit has completed successfully. 

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

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

611 DiaSources. (`pandas.DataFrame`) 

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

613 forced DiaSources. (`pandas.DataFrame`) 

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

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

616 (`pandas.DataFrame`) 

617 

618 Raises 

619 ------ 

620 RuntimeError 

621 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

622 """ 

623 if associationResults is None: 

624 associationResults = pipeBase.Struct() 

625 self._validateExposure(exposure, "Science") 

626 self._validateExposure(diffIm, "Difference") 

627 self._validateExposure(template, "Template") 

628 

629 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

630 if legacySolarSystemTable is not None and solarSystemObjectTable is None: 

631 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

632 if self.config.doReloadDiaObjects: 

633 try: 

634 preloadedDiaObjects = self.loadRefreshedDiaObjects(getRegion(exposure), preloadedDiaObjects) 

635 except Exception as e: 

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

637 "the latest diaObjects from the APDB. Processing " 

638 "will continue with only the diaObjects from " 

639 "preload.", exc_info=e) 

640 finally: 

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

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

643 ) 

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

645 self.metadata["loadDiaObjectsDuration"]) 

646 

647 else: 

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

649 

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

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

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

653 

654 if not preloadedDiaObjects.empty: 

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

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

657 buffer=self.config.imagePixelMargin) 

658 else: 

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

660 diaObjects = preloadedDiaObjects 

661 

662 # Associate DiaSources with DiaObjects 

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

664 

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

666 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

667 assocResults.associatedDiaSources, 

668 "DiaSource", 

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

670 ) 

671 

672 # Merge associated diaSources 

673 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

674 preloadedDiaSources, 

675 assocResults.associatedDiaSources, 

676 diaObjects, 

677 assocResults.newDiaObjects, 

678 diffIm 

679 ) 

680 

681 # Compute DiaObject Summary statistics from their full DiaSource 

682 # history. 

683 diaCalResult = self.diaCalculation.run( 

684 mergedDiaObjects, 

685 mergedDiaSourceHistory, 

686 updatedDiaObjectIds, 

687 self.config.validBands) 

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

689 tableName="DiaObject", skipIndex=True) 

690 

691 # Test for duplication in the updated DiaObjects. 

692 if self.testDataFrameIndex(diaCalResult.diaObjectCat): 

693 raise RuntimeError( 

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

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

696 "reported. Exiting.") 

697 if self.testDataFrameIndex(updatedDiaObjects): 

698 raise RuntimeError( 

699 "Duplicate DiaObjects (updated) created after " 

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

701 "reported. Exiting.") 

702 

703 # Forced source measurement 

704 if self.config.doRunForcedMeasurement: 

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

706 diaForcedSources = self.runForcedMeasurement( 

707 diaObjectsForced, updatedDiaObjects, exposure, diffIm, idGenerator 

708 ) 

709 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

710 else: 

711 # alertPackager needs correct columns 

712 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

713 forcedSourceHistoryThreshold = 0 

714 

715 # Write results to Alert Production Database (APDB) 

716 try: 

717 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

718 diaForcedSources) 

719 finally: 

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

721 clock="Utc") 

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

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

724 

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

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

727 associationResults.apdbMarker = pexConfig.Config() 

728 associationResults.associatedDiaSources = assocResults.associatedDiaSources 

729 associationResults.diaForcedSources = diaForcedSources 

730 associationResults.diaObjects = diaCalResult.diaObjectCat 

731 associationResults.unassociatedSsObjects = assocResults.unassociatedSsObjects 

732 associationResults.newDiaSources = assocResults.newDiaSources 

733 associationResults.marginalDiaSources = assocResults.marginalDiaSources 

734 

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

736 try: 

737 associatedSsSources = assocResults.associatedSsSources 

738 associatedSsSourcesPlusMpcorb = None 

739 if associatedSsSources is not None: 

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

741 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

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

743 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

744 for mpcorbColumn in mpcorbColumns: 

745 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

746 associationResults.associatedSsSources = associatedSsSources 

747 

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

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

750 scale=DateTime.TAI) 

751 

752 # Package alerts 

753 if self.config.doPackageAlerts: 

754 # Append new forced sources to the full history 

755 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

756 tableName="DiaForcedSource") 

757 if self.testDataFrameIndex(diaForcedSourcesFull): 

758 self.log.warning( 

759 "Duplicate DiaForcedSources created after merge with " 

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

761 "problems. Dropping duplicates.") 

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

763 diaForcedSourcesFull = diaForcedSourcesFull[ 

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

765 

766 self.alertPackager.run(assocResults.associatedDiaSources, 

767 updatedDiaObjects, 

768 preloadedDiaSources, 

769 diaForcedSourcesFull, 

770 diffIm, 

771 exposure, 

772 template, 

773 ssSrc=associatedSsSourcesPlusMpcorb, 

774 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

775 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

776 ) 

777 except Exception as e: 

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

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

780 

781 return associationResults 

782 

783 def _validateExposure(self, exposure, label): 

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

785 the APDB. 

786 

787 Parameters 

788 ---------- 

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

790 Calibrated science image that was subtracted. Corresponds to the 

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

792 label : `str` 

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

794 

795 Raises 

796 ------ 

797 ValueError 

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

799 """ 

800 errors = [] 

801 

802 photo_calib = exposure.getPhotoCalib() 

803 if photo_calib is None: 

804 errors.append( 

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

806 "in alert packaging will fail." 

807 ) 

808 else: 

809 mean = photo_calib.getCalibrationMean() 

810 if not np.isfinite(mean) or mean <= 0.0: 

811 errors.append( 

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

813 "must be a finite, positive value." 

814 ) 

815 

816 if exposure.getWcs() is None: 

817 errors.append( 

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

819 "DIA Object association will fail." 

820 ) 

821 

822 # Skip visitInfo checks for the template 

823 if label != "Template": 

824 visit_info = exposure.getInfo().getVisitInfo() 

825 if visit_info is None: 

826 errors.append( 

827 f"{label}: VisitInfo is None; observation time, pointing, " 

828 "rotation angle, and exposure duration are unavailable." 

829 ) 

830 else: 

831 obs_date = visit_info.getDate() 

832 if not obs_date.isValid(): 

833 errors.append( 

834 f"{label}: VisitInfo.date is invalid." 

835 ) 

836 

837 # boresightRaDec — written to APDB rows and alert payloads. 

838 boresight = visit_info.getBoresightRaDec() 

839 ra_deg = boresight.getRa().asDegrees() 

840 dec_deg = boresight.getDec().asDegrees() 

841 if not np.isfinite(ra_deg) or not np.isfinite(dec_deg): 

842 errors.append( 

843 f"{label}: VisitInfo.boresightRaDec contains NaN " 

844 f"(RA={ra_deg}, Dec={dec_deg}); pointing metadata " 

845 "is missing or corrupt." 

846 ) 

847 

848 rot_deg = visit_info.getBoresightRotAngle().asDegrees() 

849 if not np.isfinite(rot_deg): 

850 errors.append( 

851 f"{label}: VisitInfo.boresightRotAngle is NaN; the " 

852 "ROTPA keyword will be corrupt in all alert cutout " 

853 "FITS headers." 

854 ) 

855 

856 filter_label = exposure.getFilter() 

857 if filter_label is None: 

858 errors.append( 

859 f"{label}: filter label is None." 

860 ) 

861 else: 

862 band = filter_label.bandLabel if filter_label.hasBandLabel() else "" 

863 if not band: 

864 errors.append( 

865 f"{label}: filter label carries no band name." 

866 ) 

867 

868 bbox = exposure.getBBox() 

869 if bbox.getWidth() <= 0 or bbox.getHeight() <= 0: 

870 errors.append( 

871 f"{label}: bounding box is empty " 

872 f"({bbox.getWidth()}×{bbox.getHeight()} px); the " 

873 "readout or ISR step may have failed." 

874 ) 

875 

876 if exposure.getPsf() is None: 

877 errors.append( 

878 f"{label}: PSF model is None." 

879 ) 

880 

881 if errors: 

882 detail = "\n ".join(errors) 

883 raise ValueError( 

884 "Input exposure metadata is invalid; aborting before any " 

885 "APDB writes to prevent non-recoverable data corruption:\n " 

886 + detail 

887 ) 

888 

889 def _selectGoodDiaObjects(self, diaObjects, diaSources): 

890 """ 

891 Extract a subset of diaObjects with multiple associated diaSources that 

892 pass selection requirements. 

893 """ 

894 n_history = self.diaForcedSource.config.historyThreshold 

895 required_columns = ["reliability", "trailLength"] 

896 for col in required_columns: 

897 if col not in diaSources.columns: 

898 raise RuntimeError("Required column '%s' missing from diaSource table." % col) 

899 

900 rel_ok = diaSources["reliability"] >= self.config.forcedReliabilityThreshold 

901 trailed_ok = diaSources["trailLength"] <= self.config.forcedTrailLengthThreshold 

902 badFlags = [f for f in self.config.forcedBadFlags if f in diaSources.columns] 

903 if badFlags: 

904 self.log.info("Excluding diaSources with %s flags set when calculating" 

905 " diaObject history for forced photometry." % badFlags) 

906 # Exclude rows where ANY bad flag is True. 

907 # If bad flag columns can contain NA, treat NA as True (i.e., bad) via fillna(True). 

908 bad_any = diaSources[badFlags].fillna(True).any(axis=1)