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

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 

38 

39from astropy.table import Table 

40import numpy as np 

41import pandas as pd 

42from lsst.ap.association import ( 

43 AssociationTask, 

44 DiaForcedSourceTask, 

45 PackageAlertsTask) 

46 

47from lsst.ap.association.utils import convertDataFrameToSdmSchema, checkSdmSchemaColumns, \ 

48 readSchemaFromApdb, dropEmptyColumns, make_empty_catalog, makeEmptyForcedSourceTable 

49from lsst.daf.base import DateTime 

50from lsst.meas.base import DetectorVisitIdGeneratorConfig, \ 

51 DiaObjectCalculationTask 

52from lsst.pipe.tasks.ssoAssociation import SolarSystemAssociationTask 

53from lsst.utils.timer import timeMethod, duration_from_timeMethod 

54 

55 

56class TooManyDiaObjectsError(pipeBase.AlgorithmError): 

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

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

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

60 """ 

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

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

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

64 f" configured maximum of {threshold}") 

65 super().__init__(msg) 

66 self.nNewDiaObjects = nNewDiaObjects 

67 self.threshold = threshold 

68 

69 @property 

70 def metadata(self): 

71 return {"nNewDiaObjects": self.nNewDiaObjects, 

72 "threshold": self.threshold 

73 } 

74 

75 

76class DiaPipelineConnections( 

77 pipeBase.PipelineTaskConnections, 

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

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

80 """Butler connections for DiaPipelineTask. 

81 """ 

82 diaSourceTable = connTypes.Input( 

83 doc="Catalog of calibrated DiaSources.", 

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

85 storageClass="DataFrame", 

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

87 ) 

88 solarSystemObjectTable = connTypes.Input( 

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

90 "this detectorVisit.", 

91 name="preloaded_SsObjects", 

92 storageClass="ArrowAstropy", 

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

94 minimum=0, 

95 ) 

96 diffIm = connTypes.Input( 

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

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

99 storageClass="ExposureF", 

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

101 ) 

102 exposure = connTypes.Input( 

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

104 "image differencing.", 

105 name="{fakesType}calexp", 

106 storageClass="ExposureF", 

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

108 ) 

109 template = connTypes.Input( 

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

111 "matched.", 

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

113 storageClass="ExposureF", 

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

115 ) 

116 preloadedDiaObjects = connTypes.Input( 

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

118 name="preloaded_diaObjects", 

119 storageClass="DataFrame", 

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

121 ) 

122 preloadedDiaSources = connTypes.Input( 

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

124 name="preloaded_diaSources", 

125 storageClass="DataFrame", 

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

127 ) 

128 preloadedDiaForcedSources = connTypes.Input( 

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

130 name="preloaded_diaForcedSources", 

131 storageClass="DataFrame", 

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

133 ) 

134 apdbMarker = connTypes.Output( 

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

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

137 name="apdb_marker", 

138 storageClass="Config", 

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

140 ) 

141 associatedDiaSources = connTypes.Output( 

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

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

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

145 storageClass="ArrowAstropy", 

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

147 ) 

148 associatedSsSources = connTypes.Output( 

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

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

151 storageClass="ArrowAstropy", 

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

153 ) 

154 unassociatedSsObjects = connTypes.Output( 

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

156 name="ssUnassociatedObjects", 

157 storageClass="ArrowAstropy", 

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

159 ) 

160 

161 diaForcedSources = connTypes.Output( 

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

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

164 storageClass="ArrowAstropy", 

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

166 ) 

167 diaObjects = connTypes.Output( 

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

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

170 storageClass="ArrowAstropy", 

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

172 ) 

173 newDiaSources = connTypes.Output( 

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

175 " were used to create a new diaObject", 

176 name="{fakesType}new_dia_source", 

177 storageClass="ArrowAstropy", 

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

179 ) 

180 marginalDiaSources = connTypes.Output( 

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

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

183 name="{fakesType}marginal_new_dia_source", 

184 storageClass="ArrowAstropy", 

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

186 ) 

187 

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

189 super().__init__(config=config) 

190 

191 if not config.doWriteAssociatedSources: 

192 self.outputs.remove("associatedDiaSources") 

193 self.outputs.remove("diaForcedSources") 

194 self.outputs.remove("diaObjects") 

195 self.outputs.remove("newDiaSources") 

196 self.outputs.remove("marginalDiaSources") 

197 else: 

198 if not config.doRunForcedMeasurement: 

199 self.outputs.remove("diaForcedSources") 

200 if not config.filterUnAssociatedSources: 

201 self.outputs.remove("newDiaSources") 

202 self.outputs.remove("marginalDiaSources") 

203 if not config.doSolarSystemAssociation: 

204 self.inputs.remove("solarSystemObjectTable") 

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

206 self.outputs.remove("associatedSsSources") 

207 self.outputs.remove("unassociatedSsObjects") 

208 

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

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

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

212 of the activator. 

213 

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

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

216 

217 Parameters 

218 ---------- 

219 inputs : `dict` 

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

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

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

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

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

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

226 dictionaries are guaranteed to be temporary copies that are true 

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

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

229 outputs : `dict` 

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

231 label : `str` 

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

233 diagnostic messages). 

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

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

236 diagnostic messages). 

237 

238 Returns 

239 ------- 

240 adjusted_inputs : `dict` 

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

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

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

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

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

246 `PipelineTask.runQuantum`. 

247 adjusted_outputs : `dict` 

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

249 interpretation as ``adjusted_inputs``. 

250 

251 Raises 

252 ------ 

253 ScalarError 

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

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

256 NoWorkFound 

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

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

259 quantum should be pruned or skipped. 

260 FileNotFoundError 

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

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

263 `PrerequisiteInput` connection. 

264 """ 

265 _, refs = inputs["diffIm"] 

266 for ref in refs: 

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

268 raise ValueError( 

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

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

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

272 "the validBands list in DiaPipelineConfig and add the " 

273 "appropriate columns to the Apdb schema.") 

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

275 

276 

277class DiaPipelineConfig(pipeBase.PipelineTaskConfig, 

278 pipelineConnections=DiaPipelineConnections): 

279 """Config for DiaPipelineTask. 

280 """ 

281 coaddName = pexConfig.Field( 

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

283 dtype=str, 

284 default="deep", 

285 ) 

286 apdb_config_url = pexConfig.Field( 

287 dtype=str, 

288 default=None, 

289 optional=False, 

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

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

292 "The database must already be initialized.", 

293 ) 

294 validBands = pexConfig.ListField( 

295 dtype=str, 

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

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

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

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

300 ) 

301 associator = pexConfig.ConfigurableField( 

302 target=AssociationTask, 

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

304 ) 

305 doSolarSystemAssociation = pexConfig.Field( 

306 dtype=bool, 

307 default=True, 

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

309 ) 

310 solarSystemAssociator = pexConfig.ConfigurableField( 

311 target=SolarSystemAssociationTask, 

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

313 ) 

314 diaCalculation = pexConfig.ConfigurableField( 

315 target=DiaObjectCalculationTask, 

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

317 ) 

318 doRunForcedMeasurement = pexConfig.Field( 

319 dtype=bool, 

320 default=True, 

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

322 "reasons during the ops rehearsal. " 

323 "It is expected to be removed.", 

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

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

326 ) 

327 diaForcedSource = pexConfig.ConfigurableField( 

328 target=DiaForcedSourceTask, 

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

330 "difference images.", 

331 ) 

332 alertPackager = pexConfig.ConfigurableField( 

333 target=PackageAlertsTask, 

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

335 ) 

336 doPackageAlerts = pexConfig.Field( 

337 dtype=bool, 

338 default=False, 

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

340 "write them to disk.", 

341 ) 

342 doWriteAssociatedSources = pexConfig.Field( 

343 dtype=bool, 

344 default=True, 

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

346 "formatted following the Science Data Model.", 

347 ) 

348 imagePixelMargin = pexConfig.RangeField( 

349 dtype=int, 

350 default=10, 

351 min=0, 

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

353 "diaObjects for association.", 

354 ) 

355 filterUnAssociatedSources = pexConfig.Field( 

356 dtype=bool, 

357 default=True, 

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

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

360 "will not be affected." 

361 ) 

362 newObjectBadFlags = pexConfig.ListField( 

363 dtype=str, 

364 default=("centroid_flag", 

365 "pixelFlags_crCenter", 

366 "pixelFlags_nodataCenter", 

367 "pixelFlags_interpolatedCenter", 

368 "pixelFlags_saturatedCenter", 

369 "pixelFlags_suspectCenter", 

370 "pixelFlags_streakCenter", 

371 "glint_trail"), 

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

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

374 ) 

375 maxNewDiaObjects = pexConfig.RangeField( 

376 dtype=float, 

377 default=0, 

378 min=0, 

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

380 "Set to zero to disable.", 

381 ) 

382 newObjectSnrThreshold = pexConfig.RangeField( 

383 dtype=float, 

384 default=0, 

385 min=0, 

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

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

388 "diaObjects." 

389 "Set to zero to disable.", 

390 ) 

391 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

392 dtype=float, 

393 default=0, 

394 min=0, 

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

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

397 " low reliability scores before creating new diaObjects." 

398 "Set to zero to disable.", 

399 ) 

400 newObjectReliabilityThreshold = pexConfig.RangeField( 

401 dtype=float, 

402 default=0.1, 

403 min=0, 

404 max=1, 

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

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

407 "diaObjects." 

408 "Set to zero to disable.", 

409 ) 

410 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

411 dtype=float, 

412 default=0.1, 

413 min=0, 

414 max=1, 

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

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

417 "before creating new diaObjects." 

418 "Set to zero to disable.", 

419 ) 

420 newObjectFluxField = pexConfig.Field( 

421 dtype=str, 

422 default="apFlux", 

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

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

425 ) 

426 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

427 

428 def setDefaults(self): 

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

430 "ap_nDiaSources", 

431 "ap_meanFlux", 

432 "ap_sigmaFlux", 

433 "ap_minMaxFlux", 

434 "ap_maxSlopeFlux", 

435 "ap_meanErrFlux", 

436 "ap_meanTotFlux"] 

437 

438 

439class DiaPipelineTask(pipeBase.PipelineTask): 

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

441 (DIA) Objects and Sources. 

442 """ 

443 ConfigClass = DiaPipelineConfig 

444 _DefaultName = "diaPipe" 

445 

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

447 super().__init__(**kwargs) 

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

449 self.schema = readSchemaFromApdb(self.apdb) 

450 self.makeSubtask("associator") 

451 self.makeSubtask("diaCalculation") 

452 if self.config.doRunForcedMeasurement: 

453 self.makeSubtask("diaForcedSource") 

454 if self.config.doPackageAlerts: 

455 self.makeSubtask("alertPackager") 

456 if self.config.doSolarSystemAssociation: 

457 self.makeSubtask("solarSystemAssociator") 

458 if self.config.filterUnAssociatedSources: 

459 columns = [self.config.newObjectFluxField, 

460 self.config.newObjectFluxField + "Err", 

461 "reliability"] 

462 columns += self.config.newObjectBadFlags 

463 

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

465 if missing: 

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

467 

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

469 inputs = butlerQC.get(inputRefs) 

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

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

472 inputs["legacySolarSystemTable"] = None 

473 if not self.config.doSolarSystemAssociation: 

474 inputs["solarSystemObjectTable"] = None 

475 

476 outputs = self.run(**inputs) 

477 

478 butlerQC.put(outputs, outputRefs) 

479 

480 @timeMethod 

481 def run(self, 

482 diaSourceTable, 

483 legacySolarSystemTable, 

484 diffIm, 

485 exposure, 

486 template, 

487 preloadedDiaObjects, 

488 preloadedDiaSources, 

489 preloadedDiaForcedSources, 

490 band, 

491 idGenerator, 

492 solarSystemObjectTable=None): 

493 """Process DiaSources and DiaObjects. 

494 

495 Load previous DiaObjects and their DiaSource history. Calibrate the 

496 values in the diaSourceCat. Associate new DiaSources with previous 

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

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

499 

500 Parameters 

501 ---------- 

502 diaSourceTable : `pandas.DataFrame` 

503 Newly detected DiaSources. 

504 legacySolarSystemTable : `pandas.DataFrame` 

505 Not used 

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

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

508 were detected. 

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

510 Calibrated exposure differenced with a template to create 

511 ``diffIm``. 

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

513 Template exposure used to create diffIm. 

514 preloadedDiaObjects : `pandas.DataFrame` 

515 Previously detected DiaObjects, loaded from the APDB. 

516 preloadedDiaSources : `pandas.DataFrame` 

517 Previously detected DiaSources, loaded from the APDB. 

518 preloadedDiaForcedSources : `pandas.DataFrame` 

519 Catalog of previously detected forced DiaSources, from the APDB 

520 band : `str` 

521 The band in which the new DiaSources were detected. 

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

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

524 solarSystemObjectTable : `astropy.table.Table` 

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

526 

527 Returns 

528 ------- 

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

530 Results struct with components. 

531 

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

533 that this ccdVisit has completed successfully. 

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

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

536 DiaSources. (`pandas.DataFrame`) 

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

538 forced DiaSources. (`pandas.DataFrame`) 

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

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

541 (`pandas.DataFrame`) 

542 

543 Raises 

544 ------ 

545 RuntimeError 

546 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

547 """ 

548 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

549 if legacySolarSystemTable is not None and solarSystemObjectTable is None: 

550 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

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

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

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

554 

555 if not preloadedDiaObjects.empty: 

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

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

558 buffer=self.config.imagePixelMargin) 

559 else: 

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

561 diaObjects = preloadedDiaObjects 

562 

563 # Associate DiaSources with DiaObjects 

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

565 

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

567 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

568 assocResults.associatedDiaSources, 

569 "DiaSource", 

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

571 ) 

572 

573 # Merge associated diaSources 

574 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

575 preloadedDiaSources, 

576 assocResults.associatedDiaSources, 

577 diaObjects, 

578 assocResults.newDiaObjects, 

579 diffIm 

580 ) 

581 

582 # Compute DiaObject Summary statistics from their full DiaSource 

583 # history. 

584 diaCalResult = self.diaCalculation.run( 

585 mergedDiaObjects, 

586 mergedDiaSourceHistory, 

587 updatedDiaObjectIds, 

588 [band]) 

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

590 tableName="DiaObject") 

591 

592 # Test for duplication in the updated DiaObjects. 

593 if self.testDataFrameIndex(diaCalResult.diaObjectCat): 

594 raise RuntimeError( 

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

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

597 "reported. Exiting.") 

598 if self.testDataFrameIndex(updatedDiaObjects): 

599 raise RuntimeError( 

600 "Duplicate DiaObjects (updated) created after " 

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

602 "reported. Exiting.") 

603 

604 # Forced source measurement 

605 if self.config.doRunForcedMeasurement: 

606 diaForcedSources = self.runForcedMeasurement( 

607 diaCalResult.diaObjectCat, updatedDiaObjects, exposure, diffIm, idGenerator 

608 ) 

609 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

610 else: 

611 # alertPackager needs correct columns 

612 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

613 forcedSourceHistoryThreshold = 0 

614 

615 # Write results to Alert Production Database (APDB) 

616 try: 

617 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

618 diaForcedSources) 

619 finally: 

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

621 clock="Utc") 

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

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

624 

625 associatedSsSources = assocResults.associatedSsSources 

626 associatedSsSourcesPlusMpcorb = None 

627 if associatedSsSources is not None: 

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

629 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

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

631 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

632 for mpcorbColumn in mpcorbColumns: 

633 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

634 

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

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

637 

638 # Package alerts 

639 if self.config.doPackageAlerts: 

640 # Append new forced sources to the full history 

641 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

642 tableName="DiaForcedSource") 

643 if self.testDataFrameIndex(diaForcedSources): 

644 self.log.warning( 

645 "Duplicate DiaForcedSources created after merge with " 

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

647 "problems. Dropping duplicates.") 

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

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

650 # expected. 

651 diaForcedSourcesFull = diaForcedSourcesFull.groupby( 

652 diaForcedSourcesFull.index).first() 

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

654 diaForcedSourcesFull.set_index( 

655 ["diaObjectId", "diaForcedSourceId"], 

656 drop=False, 

657 inplace=True) 

658 

659 self.alertPackager.run(assocResults.associatedDiaSources, 

660 updatedDiaObjects, 

661 preloadedDiaSources, 

662 diaForcedSourcesFull, 

663 diffIm, 

664 exposure, 

665 template, 

666 ssSrc=associatedSsSourcesPlusMpcorb, 

667 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

668 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

669 ) 

670 

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

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

673 marker = pexConfig.Config() 

674 return pipeBase.Struct(apdbMarker=marker, 

675 associatedDiaSources=assocResults.associatedDiaSources, 

676 diaForcedSources=diaForcedSources, 

677 diaObjects=diaCalResult.diaObjectCat, 

678 associatedSsSources=associatedSsSources, 

679 unassociatedSsObjects=assocResults.unassociatedSsObjects, 

680 newDiaSources=assocResults.newDiaSources, 

681 marginalDiaSources=assocResults.marginalDiaSources 

682 ) 

683 

684 def createNewDiaObjects(self, unassociatedDiaSources): 

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

686 for unassociated DiaSources. 

687 

688 Parameters 

689 ---------- 

690 unassociatedDiaSources : `pandas.DataFrame` 

691 Set of DiaSources to create new DiaObjects from. 

692 

693 Returns 

694 ------- 

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

696 Results struct containing: 

697 

698 - diaSources : `pandas.DataFrame` 

699 DiaSource catalog with updated DiaObject ids. 

700 - newDiaObjects : `pandas.DataFrame` 

701 Newly created DiaObjects from the unassociated DiaSources. 

702 - nNewDiaObjects : `int` 

703 Number of newly created diaObjects. 

704 - marginalDiaSources : `pandas.DataFrame` 

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

706 reliability, which are excluded from the new DiaObjects. 

707 """ 

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

709 if len(unassociatedDiaSources) == 0: 

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

711 else: 

712 if self.config.filterUnAssociatedSources: 

713 results = self.filterSources(unassociatedDiaSources) 

714 unassociatedDiaSources = results.goodSources 

715 marginalDiaSources = results.badSources 

716 unassociatedDiaSources["diaObjectId"] = unassociatedDiaSources["diaSourceId"] 

717 newDiaObjects = convertDataFrameToSdmSchema(self.schema, unassociatedDiaSources, 

718 tableName="DiaObject") 

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

720 return pipeBase.Struct(diaSources=unassociatedDiaSources, 

721 newDiaObjects=newDiaObjects, 

722 nNewDiaObjects=len(newDiaObjects), 

723 marginalDiaSources=marginalDiaSources) 

724 

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

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

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

728 

729 Parameters 

730 ---------- 

731 sources : `pandas.DataFrame` 

732 Set of DiaSources to check. 

733 snrThreshold : `float`, optional 

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

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

736 lowReliabilitySnrThreshold : `float`, optional 

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

738 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

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

740 reliabilityThreshold : `float`, optional 

741 The minimum reliability score diaSource to make a new diaObject 

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

743 lowSnrReliabilityThreshold : `float`, optional 

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

745 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

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

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

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

749 flags set. 

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

751 

752 Returns 

753 ------- 

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

755 Results struct containing: 

756 

757 - goodSources : `pandas.DataFrame` 

758 Subset of the input sources that pass all checks. 

759 - badSources : `pandas.DataFrame` 

760 Subset of the input sources that fail any check. 

761 """ 

762 if snrThreshold is None: 

763 snrThreshold = self.config.newObjectSnrThreshold 

764 if lowReliabilitySnrThreshold is None: 

765 lowReliabilitySnrThreshold = self.config.newObjectLowReliabilitySnrThreshold 

766 if reliabilityThreshold is None: 

767 reliabilityThreshold = self.config.newObjectReliabilityThreshold 

768 if lowSnrReliabilityThreshold is None: 

769 lowSnrReliabilityThreshold = self.config.newObjectLowSnrReliabilityThreshold 

770 if badFlags is None: 

771 badFlags = self.config.newObjectBadFlags 

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

773 fluxField = self.config.newObjectFluxField 

774 fluxErrField = fluxField + "Err" 

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

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

777 nFlagged = np.count_nonzero(flagged) 

778 if nFlagged > 0: 

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

780 

781 if snrThreshold > 0: 

782 snr_flag = signalToNoise < snrThreshold 

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

784 " signal to noise < %f", 

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

786 flagged += snr_flag 

787 if reliabilityThreshold > 0: 

788 reliability_flag = reliability < reliabilityThreshold 

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

790 "reliability<%f", 

791 np.sum(reliability_flag), reliabilityThreshold) 

792 flagged += reliability_flag 

793 if min(lowReliabilitySnrThreshold, lowSnrReliabilityThreshold) > 0: 

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

795 lowSnrReliability_flag = ((signalToNoise < lowReliabilitySnrThreshold) 

796 & (reliability < lowSnrReliabilityThreshold)) 

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

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

799 np.sum(lowSnrReliability_flag), 

800 fluxField, 

801 lowReliabilitySnrThreshold, 

802 lowSnrReliabilityThreshold) 

803 flagged += lowSnrReliability_flag 

804 

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

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

807 else: 

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

809 if np.count_nonzero(flagged) > 0: 

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

811 else: 

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

813 return pipeBase.Struct(goodSources=goodSources, 

814 badSources=badSources 

815 ) 

816 

817 @timeMethod 

818 def associateDiaSources(self, diaSourceTable, solarSystemObjectTable, diffIm, diaObjects): 

819 """Associate DiaSources with DiaObjects. 

820 

821 Associate new DiaSources with existing DiaObjects. Create new 

822 DiaObjects fron unassociated DiaSources. Index DiaSource catalogue 

823 after associations. Append new DiaObjects and DiaSources to their 

824 previous history. Test for DiaSource and DiaObject duplications. 

825 Compute DiaObject Summary statistics from their full DiaSource 

826 history. Test for duplication in the updated DiaObjects. 

827 

828 Parameters 

829 ---------- 

830 diaSourceTable : `pandas.DataFrame` 

831 Newly detected DiaSources. 

832 solarSystemObjectTable : `astropy.table.Table` 

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

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

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

836 were detected. 

837 diaObjects : `pandas.DataFrame` 

838 Table of DiaObjects from preloaded DiaObjects. 

839 

840 Returns 

841 ------- 

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

843 Results struct containing: 

844 

845 - associatedDiaSources : `pandas.DataFrame` 

846 Associated DiaSources with DiaObjects. 

847 - newDiaObjects : `pandas.DataFrame` 

848 Table of new DiaObjects after association. 

849 - associatedSsSources : `astropy.table.Table` 

850 Table of new ssSources after association. 

851 - unassociatedSsObjects : `astropy.table.Table` 

852 Table of expected ssSources that were not associated with a 

853 diaSource. 

854 - newDiaSources : `pandas.DataFrame` 

855 Subset of `associatedDiaSources` consisting of only the 

856 unassociated diaSources that were added as new diaObjects. 

857 - marginalDiaSources : `pandas.DataFrame` 

858 Unassociated diaSources with marginal detections, which were 

859 removed from `associatedDiaSources` and were not added as new 

860 diaObjects. 

861 """ 

862 # Associate new DiaSources with existing DiaObjects. 

863 assocResults = self.associator.run(diaSourceTable, diaObjects) 

864 # TODO: work out how to mark asteroid-diaObj overlaps as contaminated 

865 toAssociate = [] 

866 if self.config.doSolarSystemAssociation and solarSystemObjectTable is not None: 

867 ssoAssocResult = self.solarSystemAssociator.run( 

868 Table.from_pandas(assocResults.unAssocDiaSources), 

869 solarSystemObjectTable, 

870 diffIm.visitInfo, 

871 diffIm.getBBox(), 

872 diffIm.wcs 

873 ) 

874 # Create new DiaObjects from unassociated diaSources. 

875 if len(ssoAssocResult.unAssocDiaSources) > 0: 

876 # If the table is empty then converting time fields to pandas 

877 # will raise an error. Pass in an empty Dataframe in that case. 

878 createResults = self.createNewDiaObjects(ssoAssocResult.unAssocDiaSources.to_pandas()) 

879 else: 

880 createResults = self.createNewDiaObjects(pd.DataFrame()) 

881 if len(ssoAssocResult.ssoAssocDiaSources) > 0: 

882 toAssociate.append(ssoAssocResult.ssoAssocDiaSources.to_pandas()) 

883 nTotalSsObjects = ssoAssocResult.nTotalSsObjects 

884 nAssociatedSsObjects = ssoAssocResult.nAssociatedSsObjects 

885 associatedSsSources = ssoAssocResult.associatedSsSources 

886 unassociatedSsObjects = ssoAssocResult.unassociatedSsObjects 

887 

888 else: 

889 # Create new DiaObjects from unassociated diaSources. 

890 createResults = self.createNewDiaObjects(assocResults.unAssocDiaSources) 

891 nTotalSsObjects = 0 

892 nAssociatedSsObjects = 0 

893 associatedSsSources = None 

894 unassociatedSsObjects = None 

895 if not assocResults.matchedDiaSources.empty: 

896 toAssociate.append(assocResults.matchedDiaSources) 

897 if not createResults.diaSources.empty: 

898 toAssociate.append(createResults.diaSources) 

899 if len(toAssociate) == 0: 

900 associatedDiaSources = make_empty_catalog(self.schema, tableName="DiaSource") 

901 else: 

902 associatedDiaSources = pd.concat(toAssociate) 

903 

904 self._add_association_meta_data(assocResults.nUpdatedDiaObjects, 

905 assocResults.nUnassociatedDiaObjects, 

906 createResults.nNewDiaObjects, 

907 nTotalSsObjects, 

908 nAssociatedSsObjects) 

909 self.log.info("%i updated and %i unassociated diaObjects. Creating %i new diaObjects" 

910 " and dropping %i marginal diaSources.", 

911 assocResults.nUpdatedDiaObjects, 

912 assocResults.nUnassociatedDiaObjects, 

913 createResults.nNewDiaObjects, 

914 len(createResults.marginalDiaSources), 

915 ) 

916 if createResults.nNewDiaObjects > self.config.maxNewDiaObjects > 0: