Coverage for python / lsst / ip / diffim / subtractImages.py: 20%

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 

22from astropy import units as u 

23from astropy.stats import gaussian_fwhm_to_sigma 

24import numpy as np 

25 

26import lsst.afw.detection as afwDetection 

27import lsst.afw.image 

28import lsst.afw.math 

29import lsst.geom 

30from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, 

31 checkMask, setSourceFootprints) 

32from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask 

33import lsst.pex.config 

34import lsst.pipe.base 

35import lsst.pex.exceptions 

36from lsst.pipe.base import connectionTypes 

37from . import MakeKernelTask, DecorrelateALKernelTask 

38from lsst.utils.timer import timeMethod 

39 

40__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask", 

41 "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask", 

42 "SimplifiedSubtractConfig", "SimplifiedSubtractTask", 

43 "InsufficientKernelSourcesError"] 

44 

45_dimensions = ("instrument", "visit", "detector") 

46_defaultTemplates = {"coaddName": "deep", "fakesType": ""} 

47 

48 

49class InsufficientKernelSourcesError(lsst.pipe.base.AlgorithmError): 

50 """Raised when there are too few sources to calculate the PSF matching 

51 kernel. 

52 """ 

53 def __init__(self, *, nSources, nRequired): 

54 msg = (f"Only {nSources} sources were selected for PSF matching," 

55 f" but {nRequired} are required.") 

56 super().__init__(msg) 

57 self.nSources = nSources 

58 self.nRequired = nRequired 

59 

60 @property 

61 def metadata(self): 

62 return {"nSources": self.nSources, 

63 "nRequired": self.nRequired 

64 } 

65 

66 

67class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections, 

68 dimensions=_dimensions, 

69 defaultTemplates=_defaultTemplates): 

70 template = connectionTypes.Input( 

71 doc="Input warped template to subtract.", 

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

73 storageClass="ExposureF", 

74 name="{fakesType}{coaddName}Diff_templateExp" 

75 ) 

76 science = connectionTypes.Input( 

77 doc="Input science exposure to subtract from.", 

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

79 storageClass="ExposureF", 

80 name="{fakesType}calexp" 

81 ) 

82 sources = connectionTypes.Input( 

83 doc="Sources measured on the science exposure; " 

84 "used to select sources for making the matching kernel.", 

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

86 storageClass="SourceCatalog", 

87 name="{fakesType}src" 

88 ) 

89 visitSummary = connectionTypes.Input( 

90 doc=("Per-visit catalog with final calibration objects. " 

91 "These catalogs use the detector id for the catalog id, " 

92 "sorted on id for fast lookup."), 

93 dimensions=("instrument", "visit"), 

94 storageClass="ExposureCatalog", 

95 name="finalVisitSummary", 

96 ) 

97 

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

99 super().__init__(config=config) 

100 if not config.doApplyExternalCalibrations: 

101 del self.visitSummary 

102 

103 

104class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections, 

105 dimensions=_dimensions, 

106 defaultTemplates=_defaultTemplates): 

107 difference = connectionTypes.Output( 

108 doc="Result of subtracting convolved template from science image.", 

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

110 storageClass="ExposureF", 

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

112 ) 

113 matchedTemplate = connectionTypes.Output( 

114 doc="Warped and PSF-matched template used to create `subtractedExposure`.", 

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

116 storageClass="ExposureF", 

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

118 ) 

119 psfMatchingKernel = connectionTypes.Output( 

120 doc="Kernel used to PSF match the science and template images.", 

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

122 storageClass="MatchingKernel", 

123 name="{fakesType}{coaddName}Diff_psfMatchKernel", 

124 ) 

125 kernelSources = connectionTypes.Output( 

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

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

128 storageClass="SourceCatalog", 

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

130 ) 

131 

132 

133class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections, 

134 dimensions=_dimensions, 

135 defaultTemplates=_defaultTemplates): 

136 scoreExposure = connectionTypes.Output( 

137 doc="The maximum likelihood image, used for the detection of diaSources.", 

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

139 storageClass="ExposureF", 

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

141 ) 

142 psfMatchingKernel = connectionTypes.Output( 

143 doc="Kernel used to PSF match the science and template images.", 

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

145 storageClass="MatchingKernel", 

146 name="{fakesType}{coaddName}Diff_psfScoreMatchKernel", 

147 ) 

148 kernelSources = connectionTypes.Output( 

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

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

151 storageClass="SourceCatalog", 

152 name="{fakesType}{coaddName}Diff_psfScoreMatchSources" 

153 ) 

154 

155 

156class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

157 pass 

158 

159 

160class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

161 inputPsfMatchingKernel = connectionTypes.Input( 

162 doc="Kernel used to PSF match the science and template images.", 

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

164 storageClass="MatchingKernel", 

165 name="{fakesType}{coaddName}Diff_psfMatchKernel", 

166 ) 

167 

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

169 super().__init__(config=config) 

170 del self.sources 

171 if config.useExistingKernel: 

172 del self.psfMatchingKernel 

173 del self.kernelSources 

174 else: 

175 del self.inputPsfMatchingKernel 

176 

177 

178class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): 

179 makeKernel = lsst.pex.config.ConfigurableField( 

180 target=MakeKernelTask, 

181 doc="Task to construct a matching kernel for convolution.", 

182 ) 

183 doDecorrelation = lsst.pex.config.Field( 

184 dtype=bool, 

185 default=True, 

186 doc="Perform diffim decorrelation to undo pixel correlation due to A&L " 

187 "kernel convolution? If True, also update the diffim PSF." 

188 ) 

189 decorrelate = lsst.pex.config.ConfigurableField( 

190 target=DecorrelateALKernelTask, 

191 doc="Task to decorrelate the image difference.", 

192 ) 

193 requiredTemplateFraction = lsst.pex.config.Field( 

194 dtype=float, 

195 default=0.1, 

196 doc="Raise NoWorkFound and do not attempt image subtraction if template covers less than this " 

197 " fraction of pixels. Setting to 0 will always attempt image subtraction." 

198 ) 

199 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field( 

200 dtype=float, 

201 default=0.2, 

202 doc="Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels." 

203 " If the fraction of pixels covered by the template is less than this value (and greater than" 

204 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated." 

205 ) 

206 doScaleVariance = lsst.pex.config.Field( 

207 dtype=bool, 

208 default=True, 

209 doc="Scale variance of the image difference?" 

210 ) 

211 scaleVariance = lsst.pex.config.ConfigurableField( 

212 target=ScaleVarianceTask, 

213 doc="Subtask to rescale the variance of the template to the statistically expected level." 

214 ) 

215 doSubtractBackground = lsst.pex.config.Field( 

216 doc="Subtract the background fit when solving the kernel? " 

217 "It is generally better to instead subtract the background in detectAndMeasure.", 

218 dtype=bool, 

219 default=False, 

220 ) 

221 doApplyExternalCalibrations = lsst.pex.config.Field( 

222 doc=( 

223 "Replace science Exposure's calibration objects with those" 

224 " in visitSummary. Ignored if `doApplyFinalizedPsf is True." 

225 ), 

226 dtype=bool, 

227 default=False, 

228 ) 

229 sourceSelector = lsst.pex.config.ConfigurableField( 

230 target=ScienceSourceSelectorTask, 

231 doc="Task to select sources to be used for PSF matching.", 

232 ) 

233 fallbackSourceSelector = lsst.pex.config.ConfigurableField( 

234 target=ScienceSourceSelectorTask, 

235 doc="Task to select sources to be used for PSF matching." 

236 "Used only if the kernel calculation fails and" 

237 "`allowKernelSourceDetection` is set. The fallback source detection" 

238 " will not include all of the same plugins as the original source " 

239 " detection, so not all of the same flags can be used.", 

240 ) 

241 detectionThreshold = lsst.pex.config.Field( 

242 dtype=float, 

243 default=10, 

244 doc="Minimum signal to noise ratio of detected sources " 

245 "to use for calculating the PSF matching kernel.", 

246 deprecated="No longer used. Will be removed after v30" 

247 ) 

248 detectionThresholdMax = lsst.pex.config.Field( 

249 dtype=float, 

250 default=500, 

251 doc="Maximum signal to noise ratio of detected sources " 

252 "to use for calculating the PSF matching kernel.", 

253 deprecated="No longer used. Will be removed after v30" 

254 ) 

255 restrictKernelEdgeSources = lsst.pex.config.Field( 

256 dtype=bool, 

257 default=True, 

258 doc="Exclude sources close to the edge from the kernel calculation?" 

259 ) 

260 maxKernelSources = lsst.pex.config.Field( 

261 dtype=int, 

262 default=1000, 

263 doc="Maximum number of sources to use for calculating the PSF matching kernel." 

264 "Set to -1 to disable." 

265 ) 

266 minKernelSources = lsst.pex.config.Field( 

267 dtype=int, 

268 default=3, 

269 doc="Minimum number of sources needed for calculating the PSF matching kernel." 

270 ) 

271 excludeMaskPlanes = lsst.pex.config.ListField( 

272 dtype=str, 

273 default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"), 

274 doc="Template mask planes to exclude when selecting sources for PSF matching.", 

275 ) 

276 badMaskPlanes = lsst.pex.config.ListField( 

277 dtype=str, 

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

279 doc="Mask planes to interpolate over." 

280 ) 

281 preserveTemplateMask = lsst.pex.config.ListField( 

282 dtype=str, 

283 default=("NO_DATA", "BAD", "HIGH_VARIANCE"), 

284 doc="Mask planes from the template to propagate to the image difference." 

285 ) 

286 renameTemplateMask = lsst.pex.config.ListField( 

287 dtype=str, 

288 default=("SAT", "INJECTED", "INJECTED_CORE",), 

289 doc="Mask planes from the template to propagate to the image difference" 

290 "with '_TEMPLATE' appended to the name." 

291 ) 

292 allowKernelSourceDetection = lsst.pex.config.Field( 

293 dtype=bool, 

294 default=False, 

295 doc="Re-run source detection for kernel candidates if an error is" 

296 " encountered while calculating the matching kernel." 

297 ) 

298 

299 def setDefaults(self): 

300 self.makeKernel.kernel.name = "AL" 

301 # Always include background fitting in the kernel fit, 

302 # even if it is not subtracted 

303 self.makeKernel.kernel.active.fitForBackground = True 

304 self.makeKernel.kernel.active.spatialKernelOrder = 1 

305 self.makeKernel.kernel.active.spatialBgOrder = 2 

306 # Shared source selector settings 

307 doSkySources = False # Do not include sky sources 

308 doSignalToNoise = True # apply signal to noise filter 

309 doUnresolved = True # apply star-galaxy separation 

310 signalToNoiseMinimum = 10 

311 signalToNoiseMaximum = 500 

312 self.sourceSelector.doIsolated = True # apply isolated star selection 

313 self.sourceSelector.doRequirePrimary = True # apply primary flag selection 

314 self.sourceSelector.doUnresolved = doUnresolved 

315 self.sourceSelector.doSkySources = doSkySources 

316 self.sourceSelector.doSignalToNoise = doSignalToNoise 

317 self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

318 self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

319 # The following two configs should not be necessary to be turned on for 

320 # PSF-matching, and the fallback kernel source selection will fail if 

321 # they are set since it does not run deblending. 

322 self.fallbackSourceSelector.doIsolated = False # Do not apply isolated star selection 

323 self.fallbackSourceSelector.doRequirePrimary = False # Do not apply primary flag selection 

324 self.fallbackSourceSelector.doUnresolved = doUnresolved 

325 self.fallbackSourceSelector.doSkySources = doSkySources 

326 self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise 

327 self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

328 self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

329 

330 

331class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig, 

332 pipelineConnections=AlardLuptonSubtractConnections): 

333 mode = lsst.pex.config.ChoiceField( 

334 dtype=str, 

335 default="convolveTemplate", 

336 allowed={"auto": "Choose which image to convolve at runtime.", 

337 "convolveScience": "Only convolve the science image.", 

338 "convolveTemplate": "Only convolve the template image."}, 

339 doc="Choose which image to convolve at runtime, or require that a specific image is convolved." 

340 ) 

341 

342 

343class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask): 

344 """Compute the image difference of a science and template image using 

345 the Alard & Lupton (1998) algorithm. 

346 """ 

347 ConfigClass = AlardLuptonSubtractConfig 

348 _DefaultName = "alardLuptonSubtract" 

349 

350 def __init__(self, **kwargs): 

351 super().__init__(**kwargs) 

352 self.makeSubtask("decorrelate") 

353 self.makeSubtask("makeKernel") 

354 self.makeSubtask("sourceSelector") 

355 self.makeSubtask("fallbackSourceSelector") 

356 if self.config.doScaleVariance: 

357 self.makeSubtask("scaleVariance") 

358 

359 self.convolutionControl = lsst.afw.math.ConvolutionControl() 

360 # Normalization is an extra, unnecessary, calculation and will result 

361 # in mis-subtraction of the images if there are calibration errors. 

362 self.convolutionControl.setDoNormalize(False) 

363 self.convolutionControl.setDoCopyEdge(True) 

364 

365 def _applyExternalCalibrations(self, exposure, visitSummary): 

366 """Replace calibrations (psf, and ApCorrMap) on this exposure with 

367 external ones.". 

368 

369 Parameters 

370 ---------- 

371 exposure : `lsst.afw.image.exposure.Exposure` 

372 Input exposure to adjust calibrations. 

373 visitSummary : `lsst.afw.table.ExposureCatalog` 

374 Exposure catalog with external calibrations to be applied. Catalog 

375 uses the detector id for the catalog id, sorted on id for fast 

376 lookup. 

377 

378 Returns 

379 ------- 

380 exposure : `lsst.afw.image.exposure.Exposure` 

381 Exposure with adjusted calibrations. 

382 """ 

383 detectorId = exposure.info.getDetector().getId() 

384 

385 row = visitSummary.find(detectorId) 

386 if row is None: 

387 self.log.warning("Detector id %s not found in external calibrations catalog; " 

388 "Using original calibrations.", detectorId) 

389 else: 

390 psf = row.getPsf() 

391 apCorrMap = row.getApCorrMap() 

392 if psf is None: 

393 self.log.warning("Detector id %s has None for psf in " 

394 "external calibrations catalog; Using original psf and aperture correction.", 

395 detectorId) 

396 elif apCorrMap is None: 

397 self.log.warning("Detector id %s has None for apCorrMap in " 

398 "external calibrations catalog; Using original psf and aperture correction.", 

399 detectorId) 

400 else: 

401 exposure.setPsf(psf) 

402 exposure.info.setApCorrMap(apCorrMap) 

403 

404 return exposure 

405 

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

407 inputs = butlerQC.get(inputRefs) 

408 

409 try: 

410 results = self.run(**inputs) 

411 except lsst.pipe.base.AlgorithmError as e: 

412 error = lsst.pipe.base.AnnotatedPartialOutputsError.annotate(e, self, log=self.log) 

413 # No partial outputs for butler to put 

414 raise error from e 

415 

416 butlerQC.put(results, outputRefs) 

417 

418 @timeMethod 

419 def run(self, template, science, sources, visitSummary=None): 

420 """PSF match, subtract, and decorrelate two images. 

421 

422 Parameters 

423 ---------- 

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

425 Template exposure, warped to match the science exposure. 

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

427 Science exposure to subtract from the template. 

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

429 Identified sources on the science exposure. This catalog is used to 

430 select sources in order to perform the AL PSF matching on stamp 

431 images around them. 

432 visitSummary : `lsst.afw.table.ExposureCatalog`, optional 

433 Exposure catalog with external calibrations to be applied. Catalog 

434 uses the detector id for the catalog id, sorted on id for fast 

435 lookup. 

436 

437 Returns 

438 ------- 

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

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

441 Result of subtracting template and science. 

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

443 Warped and PSF-matched template exposure. 

444 ``backgroundModel`` : `lsst.afw.math.Function2D` 

445 Background model that was fit while solving for the 

446 PSF-matching kernel 

447 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

448 Kernel used to PSF-match the convolved image. 

449 ``kernelSources` : `lsst.afw.table.SourceCatalog` 

450 Sources from the input catalog that were used to construct the 

451 PSF-matching kernel. 

452 """ 

453 self._prepareInputs(template, science, visitSummary=visitSummary) 

454 

455 convolveTemplate = self.chooseConvolutionMethod(template, science) 

456 

457 kernelResult = self.runMakeKernel(template, science, sources=sources, 

458 convolveTemplate=convolveTemplate, 

459 runSourceDetection=False) 

460 

461 if self.config.doSubtractBackground: 

462 backgroundModel = kernelResult.backgroundModel 

463 else: 

464 backgroundModel = None 

465 if convolveTemplate: 

466 subtractResults = self.runConvolveTemplate(template, science, kernelResult.psfMatchingKernel, 

467 backgroundModel=backgroundModel) 

468 else: 

469 subtractResults = self.runConvolveScience(template, science, kernelResult.psfMatchingKernel, 

470 backgroundModel=backgroundModel) 

471 subtractResults.kernelSources = kernelResult.kernelSources 

472 

473 metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources) 

474 

475 self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean 

476 self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev 

477 self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean 

478 self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev 

479 self.log.info("Mean, stdev of ratio of difference to science " 

480 "pixels in star footprints: %5.4f, %5.4f", 

481 self.metadata["differenceFootprintRatioMean"], 

482 self.metadata["differenceFootprintRatioStdev"]) 

483 

484 return subtractResults 

485 

486 def chooseConvolutionMethod(self, template, science): 

487 """Determine whether the template should be convolved with the PSF 

488 matching kernel. 

489 

490 Parameters 

491 ---------- 

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

493 Template exposure, warped to match the science exposure. 

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

495 Science exposure to subtract from the template. 

496 

497 Returns 

498 ------- 

499 convolveTemplate : `bool` 

500 Convolve the template to match the two images? 

501 

502 Raises 

503 ------ 

504 RuntimeError 

505 If an unsupported convolution mode is supplied. 

506 """ 

507 if self.config.mode == "auto": 

508 convolveTemplate = _shapeTest(template, 

509 science, 

510 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer, 

511 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid) 

512 if convolveTemplate: 

513 if self.sciencePsfSize < self.templatePsfSize: 

514 self.log.info("Average template PSF size is greater, " 

515 "but science PSF greater in one dimension: convolving template image.") 

516 else: 

517 self.log.info("Science PSF size is greater: convolving template image.") 

518 else: 

519 self.log.info("Template PSF size is greater: convolving science image.") 

520 elif self.config.mode == "convolveTemplate": 

521 self.log.info("`convolveTemplate` is set: convolving template image.") 

522 convolveTemplate = True 

523 elif self.config.mode == "convolveScience": 

524 self.log.info("`convolveScience` is set: convolving science image.") 

525 convolveTemplate = False 

526 else: 

527 raise RuntimeError("Cannot handle AlardLuptonSubtract mode: %s", self.config.mode) 

528 return convolveTemplate 

529 

530 def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False): 

531 """Construct the PSF-matching kernel. 

532 

533 Parameters 

534 ---------- 

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

536 Template exposure, warped to match the science exposure. 

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

538 Science exposure to subtract from the template. 

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

540 Identified sources on the science exposure. This catalog is used to 

541 select sources in order to perform the AL PSF matching on stamp 

542 images around them. 

543 Not used if ``runSourceDetection`` is set. 

544 convolveTemplate : `bool`, optional 

545 Construct the matching kernel to convolve the template? 

546 runSourceDetection : `bool`, optional 

547 Run a minimal version of source detection to determine kernel 

548 candidates? If False, a source list to select kernel candidates 

549 from must be supplied. 

550 

551 Returns 

552 ------- 

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

554 ``backgroundModel`` : `lsst.afw.math.Function2D` 

555 Background model that was fit while solving for the 

556 PSF-matching kernel 

557 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

558 Kernel used to PSF-match the convolved image. 

559 ``kernelSources` : `lsst.afw.table.SourceCatalog` 

560 Sources from the input catalog that were used to construct the 

561 PSF-matching kernel. 

562 """ 

563 

564 if convolveTemplate: 

565 reference = template 

566 target = science 

567 referenceFwhmPix = self.templatePsfSize 

568 targetFwhmPix = self.sciencePsfSize 

569 else: 

570 reference = science 

571 target = template 

572 referenceFwhmPix = self.sciencePsfSize 

573 targetFwhmPix = self.templatePsfSize 

574 try: 

575 if runSourceDetection: 

576 kernelSources = self.runKernelSourceDetection(template, science) 

577 else: 

578 kernelSources = self._sourceSelector(template, science, sources) 

579 kernelResult = self.makeKernel.run(reference, target, kernelSources, 

580 preconvolved=False, 

581 templateFwhmPix=referenceFwhmPix, 

582 scienceFwhmPix=targetFwhmPix) 

583 except (RuntimeError, lsst.pex.exceptions.Exception) as e: 

584 self.log.warning("Failed to match template. Checking coverage") 

585 # Raise NoWorkFound if template fraction is insufficient 

586 checkTemplateIsSufficient(template[science.getBBox()], science, self.log, 

587 self.config.minTemplateFractionForExpectedSuccess, 

588 exceptionMessage="Template coverage lower than expected to succeed." 

589 f" Failure is tolerable: {e}") 

590 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception 

591 raise e 

592 

593 return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel, 

594 psfMatchingKernel=kernelResult.psfMatchingKernel, 

595 kernelSources=kernelSources) 

596 

597 def runKernelSourceDetection(self, template, science): 

598 """Run detection on the science image and use the template mask plane 

599 to reject candidate sources. 

600 

601 Parameters 

602 ---------- 

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

604 Template exposure, warped to match the science exposure. 

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

606 Science exposure to subtract from the template. 

607 

608 Returns 

609 ------- 

610 kernelSources : `lsst.afw.table.SourceCatalog` 

611 Sources from the input catalog to use to construct the 

612 PSF-matching kernel. 

613 """ 

614 kernelSize = self.makeKernel.makeKernelBasisList( 

615 self.templatePsfSize, self.sciencePsfSize)[0].getWidth() 

616 candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize, 

617 candidateList=None, 

618 sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize) 

619 sources = self.makeKernel.selectKernelSources(template, science, 

620 candidateList=candidateList, 

621 preconvolved=False, 

622 templateFwhmPix=self.templatePsfSize, 

623 scienceFwhmPix=self.sciencePsfSize) 

624 

625 # return sources 

626 return self._sourceSelector(template, science, sources, fallback=True) 

627 

628 def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None): 

629 """Convolve the template image with a PSF-matching kernel and subtract 

630 from the science image. 

631 

632 Parameters 

633 ---------- 

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

635 Template exposure, warped to match the science exposure. 

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

637 Science exposure to subtract from the template. 

638 psfMatchingKernel : `lsst.afw.math.Kernel` 

639 Kernel to be used to PSF-match the science image to the template. 

640 backgroundModel : `lsst.afw.math.Function2D`, optional 

641 Background model that was fit while solving for the PSF-matching 

642 kernel. 

643 

644 Returns 

645 ------- 

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

647 

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

649 Result of subtracting template and science. 

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

651 Warped and PSF-matched template exposure. 

652 ``backgroundModel`` : `lsst.afw.math.Function2D` 

653 Background model that was fit while solving for the PSF-matching kernel 

654 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

655 Kernel used to PSF-match the template to the science image. 

656 """ 

657 self.metadata["convolvedExposure"] = "Template" 

658 

659 matchedTemplate = self._convolveExposure(template, psfMatchingKernel, 

660 self.convolutionControl, 

661 bbox=science.getBBox(), 

662 psf=science.psf, 

663 photoCalib=science.photoCalib) 

664 

665 difference = _subtractImages(science, matchedTemplate, backgroundModel=backgroundModel) 

666 correctedExposure = self.finalize(template, science, difference, 

667 psfMatchingKernel, 

668 templateMatched=True) 

669 

670 return lsst.pipe.base.Struct(difference=correctedExposure, 

671 matchedTemplate=matchedTemplate, 

672 matchedScience=science, 

673 backgroundModel=backgroundModel, 

674 psfMatchingKernel=psfMatchingKernel) 

675 

676 def runConvolveScience(self, template, science, psfMatchingKernel, backgroundModel=None): 

677 """Convolve the science image with a PSF-matching kernel and subtract 

678 the template image. 

679 

680 Parameters 

681 ---------- 

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

683 Template exposure, warped to match the science exposure. 

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

685 Science exposure to subtract from the template. 

686 psfMatchingKernel : `lsst.afw.math.Kernel` 

687 Kernel to be used to PSF-match the science image to the template. 

688 backgroundModel : `lsst.afw.math.Function2D`, optional 

689 Background model that was fit while solving for the PSF-matching 

690 kernel. 

691 

692 Returns 

693 ------- 

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

695 

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

697 Result of subtracting template and science. 

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

699 Warped template exposure. Note that in this case, the template 

700 is not PSF-matched to the science image. 

701 ``backgroundModel`` : `lsst.afw.math.Function2D` 

702 Background model that was fit while solving for the PSF-matching kernel 

703 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

704 Kernel used to PSF-match the science image to the template. 

705 """ 

706 self.metadata["convolvedExposure"] = "Science" 

707 bbox = science.getBBox() 

708 

709 kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions()) 

710 norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False) 

711 

712 matchedScience = self._convolveExposure(science, psfMatchingKernel, 

713 self.convolutionControl, 

714 psf=template.psf) 

715 

716 # Place back on native photometric scale 

717 matchedScience.maskedImage /= norm 

718 matchedTemplate = template.clone()[bbox] 

719 matchedTemplate.maskedImage /= norm 

720 matchedTemplate.setPhotoCalib(science.photoCalib) 

721 

722 if backgroundModel is not None: 

723 modelParams = backgroundModel.getParameters() 

724 # We must invert the background model if the matching kernel is solved for the science image. 

725 backgroundModel.setParameters([-p for p in modelParams]) 

726 

727 difference = _subtractImages(matchedScience, matchedTemplate, backgroundModel=backgroundModel) 

728 

729 correctedExposure = self.finalize(template, science, difference, 

730 psfMatchingKernel, 

731 templateMatched=False) 

732 

733 return lsst.pipe.base.Struct(difference=correctedExposure, 

734 matchedTemplate=matchedTemplate, 

735 matchedScience=matchedScience, 

736 backgroundModel=backgroundModel, 

737 psfMatchingKernel=psfMatchingKernel) 

738 

739 def finalize(self, template, science, difference, kernel, 

740 templateMatched=True, 

741 preConvMode=False, 

742 preConvKernel=None, 

743 spatiallyVarying=False): 

744 """Decorrelate the difference image to undo the noise correlations 

745 caused by convolution. 

746 

747 Parameters 

748 ---------- 

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

750 Template exposure, warped to match the science exposure. 

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

752 Science exposure to subtract from the template. 

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

754 Result of subtracting template and science. 

755 kernel : `lsst.afw.math.Kernel` 

756 An (optionally spatially-varying) PSF matching kernel 

757 templateMatched : `bool`, optional 

758 Was the template PSF-matched to the science image? 

759 preConvMode : `bool`, optional 

760 Was the science image preconvolved with its own PSF 

761 before PSF matching the template? 

762 preConvKernel : `lsst.afw.detection.Psf`, optional 

763 If not `None`, then the science image was pre-convolved with 

764 (the reflection of) this kernel. Must be normalized to sum to 1. 

765 spatiallyVarying : `bool`, optional 

766 Compute the decorrelation kernel spatially varying across the image? 

767 

768 Returns 

769 ------- 

770 correctedExposure : `lsst.afw.image.ExposureF` 

771 The decorrelated image difference. 

772 """ 

773 if self.config.doDecorrelation: 

774 self.log.info("Decorrelating image difference.") 

775 # We have cleared the template mask plane, so copy the mask plane of 

776 # the image difference so that we can calculate correct statistics 

777 # during decorrelation 

778 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel, 

779 templateMatched=templateMatched, 

780 preConvMode=preConvMode, 

781 preConvKernel=preConvKernel, 

782 spatiallyVarying=spatiallyVarying).correctedExposure 

783 else: 

784 self.log.info("NOT decorrelating image difference.") 

785 correctedExposure = difference 

786 return correctedExposure 

787 

788 def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None): 

789 """Calculate an exposure's limiting magnitude. 

790 

791 This method uses the photometric zeropoint together with the 

792 PSF size from the average position of the exposure. 

793 

794 Parameters 

795 ---------- 

796 exposure : `lsst.afw.image.Exposure` 

797 The target exposure to calculate the limiting magnitude for. 

798 nsigma : `float`, optional 

799 The detection threshold in sigma. 

800 fallbackPsfSize : `float`, optional 

801 PSF FWHM to use in the event the exposure PSF cannot be retrieved. 

802 

803 Returns 

804 ------- 

805 maglim : `astropy.units.Quantity` 

806 The limiting magnitude of the exposure, or np.nan. 

807 """ 

808 if exposure.photoCalib is None: 

809 return np.nan 

810 # Set maglim to nan upfront in case on an unexpected RuntimeError 

811 maglim = np.nan 

812 try: 

813 psf = exposure.getPsf() 

814 psf_shape = psf.computeShape(psf.getAveragePosition()) 

815 except (lsst.pex.exceptions.InvalidParameterError, 

816 afwDetection.InvalidPsfError, 

817 lsst.pex.exceptions.RangeError): 

818 if fallbackPsfSize is not None: 

819 self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize) 

820 psf_area = np.pi*(fallbackPsfSize/2)**2 

821 zeropoint = exposure.photoCalib.instFluxToMagnitude(1) 

822 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area)) 

823 else: 

824 self.log.info("Unable to evaluate PSF, setting maglim to nan") 

825 maglim = np.nan 

826 else: 

827 # Get a more accurate area than `psf_shape.getArea()` via moments 

828 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy()) 

829 zeropoint = exposure.photoCalib.instFluxToMagnitude(1) 

830 maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area)) 

831 finally: 

832 return maglim 

833 

834 @staticmethod 

835 def _validateExposures(template, science): 

836 """Check that the WCS of the two Exposures match, the template bbox 

837 contains the science bbox, and that the bands match. 

838 

839 Parameters 

840 ---------- 

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

842 Template exposure, warped to match the science exposure. 

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

844 Science exposure to subtract from the template. 

845 

846 Raises 

847 ------ 

848 AssertionError 

849 Raised if the WCS of the template is not equal to the science WCS, 

850 if the science image is not fully contained in the template 

851 bounding box, or if the bands do not match. 

852 """ 

853 assert template.wcs == science.wcs, \ 

854 "Template and science exposure WCS are not identical." 

855 templateBBox = template.getBBox() 

856 scienceBBox = science.getBBox() 

857 assert science.filter.bandLabel == template.filter.bandLabel, \ 

858 "Science and template exposures have different bands: %s, %s" % \ 

859 (science.filter, template.filter) 

860 

861 assert templateBBox.contains(scienceBBox), \ 

862 "Template bbox does not contain all of the science image." 

863 

864 def _convolveExposure(self, exposure, kernel, convolutionControl, 

865 bbox=None, 

866 psf=None, 

867 photoCalib=None, 

868 interpolateBadMaskPlanes=False, 

869 ): 

870 """Convolve an exposure with the given kernel. 

871 

872 Parameters 

873 ---------- 

874 exposure : `lsst.afw.Exposure` 

875 exposure to convolve. 

876 kernel : `lsst.afw.math.LinearCombinationKernel` 

877 PSF matching kernel computed in the ``makeKernel`` subtask. 

878 convolutionControl : `lsst.afw.math.ConvolutionControl` 

879 Configuration for convolve algorithm. 

880 bbox : `lsst.geom.Box2I`, optional 

881 Bounding box to trim the convolved exposure to. 

882 psf : `lsst.afw.detection.Psf`, optional 

883 Point spread function (PSF) to set for the convolved exposure. 

884 photoCalib : `lsst.afw.image.PhotoCalib`, optional 

885 Photometric calibration of the convolved exposure. 

886 interpolateBadMaskPlanes : `bool`, optional 

887 If set, interpolate over mask planes specified in 

888 ``config.badMaskPlanes`` before convolving the image. 

889 

890 Returns 

891 ------- 

892 convolvedExp : `lsst.afw.Exposure` 

893 The convolved image. 

894 """ 

895 convolvedExposure = exposure.clone() 

896 if psf is not None: 

897 convolvedExposure.setPsf(psf) 

898 if photoCalib is not None: 

899 convolvedExposure.setPhotoCalib(photoCalib) 

900 if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None: 

901 nInterp = _interpolateImage(convolvedExposure.maskedImage, 

902 self.config.badMaskPlanes) 

903 self.metadata["nInterpolated"] = nInterp