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

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 DiaPipelineConnections( 

79 pipeBase.PipelineTaskConnections, 

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

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

82 """Butler connections for DiaPipelineTask. 

83 """ 

84 diaSourceTable = connTypes.Input( 

85 doc="Catalog of calibrated DiaSources.", 

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

87 storageClass="DataFrame", 

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

89 ) 

90 solarSystemObjectTable = connTypes.Input( 

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

92 "this detectorVisit.", 

93 name="preloaded_SsObjects", 

94 storageClass="ArrowAstropy", 

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

96 minimum=0, 

97 ) 

98 diffIm = connTypes.Input( 

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

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

101 storageClass="ExposureF", 

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

103 ) 

104 exposure = connTypes.Input( 

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

106 "image differencing.", 

107 name="{fakesType}calexp", 

108 storageClass="ExposureF", 

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

110 ) 

111 template = connTypes.Input( 

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

113 "matched.", 

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

115 storageClass="ExposureF", 

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

117 ) 

118 preloadedDiaObjects = connTypes.Input( 

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

120 name="preloaded_diaObjects", 

121 storageClass="DataFrame", 

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

123 ) 

124 preloadedDiaSources = connTypes.Input( 

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

126 name="preloaded_diaSources", 

127 storageClass="DataFrame", 

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

129 ) 

130 preloadedDiaForcedSources = connTypes.Input( 

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

132 name="preloaded_diaForcedSources", 

133 storageClass="DataFrame", 

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

135 ) 

136 apdbMarker = connTypes.Output( 

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

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

139 name="apdb_marker", 

140 storageClass="Config", 

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

142 ) 

143 associatedDiaSources = connTypes.Output( 

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

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

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

147 storageClass="ArrowAstropy", 

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

149 ) 

150 associatedSsSources = connTypes.Output( 

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

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

153 storageClass="ArrowAstropy", 

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

155 ) 

156 unassociatedSsObjects = connTypes.Output( 

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

158 name="ssUnassociatedObjects", 

159 storageClass="ArrowAstropy", 

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

161 ) 

162 

163 diaForcedSources = connTypes.Output( 

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

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

166 storageClass="ArrowAstropy", 

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

168 ) 

169 diaObjects = connTypes.Output( 

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

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

172 storageClass="ArrowAstropy", 

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

174 ) 

175 newDiaSources = connTypes.Output( 

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

177 " were used to create a new diaObject", 

178 name="{fakesType}new_dia_source", 

179 storageClass="ArrowAstropy", 

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

181 ) 

182 marginalDiaSources = connTypes.Output( 

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

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

185 name="{fakesType}marginal_new_dia_source", 

186 storageClass="ArrowAstropy", 

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

188 ) 

189 

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

191 super().__init__(config=config) 

192 

193 if not config.doWriteAssociatedSources: 

194 self.outputs.remove("associatedDiaSources") 

195 self.outputs.remove("diaForcedSources") 

196 self.outputs.remove("diaObjects") 

197 self.outputs.remove("newDiaSources") 

198 self.outputs.remove("marginalDiaSources") 

199 else: 

200 if not config.doRunForcedMeasurement: 

201 self.outputs.remove("diaForcedSources") 

202 if not config.filterUnAssociatedSources: 

203 self.outputs.remove("newDiaSources") 

204 self.outputs.remove("marginalDiaSources") 

205 if not config.doSolarSystemAssociation: 

206 self.inputs.remove("solarSystemObjectTable") 

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

208 self.outputs.remove("associatedSsSources") 

209 self.outputs.remove("unassociatedSsObjects") 

210 

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

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

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

214 of the activator. 

215 

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

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

218 

219 Parameters 

220 ---------- 

221 inputs : `dict` 

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

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

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

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

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

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

228 dictionaries are guaranteed to be temporary copies that are true 

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

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

231 outputs : `dict` 

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

233 label : `str` 

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

235 diagnostic messages). 

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

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

238 diagnostic messages). 

239 

240 Returns 

241 ------- 

242 adjusted_inputs : `dict` 

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

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

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

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

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

248 `PipelineTask.runQuantum`. 

249 adjusted_outputs : `dict` 

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

251 interpretation as ``adjusted_inputs``. 

252 

253 Raises 

254 ------ 

255 ScalarError 

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

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

258 NoWorkFound 

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

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

261 quantum should be pruned or skipped. 

262 FileNotFoundError 

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

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

265 `PrerequisiteInput` connection. 

266 """ 

267 _, refs = inputs["diffIm"] 

268 for ref in refs: 

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

270 raise ValueError( 

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

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

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

274 "the validBands list in DiaPipelineConfig and add the " 

275 "appropriate columns to the Apdb schema.") 

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

277 

278 

279class DiaPipelineConfig(pipeBase.PipelineTaskConfig, 

280 pipelineConnections=DiaPipelineConnections): 

281 """Config for DiaPipelineTask. 

282 """ 

283 coaddName = pexConfig.Field( 

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

285 dtype=str, 

286 default="deep", 

287 ) 

288 apdb_config_url = pexConfig.Field( 

289 dtype=str, 

290 default=None, 

291 optional=False, 

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

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

294 "The database must already be initialized.", 

295 ) 

296 validBands = pexConfig.ListField( 

297 dtype=str, 

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

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

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

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

302 ) 

303 associator = pexConfig.ConfigurableField( 

304 target=AssociationTask, 

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

306 ) 

307 doSolarSystemAssociation = pexConfig.Field( 

308 dtype=bool, 

309 default=True, 

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

311 ) 

312 solarSystemAssociator = pexConfig.ConfigurableField( 

313 target=SolarSystemAssociationTask, 

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

315 ) 

316 diaCalculation = pexConfig.ConfigurableField( 

317 target=DiaObjectCalculationTask, 

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

319 ) 

320 doReloadDiaObjects = pexConfig.Field( 

321 dtype=bool, 

322 default=True, 

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

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

325 "are needed.", 

326 ) 

327 angleMargin = pexConfig.RangeField( 

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

329 dtype=float, 

330 default=2, 

331 min=0, 

332 ) 

333 doRunForcedMeasurement = pexConfig.Field( 

334 dtype=bool, 

335 default=True, 

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

337 "reasons during the ops rehearsal. " 

338 "It is expected to be removed.", 

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

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

341 ) 

342 diaForcedSource = pexConfig.ConfigurableField( 

343 target=DiaForcedSourceTask, 

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

345 "difference images.", 

346 ) 

347 forcedReliabilityThreshold = pexConfig.Field( 

348 dtype=float, 

349 default=0.5, 

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

351 "determining whether to calculate forced photometry for a " 

352 "diaObject.", 

353 ) 

354 forcedTrailLengthThreshold = pexConfig.Field( 

355 dtype=float, 

356 default=1.4, 

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

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

359 "diaObject.", 

360 ) 

361 forcedBadFlags = pexConfig.ListField( 

362 dtype=str, 

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

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

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

366 ) 

367 alertPackager = pexConfig.ConfigurableField( 

368 target=PackageAlertsTask, 

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

370 ) 

371 doPackageAlerts = pexConfig.Field( 

372 dtype=bool, 

373 default=False, 

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

375 "write them to disk.", 

376 ) 

377 doWriteAssociatedSources = pexConfig.Field( 

378 dtype=bool, 

379 default=True, 

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

381 "formatted following the Science Data Model.", 

382 ) 

383 imagePixelMargin = pexConfig.RangeField( 

384 dtype=int, 

385 default=10, 

386 min=0, 

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

388 "diaObjects for association.", 

389 ) 

390 filterUnAssociatedSources = pexConfig.Field( 

391 dtype=bool, 

392 default=True, 

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

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

395 "will not be affected." 

396 ) 

397 newObjectBadFlags = pexConfig.ListField( 

398 dtype=str, 

399 default=("centroid_flag", 

400 "pixelFlags_crCenter", 

401 "pixelFlags_nodataCenter", 

402 "pixelFlags_interpolatedCenter", 

403 "pixelFlags_saturatedCenter", 

404 "pixelFlags_suspectCenter", 

405 "pixelFlags_streakCenter", 

406 "glint_trail"), 

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

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

409 ) 

410 maxNewDiaObjects = pexConfig.RangeField( 

411 dtype=float, 

412 default=0, 

413 min=0, 

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

415 "Set to zero to disable.", 

416 ) 

417 newObjectSnrThreshold = pexConfig.RangeField( 

418 dtype=float, 

419 default=0, 

420 min=0, 

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

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

423 "diaObjects." 

424 "Set to zero to disable.", 

425 ) 

426 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

427 dtype=float, 

428 default=0, 

429 min=0, 

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

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

432 " low reliability scores before creating new diaObjects." 

433 "Set to zero to disable.", 

434 ) 

435 newObjectReliabilityThreshold = pexConfig.RangeField( 

436 dtype=float, 

437 default=0.1, 

438 min=0, 

439 max=1, 

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

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

442 "diaObjects." 

443 "Set to zero to disable.", 

444 ) 

445 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

446 dtype=float, 

447 default=0.1, 

448 min=0, 

449 max=1, 

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

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

452 "before creating new diaObjects." 

453 "Set to zero to disable.", 

454 ) 

455 newObjectFluxField = pexConfig.Field( 

456 dtype=str, 

457 default="apFlux", 

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

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

460 ) 

461 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

462 

463 def setDefaults(self): 

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

465 "ap_nDiaSources", 

466 "ap_meanFlux", 

467 "ap_sigmaFlux", 

468 "ap_minMaxFlux", 

469 "ap_maxSlopeFlux", 

470 "ap_meanErrFlux", 

471 "ap_meanTotFlux"] 

472 

473 

474class DiaPipelineTask(pipeBase.PipelineTask): 

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

476 (DIA) Objects and Sources. 

477 """ 

478 ConfigClass = DiaPipelineConfig 

479 _DefaultName = "diaPipe" 

480 

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

482 super().__init__(**kwargs) 

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

484 self.schema = readSchemaFromApdb(self.apdb) 

485 self.makeSubtask("associator") 

486 self.makeSubtask("diaCalculation") 

487 if self.config.doRunForcedMeasurement: 

488 self.makeSubtask("diaForcedSource") 

489 if self.config.doPackageAlerts: 

490 self.makeSubtask("alertPackager") 

491 if self.config.doSolarSystemAssociation: 

492 self.makeSubtask("solarSystemAssociator") 

493 if self.config.filterUnAssociatedSources: 

494 columns = [self.config.newObjectFluxField, 

495 self.config.newObjectFluxField + "Err", 

496 "reliability"] 

497 columns += self.config.newObjectBadFlags 

498 

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

500 if missing: 

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

502 

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

504 inputs = butlerQC.get(inputRefs) 

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

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

507 inputs["legacySolarSystemTable"] = None 

508 if not self.config.doSolarSystemAssociation: 

509 inputs["solarSystemObjectTable"] = None 

510 

511 outputs = self.run(**inputs) 

512 

513 butlerQC.put(outputs, outputRefs) 

514 

515 @timeMethod 

516 def run(self, 

517 diaSourceTable, 

518 legacySolarSystemTable, 

519 diffIm, 

520 exposure, 

521 template, 

522 preloadedDiaObjects, 

523 preloadedDiaSources, 

524 preloadedDiaForcedSources, 

525 band, 

526 idGenerator, 

527 solarSystemObjectTable=None): 

528 """Process DiaSources and DiaObjects. 

529 

530 Load previous DiaObjects and their DiaSource history. Calibrate the 

531 values in the diaSourceCat. Associate new DiaSources with previous 

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

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

534 

535 Parameters 

536 ---------- 

537 diaSourceTable : `pandas.DataFrame` 

538 Newly detected DiaSources. 

539 legacySolarSystemTable : `pandas.DataFrame` 

540 Not used 

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

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

543 were detected. 

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

545 Calibrated exposure differenced with a template to create 

546 ``diffIm``. 

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

548 Template exposure used to create diffIm. 

549 preloadedDiaObjects : `pandas.DataFrame` 

550 Previously detected DiaObjects, loaded from the APDB. 

551 preloadedDiaSources : `pandas.DataFrame` 

552 Previously detected DiaSources, loaded from the APDB. 

553 preloadedDiaForcedSources : `pandas.DataFrame` 

554 Catalog of previously detected forced DiaSources, from the APDB 

555 band : `str` 

556 The band in which the new DiaSources were detected. 

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

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

559 solarSystemObjectTable : `astropy.table.Table` 

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

561 

562 Returns 

563 ------- 

564 results : `lsst.pipe.base.Struct` 

565 Results struct with components. 

566 

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

568 that this ccdVisit has completed successfully. 

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

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

571 DiaSources. (`pandas.DataFrame`) 

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

573 forced DiaSources. (`pandas.DataFrame`) 

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

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

576 (`pandas.DataFrame`) 

577 

578 Raises 

579 ------ 

580 RuntimeError 

581 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

582 """ 

583 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

584 if legacySolarSystemTable is not None and solarSystemObjectTable is None: 

585 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

586 if self.config.doReloadDiaObjects: 

587 try: 

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

589 except Exception as e: 

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

591 "the latest diaObjects from the APDB. Processing " 

592 "will continue with only the diaObjects from " 

593 "preload.", exc_info=e) 

594 finally: 

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

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

597 ) 

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

599 self.metadata["loadDiaObjectsDuration"]) 

600 

601 else: 

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

603 

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

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

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

607 

608 if not preloadedDiaObjects.empty: 

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

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

611 buffer=self.config.imagePixelMargin) 

612 else: 

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

614 diaObjects = preloadedDiaObjects 

615 

616 # Associate DiaSources with DiaObjects 

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

618 

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

620 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

621 assocResults.associatedDiaSources, 

622 "DiaSource", 

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

624 ) 

625 

626 # Merge associated diaSources 

627 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

628 preloadedDiaSources, 

629 assocResults.associatedDiaSources, 

630 diaObjects, 

631 assocResults.newDiaObjects, 

632 diffIm 

633 ) 

634 

635 # Compute DiaObject Summary statistics from their full DiaSource 

636 # history. 

637 diaCalResult = self.diaCalculation.run( 

638 mergedDiaObjects, 

639 mergedDiaSourceHistory, 

640 updatedDiaObjectIds, 

641 [band]) 

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

643 tableName="DiaObject", skipIndex=True) 

644 

645 # Test for duplication in the updated DiaObjects. 

646 if self.testDataFrameIndex(diaCalResult.diaObjectCat): 

647 raise RuntimeError( 

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

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

650 "reported. Exiting.") 

651 if self.testDataFrameIndex(updatedDiaObjects): 

652 raise RuntimeError( 

653 "Duplicate DiaObjects (updated) created after " 

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

655 "reported. Exiting.") 

656 

657 # Forced source measurement 

658 if self.config.doRunForcedMeasurement: 

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

660 diaForcedSources = self.runForcedMeasurement( 

661 diaObjectsForced, updatedDiaObjects, exposure, diffIm, idGenerator 

662 ) 

663 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

664 else: 

665 # alertPackager needs correct columns 

666 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

667 forcedSourceHistoryThreshold = 0 

668 

669 # Write results to Alert Production Database (APDB) 

670 try: 

671 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

672 diaForcedSources) 

673 finally: 

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

675 clock="Utc") 

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

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

678 

679 associatedSsSources = assocResults.associatedSsSources 

680 associatedSsSourcesPlusMpcorb = None 

681 if associatedSsSources is not None: 

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

683 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

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

685 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

686 for mpcorbColumn in mpcorbColumns: 

687 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

688 

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

690 updatedDiaObjects['validityStartMjdTai'] = validityStart.get(system=DateTime.MJD, scale=DateTime.TAI) 

691 

692 # Package alerts 

693 if self.config.doPackageAlerts: 

694 # Append new forced sources to the full history 

695 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

696 tableName="DiaForcedSource") 

697 if self.testDataFrameIndex(diaForcedSources): 

698 self.log.warning( 

699 "Duplicate DiaForcedSources created after merge with " 

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

701 "problems. Dropping duplicates.") 

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

703 # Reset due to the index shape being slight different than 

704 # expected. 

705 diaForcedSourcesFull = diaForcedSourcesFull.groupby( 

706 diaForcedSourcesFull.index).first() 

707 diaForcedSourcesFull.reset_index(drop=True, inplace=True) 

708 diaForcedSourcesFull.set_index( 

709 ["diaObjectId", "diaForcedSourceId"], 

710 drop=False, 

711 inplace=True) 

712 

713 self.alertPackager.run(assocResults.associatedDiaSources, 

714 updatedDiaObjects, 

715 preloadedDiaSources, 

716 diaForcedSourcesFull, 

717 diffIm, 

718 exposure, 

719 template, 

720 ssSrc=associatedSsSourcesPlusMpcorb, 

721 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

722 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

723 ) 

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 marker = pexConfig.Config() 

728 return pipeBase.Struct(apdbMarker=marker, 

729 associatedDiaSources=assocResults.associatedDiaSources, 

730 diaForcedSources=diaForcedSources, 

731 diaObjects=diaCalResult.diaObjectCat, 

732 associatedSsSources=associatedSsSources, 

733 unassociatedSsObjects=assocResults.unassociatedSsObjects, 

734 newDiaSources=assocResults.newDiaSources, 

735 marginalDiaSources=assocResults.marginalDiaSources 

736 ) 

737 

738 def _selectGoodDiaObjects(self, diaObjects, diaSources): 

739 """ 

740 Extract a subset of diaObjects with multiple associated diaSources that 

741 pass selection requirements. 

742 """ 

743 n_history = self.diaForcedSource.config.historyThreshold 

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

745 for col in required_columns: 

746 if col not in diaSources.columns: 

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

748 

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

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

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

752 if badFlags: 

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

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

755 # Exclude rows where ANY bad flag is True. 

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

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

758 mask = rel_ok & trailed_ok & ~bad_any 

759 ids = diaSources.index.get_level_values("diaObjectId")[mask] 

760 

761 # Count occurrences in diaSource 

762 counts = pd.Series(ids, copy=False).value_counts() 

763 

764 # IDs that exceed threshold 

765 keep_ids = counts[counts >= n_history].index 

766 

767 # Since diaObject index is unique, direct index filtering is safe and efficient 

768 return diaObjects.loc[diaObjects.index.intersection(keep_ids)] 

769 

770 def createNewDiaObjects(self, unassociatedDiaSources): 

771 """Loop through the set of DiaSources and create new DiaObjects 

772 for unassociated DiaSources. 

773 

774 Parameters 

775 ---------- 

776 unassociatedDiaSources : `pandas.DataFrame` 

777 Set of DiaSources to create new DiaObjects from. 

778 

779 Returns 

780 ------- 

781 results : `lsst.pipe.base.Struct` 

782 Results struct containing: 

783 

784 - diaSources : `pandas.DataFrame` 

785 DiaSource catalog with updated DiaObject ids. 

786 - newDiaObjects : `pandas.DataFrame` 

787 Newly created DiaObjects from the unassociated DiaSources. 

788 - nNewDiaObjects : `int` 

789 Number of newly created diaObjects. 

790 - marginalDiaSources : `pandas.DataFrame` 

791 Unassociated diaSources with low signal-to-noise and/or 

792 reliability, which are excluded from the new DiaObjects. 

793 """ 

794 marginalDiaSources = make_empty_catalog(self.schema, tableName="DiaSource") 

795 if len(unassociatedDiaSources) == 0: 

796 newDiaObjects = make_empty_catalog(self.schema, tableName="DiaObject") 

797 else: 

798 if self.config.filterUnAssociatedSources: 

799 results = self.filterSources(unassociatedDiaSources) 

800 unassociatedDiaSources = results.goodSources 

801 marginalDiaSources = results.badSources 

802 unassociatedDiaSources["diaObjectId"] = unassociatedDiaSources["diaSourceId"] 

803 newDiaObjects = convertDataFrameToSdmSchema(self.schema, unassociatedDiaSources, 

804 tableName="DiaObject", skipIndex=True) 

805 self.metadata["nRejectedNewDiaObjects"] = len(marginalDiaSources) 

806 return pipeBase.Struct(diaSources=unassociatedDiaSources, 

807 newDiaObjects=newDiaObjects, 

808 nNewDiaObjects=len(newDiaObjects), 

809 marginalDiaSources=marginalDiaSources) 

810 

811 def filterSources(self, sources, snrThreshold=None, lowReliabilitySnrThreshold=None, 

812 reliabilityThreshold=None, lowSnrReliabilityThreshold=None, badFlags=None): 

813 """Select good sources out of a catalog. 

814 

815 Parameters 

816 ---------- 

817 sources : `pandas.DataFrame` 

818 Set of DiaSources to check. 

819 snrThreshold : `float`, optional 

820 The minimum signal to noise diaSource to make a new diaObject. 

821 Included for unit tests. Uses the task config value if not set. 

822 lowReliabilitySnrThreshold : `float`, optional 

823 Use ``lowSnrReliabilityThreshold`` as the reliability threshold for 

824 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

825 Included for unit tests. Uses the task config value if not set. 

826 reliabilityThreshold : `float`, optional 

827 The minimum reliability score diaSource to make a new diaObject 

828 Included for unit tests. Uses the task config value if not set. 

829 lowSnrReliabilityThreshold : `float`, optional 

830 Use ``lowSnrReliabilityThreshold`` as the reliability threshold for 

831 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

832 Included for unit tests. Uses the task config value if not set. 

833 badFlags : `list` of `str`, optional 

834 Do not create new diaObjects for any diaSource with any of these 

835 flags set. 

836 Included for unit tests. Uses the task config value if not set. 

837 

838 Returns 

839 ------- 

840 results : `lsst.pipe.base.Struct` 

841 Results struct containing: 

842 

843 - goodSources : `pandas.DataFrame` 

844 Subset of the input sources that pass all checks. 

845 - badSources : `pandas.DataFrame` 

846 Subset of the input sources that fail any check. 

847 """ 

848 if snrThreshold is None: 

849 snrThreshold = self.config.newObjectSnrThreshold 

850 if lowReliabilitySnrThreshold is None: 

851 lowReliabilitySnrThreshold = self.config.newObjectLowReliabilitySnrThreshold 

852 if reliabilityThreshold is None: 

853 reliabilityThreshold = self.config.newObjectReliabilityThreshold 

854 if lowSnrReliabilityThreshold is None: 

855 lowSnrReliabilityThreshold = self.config.newObjectLowSnrReliabilityThreshold 

856 if badFlags is None: 

857 badFlags = self.config.newObjectBadFlags 

858 flagged = sources[badFlags].fillna(False).any(axis=1) 

859 fluxField = self.config.newObjectFluxField 

860 fluxErrField = fluxField + "Err" 

861 signalToNoise = np.abs(np.array(sources[fluxField]/sources[fluxErrField])) 

862 reliability = np.array(sources['reliability']) 

863 nFlagged = np.count_nonzero(flagged) 

864 if nFlagged > 0: 

865 self.log.info("Not creating new diaObjects for %i unassociated diaSources due to flags", nFlagged) 

866 

867 if snrThreshold > 0: 

868 snr_flag = signalToNoise < snrThreshold 

869 self.log.info("Not creating new diaObjects for %i unassociated diaSources due to %sFlux" 

870 " signal to noise < %f", 

871 np.sum(snr_flag), fluxField, snrThreshold) 

872 flagged += snr_flag 

873 if reliabilityThreshold > 0: 

874 reliability_flag = reliability < reliabilityThreshold 

875 self.log.info("Not creating new diaObjects for %i unassociated diaSources due to " 

876 "reliability<%f", 

877 np.sum(reliability_flag), reliabilityThreshold) 

878 flagged += reliability_flag 

879 if min(lowReliabilitySnrThreshold, lowSnrReliabilityThreshold) > 0: 

880 # Only run the combined test if both thresholds are greater than zero 

881 lowSnrReliability_flag = ((signalToNoise < lowReliabilitySnrThreshold) 

882 & (reliability < lowSnrReliabilityThreshold)) 

883 self.log.info("Not creating new diaObjects for %i unassociated diaSources due to %sFlux" 

884 " signal to noise < %f combined with reliability< %f", 

885 np.sum(lowSnrReliability_flag), 

886 fluxField, 

887 lowReliabilitySnrThreshold, 

888 lowSnrReliabilityThreshold) 

889 flagged += lowSnrReliability_flag 

890 

891 if np.count_nonzero(~flagged) > 0: 

892 goodSources = sources[~flagged].copy(deep=True) 

893 else: 

894 goodSources = make_empty_catalog(self.schema, tableName="DiaSource") 

895 if np.count_nonzero(flagged) > 0: 

896 badSources = sources[flagged].copy(deep=True) 

897 else: 

898 badSources = make_empty_catalog(self.schema, tableName="DiaSource") 

899 return pipeBase.Struct(goodSources=goodSources, 

900 badSources=badSources 

901 )