Coverage for python/lsst/ip/diffim/subtractImages.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 

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, 

31 computeDifferenceImageMetrics, computePSFNoiseEquivalentArea, 

32 checkMask, setSourceFootprints) 

33from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask 

34import lsst.pex.config 

35import lsst.pipe.base 

36import lsst.pex.exceptions 

37from lsst.pipe.base import connectionTypes 

38from . import MakeKernelTask, DecorrelateALKernelTask 

39from lsst.utils.timer import timeMethod 

40 

41__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask", 

42 "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask", 

43 "SimplifiedSubtractConfig", "SimplifiedSubtractTask", 

44 "InsufficientKernelSourcesError"] 

45 

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

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

48 

49 

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

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

52 kernel. 

53 """ 

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

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

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

57 super().__init__(msg) 

58 self.nSources = nSources 

59 self.nRequired = nRequired 

60 

61 @property 

62 def metadata(self): 

63 return {"nSources": self.nSources, 

64 "nRequired": self.nRequired 

65 } 

66 

67 

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

69 dimensions=_dimensions, 

70 defaultTemplates=_defaultTemplates): 

71 template = connectionTypes.Input( 

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

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

74 storageClass="ExposureF", 

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

76 ) 

77 science = connectionTypes.Input( 

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

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

80 storageClass="ExposureF", 

81 name="{fakesType}calexp" 

82 ) 

83 sources = connectionTypes.Input( 

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

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

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

87 storageClass="SourceCatalog", 

88 name="{fakesType}src" 

89 ) 

90 visitSummary = connectionTypes.Input( 

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

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

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

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

95 storageClass="ExposureCatalog", 

96 name="finalVisitSummary", 

97 ) 

98 

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

100 super().__init__(config=config) 

101 if not config.doApplyExternalCalibrations: 

102 del self.visitSummary 

103 

104 

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

106 dimensions=_dimensions, 

107 defaultTemplates=_defaultTemplates): 

108 difference = connectionTypes.Output( 

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

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

111 storageClass="ExposureF", 

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

113 ) 

114 matchedTemplate = connectionTypes.Output( 

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

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

117 storageClass="ExposureF", 

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

119 ) 

120 psfMatchingKernel = connectionTypes.Output( 

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

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

123 storageClass="MatchingKernel", 

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

125 ) 

126 kernelSources = connectionTypes.Output( 

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

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

129 storageClass="SourceCatalog", 

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

131 ) 

132 

133 

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

135 dimensions=_dimensions, 

136 defaultTemplates=_defaultTemplates): 

137 scoreExposure = connectionTypes.Output( 

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

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

140 storageClass="ExposureF", 

141 name="{fakesType}{coaddName}Diff_scoreTempExp", 

142 ) 

143 psfMatchingKernel = connectionTypes.Output( 

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

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

146 storageClass="MatchingKernel", 

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

148 ) 

149 kernelSources = connectionTypes.Output( 

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

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

152 storageClass="SourceCatalog", 

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

154 ) 

155 

156 

157class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

158 pass 

159 

160 

161class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

162 inputPsfMatchingKernel = connectionTypes.Input( 

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

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

165 storageClass="MatchingKernel", 

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

167 ) 

168 

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

170 super().__init__(config=config) 

171 del self.sources 

172 if config.useExistingKernel: 

173 del self.psfMatchingKernel 

174 del self.kernelSources 

175 else: 

176 del self.inputPsfMatchingKernel 

177 

178 

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

180 makeKernel = lsst.pex.config.ConfigurableField( 

181 target=MakeKernelTask, 

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

183 ) 

184 doDecorrelation = lsst.pex.config.Field( 

185 dtype=bool, 

186 default=True, 

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

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

189 ) 

190 decorrelate = lsst.pex.config.ConfigurableField( 

191 target=DecorrelateALKernelTask, 

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

193 ) 

194 requiredTemplateFraction = lsst.pex.config.Field( 

195 dtype=float, 

196 default=0.1, 

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

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

199 ) 

200 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field( 

201 dtype=float, 

202 default=0.2, 

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

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

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

206 ) 

207 doScaleVariance = lsst.pex.config.Field( 

208 dtype=bool, 

209 default=True, 

210 doc="Scale variance of the science image? Note that the template variance is NOT scaled" 

211 " here. The template variance may be scaled independently in ``GetTemplateTask``." 

212 ) 

213 scaleVariance = lsst.pex.config.ConfigurableField( 

214 target=ScaleVarianceTask, 

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

216 ) 

217 doSubtractBackground = lsst.pex.config.Field( 

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

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

220 dtype=bool, 

221 default=False, 

222 ) 

223 doApplyExternalCalibrations = lsst.pex.config.Field( 

224 doc=( 

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

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

227 ), 

228 dtype=bool, 

229 default=False, 

230 ) 

231 sourceSelector = lsst.pex.config.ConfigurableField( 

232 target=ScienceSourceSelectorTask, 

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

234 ) 

235 fallbackSourceSelector = lsst.pex.config.ConfigurableField( 

236 target=ScienceSourceSelectorTask, 

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

238 "Used only if the kernel calculation fails and" 

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

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

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

242 ) 

243 detectionThreshold = lsst.pex.config.Field( 

244 dtype=float, 

245 default=10, 

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

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

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

249 ) 

250 detectionThresholdMax = lsst.pex.config.Field( 

251 dtype=float, 

252 default=500, 

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

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

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

256 ) 

257 restrictKernelEdgeSources = lsst.pex.config.Field( 

258 dtype=bool, 

259 default=True, 

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

261 ) 

262 maxKernelSources = lsst.pex.config.Field( 

263 dtype=int, 

264 default=1000, 

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

266 "Set to -1 to disable." 

267 ) 

268 minKernelSources = lsst.pex.config.Field( 

269 dtype=int, 

270 default=3, 

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

272 ) 

273 excludeMaskPlanes = lsst.pex.config.ListField( 

274 dtype=str, 

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

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

277 ) 

278 badMaskPlanes = lsst.pex.config.ListField( 

279 dtype=str, 

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

281 doc="Mask planes to interpolate over." 

282 ) 

283 preserveTemplateMask = lsst.pex.config.ListField( 

284 dtype=str, 

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

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

287 ) 

288 renameTemplateMask = lsst.pex.config.ListField( 

289 dtype=str, 

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

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

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

293 ) 

294 preserveMaskPlanes = lsst.pex.config.ListField( 

295 dtype=str, 

296 default=("INJECTED", "INJECTED_CORE", "INJECTED_TEMPLATE", "INJECTED_CORE_TEMPLATE"), 

297 doc="Mask planes to preserve without dilation when convolving the image.", 

298 ) 

299 allowKernelSourceDetection = lsst.pex.config.Field( 

300 dtype=bool, 

301 default=False, 

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

303 " encountered while calculating the matching kernel." 

304 ) 

305 

306 def setDefaults(self): 

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

308 # Always include background fitting in the kernel fit, 

309 # even if it is not subtracted 

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

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

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

313 # Shared source selector settings 

314 doSkySources = False # Do not include sky sources 

315 doSignalToNoise = True # apply signal to noise filter 

316 doUnresolved = True # apply star-galaxy separation 

317 signalToNoiseMinimum = 10 

318 signalToNoiseMaximum = 500 

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

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

321 self.sourceSelector.doUnresolved = doUnresolved 

322 self.sourceSelector.doSkySources = doSkySources 

323 self.sourceSelector.doSignalToNoise = doSignalToNoise 

324 self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

325 self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

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

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

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

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

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

331 self.fallbackSourceSelector.doUnresolved = doUnresolved 

332 self.fallbackSourceSelector.doSkySources = doSkySources 

333 self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise 

334 self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

335 self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

336 

337 

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

339 pipelineConnections=AlardLuptonSubtractConnections): 

340 mode = lsst.pex.config.ChoiceField( 

341 dtype=str, 

342 default="convolveTemplate", 

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

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

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

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

347 ) 

348 

349 

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

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

352 the Alard & Lupton (1998) algorithm. 

353 """ 

354 ConfigClass = AlardLuptonSubtractConfig 

355 _DefaultName = "alardLuptonSubtract" 

356 usePreconvolution = False 

357 """Whether this task preconvolves the science image with its own PSF 

358 before kernel-matching. Subclasses that preconvolve override this to 

359 `True`.""" 

360 

361 def __init__(self, **kwargs): 

362 super().__init__(**kwargs) 

363 self.makeSubtask("decorrelate") 

364 self.makeSubtask("makeKernel") 

365 self.makeSubtask("sourceSelector") 

366 self.makeSubtask("fallbackSourceSelector") 

367 if self.config.doScaleVariance: 

368 self.makeSubtask("scaleVariance") 

369 

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

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

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

373 self.convolutionControl.setDoNormalize(False) 

374 self.convolutionControl.setDoCopyEdge(True) 

375 

376 def _applyExternalCalibrations(self, exposure, visitSummary): 

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

378 external ones.". 

379 

380 Parameters 

381 ---------- 

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

383 Input exposure to adjust calibrations. 

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

385 Exposure catalog with external calibrations to be applied. Catalog 

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

387 lookup. 

388 

389 Returns 

390 ------- 

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

392 Exposure with adjusted calibrations. 

393 """ 

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

395 

396 row = visitSummary.find(detectorId) 

397 if row is None: 

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

399 "Using original calibrations.", detectorId) 

400 else: 

401 psf = row.getPsf() 

402 apCorrMap = row.getApCorrMap() 

403 if psf is None: 

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

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

406 detectorId) 

407 elif apCorrMap is None: 

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

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

410 detectorId) 

411 else: 

412 exposure.setPsf(psf) 

413 exposure.info.setApCorrMap(apCorrMap) 

414 

415 return exposure 

416 

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

418 inputs = butlerQC.get(inputRefs) 

419 

420 try: 

421 results = self.run(**inputs) 

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

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

424 # No partial outputs for butler to put 

425 raise error from e 

426 

427 butlerQC.put(results, outputRefs) 

428 

429 @timeMethod 

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

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

432 

433 Parameters 

434 ---------- 

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

436 Template exposure, warped to match the science exposure. 

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

438 Science exposure to subtract from the template. 

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

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

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

442 images around them. 

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

444 Exposure catalog with external calibrations to be applied. Catalog 

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

446 lookup. 

447 

448 Returns 

449 ------- 

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

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

452 Result of subtracting template and science. 

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

454 Warped and PSF-matched template exposure. 

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

456 Background model that was fit while solving for the 

457 PSF-matching kernel 

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

459 Kernel used to PSF-match the convolved image. 

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

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

462 PSF-matching kernel. 

463 """ 

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

465 

466 convolveTemplate = self.chooseConvolutionMethod(template, science) 

467 self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize 

468 

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

470 convolveTemplate=convolveTemplate, 

471 runSourceDetection=False) 

472 

473 if self.config.doSubtractBackground: 

474 backgroundModel = kernelResult.backgroundModel 

475 else: 

476 backgroundModel = None 

477 if convolveTemplate: 

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

479 backgroundModel=backgroundModel) 

480 else: 

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

482 backgroundModel=backgroundModel) 

483 subtractResults.kernelSources = kernelResult.kernelSources 

484 

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

486 

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

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

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

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

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

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

493 self.metadata["differenceFootprintRatioMean"], 

494 self.metadata["differenceFootprintRatioStdev"]) 

495 

496 return subtractResults 

497 

498 def chooseConvolutionMethod(self, template, science): 

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

500 matching kernel. 

501 

502 Parameters 

503 ---------- 

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

505 Template exposure, warped to match the science exposure. 

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

507 Science exposure to subtract from the template. 

508 

509 Returns 

510 ------- 

511 convolveTemplate : `bool` 

512 Convolve the template to match the two images? 

513 

514 Raises 

515 ------ 

516 RuntimeError 

517 If an unsupported convolution mode is supplied. 

518 """ 

519 if self.usePreconvolution: 519 ↛ 520line 519 didn't jump to line 520 because the condition on line 519 was never true

520 raise RuntimeError("Choosing a convolution method is incompatible with preconvolution!") 

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

522 convolveTemplate = _shapeTest(template, 

523 science, 

524 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer, 

525 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid) 

526 if convolveTemplate: 

527 if self.sciencePsfSize < self.templatePsfSize: 527 ↛ 528line 527 didn't jump to line 528 because the condition on line 527 was never true

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

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

530 else: 

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

532 else: 

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

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

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

536 convolveTemplate = True 

537 elif self.config.mode == "convolveScience": 537 ↛ 541line 537 didn't jump to line 541 because the condition on line 537 was always true

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

539 convolveTemplate = False 

540 else: 

541 raise RuntimeError(f"Cannot handle AlardLuptonSubtract mode: {self.config.mode}") 

542 return convolveTemplate 

543 

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

545 """Construct the PSF-matching kernel. Not used for preconvolution. 

546 

547 Parameters 

548 ---------- 

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

550 Template exposure, warped to match the science exposure. 

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

552 Science exposure to subtract from the template. 

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

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

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

556 images around them. 

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

558 convolveTemplate : `bool`, optional 

559 Construct the matching kernel to convolve the template? 

560 runSourceDetection : `bool`, optional 

561 Run a minimal version of source detection to determine kernel 

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

563 from must be supplied. 

564 

565 Returns 

566 ------- 

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

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

569 Background model that was fit while solving for the 

570 PSF-matching kernel 

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

572 Kernel used to PSF-match the convolved image. 

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

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

575 PSF-matching kernel. 

576 """ 

577 if self.usePreconvolution: 577 ↛ 578line 577 didn't jump to line 578 because the condition on line 577 was never true

578 raise RuntimeError("Incorrect matching kernel calculation configured. " 

579 "`runMakeKernel` can't be called if `usePreconvolution` is set.") 

580 if convolveTemplate: 

581 reference = template 

582 target = science 

583 referenceFwhmPix = self.templatePsfSize 

584 targetFwhmPix = self.sciencePsfSize 

585 else: 

586 reference = science 

587 target = template 

588 referenceFwhmPix = self.sciencePsfSize 

589 targetFwhmPix = self.templatePsfSize 

590 try: 

591 if runSourceDetection: 

592 kernelSources = self.runKernelSourceDetection(template, science) 

593 else: 

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

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

596 preconvolved=False, 

597 templateFwhmPix=referenceFwhmPix, 

598 scienceFwhmPix=targetFwhmPix) 

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

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

601 # Raise NoWorkFound if template fraction is insufficient 

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

603 self.config.minTemplateFractionForExpectedSuccess, 

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

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

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

607 raise e 

608 

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

610 psfMatchingKernel=kernelResult.psfMatchingKernel, 

611 kernelSources=kernelSources) 

612 

613 def runKernelSourceDetection(self, template, science): 

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

615 to reject candidate sources. 

616 

617 Parameters 

618 ---------- 

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

620 Template exposure, warped to match the science exposure. 

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

622 Science exposure to subtract from the template. 

623 

624 Returns 

625 ------- 

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

627 Sources from the input catalog to use to construct the 

628 PSF-matching kernel. 

629 """ 

630 kernelSize = self.makeKernel.makeKernelBasisList( 

631 self.templatePsfSize, self.matchedPsfSize)[0].getWidth() 

632 sources = self.makeKernel.makeCandidateList(template, science, kernelSize, 

633 candidateList=None, 

634 sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize) 

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

636 

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

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

639 from the science image. 

640 

641 Parameters 

642 ---------- 

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

644 Template exposure, warped to match the science exposure. 

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

646 Science exposure to subtract from the template. 

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

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

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

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

651 kernel. 

652 

653 Returns 

654 ------- 

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

656 

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

658 Result of subtracting template and science. 

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

660 Warped and PSF-matched template exposure. 

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

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

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

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

665 """ 

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

667 

668 matchedTemplate = self._convolveExposure(template, psfMatchingKernel, 

669 self.convolutionControl, 

670 bbox=science.getBBox(), 

671 psf=science.psf, 

672 photoCalib=science.photoCalib) 

673 

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

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

676 psfMatchingKernel, 

677 templateMatched=True) 

678 

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

680 matchedTemplate=matchedTemplate, 

681 matchedScience=science, 

682 backgroundModel=backgroundModel, 

683 psfMatchingKernel=psfMatchingKernel) 

684 

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

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

687 the template image. 

688 

689 Parameters 

690 ---------- 

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

692 Template exposure, warped to match the science exposure. 

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

694 Science exposure to subtract from the template. 

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

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

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

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

699 kernel. 

700 

701 Returns 

702 ------- 

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

704 

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

706 Result of subtracting template and science. 

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

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

709 is not PSF-matched to the science image. 

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

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

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

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

714 """ 

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

716 bbox = science.getBBox() 

717 

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

719 xcen, ycen = bbox.getCenter() 

720 norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False, x=xcen, y=ycen) 

721 

722 matchedScience = self._convolveExposure(science, psfMatchingKernel, 

723 self.convolutionControl, 

724 psf=template.psf) 

725 

726 # Place back on native photometric scale 

727 matchedScience.maskedImage /= norm 

728 matchedTemplate = template.clone()[bbox] 

729 matchedTemplate.setPhotoCalib(science.photoCalib) 

730 

731 if backgroundModel is not None: 

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

733 invertedBackground = backgroundModel.clone() 

734 invertedBackground.setParameters([-p for p in backgroundModel.getParameters()]) 

735 backgroundModel = invertedBackground 

736 

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

738 

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

740 psfMatchingKernel, 

741 templateMatched=False) 

742 

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

744 matchedTemplate=matchedTemplate, 

745 matchedScience=matchedScience, 

746 backgroundModel=backgroundModel, 

747 psfMatchingKernel=psfMatchingKernel) 

748 

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

750 templateMatched=True, 

751 preConvMode=False, 

752 preConvKernel=None, 

753 spatiallyVarying=False): 

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

755 caused by convolution. 

756 

757 Parameters 

758 ---------- 

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

760 Template exposure, warped to match the science exposure. 

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

762 Science exposure to subtract from the template. 

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

764 Result of subtracting template and science. 

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

766 An (optionally spatially-varying) PSF matching kernel 

767 templateMatched : `bool`, optional 

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

769 preConvMode : `bool`, optional 

770 Was the science image preconvolved with its own PSF 

771 before PSF matching the template? 

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

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

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

775 spatiallyVarying : `bool`, optional 

776 Compute the decorrelation kernel spatially varying across the image? 

777 

778 Returns 

779 ------- 

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

781 The decorrelated image difference. 

782 """ 

783 if self.config.doDecorrelation: 

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

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

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

787 # during decorrelation 

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

789 templateMatched=templateMatched, 

790 preConvMode=preConvMode, 

791 preConvKernel=preConvKernel, 

792 spatiallyVarying=spatiallyVarying).correctedExposure 

793 else: 

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

795 correctedExposure = difference 

796 return correctedExposure 

797 

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

799 """Calculate the point-source limiting magnitude of an exposure. 

800 

801 The limit is the flux of a point source detected at ``nsigma`` in 

802 the median pixel noise, integrated over the effective area of the 

803 PSF, so it does not depend on the units of the pixels. 

804 

805 Parameters 

806 ---------- 

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

808 The exposure to calculate the limiting magnitude for. 

809 nsigma : `float`, optional 

810 The detection threshold in sigma. 

811 fallbackPsfSize : `float`, optional 

812 PSF FWHM in pixels to use if the exposure PSF cannot be 

813 evaluated. 

814 

815 Returns 

816 ------- 

817 maglim : `float` 

818 The limiting AB magnitude, or `numpy.nan` if the exposure has no 

819 photometric calibration, no usable pixels, or no PSF. 

820 """ 

821 if exposure.photoCalib is None: 821 ↛ 822line 821 didn't jump to line 822 because the condition on line 821 was never true

822 return np.nan 

823 try: 

824 psfArea = computePSFNoiseEquivalentArea(exposure.getPsf()) 

825 except (lsst.pex.exceptions.InvalidParameterError, 

826 afwDetection.InvalidPsfError, 

827 lsst.pex.exceptions.RangeError): 

828 if fallbackPsfSize is None: 

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

830 return np.nan 

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

832 # Effective area of a Gaussian PSF with this FWHM. 

833 psfArea = 4*np.pi*(fallbackPsfSize/np.sqrt(8*np.log(2)))**2 

834 mask = exposure.mask 

835 good = (mask.array & mask.getPlaneBitMask(self.config.badMaskPlanes)) == 0 

836 variance = exposure.variance.array[good] 

837 variance = variance[np.isfinite(variance) & (variance > 0)] 

838 if variance.size == 0: 

839 return np.nan 

840 fluxLim = nsigma*np.sqrt(np.median(variance)*psfArea) 

841 return exposure.photoCalib.instFluxToMagnitude(fluxLim) 

842 

843 @staticmethod 

844 def _validateExposures(template, science): 

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

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

847 

848 Parameters 

849 ---------- 

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

851 Template exposure, warped to match the science exposure. 

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

853 Science exposure to subtract from the template. 

854 

855 Raises 

856 ------ 

857 AssertionError 

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

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

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

861 """ 

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

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

864 templateBBox = template.getBBox() 

865 scienceBBox = science.getBBox() 

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

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

868 (science.filter, template.filter) 

869 

870 assert templateBBox.contains(scienceBBox), \ 

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

872 

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

874 bbox=None, 

875 psf=None, 

876 photoCalib=None, 

877 interpolateBadMaskPlanes=False, 

878 ): 

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

880 

881 Parameters 

882 ---------- 

883 exposure : `lsst.afw.Exposure` 

884 exposure to convolve. 

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

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

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

888 Configuration for convolve algorithm. 

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

890 Bounding box to trim the convolved exposure to. 

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

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

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

894 Photometric calibration of the convolved exposure. 

895 interpolateBadMaskPlanes : `bool`, optional 

896 If set, interpolate over mask planes specified in 

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

898 

899 Returns 

900 ------- 

901 convolvedExp : `lsst.afw.Exposure` 

902 The convolved image.