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, 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_edge", 

335 "base_PixelFlags_flag_interpolatedCenterAll", 

336 "base_PixelFlags_flag_badCenter", 

337 "base_PixelFlags_flag_edgeCenter", 

338 "base_PixelFlags_flag_nodataCenter", 

339 "base_PixelFlags_flag_saturatedCenter", 

340 "base_PixelFlags_flag_saturated_templateCenter", 

341 "base_PixelFlags_flag_spikeCenter", 

342 ), 

343 ) 

344 clearMaskPlanes = lsst.pex.config.ListField( 

345 dtype=str, 

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

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

348 ) 

349 raiseOnBadSubtractionRatio = pexConfig.Field( 

350 dtype=bool, 

351 default=True, 

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

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

354 " image exceeds ``badSubtractionRatioThreshold``", 

355 ) 

356 badSubtractionRatioThreshold = pexConfig.Field( 

357 dtype=float, 

358 default=0.2, 

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

360 " the same footprints on the science image." 

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

362 ) 

363 badSubtractionVariationThreshold = pexConfig.Field( 

364 dtype=float, 

365 default=0.4, 

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

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

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

369 ) 

370 raiseOnNoDiaSources = pexConfig.Field( 

371 dtype=bool, 

372 default=True, 

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

374 ) 

375 run_sattle = pexConfig.Field( 

376 dtype=bool, 

377 default=False, 

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

379 "to the sattle service." 

380 ) 

381 sattle_historical = pexConfig.Field( 

382 dtype=bool, 

383 default=False, 

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

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

386 "closest in time to the exposure." 

387 ) 

388 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

389 

390 def setDefaults(self): 

391 # Background subtraction 

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

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

394 # faint diaSources 

395 self.subtractInitialBackground.binSize = 8 

396 self.subtractInitialBackground.useApprox = False 

397 self.subtractInitialBackground.statisticsProperty = "MEDIAN" 

398 self.subtractInitialBackground.doFilterSuperPixels = True 

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

400 "EDGE", 

401 "DETECTED", 

402 "DETECTED_NEGATIVE", 

403 "NO_DATA", 

404 ] 

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

406 # over-subtraction of bright objects. 

407 self.subtractFinalBackground.binSize = 40 

408 self.subtractFinalBackground.useApprox = False 

409 self.subtractFinalBackground.statisticsProperty = "MEDIAN" 

410 self.subtractFinalBackground.doFilterSuperPixels = True 

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

412 "EDGE", 

413 "DETECTED", 

414 "DETECTED_NEGATIVE", 

415 "NO_DATA", 

416 ] 

417 # Cosmic ray detection 

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

419 # DiaSource Detection 

420 self.detection.thresholdPolarity = "both" 

421 self.detection.thresholdValue = 5.0 

422 self.detection.reEstimateBackground = False 

423 self.detection.thresholdType = "pixel_stdev" 

424 self.detection.excludeMaskPlanes = [] 

425 

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

427 self.streakDetection.thresholdType = self.detection.thresholdType 

428 self.streakDetection.reEstimateBackground = False 

429 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes 

430 self.streakDetection.thresholdValue = self.detection.thresholdValue 

431 # Only detect positive streaks 

432 self.streakDetection.thresholdPolarity = "positive" 

433 # Do not grow detected mask for streaks 

434 self.streakDetection.nSigmaToGrow = 0 

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

436 # detected mask is set. 

437 self.maskStreaks.onlyMaskDetected = False 

438 # Restrict streak masking from growing too large 

439 self.maskStreaks.maxStreakWidth = 100 

440 # Restrict the number of iterations allowed for fitting streaks 

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

442 self.maskStreaks.maxFitIter = 10 

443 # Only mask to 2 sigma in width 

444 self.maskStreaks.nSigmaMask = 2 

445 # Threshold for including streaks after the Hough Transform. 

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

447 self.maskStreaks.absMinimumKernelHeight = 2 

448 

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

450 "base_LocalPhotoCalib", 

451 "base_LocalWcs", 

452 "ext_shapeHSM_HsmSourceMoments", 

453 "ext_shapeHSM_HsmPsfMoments", 

454 "base_ClassificationSizeExtendedness", 

455 ] 

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

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

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

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

460 self.forcedMeasurement.copyColumns = { 

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

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

463 self.forcedMeasurement.slots.shape = None 

464 

465 # Keep track of which footprints contain streaks 

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

467 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SATURATED_TEMPLATE", "SPIKE"] 

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

469 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SATURATED_TEMPLATE", "SPIKE"] 

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

471 

472 def validate(self): 

473 super().validate() 

474 

475 if self.run_sattle: 

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

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

478 "Sattle requested but SATTLE_URI_BASE " 

479 "environment variable not set.") 

480 

481 

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

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

484 """ 

485 ConfigClass = DetectAndMeasureConfig 

486 _DefaultName = "detectAndMeasure" 

487 

488 def __init__(self, **kwargs): 

489 super().__init__(**kwargs) 

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

491 

492 self.algMetadata = dafBase.PropertyList() 

493 if self.config.doSubtractBackground: 

494 self.makeSubtask("subtractInitialBackground") 

495 self.makeSubtask("subtractFinalBackground") 

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

497 if self.config.doDeblend: 

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

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

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

501 algMetadata=self.algMetadata) 

502 if self.config.doApCorr: 

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

504 if self.config.doForcedMeasurement: 

505 self.schema.addField( 

506 "ip_diffim_forced_PsfFlux_instFlux", "D", 

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

508 units="count") 

509 self.schema.addField( 

510 "ip_diffim_forced_PsfFlux_instFluxErr", "D", 

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

512 units="count") 

513 self.schema.addField( 

514 "ip_diffim_forced_PsfFlux_area", "F", 

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

516 units="pixel") 

517 self.schema.addField( 

518 "ip_diffim_forced_PsfFlux_flag", "Flag", 

519 "Forced PSF flux general failure flag.") 

520 self.schema.addField( 

521 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag", 

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

523 self.schema.addField( 

524 "ip_diffim_forced_PsfFlux_flag_edge", "Flag", 

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

526 self.schema.addField( 

527 "ip_diffim_forced_template_PsfFlux_instFlux", "D", 

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

529 units="count") 

530 self.schema.addField( 

531 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D", 

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

533 units="count") 

534 self.schema.addField( 

535 "ip_diffim_forced_template_PsfFlux_area", "F", 

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

537 units="pixel") 

538 self.schema.addField( 

539 "ip_diffim_forced_template_PsfFlux_flag", "Flag", 

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

541 self.schema.addField( 

542 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag", 

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

544 self.schema.addField( 

545 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag", 

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

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

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

549 

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

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

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

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

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

555 if self.config.doMaskStreaks: 

556 self.makeSubtask("maskStreaks") 

557 self.makeSubtask("streakDetection") 

558 self.makeSubtask("findGlints") 

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

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

561 

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

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

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

565 

566 # Check that the schema and config are consistent 

567 for flag in self.config.badSourceFlags: 

568 if flag not in self.schema: 

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

570 

571 # initialize InitOutputs 

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

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

574 

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

576 inputRefs: pipeBase.InputQuantizedConnection, 

577 outputRefs: pipeBase.OutputQuantizedConnection): 

578 inputs = butlerQC.get(inputRefs) 

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

580 idFactory = idGenerator.make_table_id_factory() 

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

582 # exist, even as None. 

583 measurementResults = pipeBase.Struct( 

584 subtractedMeasuredExposure=None, 

585 diaSources=None, 

586 maskedStreaks=None, 

587 differenceBackground=None, 

588 ) 

589 try: 

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

591 except pipeBase.AlgorithmError as e: 

592 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

593 e, 

594 self, 

595 measurementResults.subtractedMeasuredExposure, 

596 measurementResults.diaSources, 

597 measurementResults.maskedStreaks, 

598 log=self.log 

599 ) 

600 butlerQC.put(measurementResults, outputRefs) 

601 raise error from e 

602 butlerQC.put(measurementResults, outputRefs) 

603 

604 @timeMethod 

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

606 idFactory=None, measurementResults=None): 

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

608 

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

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

611 Close positive and negative detections will optionally be merged into 

612 dipole diaSources. 

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

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

615 detections. 

616 

617 Parameters 

618 ---------- 

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

620 Science exposure that the template was subtracted from. 

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

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

623 difference image. 

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

625 Result of subtracting template from the science image. 

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

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

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

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

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

631 deblending and merging positive/negative peaks. 

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

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

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

635 this is also returned. 

636 

637 Returns 

638 ------- 

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

640 

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

642 Subtracted exposure with detection mask applied. 

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

644 The catalog of detected sources. 

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

646 Background that was subtracted from the difference image. 

647 """ 

648 if measurementResults is None: 

649 measurementResults = pipeBase.Struct() 

650 if idFactory is None: 

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

652 

653 if self.config.doSubtractBackground: 

654 # Run background subtraction before clearing the mask planes 

655 detectionExposure = difference.clone() 

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

657 else: 

658 detectionExposure = difference 

659 background = afwMath.BackgroundList() 

660 

661 self._prepareInputs(detectionExposure) 

662 

663 if self.config.doFindCosmicRays and not self.config.doSubtractBackground: 

664 # Detect and interpolate over any remaining cosmic rays 

665 # This only needs to run once, right before the final detection. 

666 self.findAndMaskCosmicRays(detectionExposure) 

667 

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

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

670 # know about past ids). 

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

672 results = self.detection.run( 

673 table=table, 

674 exposure=detectionExposure, 

675 doSmooth=True, 

676 background=background 

677 ) 

678 

679 if self.config.doSubtractBackground: 

680 # Run background subtraction again after detecting peaks 

681 # but before measurement 

682 # First update the mask using the detection image 

683 difference.setMask(detectionExposure.mask) 

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

685 

686 if self.config.doFindCosmicRays: 

687 # Detect and interpolate over any remaining cosmic rays 

688 self.findAndMaskCosmicRays(difference) 

689 

690 # Re-run detection to get final footprints 

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

692 results = self.detection.run( 

693 table=table, 

694 exposure=difference, 

695 doSmooth=True, 

696 background=background 

697 ) 

698 measurementResults.differenceBackground = background 

699 

700 if self.config.doDeblend: 

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

702 results.positive, 

703 results.negative) 

704 

705 else: 

706 positives = afwTable.SourceCatalog(self.schema) 

707 results.positive.makeSources(positives) 

708 negatives = afwTable.SourceCatalog(self.schema) 

709 results.negative.makeSources(negatives) 

710 sources = results.sources 

711 

712 self.processResults(science, matchedTemplate, difference, 

713 sources, idFactory, kernelSources, 

714 positives=positives, 

715 negatives=negatives, 

716 measurementResults=measurementResults) 

717 return measurementResults 

718 

719 def _prepareInputs(self, difference): 

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

721 

722 Parameters 

723 ---------- 

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

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

726 The mask plane will be modified in place. 

727 

728 Raises 

729 ------ 

730 lsst.pipe.base.UpstreamFailureNoWorkFound 

731 If the PSF is not usable for measurement. 

732 """ 

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

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

735 if np.isnan(sigma): 

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

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

738 mask = difference.mask 

739 for mp in self.config.clearMaskPlanes: 

740 if mp not in mask.getMaskPlaneDict(): 

741 mask.addMaskPlane(mp) 

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

743 

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

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

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

747 

748 Parameters 

749 ---------- 

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

751 Science exposure that the template was subtracted from. 

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

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

754 difference image. 

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

756 Result of subtracting template from the science image. 

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

758 Detected sources on the difference exposure. 

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

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

761 difference image. 

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

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

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

765 Positive polarity footprints. 

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

767 Negative polarity footprints. 

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

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

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

771 this is also returned. 

772 

773 Returns 

774 ------- 

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

776 

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

778 Subtracted exposure with detection mask applied. 

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

780 The catalog of detected sources. 

781 """ 

782 if measurementResults is None: 

783 measurementResults = pipeBase.Struct() 

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

785 if self.config.doMerge: 

786 # preserve peak schema, if there are any footprints 

787 if len(positives) > 0: 

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

789 elif len(negatives) > 0: 

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

791 else: 

792 peakSchema = afwDetection.PeakTable.makeMinimalSchema() 

793 mergeList = afwDetection.FootprintMergeList(self.schema, 

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

795 initialDiaSources = afwTable.SourceCatalog(self.schema) 

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

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

798 # deblended positive sources. 

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

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

801 mergeList.getFinalSources(initialDiaSources) 

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

803 # footprint set. 

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

805 ~initialDiaSources["merge_footprint_positive"] 

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

807 else: 

808 initialDiaSources = sources 

809 

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

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

812 for source in initialDiaSources: 

813 source.setId(idFactory()) 

814 # Ensure sources added after this get correct ids. 

815 initialDiaSources.getTable().setIdFactory(idFactory) 

816 initialDiaSources.setMetadata(self.algMetadata) 

817 

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

819 

820 if self.config.doMaskStreaks: 

821 streakInfo = self._runStreakMasking(difference) 

822 

823 if self.config.doSkySources: 

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

825 

826 if not initialDiaSources.isContiguous(): 

827 initialDiaSources = initialDiaSources.copy(deep=True) 

828 

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

830 

831 # Remove unphysical diaSources per config.badSourceFlags 

832 diaSources = self._removeBadSources(initialDiaSources) 

833 

834 if self.config.run_sattle: 

835 diaSources = self.filterSatellites(diaSources, science) 

836 

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

838 diaSources, trail_parameters = self._find_glint_trails(diaSources) 

839 if self.config.writeGlintInfo: 

840 measurementResults.mergeItems(trail_parameters, 'glintTrailInfo') 

841 

842 if self.config.doForcedMeasurement: 

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

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

845 template=True) 

846 

847 # Clear the image plane for regions with NO_DATA. 

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

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

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

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