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_firstLastDiaSourceTime", 

520 "ap_nDiaSources", 

521 "ap_meanFlux", 

522 "ap_sigmaFlux", 

523 "ap_minMaxFlux", 

524 "ap_maxSlopeFlux", 

525 "ap_meanErrFlux", 

526 "ap_meanTotFlux"] 

527 

528 

529class DiaPipelineTask(pipeBase.PipelineTask): 

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

531 (DIA) Objects and Sources. 

532 """ 

533 ConfigClass = DiaPipelineConfig 

534 _DefaultName = "diaPipe" 

535 

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

537 super().__init__(**kwargs) 

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

539 self.schema = readSchemaFromApdb(self.apdb) 

540 self.makeSubtask("associator") 

541 self.makeSubtask("diaCalculation") 

542 if self.config.doRunForcedMeasurement: 

543 self.makeSubtask("diaForcedSource") 

544 if self.config.doPackageAlerts: 

545 self.makeSubtask("alertPackager") 

546 if self.config.doSolarSystemAssociation: 

547 self.makeSubtask("solarSystemAssociator") 

548 if self.config.filterUnAssociatedSources: 

549 columns = [self.config.newObjectFluxField, 

550 self.config.newObjectFluxField + "Err", 

551 "reliability"] 

552 columns += self.config.newObjectBadFlags 

553 

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

555 if missing: 

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

557 

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

559 inputs = butlerQC.get(inputRefs) 

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

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

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

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

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

565 inputs["legacySolarSystemTable"] = None 

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

567 inputs["solarSystemObjectTable"] = None 

568 

569 associationResults = pipeBase.Struct( 

570 apdbMarker=None, 

571 associatedDiaSources=None, 

572 diaForcedSources=None, 

573 diaObjects=None, 

574 associatedSsSources=None, 

575 unassociatedSsObjects=None, 

576 newDiaSources=None, 

577 marginalDiaSources=None, 

578 ) 

579 try: 

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

581 except pipeBase.AlgorithmError as e: 

582 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

583 e, 

584 self, 

585 log=self.log 

586 ) 

587 butlerQC.put(associationResults, outputRefs) 

588 raise error from e 

589 butlerQC.put(associationResults, outputRefs) 

590 

591 @timeMethod 

592 def run(self, 

593 diaSourceTable, 

594 legacySolarSystemTable, 

595 diffIm, 

596 exposure, 

597 template, 

598 preloadedDiaObjects=None, 

599 preloadedDiaSources=None, 

600 preloadedDiaForcedSources=None, 

601 band=None, 

602 idGenerator=None, 

603 solarSystemObjectTable=None, 

604 associationResults=None): 

605 """Process DiaSources and DiaObjects. 

606 

607 Load previous DiaObjects and their DiaSource history. Calibrate the 

608 values in the diaSourceCat. Associate new DiaSources with previous 

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

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

611 

612 Parameters 

613 ---------- 

614 diaSourceTable : `pandas.DataFrame` 

615 Newly detected DiaSources. 

616 legacySolarSystemTable : `pandas.DataFrame` 

617 Not used 

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

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

620 were detected. 

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

622 Calibrated exposure differenced with a template to create 

623 ``diffIm``. 

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

625 Template exposure used to create diffIm. 

626 preloadedDiaObjects : `pandas.DataFrame`, optional 

627 Previously detected DiaObjects, loaded from the APDB. 

628 `None` if ``doReloadAllApdbCatalogs`` is set. 

629 preloadedDiaSources : `pandas.DataFrame`, optional 

630 Previously detected DiaSources, loaded from the APDB. 

631 `None` if ``doReloadAllApdbCatalogs`` is set. 

632 preloadedDiaForcedSources : `pandas.DataFrame`, optional 

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

634 `None` if ``doReloadAllApdbCatalogs`` is set. 

635 band : `str` 

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

637 despite the `None` default. 

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

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

640 Required, despite the `None` default. 

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

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

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

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

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

646 this is also returned. 

647 

648 Returns 

649 ------- 

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

651 Results struct with components. 

652 

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

654 that this ccdVisit has completed successfully. 

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

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

657 DiaSources. (`pandas.DataFrame`) 

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

659 forced DiaSources. (`pandas.DataFrame`) 

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

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

662 (`pandas.DataFrame`) 

663 

664 Raises 

665 ------ 

666 RuntimeError 

667 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

668 ValueError 

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

670 """ 

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

672 # precede them, can be omitted. 

673 if band is None: 

674 raise ValueError("band is required.") 

675 if idGenerator is None: 

676 raise ValueError("idGenerator is required.") 

677 

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

679 associationResults = pipeBase.Struct() 

680 self._validateExposure(exposure, "Science") 

681 self._validateExposure(diffIm, "Difference") 

682 self._validateExposure(template, "Template") 

683 

684 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

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

686 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

687 region = getRegion(exposure) 

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

689 try: 

690 preloadedDiaObjects = self.loadRefreshedDiaObjects(region, preloadedDiaObjects) 

691 except Exception as e: 

692 if self.config.doReloadAllApdbCatalogs: 

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

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

695 # no preloaded diaObjects to fall back on. 

696 raise 

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

698 "the latest diaObjects from the APDB. Processing " 

699 "will continue with only the diaObjects from " 

700 "preload.", exc_info=e) 

701 finally: 

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

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

704 ) 

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

706 self.metadata["loadDiaObjectsDuration"]) 

707 

708 else: 

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

710 

711 if self.config.doReloadAllApdbCatalogs: 

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

713 # LoadDiaCatalogsTask to keep the metrics upload consistent. 

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

715 try: 

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

717 finally: 

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

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

720 ) 

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

722 self.metadata["loadDiaSourcesDuration"]) 

723 try: 

724 preloadedDiaForcedSources = self.loadRefreshedDiaForcedSources( 

725 region, preloadedDiaObjects, visitTime) 

726 finally: 

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

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

729 ) 

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

731 self.metadata["loadDiaForcedSourcesDuration"]) 

732 else: 

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

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

735 

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

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

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

739 

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

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

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

743 buffer=self.config.imagePixelMargin) 

744 else: 

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

746 diaObjects = preloadedDiaObjects 

747 

748 # Associate DiaSources with DiaObjects 

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

750 

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

752 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

753 assocResults.associatedDiaSources, 

754 "DiaSource", 

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

756 ) 

757 

758 # Merge associated diaSources 

759 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

760 preloadedDiaSources, 

761 assocResults.associatedDiaSources, 

762 diaObjects, 

763 assocResults.newDiaObjects, 

764 diffIm 

765 ) 

766 

767 # Compute DiaObject Summary statistics from their full DiaSource 

768 # history. 

769 diaCalResult = self.diaCalculation.run( 

770 mergedDiaObjects, 

771 mergedDiaSourceHistory, 

772 updatedDiaObjectIds, 

773 self.config.validBands) 

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

775 tableName="DiaObject", skipIndex=True) 

776 

777 # Test for duplication in the updated DiaObjects. 

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

779 raise RuntimeError( 

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

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

782 "reported. Exiting.") 

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

784 raise RuntimeError( 

785 "Duplicate DiaObjects (updated) created after " 

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

787 "reported. Exiting.") 

788 

789 # Forced source measurement 

790 if self.config.doRunForcedMeasurement: 

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

792 diaForcedSources = self.runForcedMeasurement( 

793 diaObjectsForced, updatedDiaObjects, exposure, diffIm, idGenerator 

794 ) 

795 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

796 else: 

797 # alertPackager needs correct columns 

798 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

799 forcedSourceHistoryThreshold = 0 

800 

801 # Write results to Alert Production Database (APDB) 

802 try: 

803 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

804 diaForcedSources) 

805 finally: 

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

807 clock="Utc") 

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

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

810 

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

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

813 associationResults.apdbMarker = pexConfig.Config() 

814 associationResults.associatedDiaSources = assocResults.associatedDiaSources 

815 associationResults.diaForcedSources = diaForcedSources 

816 associationResults.diaObjects = diaCalResult.diaObjectCat 

817 associationResults.unassociatedSsObjects = assocResults.unassociatedSsObjects 

818 associationResults.newDiaSources = assocResults.newDiaSources 

819 associationResults.marginalDiaSources = assocResults.marginalDiaSources 

820 

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

822 try: 

823 associatedSsSources = assocResults.associatedSsSources 

824 associatedSsSourcesPlusMpcorb = None 

825 if associatedSsSources is not None: 

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

827 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

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

829 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

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

831 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

832 associationResults.associatedSsSources = associatedSsSources 

833 

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

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

836 scale=DateTime.TAI) 

837 

838 # Package alerts 

839 if self.config.doPackageAlerts: 

840 # Append new forced sources to the full history 

841 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

842 tableName="DiaForcedSource") 

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

844 self.log.warning( 

845 "Duplicate DiaForcedSources created after merge with " 

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

847 "problems. Dropping duplicates.") 

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

849 diaForcedSourcesFull = diaForcedSourcesFull[ 

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

851 

852 self.alertPackager.run(assocResults.associatedDiaSources, 

853 updatedDiaObjects, 

854 preloadedDiaSources, 

855 diaForcedSourcesFull, 

856 diffIm, 

857 exposure, 

858 template, 

859 ssSrc=associatedSsSourcesPlusMpcorb, 

860 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

861 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

862 ) 

863 except Exception as e: 

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

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

866 

867 return associationResults 

868 

869 def _validateExposure(self, exposure, label): 

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

871 the APDB. 

872 

873 Parameters 

874 ---------- 

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

876 Calibrated science image that was subtracted. Corresponds to the 

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

878 label : `str` 

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

880 

881 Raises 

882 ------ 

883 ValueError 

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

885 """ 

886 errors = [] 

887 

888 photo_calib = exposure.getPhotoCalib() 

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

890 errors.append( 

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

892 "in alert packaging will fail." 

893 ) 

894 else: 

895 mean = photo_calib.getCalibrationMean() 

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

897 errors.append( 

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

899 "must be a finite, positive value." 

900 ) 

901 

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

903 errors.append( 

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

905 "DIA Object association will fail." 

906 ) 

907 

908 # Skip visitInfo checks for the template 

909 if label != "Template":