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

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, FindCosmicRaysConfig, findCosmicRays 

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, LengthError 

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 TooManyCosmicRays(pipeBase.AlgorithmError): 

85 """Raised if the cosmic ray task fails with too many cosmics. 

86 

87 Parameters 

88 ---------- 

89 maxCosmicRays : `int` 

90 Maximum number of cosmic rays allowed. 

91 """ 

92 def __init__(self, maxCosmicRays, **kwargs): 

93 msg = f"Cosmic ray task found more than {maxCosmicRays} cosmics." 

94 self.msg = msg 

95 self._metadata = kwargs 

96 super().__init__(msg, **kwargs) 

97 self._metadata["maxCosmicRays"] = maxCosmicRays 

98 

99 def __str__(self): 

100 # Exception doesn't handle **kwargs, so we need a custom str. 

101 return f"{self.msg}: {self.metadata}" 

102 

103 @property 

104 def metadata(self): 

105 for key, value in self._metadata.items(): 

106 if not (isinstance(value, int) or isinstance(value, float) or isinstance(value, str)): 

107 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.") 

108 return self._metadata 

109 

110 

111class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, 

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

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

114 "warpTypeSuffix": "", 

115 "fakesType": ""}): 

116 science = pipeBase.connectionTypes.Input( 

117 doc="Input science exposure.", 

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

119 storageClass="ExposureF", 

120 name="{fakesType}calexp" 

121 ) 

122 matchedTemplate = pipeBase.connectionTypes.Input( 

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

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

125 storageClass="ExposureF", 

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

127 ) 

128 difference = pipeBase.connectionTypes.Input( 

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

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

131 storageClass="ExposureF", 

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

133 ) 

134 kernelSources = pipeBase.connectionTypes.Input( 

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

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

137 storageClass="SourceCatalog", 

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

139 ) 

140 outputSchema = pipeBase.connectionTypes.InitOutput( 

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

142 storageClass="SourceCatalog", 

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

144 ) 

145 diaSources = pipeBase.connectionTypes.Output( 

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

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

148 storageClass="SourceCatalog", 

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

150 ) 

151 subtractedMeasuredExposure = pipeBase.connectionTypes.Output( 

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

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

154 storageClass="ExposureF", 

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

156 ) 

157 differenceBackground = pipeBase.connectionTypes.Output( 

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

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

160 storageClass="Background", 

161 name="difference_background", 

162 ) 

163 maskedStreaks = pipeBase.connectionTypes.Output( 

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

165 storageClass="ArrowNumpyDict", 

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

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

168 ) 

169 glintTrailInfo = pipeBase.connectionTypes.Output( 

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

171 storageClass="ArrowNumpyDict", 

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

173 name="trailed_glints", 

174 ) 

175 

176 def __init__(self, *, config): 

177 super().__init__(config=config) 

178 if not (self.config.doCalculateResidualMetics): 

179 self.inputs.remove("kernelSources") 

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

181 self.outputs.remove("maskedStreaks") 

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

183 self.outputs.remove("differenceBackground") 

184 if not (self.config.writeGlintInfo): 

185 self.outputs.remove("glintTrailInfo") 

186 

187 

188class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, 

189 pipelineConnections=DetectAndMeasureConnections): 

190 """Config for DetectAndMeasureTask 

191 """ 

192 doMerge = pexConfig.Field( 

193 dtype=bool, 

194 default=True, 

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

196 "set by growFootprint" 

197 ) 

198 doForcedMeasurement = pexConfig.Field( 

199 dtype=bool, 

200 default=True, 

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

202 ) 

203 doAddMetrics = pexConfig.Field( 

204 dtype=bool, 

205 default=False, 

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

207 ) 

208 doSubtractBackground = pexConfig.Field( 

209 dtype=bool, 

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

211 default=True, 

212 ) 

213 doWriteBackground = pexConfig.Field( 

214 dtype=bool, 

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

216 default=False, 

217 ) 

218 doCalculateResidualMetics = pexConfig.Field( 

219 dtype=bool, 

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

221 "metadata?", 

222 default=True, 

223 ) 

224 subtractInitialBackground = pexConfig.ConfigurableField( 

225 target=lsst.meas.algorithms.SubtractBackgroundTask, 

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

227 ) 

228 subtractFinalBackground = pexConfig.ConfigurableField( 

229 target=lsst.meas.algorithms.SubtractBackgroundTask, 

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

231 ) 

232 doFindCosmicRays = pexConfig.Field( 

233 dtype=bool, 

234 doc="Detect and mask cosmic rays on the difference image?" 

235 "CRs can be interpolated over by setting cosmicray.keepCRs=False", 

236 default=True, 

237 ) 

238 cosmicray = pexConfig.ConfigField( 

239 dtype=FindCosmicRaysConfig, 

240 doc="Options for finding and masking cosmic rays", 

241 ) 

242 detection = pexConfig.ConfigurableField( 

243 target=SourceDetectionTask, 

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

245 ) 

246 streakDetection = pexConfig.ConfigurableField( 

247 target=SourceDetectionTask, 

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

249 ) 

250 doDeblend = pexConfig.Field( 

251 dtype=bool, 

252 default=False, 

253 doc="Deblend DIASources after detection?" 

254 ) 

255 deblend = pexConfig.ConfigurableField( 

256 target=lsst.meas.deblender.SourceDeblendTask, 

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

258 ) 

259 measurement = pexConfig.ConfigurableField( 

260 target=DipoleFitTask, 

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

262 ) 

263 doApCorr = lsst.pex.config.Field( 

264 dtype=bool, 

265 default=True, 

266 doc="Run subtask to apply aperture corrections" 

267 ) 

268 applyApCorr = lsst.pex.config.ConfigurableField( 

269 target=ApplyApCorrTask, 

270 doc="Task to apply aperture corrections" 

271 ) 

272 forcedMeasurement = pexConfig.ConfigurableField( 

273 target=ForcedMeasurementTask, 

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

275 ) 

276 growFootprint = pexConfig.Field( 

277 dtype=int, 

278 default=2, 

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

280 ) 

281 diaSourceMatchRadius = pexConfig.Field( 

282 dtype=float, 

283 default=0.5, 

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

285 ) 

286 doSkySources = pexConfig.Field( 

287 dtype=bool, 

288 default=False, 

289 doc="Generate sky sources?", 

290 ) 

291 skySources = pexConfig.ConfigurableField( 

292 target=SkyObjectsTask, 

293 doc="Generate sky sources", 

294 ) 

295 doMaskStreaks = pexConfig.Field( 

296 dtype=bool, 

297 default=True, 

298 doc="Turn on streak masking", 

299 ) 

300 maskStreaks = pexConfig.ConfigurableField( 

301 target=MaskStreaksTask, 

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

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

304 ) 

305 streakBinFactor = pexConfig.Field( 

306 dtype=int, 

307 default=4, 

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

309 "Only used if doMaskStreaks is True.", 

310 ) 

311 writeStreakInfo = pexConfig.Field( 

312 dtype=bool, 

313 default=False, 

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

315 "development work." 

316 ) 

317 findGlints = pexConfig.ConfigurableField( 

318 target=FindGlintTrailsTask, 

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

320 ) 

321 writeGlintInfo = pexConfig.Field( 

322 dtype=bool, 

323 default=True, 

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

325 ) 

326 setPrimaryFlags = pexConfig.ConfigurableField( 

327 target=SetPrimaryFlagsTask, 

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

329 ) 

330 badSourceFlags = lsst.pex.config.ListField( 

331 dtype=str, 

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

333 default=("base_PixelFlags_flag_offimage", 

334 "base_PixelFlags_flag_interpolatedCenterAll", 

335 "base_PixelFlags_flag_badCenter", 

336 "base_PixelFlags_flag_edgeCenter", 

337 "base_PixelFlags_flag_nodataCenter", 

338 "base_PixelFlags_flag_saturatedCenter", 

339 "base_PixelFlags_flag_sat_templateCenter", 

340 "base_PixelFlags_flag_spikeCenter", 

341 ), 

342 ) 

343 clearMaskPlanes = lsst.pex.config.ListField( 

344 dtype=str, 

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

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

347 ) 

348 doRejectBadMaskPlaneDetections = pexConfig.Field( 

349 dtype=bool, 

350 default=True, 

351 doc="Reject any peaks detected on ``badMaskPlanes`` before measurement." 

352 "These should all be rejected downstream by ``badSourceFlags``, " 

353 "but filtering earlier can save time." 

354 ) 

355 badMaskPlanes = lsst.pex.config.ListField( 

356 dtype=str, 

357 doc="Detections whose footprint peak lies on a pixel with any of these" 

358 " mask planes set will be rejected before measurement." 

359 " Any missing mask planes will be silently ignored.", 

360 default=("NO_DATA", "BAD", "SAT", "EDGE"), 

361 ) 

362 raiseOnBadSubtractionRatio = pexConfig.Field( 

363 dtype=bool, 

364 default=True, 

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

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

367 " image exceeds ``badSubtractionRatioThreshold``", 

368 ) 

369 badSubtractionRatioThreshold = pexConfig.Field( 

370 dtype=float, 

371 default=0.2, 

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

373 " the same footprints on the science image." 

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

375 ) 

376 badSubtractionVariationThreshold = pexConfig.Field( 

377 dtype=float, 

378 default=0.4, 

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

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

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

382 ) 

383 raiseOnNoDiaSources = pexConfig.Field( 

384 dtype=bool, 

385 default=True, 

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

387 ) 

388 run_sattle = pexConfig.Field( 

389 dtype=bool, 

390 default=False, 

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

392 "to the sattle service." 

393 ) 

394 sattle_historical = pexConfig.Field( 

395 dtype=bool, 

396 default=False, 

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

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

399 "closest in time to the exposure." 

400 ) 

401 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

402 

403 def setDefaults(self): 

404 # Background subtraction 

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

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

407 # faint diaSources 

408 self.subtractInitialBackground.binSize = 8 

409 self.subtractInitialBackground.useApprox = False 

410 self.subtractInitialBackground.statisticsProperty = "MEDIAN" 

411 self.subtractInitialBackground.doFilterSuperPixels = True 

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

413 "EDGE", 

414 "DETECTED", 

415 "DETECTED_NEGATIVE", 

416 "NO_DATA", 

417 ] 

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

419 # over-subtraction of bright objects. 

420 self.subtractFinalBackground.binSize = 40 

421 self.subtractFinalBackground.useApprox = False 

422 self.subtractFinalBackground.statisticsProperty = "MEDIAN" 

423 self.subtractFinalBackground.doFilterSuperPixels = True 

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

425 "EDGE", 

426 "DETECTED", 

427 "DETECTED_NEGATIVE", 

428 "NO_DATA", 

429 ] 

430 # Cosmic ray detection 

431 self.cosmicray.keepCRs = True # do not interpolate over detected CRs 

432 # DiaSource Detection 

433 self.detection.thresholdPolarity = "both" 

434 self.detection.thresholdValue = 5.0 

435 self.detection.reEstimateBackground = False 

436 self.detection.thresholdType = "pixel_stdev" 

437 self.detection.excludeMaskPlanes = [] 

438 

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

440 self.streakDetection.thresholdType = self.detection.thresholdType 

441 self.streakDetection.reEstimateBackground = False 

442 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes 

443 self.streakDetection.thresholdValue = self.detection.thresholdValue 

444 # Only detect positive streaks 

445 self.streakDetection.thresholdPolarity = "positive" 

446 # Do not grow detected mask for streaks 

447 self.streakDetection.nSigmaToGrow = 0 

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

449 # detected mask is set. 

450 self.maskStreaks.onlyMaskDetected = False 

451 # Restrict streak masking from growing too large 

452 self.maskStreaks.maxStreakWidth = 100 

453 # Restrict the number of iterations allowed for fitting streaks 

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

455 self.maskStreaks.maxFitIter = 10 

456 # Only mask to 2 sigma in width 

457 self.maskStreaks.nSigmaMask = 2 

458 # Threshold for including streaks after the Hough Transform. 

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

460 self.maskStreaks.absMinimumKernelHeight = 2 

461 

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

463 "base_LocalPhotoCalib", 

464 "base_LocalWcs", 

465 "ext_shapeHSM_HsmSourceMoments", 

466 "ext_shapeHSM_HsmPsfMoments", 

467 "base_ClassificationSizeExtendedness", 

468 ] 

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

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

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

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

473 self.forcedMeasurement.copyColumns = { 

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

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

476 self.forcedMeasurement.slots.shape = None 

477 

478 # Keep track of which footprints contain streaks 

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

480 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"] 

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

482 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"] 

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

484 

485 def validate(self): 

486 super().validate() 

487 

488 if self.run_sattle: 

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

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

491 "Sattle requested but SATTLE_URI_BASE " 

492 "environment variable not set.") 

493 

494 

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

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

497 """ 

498 ConfigClass = DetectAndMeasureConfig 

499 _DefaultName = "detectAndMeasure" 

500 

501 def __init__(self, **kwargs): 

502 super().__init__(**kwargs) 

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

504 

505 self.algMetadata = dafBase.PropertyList() 

506 if self.config.doSubtractBackground: 

507 self.makeSubtask("subtractInitialBackground") 

508 self.makeSubtask("subtractFinalBackground") 

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

510 if self.config.doDeblend: 

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

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

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

514 algMetadata=self.algMetadata) 

515 if self.config.doApCorr: 515 ↛ 517line 515 didn't jump to line 517 because the condition on line 515 was always true

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

517 if self.config.doForcedMeasurement: 

518 self.schema.addField( 

519 "ip_diffim_forced_PsfFlux_instFlux", "D", 

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

521 units="count") 

522 self.schema.addField( 

523 "ip_diffim_forced_PsfFlux_instFluxErr", "D", 

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

525 units="count") 

526 self.schema.addField( 

527 "ip_diffim_forced_PsfFlux_area", "F", 

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

529 units="pixel") 

530 self.schema.addField( 

531 "ip_diffim_forced_PsfFlux_flag", "Flag", 

532 "Forced PSF flux general failure flag.") 

533 self.schema.addField( 

534 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag", 

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

536 self.schema.addField( 

537 "ip_diffim_forced_PsfFlux_flag_edge", "Flag", 

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

539 self.schema.addField( 

540 "ip_diffim_forced_template_PsfFlux_instFlux", "D", 

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

542 units="count") 

543 self.schema.addField( 

544 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D", 

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

546 units="count") 

547 self.schema.addField( 

548 "ip_diffim_forced_template_PsfFlux_area", "F", 

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

550 units="pixel") 

551 self.schema.addField( 

552 "ip_diffim_forced_template_PsfFlux_flag", "Flag", 

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

554 self.schema.addField( 

555 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag", 

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

557 self.schema.addField( 

558 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag", 

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

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

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

562 

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

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

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

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

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

568 if self.config.doMaskStreaks: 568 ↛ 571line 568 didn't jump to line 571 because the condition on line 568 was always true

569 self.makeSubtask("maskStreaks") 

570 self.makeSubtask("streakDetection") 

571 self.makeSubtask("findGlints") 

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

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

574 self.schema.addField("reliabilityVersion", type=str, size=7, doc="Version of the reliability model") 

575 

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

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

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

579 

580 # Check that the schema and config are consistent 

581 for flag in self.config.badSourceFlags: 

582 if flag not in self.schema: 

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

584 

585 # initialize InitOutputs 

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

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

588 

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

590 inputRefs: pipeBase.InputQuantizedConnection, 

591 outputRefs: pipeBase.OutputQuantizedConnection): 

592 inputs = butlerQC.get(inputRefs) 

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

594 idFactory = idGenerator.make_table_id_factory() 

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

596 # exist, even as None. 

597 measurementResults = pipeBase.Struct( 

598 subtractedMeasuredExposure=None, 

599 diaSources=None, 

600 maskedStreaks=None, 

601 differenceBackground=None, 

602 ) 

603 try: 

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

605 except pipeBase.AlgorithmError as e: 

606 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

607 e, 

608 self, 

609 measurementResults.subtractedMeasuredExposure, 

610 measurementResults.diaSources, 

611 measurementResults.maskedStreaks, 

612 log=self.log 

613 ) 

614 butlerQC.put(measurementResults, outputRefs) 

615 raise error from e 

616 butlerQC.put(measurementResults, outputRefs) 

617 

618 @timeMethod 

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

620 idFactory=None, measurementResults=None): 

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

622 

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

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

625 Close positive and negative detections will optionally be merged into 

626 dipole diaSources. 

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

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

629 detections. 

630 

631 Parameters 

632 ---------- 

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

634 Science exposure that the template was subtracted from. 

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

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

637 difference image. 

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

639 Result of subtracting template from the science image. 

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

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

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

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

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

645 deblending and merging positive/negative peaks. 

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

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

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

649 this is also returned. 

650 

651 Returns 

652 ------- 

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

654 

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

656 Subtracted exposure with detection mask applied. 

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

658 The catalog of detected sources. 

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

660 Background that was subtracted from the difference image. 

661 """ 

662 if measurementResults is None: 662 ↛ 664line 662 didn't jump to line 664 because the condition on line 662 was always true

663 measurementResults = pipeBase.Struct() 

664 if idFactory is None: 

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

666 

667 # Check image properties and clear detection mask planes. 

668 self._prepareInputs(difference) 

669 if self.config.doSubtractBackground: 

670 background = self._fitAndSubtractBackground(differenceExposure=difference) 

671 else: 

672 background = afwMath.BackgroundList() 

673 

674 if self.config.doFindCosmicRays: 674 ↛ 680line 674 didn't jump to line 680 because the condition on line 674 was always true

675 self.findAndMaskCosmicRays(difference) 

676 

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

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

679 # know about past ids). 

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

681 results = self.detection.run( 

682 table=table, 

683 exposure=difference, 

684 doSmooth=True, 

685 background=background, 

686 clearMask=True, 

687 ) 

688 measurementResults.differenceBackground = background 

689 

690 if self.config.doDeblend: 

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

692 results.positive, 

693 results.negative) 

694 

695 else: 

696 positives = afwTable.SourceCatalog(self.schema) 

697 results.positive.makeSources(positives) 

698 negatives = afwTable.SourceCatalog(self.schema) 

699 results.negative.makeSources(negatives) 

700 sources = results.sources 

701 

702 self.processResults(science, matchedTemplate, difference, 

703 sources, idFactory, kernelSources, 

704 positives=positives, 

705 negatives=negatives, 

706 measurementResults=measurementResults) 

707 return measurementResults 

708 

709 def _prepareInputs(self, difference): 

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

711 

712 Parameters 

713 ---------- 

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

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

716 The mask plane will be modified in place. 

717 

718 Raises 

719 ------ 

720 lsst.pipe.base.UpstreamFailureNoWorkFound 

721 If the PSF is not usable for measurement. 

722 """ 

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

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

725 if np.isnan(sigma): 

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

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

728 mask = difference.mask 

729 for mp in self.config.clearMaskPlanes: 

730 if mp not in mask.getMaskPlaneDict(): 

731 mask.addMaskPlane(mp) 

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

733 

734 def _fitAndSubtractBackground(self, differenceExposure, scoreExposure=None): 

735 """Fit the final background to ``differenceExposure``. 

736 

737 A rough background is first fit on the difference and an internal 

738 first-pass detection is run on the (background-subtracted) clone; the 

739 resulting detection mask is then transferred to ``differenceExposure`` 

740 so the final background fit can mask out real sources. This two-pass 

741 scheme is what suppresses spurious detections from glints and 

742 extended structures: the small-binsize initial fit aggressively 

743 over-subtracts those features, and the resulting peak mask lets the 

744 large-binsize final fit produce a clean background. When 

745 ``scoreExposure`` is supplied (score-image variant), the final 

746 background is also subtracted from its image so the caller's final 

747 detection runs on a properly background-subtracted score. 

748 

749 The first-pass detection's `Struct` is discarded; only the mask it 

750 leaves on the clone is used downstream. 

751 

752 Parameters 

753 ---------- 

754 differenceExposure : `lsst.afw.image.ExposureF` 

755 Difference image. The final background is subtracted from its 

756 image in place but its mask is not modified. 

757 scoreExposure : `lsst.afw.image.ExposureF`, optional 

758 Maximum-likelihood score image. When supplied, the final 

759 background is also subtracted from its image. 

760 

761 Returns 

762 ------- 

763 background : `lsst.afw.math.BackgroundList` 

764 The fit final background. 

765 """ 

766 if scoreExposure is None: 

767 detectionExposure = differenceExposure.clone() 

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

769 doSmooth = True 

770 else: 

771 # Use a clone of differenceExposure because the background is 

772 # subtracted in place. 

773 background = self.subtractInitialBackground.run(differenceExposure.clone()).background 

774 detectionExposure = scoreExposure.clone() 

775 detectionExposure.image -= background.getImage() 

776 doSmooth = False 

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

778 self.detection.run( 

779 table=table, 

780 exposure=detectionExposure, 

781 doSmooth=doSmooth, 

782 background=background, 

783 clearMask=True, 

784 ) 

785 # Use the temporary detection mask for the final background subtraction. 

786 # The detection mask planes will be cleared before the final detection 

787 # step, so it is OK if they get set for differenceExposure. 

788 detectedBit = differenceExposure.mask.getPlaneBitMask(["DETECTED"]) 

789 detectedPix = detectionExposure.mask.array & detectedBit > 0 

790 detectedNegativeBit = differenceExposure.mask.getPlaneBitMask(["DETECTED_NEGATIVE"]) 

791 detectedNegativePix = detectionExposure.mask.array & detectedNegativeBit > 0 

792 differenceExposure.mask.array[detectedPix] |= detectedBit 

793 differenceExposure.mask.array[detectedNegativePix] |= detectedNegativeBit 

794 background = self.subtractFinalBackground.run(differenceExposure).background 

795 if scoreExposure is not None: 

796 # The preconvolution kernel is normalized to 1, so the same 

797 # background level applies to the difference and score images. 

798 scoreExposure.image -= background.getImage() 

799 return background 

800 

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

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

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

804 

805 Parameters 

806 ---------- 

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

808 Science exposure that the template was subtracted from. 

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

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

811 difference image. 

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

813 Result of subtracting template from the science image. 

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

815 Detected sources on the difference exposure. 

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

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

818 difference image. 

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

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

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

822 Positive polarity footprints. 

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

824 Negative polarity footprints. 

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

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

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

828 this is also returned. 

829 

830 Returns 

831 ------- 

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

833 

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

835 Subtracted exposure with detection mask applied. 

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

837 The catalog of detected sources. 

838 """ 

839 if measurementResults is None: 839 ↛ 840line 839 didn't jump to line 840 because the condition on line 839 was never true

840 measurementResults = pipeBase.Struct() 

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

842 if self.config.doMerge: 

843 # preserve peak schema, if there are any footprints 

844 if len(positives) > 0: 

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

846 elif len(negatives) > 0: 846 ↛ 849line 846 didn't jump to line 849 because the condition on line 846 was always true

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

848 else: 

849 peakSchema = afwDetection.PeakTable.makeMinimalSchema() 

850 mergeList = afwDetection.FootprintMergeList(self.schema, 

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

852 initialDiaSources = afwTable.SourceCatalog(self.schema) 

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

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

855 # deblended positive sources. 

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

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

858 mergeList.getFinalSources(initialDiaSources)