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 convertDataFrameToSdmSchema, checkSdmSchemaColumns, \ 

49 readSchemaFromApdb, dropEmptyColumns, make_empty_catalog, makeEmptyForcedSourceTable, getRegion, \ 

50 paddedRegion 

51from lsst.daf.base import DateTime 

52from lsst.meas.base import DetectorVisitIdGeneratorConfig, \ 

53 DiaObjectCalculationTask 

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=False, 

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 alertPackager = pexConfig.ConfigurableField( 

348 target=PackageAlertsTask, 

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

350 ) 

351 doPackageAlerts = pexConfig.Field( 

352 dtype=bool, 

353 default=False, 

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

355 "write them to disk.", 

356 ) 

357 doWriteAssociatedSources = pexConfig.Field( 

358 dtype=bool, 

359 default=True, 

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

361 "formatted following the Science Data Model.", 

362 ) 

363 imagePixelMargin = pexConfig.RangeField( 

364 dtype=int, 

365 default=10, 

366 min=0, 

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

368 "diaObjects for association.", 

369 ) 

370 filterUnAssociatedSources = pexConfig.Field( 

371 dtype=bool, 

372 default=True, 

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

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

375 "will not be affected." 

376 ) 

377 newObjectBadFlags = pexConfig.ListField( 

378 dtype=str, 

379 default=("centroid_flag", 

380 "pixelFlags_crCenter", 

381 "pixelFlags_nodataCenter", 

382 "pixelFlags_interpolatedCenter", 

383 "pixelFlags_saturatedCenter", 

384 "pixelFlags_suspectCenter", 

385 "pixelFlags_streakCenter", 

386 "glint_trail"), 

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

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

389 ) 

390 maxNewDiaObjects = pexConfig.RangeField( 

391 dtype=float, 

392 default=0, 

393 min=0, 

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

395 "Set to zero to disable.", 

396 ) 

397 newObjectSnrThreshold = pexConfig.RangeField( 

398 dtype=float, 

399 default=0, 

400 min=0, 

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

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

403 "diaObjects." 

404 "Set to zero to disable.", 

405 ) 

406 newObjectLowReliabilitySnrThreshold = pexConfig.RangeField( 

407 dtype=float, 

408 default=0, 

409 min=0, 

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

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

412 " low reliability scores before creating new diaObjects." 

413 "Set to zero to disable.", 

414 ) 

415 newObjectReliabilityThreshold = pexConfig.RangeField( 

416 dtype=float, 

417 default=0.1, 

418 min=0, 

419 max=1, 

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

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

422 "diaObjects." 

423 "Set to zero to disable.", 

424 ) 

425 newObjectLowSnrReliabilityThreshold = pexConfig.RangeField( 

426 dtype=float, 

427 default=0.1, 

428 min=0, 

429 max=1, 

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

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

432 "before creating new diaObjects." 

433 "Set to zero to disable.", 

434 ) 

435 newObjectFluxField = pexConfig.Field( 

436 dtype=str, 

437 default="apFlux", 

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

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

440 ) 

441 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

442 

443 def setDefaults(self): 

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

445 "ap_nDiaSources", 

446 "ap_meanFlux", 

447 "ap_sigmaFlux", 

448 "ap_minMaxFlux", 

449 "ap_maxSlopeFlux", 

450 "ap_meanErrFlux", 

451 "ap_meanTotFlux"] 

452 

453 

454class DiaPipelineTask(pipeBase.PipelineTask): 

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

456 (DIA) Objects and Sources. 

457 """ 

458 ConfigClass = DiaPipelineConfig 

459 _DefaultName = "diaPipe" 

460 

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

462 super().__init__(**kwargs) 

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

464 self.schema = readSchemaFromApdb(self.apdb) 

465 self.makeSubtask("associator") 

466 self.makeSubtask("diaCalculation") 

467 if self.config.doRunForcedMeasurement: 

468 self.makeSubtask("diaForcedSource") 

469 if self.config.doPackageAlerts: 

470 self.makeSubtask("alertPackager") 

471 if self.config.doSolarSystemAssociation: 

472 self.makeSubtask("solarSystemAssociator") 

473 if self.config.filterUnAssociatedSources: 

474 columns = [self.config.newObjectFluxField, 

475 self.config.newObjectFluxField + "Err", 

476 "reliability"] 

477 columns += self.config.newObjectBadFlags 

478 

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

480 if missing: 

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

482 

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

484 inputs = butlerQC.get(inputRefs) 

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

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

487 inputs["legacySolarSystemTable"] = None 

488 if not self.config.doSolarSystemAssociation: 

489 inputs["solarSystemObjectTable"] = None 

490 

491 outputs = self.run(**inputs) 

492 

493 butlerQC.put(outputs, outputRefs) 

494 

495 @timeMethod 

496 def run(self, 

497 diaSourceTable, 

498 legacySolarSystemTable, 

499 diffIm, 

500 exposure, 

501 template, 

502 preloadedDiaObjects, 

503 preloadedDiaSources, 

504 preloadedDiaForcedSources, 

505 band, 

506 idGenerator, 

507 solarSystemObjectTable=None): 

508 """Process DiaSources and DiaObjects. 

509 

510 Load previous DiaObjects and their DiaSource history. Calibrate the 

511 values in the diaSourceCat. Associate new DiaSources with previous 

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

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

514 

515 Parameters 

516 ---------- 

517 diaSourceTable : `pandas.DataFrame` 

518 Newly detected DiaSources. 

519 legacySolarSystemTable : `pandas.DataFrame` 

520 Not used 

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

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

523 were detected. 

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

525 Calibrated exposure differenced with a template to create 

526 ``diffIm``. 

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

528 Template exposure used to create diffIm. 

529 preloadedDiaObjects : `pandas.DataFrame` 

530 Previously detected DiaObjects, loaded from the APDB. 

531 preloadedDiaSources : `pandas.DataFrame` 

532 Previously detected DiaSources, loaded from the APDB. 

533 preloadedDiaForcedSources : `pandas.DataFrame` 

534 Catalog of previously detected forced DiaSources, from the APDB 

535 band : `str` 

536 The band in which the new DiaSources were detected. 

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

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

539 solarSystemObjectTable : `astropy.table.Table` 

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

541 

542 Returns 

543 ------- 

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

545 Results struct with components. 

546 

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

548 that this ccdVisit has completed successfully. 

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

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

551 DiaSources. (`pandas.DataFrame`) 

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

553 forced DiaSources. (`pandas.DataFrame`) 

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

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

556 (`pandas.DataFrame`) 

557 

558 Raises 

559 ------ 

560 RuntimeError 

561 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

562 """ 

563 # Accept either legacySolarSystemTable or optional solarSystemObjectTable. 

564 if legacySolarSystemTable is not None and solarSystemObjectTable is None: 

565 solarSystemObjectTable = Table.from_pandas(legacySolarSystemTable) 

566 if self.config.doReloadDiaObjects: 

567 try: 

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

569 except Exception as e: 

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

571 "the latest diaObjects from the APDB. Processing " 

572 "will continue with only the diaObjects from " 

573 "preload. Error:", e) 

574 finally: 

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

576 self.metadata, "loadRefreshedDiaObjects", clock="Utc") 

577 self.log.verbose("DiaObjects: Took %.4f seconds", self.metadata["loadDiaObjectsDuration"]) 

578 

579 else: 

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

581 

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

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

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

585 

586 if not preloadedDiaObjects.empty: 

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

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

589 buffer=self.config.imagePixelMargin) 

590 else: 

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

592 diaObjects = preloadedDiaObjects 

593 

594 # Associate DiaSources with DiaObjects 

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

596 

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

598 standardizedAssociatedDiaSources = self.standardizeDataFrame( 

599 assocResults.associatedDiaSources, 

600 "DiaSource", 

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

602 ) 

603 

604 # Merge associated diaSources 

605 mergedDiaSourceHistory, mergedDiaObjects, updatedDiaObjectIds = self.mergeAssociatedCatalogs( 

606 preloadedDiaSources, 

607 assocResults.associatedDiaSources, 

608 diaObjects, 

609 assocResults.newDiaObjects, 

610 diffIm 

611 ) 

612 

613 # Compute DiaObject Summary statistics from their full DiaSource 

614 # history. 

615 diaCalResult = self.diaCalculation.run( 

616 mergedDiaObjects, 

617 mergedDiaSourceHistory, 

618 updatedDiaObjectIds, 

619 [band]) 

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

621 tableName="DiaObject") 

622 

623 # Test for duplication in the updated DiaObjects. 

624 if self.testDataFrameIndex(diaCalResult.diaObjectCat): 

625 raise RuntimeError( 

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

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

628 "reported. Exiting.") 

629 if self.testDataFrameIndex(updatedDiaObjects): 

630 raise RuntimeError( 

631 "Duplicate DiaObjects (updated) created after " 

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

633 "reported. Exiting.") 

634 

635 # Forced source measurement 

636 if self.config.doRunForcedMeasurement: 

637 diaForcedSources = self.runForcedMeasurement( 

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

639 ) 

640 forcedSourceHistoryThreshold = self.diaForcedSource.config.historyThreshold 

641 else: 

642 # alertPackager needs correct columns 

643 diaForcedSources = makeEmptyForcedSourceTable(self.schema) 

644 forcedSourceHistoryThreshold = 0 

645 

646 # Write results to Alert Production Database (APDB) 

647 try: 

648 validityStart = self.writeToApdb(updatedDiaObjects, standardizedAssociatedDiaSources, 

649 diaForcedSources) 

650 finally: 

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

652 clock="Utc") 

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

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

655 

656 associatedSsSources = assocResults.associatedSsSources 

657 associatedSsSourcesPlusMpcorb = None 

658 if associatedSsSources is not None: 

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

660 associatedSsSourceMpcorb = associatedSsSources[mpcorbColumns].copy() 

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

662 associatedSsSourcesPlusMpcorb = associatedSsSources.copy() 

663 for mpcorbColumn in mpcorbColumns: 

664 associatedSsSourcesPlusMpcorb[mpcorbColumn] = associatedSsSourceMpcorb[mpcorbColumn] 

665 

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

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

668 

669 # Package alerts 

670 if self.config.doPackageAlerts: 

671 # Append new forced sources to the full history 

672 diaForcedSourcesFull = self.mergeCatalogs(preloadedDiaForcedSources, diaForcedSources, 

673 tableName="DiaForcedSource") 

674 if self.testDataFrameIndex(diaForcedSources): 

675 self.log.warning( 

676 "Duplicate DiaForcedSources created after merge with " 

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

678 "problems. Dropping duplicates.") 

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

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

681 # expected. 

682 diaForcedSourcesFull = diaForcedSourcesFull.groupby( 

683 diaForcedSourcesFull.index).first() 

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

685 diaForcedSourcesFull.set_index( 

686 ["diaObjectId", "diaForcedSourceId"], 

687 drop=False, 

688 inplace=True) 

689 

690 self.alertPackager.run(assocResults.associatedDiaSources, 

691 updatedDiaObjects, 

692 preloadedDiaSources, 

693 diaForcedSourcesFull, 

694 diffIm, 

695 exposure, 

696 template, 

697 ssSrc=associatedSsSourcesPlusMpcorb, 

698 doRunForcedMeasurement=self.config.doRunForcedMeasurement, 

699 forcedSourceHistoryThreshold=forcedSourceHistoryThreshold, 

700 ) 

701 

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

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

704 marker = pexConfig.Config() 

705 return pipeBase.Struct(apdbMarker=marker, 

706 associatedDiaSources=assocResults.associatedDiaSources, 

707 diaForcedSources=diaForcedSources, 

708 diaObjects=diaCalResult.diaObjectCat, 

709 associatedSsSources=associatedSsSources, 

710 unassociatedSsObjects=assocResults.unassociatedSsObjects, 

711 newDiaSources=assocResults.newDiaSources, 

712 marginalDiaSources=assocResults.marginalDiaSources 

713 ) 

714 

715 def createNewDiaObjects(self, unassociatedDiaSources): 

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

717 for unassociated DiaSources. 

718 

719 Parameters 

720 ---------- 

721 unassociatedDiaSources : `pandas.DataFrame` 

722 Set of DiaSources to create new DiaObjects from. 

723 

724 Returns 

725 ------- 

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

727 Results struct containing: 

728 

729 - diaSources : `pandas.DataFrame` 

730 DiaSource catalog with updated DiaObject ids. 

731 - newDiaObjects : `pandas.DataFrame` 

732 Newly created DiaObjects from the unassociated DiaSources. 

733 - nNewDiaObjects : `int` 

734 Number of newly created diaObjects. 

735 - marginalDiaSources : `pandas.DataFrame` 

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

737 reliability, which are excluded from the new DiaObjects. 

738 """ 

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

740 if len(unassociatedDiaSources) == 0: 

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

742 else: 

743 if self.config.filterUnAssociatedSources: 

744 results = self.filterSources(unassociatedDiaSources) 

745 unassociatedDiaSources = results.goodSources 

746 marginalDiaSources = results.badSources 

747 unassociatedDiaSources["diaObjectId"] = unassociatedDiaSources["diaSourceId"] 

748 newDiaObjects = convertDataFrameToSdmSchema(self.schema, unassociatedDiaSources, 

749 tableName="DiaObject") 

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

751 return pipeBase.Struct(diaSources=unassociatedDiaSources, 

752 newDiaObjects=newDiaObjects, 

753 nNewDiaObjects=len(newDiaObjects), 

754 marginalDiaSources=marginalDiaSources) 

755 

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

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

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

759 

760 Parameters 

761 ---------- 

762 sources : `pandas.DataFrame` 

763 Set of DiaSources to check. 

764 snrThreshold : `float`, optional 

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

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

767 lowReliabilitySnrThreshold : `float`, optional 

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

769 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

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

771 reliabilityThreshold : `float`, optional 

772 The minimum reliability score diaSource to make a new diaObject 

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

774 lowSnrReliabilityThreshold : `float`, optional 

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

776 diaSources with ``snrThreshold`` < SNR < ``lowReliabilitySnrThreshold`` 

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

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

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

780 flags set. 

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

782 

783 Returns 

784 ------- 

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

786 Results struct containing: 

787 

788 - goodSources : `pandas.DataFrame` 

789 Subset of the input sources that pass all checks. 

790 - badSources : `pandas.DataFrame` 

791 Subset of the input sources that fail any check. 

792 """ 

793 if snrThreshold is None: 

794 snrThreshold = self.config.newObjectSnrThreshold 

795 if lowReliabilitySnrThreshold is None: 

796 lowReliabilitySnrThreshold = self.config.newObjectLowReliabilitySnrThreshold 

797 if reliabilityThreshold is None: 

798 reliabilityThreshold = self.config.newObjectReliabilityThreshold 

799 if lowSnrReliabilityThreshold is None: 

800 lowSnrReliabilityThreshold = self.config.newObjectLowSnrReliabilityThreshold 

801 if badFlags is None: 

802 badFlags = self.config.newObjectBadFlags 

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

804 fluxField = self.config.newObjectFluxField 

805 fluxErrField = fluxField + "Err" 

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

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

808 nFlagged = np.count_nonzero(flagged) 

809 if nFlagged > 0: 

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

811 

812 if snrThreshold > 0: 

813 snr_flag = signalToNoise < snrThreshold 

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

815 " signal to noise < %f", 

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

817 flagged += snr_flag 

818 if reliabilityThreshold > 0: 

819 reliability_flag = reliability < reliabilityThreshold 

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

821 "reliability<%f", 

822 np.sum(reliability_flag), reliabilityThreshold) 

823 flagged += reliability_flag 

824 if min(lowReliabilitySnrThreshold, lowSnrReliabilityThreshold) > 0: 

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

826 lowSnrReliability_flag = ((signalToNoise < lowReliabilitySnrThreshold) 

827 & (reliability < lowSnrReliabilityThreshold)) 

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

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

830 np.sum(lowSnrReliability_flag), 

831 fluxField, 

832 lowReliabilitySnrThreshold, 

833 lowSnrReliabilityThreshold) 

834 flagged += lowSnrReliability_flag 

835 

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

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

838 else: 

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

840 if np.count_nonzero(flagged) > 0: 

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

842 else: 

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

844 return pipeBase.Struct(goodSources=goodSources, 

845 badSources=badSources 

846 ) 

847 

848 @timeMethod 

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

850 """Associate DiaSources with DiaObjects. 

851 

852 Associate new DiaSources with existing DiaObjects. Create new 

853 DiaObjects fron unassociated DiaSources. Index DiaSource catalogue 

854 after associations. Append new DiaObjects and DiaSources to their 

855 previous history. Test for DiaSource and DiaObject duplications. 

856 Compute DiaObject Summary statistics from their full DiaSource 

857 history. Test for duplication in the updated DiaObjects. 

858 

859 Parameters 

860 ---------- 

861 diaSourceTable : `pandas.DataFrame` 

862 Newly detected DiaSources. 

863 solarSystemObjectTable : `astropy.table.Table` 

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

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

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

867 were detected. 

868 diaObjects : `pandas.DataFrame` 

869 Table of DiaObjects from preloaded DiaObjects. 

870 

871 Returns 

872 ------- 

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

874 Results struct containing: 

875 

876 - associatedDiaSources : `pandas.DataFrame` 

877 Associated DiaSources with DiaObjects. 

878 - newDiaObjects : `pandas.DataFrame` 

879 Table of new DiaObjects after association. 

880 - associatedSsSources : `astropy.table.Table` 

881 Table of new ssSources after association. 

882 - unassociatedSsObjects : `astropy.table.Table` 

883 Table of expected ssSources that were not associated with a 

884 diaSource. 

885 - newDiaSources : `pandas.DataFrame` 

886 Subset of `associatedDiaSources` consisting of only the 

887 unassociated diaSources that were added as new diaObjects. 

888 - marginalDiaSources : `pandas.DataFrame` 

889 Unassociated diaSources with marginal detections, which were 

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

891 diaObjects. 

892 """ 

893 # Associate new DiaSources with existing DiaObjects. 

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

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

896 toAssociate = [] 

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

898 ssoAssocResult = self.solarSystemAssociator.run( 

899 Table.from_pandas(assocResults.unAssocDiaSources), 

900 solarSystemObjectTable, 

901 diffIm.visitInfo, 

902 diffIm.getBBox(), 

903 diffIm.wcs 

904 ) 

905 # Create new DiaObjects from unassociated diaSources. 

906 if len(ssoAssocResult.unAssocDiaSources) > 0: 

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

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

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

910 else: 

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

912 if len(ssoAssocResult.ssoAssocDiaSources) > 0: 

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

914 nTotalSsObjects = ssoAssocResult.nTotalSsObjects 

915 nAssociatedSsObjects = ssoAssocResult.nAssociatedSsObjects 

916 associatedSsSources = ssoAssocResult.associatedSsSources 

917 unassociatedSsObjects = ssoAssocResult.unassociatedSsObjects 

918