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

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.loadDiaCatalogs import loadDiaObjectsFromApdb, loadDiaSourcesFromApdb, \ 

49 loadDiaForcedSourcesFromApdb 

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

51from lsst.daf.base import DateTime 

52from lsst.meas.base import DetectorVisitIdGeneratorConfig, \ 

53 DiaObjectCalculationTask 

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

55 dropEmptyColumns, make_empty_catalog 

56from lsst.pipe.tasks.ssoAssociation import SolarSystemAssociationTask 

57from lsst.utils.timer import timeMethod, duration_from_timeMethod 

58 

59 

60class TooManyDiaObjectsError(pipeBase.AlgorithmError): 

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

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

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

64 """ 

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

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

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

68 f" configured maximum of {threshold}") 

69 super().__init__(msg) 

70 self.nNewDiaObjects = nNewDiaObjects 

71 self.threshold = threshold 

72 

73 @property 

74 def metadata(self): 

75 return {"nNewDiaObjects": self.nNewDiaObjects, 

76 "threshold": self.threshold 

77 } 

78 

79 

80class PostApdbUpdateError(pipeBase.AlgorithmError): 

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

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

83 can't be retried. 

84 """ 

85 def __init__(self, *, errorMsg): 

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

87 "This image cannot be retried! " 

88 "The original error was:\n" 

89 f"{errorMsg}") 

90 super().__init__(msg) 

91 

92 @property 

93 def metadata(self): 

94 return {} 

95 

96 

97class DiaPipelineConnections( 

98 pipeBase.PipelineTaskConnections, 

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

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

101 """Butler connections for DiaPipelineTask. 

102 """ 

103 diaSourceTable = connTypes.Input( 

104 doc="Catalog of calibrated DiaSources.", 

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

106 storageClass="DataFrame", 

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

108 ) 

109 solarSystemObjectTable = connTypes.Input( 

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

111 "this detectorVisit.", 

112 name="preloaded_SsObjects", 

113 storageClass="ArrowAstropy", 

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

115 minimum=0, 

116 ) 

117 diffIm = connTypes.Input( 

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

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

120 storageClass="ExposureF", 

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

122 ) 

123 exposure = connTypes.Input( 

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

125 "image differencing.", 

126 name="{fakesType}calexp", 

127 storageClass="ExposureF", 

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

129 ) 

130 template = connTypes.Input( 

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

132 "matched.", 

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

134 storageClass="ExposureF", 

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

136 ) 

137 preloadedDiaObjects = connTypes.Input( 

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

139 name="preloaded_diaObjects", 

140 storageClass="DataFrame", 

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

142 ) 

143 preloadedDiaSources = connTypes.Input( 

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

145 name="preloaded_diaSources", 

146 storageClass="DataFrame", 

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

148 ) 

149 preloadedDiaForcedSources = connTypes.Input( 

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

151 name="preloaded_diaForcedSources", 

152 storageClass="DataFrame", 

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

154 ) 

155 apdbMarker = connTypes.Output( 

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

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

158 name="apdb_marker", 

159 storageClass="Config", 

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

161 ) 

162 associatedDiaSources = connTypes.Output( 

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

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

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

166 storageClass="ArrowAstropy", 

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

168 ) 

169 associatedSsSources = connTypes.Output( 

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

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

172 storageClass="ArrowAstropy", 

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

174 ) 

175 unassociatedSsObjects = connTypes.Output( 

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

177 name="ssUnassociatedObjects", 

178 storageClass="ArrowAstropy", 

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

180 ) 

181 

182 diaForcedSources = connTypes.Output( 

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

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

185 storageClass="ArrowAstropy", 

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

187 ) 

188 diaObjects = connTypes.Output( 

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

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

191 storageClass="ArrowAstropy", 

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

193 ) 

194 newDiaSources = connTypes.Output( 

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

196 " were used to create a new diaObject", 

197 name="{fakesType}new_dia_source", 

198 storageClass="ArrowAstropy", 

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

200 ) 

201 marginalDiaSources = connTypes.Output( 

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

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

204 name="{fakesType}marginal_new_dia_source", 

205 storageClass="ArrowAstropy", 

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

207 ) 

208 

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

210 super().__init__(config=config) 

211 

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

213 self.outputs.remove("associatedDiaSources") 

214 self.outputs.remove("diaForcedSources") 

215 self.outputs.remove("diaObjects") 

216 self.outputs.remove("newDiaSources") 

217 self.outputs.remove("marginalDiaSources") 

218 else: 

219 if not config.doRunForcedMeasurement: 

220 self.outputs.remove("diaForcedSources") 

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

222 self.outputs.remove("newDiaSources") 

223 self.outputs.remove("marginalDiaSources") 

224 if not config.doSolarSystemAssociation: 

225 self.inputs.remove("solarSystemObjectTable") 

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

227 self.outputs.remove("associatedSsSources") 

228 self.outputs.remove("unassociatedSsObjects") 

229 if config.doReloadAllApdbCatalogs: 

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

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

232 self.inputs.remove("preloadedDiaObjects") 

233 self.inputs.remove("preloadedDiaSources") 

234 self.inputs.remove("preloadedDiaForcedSources") 

235 

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

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

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

239 of the activator. 

240 

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

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

243 

244 Parameters 

245 ---------- 

246 inputs : `dict` 

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

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

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

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

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

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

253 dictionaries are guaranteed to be temporary copies that are true 

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

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

256 outputs : `dict` 

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

258 label : `str` 

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

260 diagnostic messages). 

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

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

263 diagnostic messages). 

264 

265 Returns 

266 ------- 

267 adjusted_inputs : `dict` 

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

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

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

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

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

273 `PipelineTask.runQuantum`. 

274 adjusted_outputs : `dict` 

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

276 interpretation as ``adjusted_inputs``. 

277 

278 Raises 

279 ------ 

280 ScalarError 

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

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

283 NoWorkFound 

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

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

286 quantum should be pruned or skipped. 

287 FileNotFoundError 

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

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

290 `PrerequisiteInput` connection. 

291 """ 

292 _, refs = inputs["diffIm"] 

293 for ref in refs: 

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

295 raise ValueError( 

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

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

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

299 "the validBands list in DiaPipelineConfig and add the " 

300 "appropriate columns to the Apdb schema.") 

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

302 

303 

304class DiaPipelineConfig(pipeBase.PipelineTaskConfig, 

305 pipelineConnections=DiaPipelineConnections): 

306 """Config for DiaPipelineTask. 

307 """ 

308 coaddName = pexConfig.Field( 

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

310 dtype=str, 

311 default="deep", 

312 ) 

313 apdb_config_url = pexConfig.Field( 

314 dtype=str, 

315 default=None, 

316 optional=False, 

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

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

319 "The database must already be initialized.", 

320 ) 

321 validBands = pexConfig.ListField( 

322 dtype=str, 

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

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

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

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

327 ) 

328 associator = pexConfig.ConfigurableField( 

329 target=AssociationTask, 

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

331 ) 

332 doSolarSystemAssociation = pexConfig.Field( 

333 dtype=bool, 

334 default=True, 

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

336 ) 

337 solarSystemAssociator = pexConfig.ConfigurableField( 

338 target=SolarSystemAssociationTask, 

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

340 ) 

341 diaCalculation = pexConfig.ConfigurableField( 

342 target=DiaObjectCalculationTask, 

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

344 ) 

345 doReloadDiaObjects = pexConfig.Field( 

346 dtype=bool, 

347 default=True, 

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

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

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

351 "is set.", 

352 ) 

353 doReloadAllApdbCatalogs = pexConfig.Field( 

354 dtype=bool, 

355 default=False, 

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

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

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

359 ) 

360 angleMargin = pexConfig.RangeField( 

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

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

363 "this task is after astrometric calibration.", 

364 dtype=float, 

365 default=2, 

366 min=0, 

367 ) 

368 doRunForcedMeasurement = pexConfig.Field( 

369 dtype=bool, 

370 default=True, 

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

372 "reasons during the ops rehearsal. " 

373 "It is expected to be removed.", 

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

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

376 ) 

377 diaForcedSource = pexConfig.ConfigurableField( 

378 target=DiaForcedSourceTask, 

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

380 "difference images.", 

381 ) 

382 forcedReliabilityThreshold = pexConfig.Field( 

383 dtype=float, 

384 default=0.5, 

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

386 "determining whether to calculate forced photometry for a " 

387 "diaObject.", 

388 ) 

389 forcedTrailLengthThreshold = pexConfig.Field( 

390 dtype=float, 

391 default=1.4, 

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

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

394 "diaObject.", 

395 ) 

396 forcedBadFlags = pexConfig.ListField( 

397 dtype=str, 

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

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

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

401 ) 

402 alertPackager = pexConfig.ConfigurableField( 

403 target=PackageAlertsTask, 

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

405 ) 

406 doPackageAlerts = pexConfig.Field( 

407 dtype=bool, 

408 default=False, 

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

410 "write them to disk.", 

411 ) 

412 doWriteAssociatedSources = pexConfig.Field( 

413 dtype=bool, 

414 default=True, 

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

416 "formatted following the Science Data Model.", 

417 ) 

418 imagePixelMargin = pexConfig.RangeField( 

419 dtype=int, 

420 default=10, 

421 min=0, 

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

423 "diaObjects for association.", 

424 ) 

425 filterUnAssociatedSources = pexConfig.Field( 

426 dtype=bool, 

427 default=True, 

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

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

430 "will not be affected." 

431 ) 

432 newObjectBadFlags = pexConfig.ListField( 

433 dtype=str, 

434 default=("centroid_flag", 

435 "pixelFlags_crCenter", 

436 "pixelFlags_nodataCenter", 

437 "pixelFlags_interpolatedCenter", 

438 "pixelFlags_saturatedCenter", 

439 "pixelFlags_suspectCenter", 

440 "pixelFlags_streakCenter", 

441 "glint_trail"), 

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

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

444 ) 

445 maxNewDiaObjects = pexConfig.RangeField( 

446 dtype=float, 

447 default=0, 

448 min=0, 

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

450 "Set to zero to disable.", 

451 ) 

452 newObjectSnrThreshold = pexConfig.RangeField( 

453 dtype=float, 

454 default=0, 

455 min=0, 

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

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

458 "diaObjects." 

459 "Set to zero to disable.", 

460 ) 

461 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

462 dtype=float, 

463 default=0, 

464 min=0, 

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

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

467 " low reliability scores before creating new diaObjects." 

468 "Set to zero to disable.", 

469 ) 

470 newObjectReliabilityThreshold = pexConfig.RangeField( 

471 dtype=float, 

472 default=0.1, 

473 min=0, 

474 max=1, 

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

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

477 "diaObjects." 

478 "Set to zero to disable.", 

479 ) 

480 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

481 dtype=float, 

482 default=0.1, 

483 min=0, 

484 max=1, 

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

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

487 "before creating new diaObjects." 

488 "Set to zero to disable.", 

489 ) 

490 newObjectFluxField = pexConfig.Field( 

491 dtype=str, 

492 default="apFlux", 

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

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

495 ) 

496 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

497 

498 def setDefaults(self): 

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

500 "ap_nDiaSources", 

501 "ap_meanFlux", 

502 "ap_sigmaFlux", 

503 "ap_minMaxFlux", 

504 "ap_maxSlopeFlux", 

505 "ap_meanErrFlux", 

506 "ap_meanTotFlux"] 

507 

508 

509class DiaPipelineTask(pipeBase.PipelineTask): 

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

511 (DIA) Objects and Sources. 

512 """ 

513 ConfigClass = DiaPipelineConfig 

514 _DefaultName = "diaPipe" 

515 

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

517 super().__init__(**kwargs) 

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

519 self.schema = readSchemaFromApdb(self.apdb) 

520 self.makeSubtask("associator") 

521 self.makeSubtask("diaCalculation") 

522 if self.config.doRunForcedMeasurement: 

523 self.makeSubtask("diaForcedSource") 

524 if self.config.doPackageAlerts: 

525 self.makeSubtask("alertPackager") 

526 if self.config.doSolarSystemAssociation: 

527 self.makeSubtask("solarSystemAssociator") 

528 if self.config.filterUnAssociatedSources: 

529 columns = [self.config.newObjectFluxField, 

530 self.config.newObjectFluxField + "Err", 

531 "reliability"] 

532 columns += self.config.newObjectBadFlags 

533 

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

535 if missing: 

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

537 

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

539 inputs = butlerQC.get(inputRefs) 

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

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

542 inputs["legacySolarSystemTable"] = None 

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

544 inputs["solarSystemObjectTable"] = None 

545 

546 associationResults = pipeBase.Struct( 

547 apdbMarker=None, 

548 associatedDiaSources=None, 

549 diaForcedSources=None, 

550 diaObjects=None, 

551 associatedSsSources=None, 

552 unassociatedSsObjects=None, 

553 newDiaSources=None, 

554 marginalDiaSources=None, 

555 ) 

556 try: 

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

558 except pipeBase.AlgorithmError as e: 

559 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

560 e, 

561 self, 

562 log=self.log 

563 ) 

564 butlerQC.put(associationResults, outputRefs) 

565 raise error from e 

566 butlerQC.put(associationResults, outputRefs) 

567 

568 @timeMethod 

569 def run(self, 

570 diaSourceTable, 

571 legacySolarSystemTable, 

572 diffIm, 

573 exposure, 

574 template, 

575 preloadedDiaObjects=None, 

576 preloadedDiaSources=None, 

577 preloadedDiaForcedSources=None, 

578 band=None, 

579 idGenerator=None, 

580 solarSystemObjectTable=None, 

581 associationResults=None): 

582 """Process DiaSources and DiaObjects. 

583 

584 Load previous DiaObjects and their DiaSource history. Calibrate the 

585 values in the diaSourceCat. Associate new DiaSources with previous 

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

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

588 

589 Parameters 

590 ---------- 

591 diaSourceTable : `pandas.DataFrame` 

592 Newly detected DiaSources. 

593 legacySolarSystemTable : `pandas.DataFrame` 

594 Not used 

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

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

597 were detected. 

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

599 Calibrated exposure differenced with a template to create 

600 ``diffIm``. 

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

602 Template exposure used to create diffIm. 

603 preloadedDiaObjects : `pandas.DataFrame`, optional 

604 Previously detected DiaObjects, loaded from the APDB. 

605 `None` if ``doReloadAllApdbCatalogs`` is set. 

606 preloadedDiaSources : `pandas.DataFrame`, optional 

607 Previously detected DiaSources, loaded from the APDB. 

608 `None` if ``doReloadAllApdbCatalogs`` is set. 

609 preloadedDiaForcedSources : `pandas.DataFrame`, optional 

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

611 `None` if ``doReloadAllApdbCatalogs`` is set. 

612 band : `str` 

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

614 despite the `None` default. 

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

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

617 Required, despite the `None` default. 

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

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

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

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

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

623 this is also returned. 

624 

625 Returns 

626 ------- 

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

628 Results struct with components. 

629 

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

631 that this ccdVisit has completed successfully. 

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

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

634 DiaSources. (`pandas.DataFrame`) 

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

636 forced DiaSources. (`pandas.DataFrame`) 

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

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

639 (`pandas.DataFrame`) 

640 

641 Raises 

642 ------ 

643 RuntimeError 

644 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

645 ValueError 

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

647 """ 

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

649 # precede them, can be omitted. 

650 if band is None: 

651 raise ValueError("band is required.") 

652 if idGenerator is None: 

653 raise ValueError("idGenerator is required.") 

654 

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

656 associationResults = pipeBase.Struct() 

657 self._validateExposure(exposure, "Science") 

658 self._validateExposure(diffIm, "Difference") 

659 self._validateExposure(template, "Template") 

660 

661 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

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

663 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

664 region = getRegion(exposure) 

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

666 try: 

667 preloadedDiaObjects = self.loadRefreshedDiaObjects(region, preloadedDiaObjects) 

668 except Exception as e: 

669 if self.config.doReloadAllApdbCatalogs: 

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

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

672 # no preloaded diaObjects to fall back on. 

673 raise 

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

675 "the latest diaObjects from the APDB. Processing " 

676 "will continue with only the diaObjects from " 

677 "preload.", exc_info=e) 

678 finally: 

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

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

681 ) 

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

683 self.metadata["loadDiaObjectsDuration"]) 

684 

685 else: 

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

687 

688 if self.config.doReloadAllApdbCatalogs: 

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

690 # LoadDiaCatalogsTask to keep the metrics upload consistent. 

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

692 try: 

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

694 finally: 

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

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

697 ) 

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

699 self.metadata["loadDiaSourcesDuration"]) 

700 try: 

701 preloadedDiaForcedSources = self.loadRefreshedDiaForcedSources( 

702 region, preloadedDiaObjects, visitTime) 

703 finally: 

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

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

706 ) 

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

708 self.metadata["loadDiaForcedSourcesDuration"]) 

709 else: 

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

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

712 

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

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

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

716 

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

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

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

720 buffer=self.config.imagePixelMargin) 

721 else: 

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

723 diaObjects = preloadedDiaObjects 

724 

725 # Associate DiaSources with DiaObjects 

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

727 

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

729 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

730 assocResults.associatedDiaSources, 

731 "DiaSource", 

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

733 ) 

734 

735 # Merge associated diaSources 

736 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

737 preloadedDiaSources, 

738 assocResults.associatedDiaSources, 

739 diaObjects, 

740 assocResults.newDiaObjects, 

741 diffIm 

742 ) 

743 

744 # Compute DiaObject Summary statistics from their full DiaSource 

745 # history. 

746 diaCalResult = self.diaCalculation.run( 

747 mergedDiaObjects, 

748 mergedDiaSourceHistory, 

749 updatedDiaObjectIds, 

750 self.config.validBands) 

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

752 tableName="DiaObject", skipIndex=True) 

753 

754 # Test for duplication in the updated DiaObjects. 

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

756 raise RuntimeError( 

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

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

759 "reported. Exiting.") 

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

761 raise RuntimeError( 

762 "Duplicate DiaObjects (updated) created after " 

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

764 "reported. Exiting.") 

765 

766 # Forced source measurement 

767 if self.config.doRunForcedMeasurement: 

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

769 diaForcedSources = self.runForcedMeasurement( 

770 diaObjectsForced, updatedDiaObjects, exposure, diffIm, idGenerator 

771 ) 

772 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

773 else: 

774 # alertPackager needs correct columns 

775 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

776 forcedSourceHistoryThreshold = 0 

777 

778 # Write results to Alert Production Database (APDB) 

779 try: 

780 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

781 diaForcedSources) 

782 finally: 

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

784 clock="Utc") 

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

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

787 

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

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

790 associationResults.apdbMarker = pexConfig.Config() 

791 associationResults.associatedDiaSources = assocResults.associatedDiaSources 

792 associationResults.diaForcedSources = diaForcedSources 

793 associationResults.diaObjects = diaCalResult.diaObjectCat 

794 associationResults.unassociatedSsObjects = assocResults.unassociatedSsObjects 

795 associationResults.newDiaSources = assocResults.newDiaSources 

796 associationResults.marginalDiaSources = assocResults.marginalDiaSources 

797 

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

799 try: 

800 associatedSsSources = assocResults.associatedSsSources 

801 associatedSsSourcesPlusMpcorb = None 

802 if associatedSsSources is not None: 

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

804 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

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

806 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

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

808 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

809 associationResults.associatedSsSources = associatedSsSources 

810 

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

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

813 scale=DateTime.TAI) 

814 

815 # Package alerts 

816 if self.config.doPackageAlerts: 

817 # Append new forced sources to the full history 

818 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

819 tableName="DiaForcedSource") 

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

821 self.log.warning( 

822 "Duplicate DiaForcedSources created after merge with " 

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

824 "problems. Dropping duplicates.") 

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

826 diaForcedSourcesFull = diaForcedSourcesFull[ 

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

828 

829 self.alertPackager.run(assocResults.associatedDiaSources, 

830 updatedDiaObjects, 

831 preloadedDiaSources, 

832 diaForcedSourcesFull, 

833 diffIm, 

834 exposure, 

835 template, 

836 ssSrc=associatedSsSourcesPlusMpcorb, 

837 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

838 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

839 ) 

840 except Exception as e: 

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

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

843 

844 return associationResults 

845 

846 def _validateExposure(self, exposure, label): 

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

848 the APDB. 

849 

850 Parameters 

851 ---------- 

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

853 Calibrated science image that was subtracted. Corresponds to the 

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

855 label : `str` 

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

857 

858 Raises 

859 ------ 

860 ValueError 

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

862 """ 

863 errors = [] 

864 

865 photo_calib = exposure.getPhotoCalib() 

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

867 errors.append( 

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

869 "in alert packaging will fail." 

870 ) 

871 else: 

872 mean = photo_calib.getCalibrationMean() 

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

874 errors.append( 

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

876 "must be a finite, positive value." 

877 ) 

878 

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

880 errors.append( 

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

882 "DIA Object association will fail." 

883 ) 

884 

885 # Skip visitInfo checks for the template 

886 if label != "Template": 

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

888 if visit_info 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}: VisitInfo is None; observation time, pointing, " 

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

892 ) 

893 else: 

894 obs_date = visit_info.getDate() 

895 if not obs_date.isValid(): 895 ↛ 896line 895 didn't jump to line 896 because the condition on line 895 was never true

896 errors.append( 

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

898 ) 

899 

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

901 boresight = visit_info.getBoresightRaDec() 

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

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

904 if not np.isfinite(ra_deg) or not np.isfinite(dec_deg): 904 ↛ 905line 904 didn't jump to line 905 because the condition on line 904 was never true

905 errors.append( 

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

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

908 "is missing or corrupt."