Coverage for python/lsst/ip/diffim/detectAndMeasure.py: 18%

1# This file is part of ip_diffim. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import numpy as np 

23import requests 

24import os 

25 

26import lsst.afw.detection as afwDetection 

27import lsst.afw.image as afwImage 

28import lsst.afw.math as afwMath 

29import lsst.afw.table as afwTable 

30import lsst.daf.base as dafBase 

31import lsst.geom 

32from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics, 

33 populate_sattle_visit_cache) 

34from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask 

35from lsst.meas.algorithms import FindGlintTrailsTask 

36from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig 

37import lsst.meas.deblender 

38import lsst.meas.extensions.trailedSources # noqa: F401 

39import lsst.meas.extensions.shapeHSM 

40import lsst.pex.config as pexConfig 

41from lsst.pex.exceptions import InvalidParameterError 

42import lsst.pipe.base as pipeBase 

43import lsst.utils 

44from lsst.utils.timer import timeMethod 

45 

46from . import DipoleFitTask 

47 

48__all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", 

49 "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] 

50 

51 

52class BadSubtractionError(pipeBase.AlgorithmError): 

53 """Raised when the residuals in footprints of stars used to compute the 

54 psf-matching kernel exceeds the configured maximum. 

55 """ 

56 def __init__(self, *, ratio, threshold): 

57 msg = ("The ratio of residual power in source footprints on the" 

58 " difference image to the power in the footprints on the" 

59 f" science image was {ratio}, which exceeds the maximum" 

60 f" threshold of {threshold}") 

61 super().__init__(msg) 

62 self.ratio = ratio 

63 self.threshold = threshold 

64 

65 @property 

66 def metadata(self): 

67 return {"ratio": self.ratio, 

68 "threshold": self.threshold 

69 } 

70 

71 

72class NoDiaSourcesError(pipeBase.AlgorithmError): 

73 """Raised when there are no diaSources detected on an image difference. 

74 """ 

75 def __init__(self): 

76 msg = ("No diaSources detected!") 

77 super().__init__(msg) 

78 

79 @property 

80 def metadata(self): 

81 return {} 

82 

83 

84class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, 

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

86 defaultTemplates={"coaddName": "deep", 

87 "warpTypeSuffix": "", 

88 "fakesType": ""}): 

89 science = pipeBase.connectionTypes.Input( 

90 doc="Input science exposure.", 

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

92 storageClass="ExposureF", 

93 name="{fakesType}calexp" 

94 ) 

95 matchedTemplate = pipeBase.connectionTypes.Input( 

96 doc="Warped and PSF-matched template used to create the difference image.", 

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

98 storageClass="ExposureF", 

99 name="{fakesType}{coaddName}Diff_matchedExp", 

100 ) 

101 difference = pipeBase.connectionTypes.Input( 

102 doc="Result of subtracting template from science.", 

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

104 storageClass="ExposureF", 

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

106 ) 

107 kernelSources = pipeBase.connectionTypes.Input( 

108 doc="Final selection of sources used for psf matching.", 

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

110 storageClass="SourceCatalog", 

111 name="{fakesType}{coaddName}Diff_psfMatchSources" 

112 ) 

113 outputSchema = pipeBase.connectionTypes.InitOutput( 

114 doc="Schema (as an example catalog) for output DIASource catalog.", 

115 storageClass="SourceCatalog", 

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

117 ) 

118 diaSources = pipeBase.connectionTypes.Output( 

119 doc="Detected diaSources on the difference image.", 

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

121 storageClass="SourceCatalog", 

122 name="{fakesType}{coaddName}Diff_diaSrc", 

123 ) 

124 subtractedMeasuredExposure = pipeBase.connectionTypes.Output( 

125 doc="Difference image with detection mask plane filled in.", 

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

127 storageClass="ExposureF", 

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

129 ) 

130 differenceBackground = pipeBase.connectionTypes.Output( 

131 doc="Background model that was subtracted from the difference image.", 

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

133 storageClass="Background", 

134 name="difference_background", 

135 ) 

136 maskedStreaks = pipeBase.connectionTypes.Output( 

137 doc='Catalog of streak fit parameters for the difference image.', 

138 storageClass="ArrowNumpyDict", 

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

140 name="{fakesType}{coaddName}Diff_streaks", 

141 ) 

142 glintTrailInfo = pipeBase.connectionTypes.Output( 

143 doc='Dict of fit parameters for glint trails in the catalog.', 

144 storageClass="ArrowNumpyDict", 

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

146 name="trailed_glints", 

147 ) 

148 

149 def __init__(self, *, config): 

150 super().__init__(config=config) 

151 if not (self.config.doCalculateResidualMetics): 

152 self.inputs.remove("kernelSources") 

153 if not (self.config.writeStreakInfo and self.config.doMaskStreaks): 

154 self.outputs.remove("maskedStreaks") 

155 if not (self.config.doSubtractBackground and self.config.doWriteBackground): 

156 self.outputs.remove("differenceBackground") 

157 if not (self.config.writeGlintInfo): 

158 self.outputs.remove("glintTrailInfo") 

159 

160 

161class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, 

162 pipelineConnections=DetectAndMeasureConnections): 

163 """Config for DetectAndMeasureTask 

164 """ 

165 doMerge = pexConfig.Field( 

166 dtype=bool, 

167 default=True, 

168 doc="Merge positive and negative diaSources with grow radius " 

169 "set by growFootprint" 

170 ) 

171 doForcedMeasurement = pexConfig.Field( 

172 dtype=bool, 

173 default=True, 

174 doc="Force photometer diaSource locations on PVI?" 

175 ) 

176 doAddMetrics = pexConfig.Field( 

177 dtype=bool, 

178 default=False, 

179 doc="Add columns to the source table to hold analysis metrics?" 

180 ) 

181 doSubtractBackground = pexConfig.Field( 

182 dtype=bool, 

183 doc="Subtract a background model from the image before detection?", 

184 default=True, 

185 ) 

186 doWriteBackground = pexConfig.Field( 

187 dtype=bool, 

188 doc="Persist the fitted background model?", 

189 default=False, 

190 ) 

191 doCalculateResidualMetics = pexConfig.Field( 

192 dtype=bool, 

193 doc="Calculate metrics to assess image subtraction quality for the task" 

194 "metadata?", 

195 default=True, 

196 ) 

197 subtractInitialBackground = pexConfig.ConfigurableField( 

198 target=lsst.meas.algorithms.SubtractBackgroundTask, 

199 doc="Task to perform intial background subtraction, before first detection pass.", 

200 ) 

201 subtractFinalBackground = pexConfig.ConfigurableField( 

202 target=lsst.meas.algorithms.SubtractBackgroundTask, 

203 doc="Task to perform final background subtraction, after first detection pass.", 

204 ) 

205 detection = pexConfig.ConfigurableField( 

206 target=SourceDetectionTask, 

207 doc="Final source detection for diaSource measurement", 

208 ) 

209 streakDetection = pexConfig.ConfigurableField( 

210 target=SourceDetectionTask, 

211 doc="Separate source detection used only for streak masking", 

212 ) 

213 doDeblend = pexConfig.Field( 

214 dtype=bool, 

215 default=False, 

216 doc="Deblend DIASources after detection?" 

217 ) 

218 deblend = pexConfig.ConfigurableField( 

219 target=lsst.meas.deblender.SourceDeblendTask, 

220 doc="Task to split blended sources into their components." 

221 ) 

222 measurement = pexConfig.ConfigurableField( 

223 target=DipoleFitTask, 

224 doc="Task to measure sources on the difference image.", 

225 ) 

226 doApCorr = lsst.pex.config.Field( 

227 dtype=bool, 

228 default=True, 

229 doc="Run subtask to apply aperture corrections" 

230 ) 

231 applyApCorr = lsst.pex.config.ConfigurableField( 

232 target=ApplyApCorrTask, 

233 doc="Task to apply aperture corrections" 

234 ) 

235 forcedMeasurement = pexConfig.ConfigurableField( 

236 target=ForcedMeasurementTask, 

237 doc="Task to force photometer science image at diaSource locations.", 

238 ) 

239 growFootprint = pexConfig.Field( 

240 dtype=int, 

241 default=2, 

242 doc="Grow positive and negative footprints by this many pixels before merging" 

243 ) 

244 diaSourceMatchRadius = pexConfig.Field( 

245 dtype=float, 

246 default=0.5, 

247 doc="Match radius (in arcseconds) for DiaSource to Source association" 

248 ) 

249 doSkySources = pexConfig.Field( 

250 dtype=bool, 

251 default=False, 

252 doc="Generate sky sources?", 

253 ) 

254 skySources = pexConfig.ConfigurableField( 

255 target=SkyObjectsTask, 

256 doc="Generate sky sources", 

257 ) 

258 doMaskStreaks = pexConfig.Field( 

259 dtype=bool, 

260 default=True, 

261 doc="Turn on streak masking", 

262 ) 

263 maskStreaks = pexConfig.ConfigurableField( 

264 target=MaskStreaksTask, 

265 doc="Subtask for masking streaks. Only used if doMaskStreaks is True. " 

266 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.", 

267 ) 

268 streakBinFactor = pexConfig.Field( 

269 dtype=int, 

270 default=4, 

271 doc="Bin scale factor to use when rerunning detection for masking streaks. " 

272 "Only used if doMaskStreaks is True.", 

273 ) 

274 writeStreakInfo = pexConfig.Field( 

275 dtype=bool, 

276 default=False, 

277 doc="Record the parameters of any detected streaks. For LSST, this should be turned off except for " 

278 "development work." 

279 ) 

280 findGlints = pexConfig.ConfigurableField( 

281 target=FindGlintTrailsTask, 

282 doc="Subtask for finding glint trails, usually caused by satellites or debris." 

283 ) 

284 writeGlintInfo = pexConfig.Field( 

285 dtype=bool, 

286 default=True, 

287 doc="Record the parameters of any detected glint trails." 

288 ) 

289 setPrimaryFlags = pexConfig.ConfigurableField( 

290 target=SetPrimaryFlagsTask, 

291 doc="Task to add isPrimary and deblending-related flags to the catalog." 

292 ) 

293 badSourceFlags = lsst.pex.config.ListField( 

294 dtype=str, 

295 doc="Sources with any of these flags set are removed before writing the output catalog.", 

296 default=("base_PixelFlags_flag_offimage", 

297 "base_PixelFlags_flag_edge", 

298 "base_PixelFlags_flag_interpolatedCenterAll", 

299 "base_PixelFlags_flag_badCenterAll", 

300 "base_PixelFlags_flag_edgeCenterAll", 

301 "base_PixelFlags_flag_nodataCenterAll", 

302 "base_PixelFlags_flag_saturatedCenterAll", 

303 "base_PixelFlags_flag_saturated_templateCenterAll", 

304 ), 

305 ) 

306 clearMaskPlanes = lsst.pex.config.ListField( 

307 dtype=str, 

308 doc="Mask planes to clear before running detection.", 

309 default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"), 

310 ) 

311 raiseOnBadSubtractionRatio = pexConfig.Field( 

312 dtype=bool, 

313 default=True, 

314 doc="Raise an error if the ratio of power in detected footprints" 

315 " on the difference image to the power in footprints on the science" 

316 " image exceeds ``badSubtractionRatioThreshold``", 

317 ) 

318 badSubtractionRatioThreshold = pexConfig.Field( 

319 dtype=float, 

320 default=0.2, 

321 doc="Maximum ratio of power in footprints on the difference image to" 

322 " the same footprints on the science image." 

323 "Only used if ``raiseOnBadSubtractionRatio`` is set", 

324 ) 

325 badSubtractionVariationThreshold = pexConfig.Field( 

326 dtype=float, 

327 default=0.4, 

328 doc="Maximum standard deviation of the ratio of power in footprints on" 

329 " the difference image to the same footprints on the science image." 

330 "Only used if ``raiseOnBadSubtractionRatio`` is set", 

331 ) 

332 raiseOnNoDiaSources = pexConfig.Field( 

333 dtype=bool, 

334 default=True, 

335 doc="Raise an algorithm error if no diaSources are detected.", 

336 ) 

337 run_sattle = pexConfig.Field( 

338 dtype=bool, 

339 default=False, 

340 doc="If true, dia source bounding boxes will be sent for verification" 

341 "to the sattle service." 

342 ) 

343 sattle_historical = pexConfig.Field( 

344 dtype=bool, 

345 default=False, 

346 doc="If re-running a pipeline that requires sattle, this should be set " 

347 "to True. This will populate sattle's cache with the historic data " 

348 "closest in time to the exposure." 

349 ) 

350 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

351 

352 def setDefaults(self): 

353 # Background subtraction 

354 # Use a small binsize for the first pass to reduce detections on glints 

355 # and extended structures. Should not affect the detectability of 

356 # faint diaSources 

357 self.subtractInitialBackground.binSize = 8 

358 self.subtractInitialBackground.useApprox = False 

359 self.subtractInitialBackground.statisticsProperty = "MEDIAN" 

360 self.subtractInitialBackground.doFilterSuperPixels = True 

361 self.subtractInitialBackground.ignoredPixelMask = ["BAD", 

362 "EDGE", 

363 "DETECTED", 

364 "DETECTED_NEGATIVE", 

365 "NO_DATA", 

366 ] 

367 # Use a larger binsize for the final background subtraction, to reduce 

368 # over-subtraction of bright objects. 

369 self.subtractFinalBackground.binSize = 40 

370 self.subtractFinalBackground.useApprox = False 

371 self.subtractFinalBackground.statisticsProperty = "MEDIAN" 

372 self.subtractFinalBackground.doFilterSuperPixels = True 

373 self.subtractFinalBackground.ignoredPixelMask = ["BAD", 

374 "EDGE", 

375 "DETECTED", 

376 "DETECTED_NEGATIVE", 

377 "NO_DATA", 

378 ] 

379 # DiaSource Detection 

380 self.detection.thresholdPolarity = "both" 

381 self.detection.thresholdValue = 5.0 

382 self.detection.reEstimateBackground = False 

383 self.detection.thresholdType = "pixel_stdev" 

384 self.detection.excludeMaskPlanes = [] 

385 

386 # Copy configs for binned streak detection from the base detection task 

387 self.streakDetection.thresholdType = self.detection.thresholdType 

388 self.streakDetection.reEstimateBackground = False 

389 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes 

390 self.streakDetection.thresholdValue = self.detection.thresholdValue 

391 # Only detect positive streaks 

392 self.streakDetection.thresholdPolarity = "positive" 

393 # Do not grow detected mask for streaks 

394 self.streakDetection.nSigmaToGrow = 0 

395 # Set the streak mask along the entire fit line, not only where the 

396 # detected mask is set. 

397 self.maskStreaks.onlyMaskDetected = False 

398 # Restrict streak masking from growing too large 

399 self.maskStreaks.maxStreakWidth = 100 

400 # Restrict the number of iterations allowed for fitting streaks 

401 # When the fit is good it should solve quickly, and exit a bad fit quickly 

402 self.maskStreaks.maxFitIter = 10 

403 # Only mask to 2 sigma in width 

404 self.maskStreaks.nSigmaMask = 2 

405 # Threshold for including streaks after the Hough Transform. 

406 # A lower value will detect more features that are less linear. 

407 self.maskStreaks.absMinimumKernelHeight = 2 

408 

409 self.measurement.plugins.names |= ["ext_trailedSources_Naive", 

410 "base_LocalPhotoCalib", 

411 "base_LocalWcs", 

412 "ext_shapeHSM_HsmSourceMoments", 

413 "ext_shapeHSM_HsmPsfMoments", 

414 "base_ClassificationSizeExtendedness", 

415 ] 

416 self.measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments" 

417 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments" 

418 self.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 5.0 

419 self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"] 

420 self.forcedMeasurement.copyColumns = { 

421 "id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"} 

422 self.forcedMeasurement.slots.centroid = "base_TransformedCentroid" 

423 self.forcedMeasurement.slots.shape = None 

424 

425 # Keep track of which footprints contain streaks 

426 self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [ 

427 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SATURATED_TEMPLATE"] 

428 self.measurement.plugins["base_PixelFlags"].masksFpCenter = [ 

429 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SATURATED_TEMPLATE"] 

430 self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"] 

431 

432 def validate(self): 

433 super().validate() 

434 

435 if self.run_sattle: 

436 if not os.getenv("SATTLE_URI_BASE"): 

437 raise pexConfig.FieldValidationError(DetectAndMeasureConfig.run_sattle, self, 

438 "Sattle requested but SATTLE_URI_BASE " 

439 "environment variable not set.") 

440 

441 

442class DetectAndMeasureTask(lsst.pipe.base.PipelineTask): 

443 """Detect and measure sources on a difference image. 

444 """ 

445 ConfigClass = DetectAndMeasureConfig 

446 _DefaultName = "detectAndMeasure" 

447 

448 def __init__(self, **kwargs): 

449 super().__init__(**kwargs) 

450 self.schema = afwTable.SourceTable.makeMinimalSchema() 

451 

452 self.algMetadata = dafBase.PropertyList() 

453 if self.config.doSubtractBackground: 

454 self.makeSubtask("subtractInitialBackground") 

455 self.makeSubtask("subtractFinalBackground") 

456 self.makeSubtask("detection", schema=self.schema) 

457 if self.config.doDeblend: 

458 self.makeSubtask("deblend", schema=self.schema) 

459 self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True) 

460 self.makeSubtask("measurement", schema=self.schema, 

461 algMetadata=self.algMetadata) 

462 if self.config.doApCorr: 

463 self.makeSubtask("applyApCorr", schema=self.measurement.schema) 

464 if self.config.doForcedMeasurement: 

465 self.schema.addField( 

466 "ip_diffim_forced_PsfFlux_instFlux", "D", 

467 "Forced PSF flux measured on the direct image.", 

468 units="count") 

469 self.schema.addField( 

470 "ip_diffim_forced_PsfFlux_instFluxErr", "D", 

471 "Forced PSF flux error measured on the direct image.", 

472 units="count") 

473 self.schema.addField( 

474 "ip_diffim_forced_PsfFlux_area", "F", 

475 "Forced PSF flux effective area of PSF.", 

476 units="pixel") 

477 self.schema.addField( 

478 "ip_diffim_forced_PsfFlux_flag", "Flag", 

479 "Forced PSF flux general failure flag.") 

480 self.schema.addField( 

481 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag", 

482 "Forced PSF flux not enough non-rejected pixels in data to attempt the fit.") 

483 self.schema.addField( 

484 "ip_diffim_forced_PsfFlux_flag_edge", "Flag", 

485 "Forced PSF flux object was too close to the edge of the image to use the full PSF model.") 

486 self.schema.addField( 

487 "ip_diffim_forced_template_PsfFlux_instFlux", "D", 

488 "Forced PSF flux measured on the template image.", 

489 units="count") 

490 self.schema.addField( 

491 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D", 

492 "Forced PSF flux error measured on the template image.", 

493 units="count") 

494 self.schema.addField( 

495 "ip_diffim_forced_template_PsfFlux_area", "F", 

496 "Forced template PSF flux effective area of PSF.", 

497 units="pixel") 

498 self.schema.addField( 

499 "ip_diffim_forced_template_PsfFlux_flag", "Flag", 

500 "Forced template PSF flux general failure flag.") 

501 self.schema.addField( 

502 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag", 

503 "Forced template PSF flux not enough non-rejected pixels in data to attempt the fit.") 

504 self.schema.addField( 

505 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag", 

506 """Forced template PSF flux object was too close to the edge of the image """ 

507 """to use the full PSF model.""") 

508 self.makeSubtask("forcedMeasurement", refSchema=self.schema) 

509 

510 self.schema.addField("refMatchId", "L", "unique id of reference catalog match") 

511 self.schema.addField("srcMatchId", "L", "unique id of source match") 

512 # Create the sky source task for use by metrics, 

513 # even if sky sources are not added to the diaSource catalog 

514 self.makeSubtask("skySources", schema=self.schema) 

515 if self.config.doMaskStreaks: 

516 self.makeSubtask("maskStreaks") 

517 self.makeSubtask("streakDetection") 

518 self.makeSubtask("findGlints") 

519 self.schema.addField("glint_trail", "Flag", "DiaSource is part of a glint trail.") 

520 self.schema.addField("reliability", type="F", doc="Reliability score of the DiaSource") 

521 

522 # To get the "merge_*" fields in the schema; have to re-initialize 

523 # this later, once we have a peak schema post-detection. 

524 lsst.afw.detection.FootprintMergeList(self.schema, ["positive", "negative"]) 

525 

526 # Check that the schema and config are consistent 

527 for flag in self.config.badSourceFlags: 

528 if flag not in self.schema: 

529 raise pipeBase.InvalidQuantumError("Field %s not in schema" % flag) 

530 

531 # initialize InitOutputs 

532 self.outputSchema = afwTable.SourceCatalog(self.schema) 

533 self.outputSchema.getTable().setMetadata(self.algMetadata) 

534 

535 def runQuantum(self, butlerQC: pipeBase.QuantumContext, 

536 inputRefs: pipeBase.InputQuantizedConnection, 

537 outputRefs: pipeBase.OutputQuantizedConnection): 

538 inputs = butlerQC.get(inputRefs) 

539 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) 

540 idFactory = idGenerator.make_table_id_factory() 

541 # Specify the fields that `annotate` needs below, to ensure they 

542 # exist, even as None. 

543 measurementResults = pipeBase.Struct( 

544 subtractedMeasuredExposure=None, 

545 diaSources=None, 

546 maskedStreaks=None, 

547 differenceBackground=None, 

548 ) 

549 try: 

550 self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults) 

551 except pipeBase.AlgorithmError as e: 

552 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

553 e, 

554 self, 

555 measurementResults.subtractedMeasuredExposure, 

556 measurementResults.diaSources, 

557 measurementResults.maskedStreaks, 

558 log=self.log 

559 ) 

560 butlerQC.put(measurementResults, outputRefs) 

561 raise error from e 

562 butlerQC.put(measurementResults, outputRefs) 

563 

564 @timeMethod 

565 def run(self, science, matchedTemplate, difference, kernelSources=None, 

566 idFactory=None, measurementResults=None): 

567 """Detect and measure sources on a difference image. 

568 

569 The difference image will be convolved with a gaussian approximation of 

570 the PSF to form a maximum likelihood image for detection. 

571 Close positive and negative detections will optionally be merged into 

572 dipole diaSources. 

573 Sky sources, or forced detections in background regions, will optionally 

574 be added, and the configured measurement algorithm will be run on all 

575 detections. 

576 

577 Parameters 

578 ---------- 

579 science : `lsst.afw.image.ExposureF` 

580 Science exposure that the template was subtracted from. 

581 matchedTemplate : `lsst.afw.image.ExposureF` 

582 Warped and PSF-matched template that was used produce the 

583 difference image. 

584 difference : `lsst.afw.image.ExposureF` 

585 Result of subtracting template from the science image. 

586 kernelSources : `lsst.afw.table.SourceCatalog`, optional 

587 Final selection of sources that was used for psf matching. 

588 idFactory : `lsst.afw.table.IdFactory`, optional 

589 Generator object used to assign ids to detected sources in the 

590 difference image. Ids from this generator are not set until after 

591 deblending and merging positive/negative peaks. 

592 measurementResults : `lsst.pipe.base.Struct`, optional 

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

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

595 this is also returned. 

596 

597 Returns 

598 ------- 

599 measurementResults : `lsst.pipe.base.Struct` 

600 

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

602 Subtracted exposure with detection mask applied. 

603 ``diaSources`` : `lsst.afw.table.SourceCatalog` 

604 The catalog of detected sources. 

605 ``differenceBackground`` : `lsst.afw.math.BackgroundList` 

606 Background that was subtracted from the difference image. 

607 """ 

608 if measurementResults is None: 

609 measurementResults = pipeBase.Struct() 

610 if idFactory is None: 

611 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory() 

612 

613 if self.config.doSubtractBackground: 

614 # Run background subtraction before clearing the mask planes 

615 detectionExposure = difference.clone() 

616 background = self.subtractInitialBackground.run(detectionExposure).background 

617 else: 

618 detectionExposure = difference 

619 background = afwMath.BackgroundList() 

620 

621 self._prepareInputs(detectionExposure) 

622 

623 # Don't use the idFactory until after deblend+merge, so that we aren't 

624 # generating ids that just get thrown away (footprint merge doesn't 

625 # know about past ids). 

626 table = afwTable.SourceTable.make(self.schema) 

627 results = self.detection.run( 

628 table=table, 

629 exposure=detectionExposure, 

630 doSmooth=True, 

631 background=background 

632 ) 

633 

634 if self.config.doSubtractBackground: 

635 # Run background subtraction again after detecting peaks 

636 # but before measurement 

637 # First update the mask using the detection image 

638 difference.setMask(detectionExposure.mask) 

639 background = self.subtractFinalBackground.run(difference).background 

640 

641 # Re-run detection to get final footprints 

642 table = afwTable.SourceTable.make(self.schema) 

643 results = self.detection.run( 

644 table=table, 

645 exposure=difference, 

646 doSmooth=True, 

647 background=background 

648 ) 

649 measurementResults.differenceBackground = background 

650 

651 if self.config.doDeblend: 

652 sources, positives, negatives = self._deblend(difference, 

653 results.positive, 

654 results.negative) 

655 

656 else: 

657 positives = afwTable.SourceCatalog(self.schema) 

658 results.positive.makeSources(positives) 

659 negatives = afwTable.SourceCatalog(self.schema) 

660 results.negative.makeSources(negatives) 

661 sources = results.sources 

662 

663 self.processResults(science, matchedTemplate, difference, 

664 sources, idFactory, kernelSources, 

665 positives=positives, 

666 negatives=negatives, 

667 measurementResults=measurementResults) 

668 return measurementResults 

669 

670 def _prepareInputs(self, difference): 

671 """Ensure that we start with an empty detection and deblended mask. 

672 

673 Parameters 

674 ---------- 

675 difference : `lsst.afw.image.ExposureF` 

676 The difference image that will be used for detecting diaSources. 

677 The mask plane will be modified in place. 

678 

679 Raises 

680 ------ 

681 lsst.pipe.base.UpstreamFailureNoWorkFound 

682 If the PSF is not usable for measurement. 

683 """ 

684 # Check that we have a valid PSF now before we do more work 

685 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius() 

686 if np.isnan(sigma): 

687 raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.") 

688 # Ensure that we start with an empty detection and deblended mask. 

689 mask = difference.mask 

690 for mp in self.config.clearMaskPlanes: 

691 if mp not in mask.getMaskPlaneDict(): 

692 mask.addMaskPlane(mp) 

693 mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes) 

694 

695 def processResults(self, science, matchedTemplate, difference, sources, idFactory, 

696 kernelSources=None, positives=None, negatives=None, measurementResults=None): 

697 """Measure and process the results of source detection. 

698 

699 Parameters 

700 ---------- 

701 science : `lsst.afw.image.ExposureF` 

702 Science exposure that the template was subtracted from. 

703 matchedTemplate : `lsst.afw.image.ExposureF` 

704 Warped and PSF-matched template that was used produce the 

705 difference image. 

706 difference : `lsst.afw.image.ExposureF` 

707 Result of subtracting template from the science image. 

708 sources : `lsst.afw.table.SourceCatalog` 

709 Detected sources on the difference exposure. 

710 idFactory : `lsst.afw.table.IdFactory` 

711 Generator object used to assign ids to detected sources in the 

712 difference image. 

713 kernelSources : `lsst.afw.table.SourceCatalog`, optional 

714 Final selection of sources that was used for psf matching. 

715 positives : `lsst.afw.table.SourceCatalog`, optional 

716 Positive polarity footprints. 

717 negatives : `lsst.afw.table.SourceCatalog`, optional 

718 Negative polarity footprints. 

719 measurementResults : `lsst.pipe.base.Struct`, optional 

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

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

722 this is also returned. 

723 

724 Returns 

725 ------- 

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

727 

728 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF` 

729 Subtracted exposure with detection mask applied. 

730 ``diaSources`` : `lsst.afw.table.SourceCatalog` 

731 The catalog of detected sources. 

732 """ 

733 if measurementResults is None: 

734 measurementResults = pipeBase.Struct() 

735 self.metadata["nUnmergedDiaSources"] = len(sources) 

736 if self.config.doMerge: 

737 # preserve peak schema, if there are any footprints 

738 if len(positives) > 0: 

739 peakSchema = positives[0].getFootprint().peaks.schema 

740 elif len(negatives) > 0: 

741 peakSchema = negatives[0].getFootprint().peaks.schema 

742 else: 

743 peakSchema = afwDetection.PeakTable.makeMinimalSchema() 

744 mergeList = afwDetection.FootprintMergeList(self.schema, 

745 ["positive", "negative"], peakSchema) 

746 initialDiaSources = afwTable.SourceCatalog(self.schema) 

747 # Start with positive, as FootprintMergeList will self-merge the 

748 # subsequent added catalogs, and we want to try to preserve 

749 # deblended positive sources. 

750 mergeList.addCatalog(initialDiaSources.table, positives, "positive", minNewPeakDist=0) 

751 mergeList.addCatalog(initialDiaSources.table, negatives, "negative", minNewPeakDist=0) 

752 mergeList.getFinalSources(initialDiaSources) 

753 # Flag as negative those sources that *only* came from the negative 

754 # footprint set. 

755 initialDiaSources["is_negative"] = initialDiaSources["merge_footprint_negative"] & \ 

756 ~initialDiaSources["merge_footprint_positive"] 

757 self.log.info("Merging detections into %d sources", len(initialDiaSources)) 

758 else: 

759 initialDiaSources = sources 

760 

761 # Assign source ids at the end: deblend/merge mean that we don't keep 

762 # track of parents and children, we only care about the final ids. 

763 for source in initialDiaSources: 

764 source.setId(idFactory()) 

765 # Ensure sources added after this get correct ids. 

766 initialDiaSources.getTable().setIdFactory(idFactory) 

767 initialDiaSources.setMetadata(self.algMetadata) 

768 

769 self.metadata["nMergedDiaSources"] = len(initialDiaSources) 

770 

771 if self.config.doMaskStreaks: 

772 streakInfo = self._runStreakMasking(difference) 

773 

774 if self.config.doSkySources: 

775 self.addSkySources(initialDiaSources, difference.mask, difference.info.id) 

776 

777 if not initialDiaSources.isContiguous(): 

778 initialDiaSources = initialDiaSources.copy(deep=True) 

779 

780 self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate) 

781 

782 # Remove unphysical diaSources per config.badSourceFlags 

783 diaSources = self._removeBadSources(initialDiaSources) 

784 

785 if self.config.run_sattle: 

786 diaSources = self.filterSatellites(diaSources, science) 

787 

788 # Flag diaSources in glint trails, but do not remove them 

789 diaSources, trail_parameters = self._find_glint_trails(diaSources) 

790 if self.config.writeGlintInfo: 

791 measurementResults.mergeItems(trail_parameters, 'glintTrailInfo') 

792 

793 if self.config.doForcedMeasurement: 

794 self.measureForcedSources(diaSources, science, difference.getWcs()) 

795 self.measureForcedSources(diaSources, matchedTemplate, difference.getWcs(), 

796 template=True) 

797 

798 # Clear the image plane for regions with NO_DATA. 

799 # These regions are most often caused by insufficient template coverage. 

800 # Do this for the final difference image after detection and measurement 

801 # since the subtasks should all be configured to handle NO_DATA properly 

802 difference.image.array[difference.mask.array & difference.mask.getPlaneBitMask('NO_DATA') > 0] = 0 

803 

804 measurementResults.subtractedMeasuredExposure = difference 

805 

806 if self.config.doMaskStreaks and self.config.writeStreakInfo: 

807 measurementResults.maskedStreaks = streakInfo.maskedStreaks 

808 

809 if kernelSources is not None: 

810 self.calculateMetrics(science, difference, diaSources, kernelSources) 

811 

812 if np.count_nonzero(~diaSources["sky_source"]) > 0: 

813 measurementResults.diaSources = diaSources 

814 elif self.config.raiseOnNoDiaSources: 

815 raise NoDiaSourcesError() 

816 elif len(diaSources) > 0: 

817 # This option allows returning sky sources, 

818 # even if there are no diaSources 

819 measurementResults.diaSources = diaSources 

820 self.log.info("Measured %d diaSources and %d sky sources", 

821 np.count_nonzero(~diaSources["sky_source"]), 

822 np.count_nonzero(diaSources["sky_source"]) 

823 ) 

824 return measurementResults 

825 

826 def _deblend(self, difference, positiveFootprints, negativeFootprints): 

827 """Deblend the positive and negative footprints and return a catalog 

828 containing just the children, and the deblended footprints. 

829 

830 Parameters 

831 ---------- 

832 difference : `lsst.afw.image.Exposure` 

833 Result of subtracting template from the science image. 

834 positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet` 

835 Positive and negative polarity footprints measured on 

836 ``difference`` to be deblended separately. 

837 

838 Returns 

839 ------- 

840 sources : `lsst.afw.table.SourceCatalog` 

841 Positive and negative deblended children. 

842 positives, negatives : `lsst.afw.table.SourceCatalog` 

843 Deblended positive and negative polarity sources with footprints 

844 detected on ``difference``. 

845 """ 

846 

847 def deblend(footprints, negative=False): 

848 """Deblend a positive or negative footprint set, 

849 and return the deblended children. 

850 

851 Parameters 

852 ---------- 

853 footprints : `lsst.afw.detection.FootprintSet` 

854 negative : `bool` 

855 Set True if the footprints contain negative fluxes 

856 

857 Returns 

858 ------- 

859 sources : `lsst.afw.table.SourceCatalog` 

860 """ 

861 sources = afwTable.SourceCatalog(self.schema) 

862 footprints.makeSources(sources)