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 

23import numpy as np 

24 

25import lsst.afw.detection as afwDetection 

26import lsst.afw.image 

27import lsst.afw.math 

28import lsst.geom 

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

30 checkMask, setSourceFootprints) 

31from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask 

32import lsst.pex.config 

33import lsst.pipe.base 

34import lsst.pex.exceptions 

35from lsst.pipe.base import connectionTypes 

36from . import MakeKernelTask, DecorrelateALKernelTask 

37from lsst.utils.timer import timeMethod 

38 

39__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask", 

40 "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask", 

41 "SimplifiedSubtractConfig", "SimplifiedSubtractTask", 

42 "InsufficientKernelSourcesError"] 

43 

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

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

46 

47 

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

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

50 kernel. 

51 """ 

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

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

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

55 super().__init__(msg) 

56 self.nSources = nSources 

57 self.nRequired = nRequired 

58 

59 @property 

60 def metadata(self): 

61 return {"nSources": self.nSources, 

62 "nRequired": self.nRequired 

63 } 

64 

65 

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

67 dimensions=_dimensions, 

68 defaultTemplates=_defaultTemplates): 

69 template = connectionTypes.Input( 

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

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

72 storageClass="ExposureF", 

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

74 ) 

75 science = connectionTypes.Input( 

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

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

78 storageClass="ExposureF", 

79 name="{fakesType}calexp" 

80 ) 

81 sources = connectionTypes.Input( 

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

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

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

85 storageClass="SourceCatalog", 

86 name="{fakesType}src" 

87 ) 

88 visitSummary = connectionTypes.Input( 

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

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

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

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

93 storageClass="ExposureCatalog", 

94 name="finalVisitSummary", 

95 ) 

96 

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

98 super().__init__(config=config) 

99 if not config.doApplyExternalCalibrations: 

100 del self.visitSummary 

101 

102 

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

104 dimensions=_dimensions, 

105 defaultTemplates=_defaultTemplates): 

106 difference = connectionTypes.Output( 

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

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

109 storageClass="ExposureF", 

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

111 ) 

112 matchedTemplate = connectionTypes.Output( 

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

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

115 storageClass="ExposureF", 

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

117 ) 

118 psfMatchingKernel = connectionTypes.Output( 

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

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

121 storageClass="MatchingKernel", 

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

123 ) 

124 kernelSources = connectionTypes.Output( 

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

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

127 storageClass="SourceCatalog", 

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

129 ) 

130 

131 

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

133 dimensions=_dimensions, 

134 defaultTemplates=_defaultTemplates): 

135 scoreExposure = connectionTypes.Output( 

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

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

138 storageClass="ExposureF", 

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

140 ) 

141 psfMatchingKernel = connectionTypes.Output( 

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

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

144 storageClass="MatchingKernel", 

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

146 ) 

147 kernelSources = connectionTypes.Output( 

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

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

150 storageClass="SourceCatalog", 

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

152 ) 

153 

154 

155class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

156 pass 

157 

158 

159class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

160 inputPsfMatchingKernel = connectionTypes.Input( 

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

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

163 storageClass="MatchingKernel", 

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

165 ) 

166 

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

168 super().__init__(config=config) 

169 del self.sources 

170 if config.useExistingKernel: 

171 del self.psfMatchingKernel 

172 del self.kernelSources 

173 else: 

174 del self.inputPsfMatchingKernel 

175 

176 

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

178 makeKernel = lsst.pex.config.ConfigurableField( 

179 target=MakeKernelTask, 

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

181 ) 

182 doDecorrelation = lsst.pex.config.Field( 

183 dtype=bool, 

184 default=True, 

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

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

187 ) 

188 decorrelate = lsst.pex.config.ConfigurableField( 

189 target=DecorrelateALKernelTask, 

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

191 ) 

192 requiredTemplateFraction = lsst.pex.config.Field( 

193 dtype=float, 

194 default=0.1, 

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

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

197 ) 

198 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field( 

199 dtype=float, 

200 default=0.2, 

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

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

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

204 ) 

205 doScaleVariance = lsst.pex.config.Field( 

206 dtype=bool, 

207 default=True, 

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

209 ) 

210 scaleVariance = lsst.pex.config.ConfigurableField( 

211 target=ScaleVarianceTask, 

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

213 ) 

214 doSubtractBackground = lsst.pex.config.Field( 

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

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

217 dtype=bool, 

218 default=False, 

219 ) 

220 doApplyExternalCalibrations = lsst.pex.config.Field( 

221 doc=( 

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

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

224 ), 

225 dtype=bool, 

226 default=False, 

227 ) 

228 sourceSelector = lsst.pex.config.ConfigurableField( 

229 target=ScienceSourceSelectorTask, 

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

231 ) 

232 fallbackSourceSelector = lsst.pex.config.ConfigurableField( 

233 target=ScienceSourceSelectorTask, 

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

235 "Used only if the kernel calculation fails and" 

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

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

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

239 ) 

240 detectionThreshold = lsst.pex.config.Field( 

241 dtype=float, 

242 default=10, 

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

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

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

246 ) 

247 detectionThresholdMax = lsst.pex.config.Field( 

248 dtype=float, 

249 default=500, 

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

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

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

253 ) 

254 restrictKernelEdgeSources = lsst.pex.config.Field( 

255 dtype=bool, 

256 default=True, 

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

258 ) 

259 maxKernelSources = lsst.pex.config.Field( 

260 dtype=int, 

261 default=1000, 

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

263 "Set to -1 to disable." 

264 ) 

265 minKernelSources = lsst.pex.config.Field( 

266 dtype=int, 

267 default=3, 

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

269 ) 

270 excludeMaskPlanes = lsst.pex.config.ListField( 

271 dtype=str, 

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

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

274 ) 

275 badMaskPlanes = lsst.pex.config.ListField( 

276 dtype=str, 

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

278 doc="Mask planes to interpolate over." 

279 ) 

280 preserveTemplateMask = lsst.pex.config.ListField( 

281 dtype=str, 

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

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

284 ) 

285 renameTemplateMask = lsst.pex.config.ListField( 

286 dtype=str, 

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

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

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

290 ) 

291 allowKernelSourceDetection = lsst.pex.config.Field( 

292 dtype=bool, 

293 default=False, 

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

295 " encountered while calculating the matching kernel." 

296 ) 

297 

298 def setDefaults(self): 

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

300 # Always include background fitting in the kernel fit, 

301 # even if it is not subtracted 

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

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

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

305 # Shared source selector settings 

306 doSkySources = False # Do not include sky sources 

307 doSignalToNoise = True # apply signal to noise filter 

308 doUnresolved = True # apply star-galaxy separation 

309 signalToNoiseMinimum = 10 

310 signalToNoiseMaximum = 500 

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

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

313 self.sourceSelector.doUnresolved = doUnresolved 

314 self.sourceSelector.doSkySources = doSkySources 

315 self.sourceSelector.doSignalToNoise = doSignalToNoise 

316 self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

317 self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

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

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

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

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

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

323 self.fallbackSourceSelector.doUnresolved = doUnresolved 

324 self.fallbackSourceSelector.doSkySources = doSkySources 

325 self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise 

326 self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

327 self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

328 

329 

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

331 pipelineConnections=AlardLuptonSubtractConnections): 

332 mode = lsst.pex.config.ChoiceField( 

333 dtype=str, 

334 default="convolveTemplate", 

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

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

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

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

339 ) 

340 

341 

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

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

344 the Alard & Lupton (1998) algorithm. 

345 """ 

346 ConfigClass = AlardLuptonSubtractConfig 

347 _DefaultName = "alardLuptonSubtract" 

348 

349 def __init__(self, **kwargs): 

350 super().__init__(**kwargs) 

351 self.makeSubtask("decorrelate") 

352 self.makeSubtask("makeKernel") 

353 self.makeSubtask("sourceSelector") 

354 self.makeSubtask("fallbackSourceSelector") 

355 if self.config.doScaleVariance: 

356 self.makeSubtask("scaleVariance") 

357 

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

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

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

361 self.convolutionControl.setDoNormalize(False) 

362 self.convolutionControl.setDoCopyEdge(True) 

363 

364 def _applyExternalCalibrations(self, exposure, visitSummary): 

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

366 external ones.". 

367 

368 Parameters 

369 ---------- 

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

371 Input exposure to adjust calibrations. 

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

373 Exposure catalog with external calibrations to be applied. Catalog 

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

375 lookup. 

376 

377 Returns 

378 ------- 

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

380 Exposure with adjusted calibrations. 

381 """ 

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

383 

384 row = visitSummary.find(detectorId) 

385 if row is None: 

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

387 "Using original calibrations.", detectorId) 

388 else: 

389 psf = row.getPsf() 

390 apCorrMap = row.getApCorrMap() 

391 if psf is None: 

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

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

394 detectorId) 

395 elif apCorrMap is None: 

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

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

398 detectorId) 

399 else: 

400 exposure.setPsf(psf) 

401 exposure.info.setApCorrMap(apCorrMap) 

402 

403 return exposure 

404 

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

406 inputs = butlerQC.get(inputRefs) 

407 

408 try: 

409 results = self.run(**inputs) 

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

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

412 # No partial outputs for butler to put 

413 raise error from e 

414 

415 butlerQC.put(results, outputRefs) 

416 

417 @timeMethod 

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

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

420 

421 Parameters 

422 ---------- 

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

424 Template exposure, warped to match the science exposure. 

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

426 Science exposure to subtract from the template. 

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

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

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

430 images around them. 

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

432 Exposure catalog with external calibrations to be applied. Catalog 

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

434 lookup. 

435 

436 Returns 

437 ------- 

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

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

440 Result of subtracting template and science. 

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

442 Warped and PSF-matched template exposure. 

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

444 Background model that was fit while solving for the 

445 PSF-matching kernel 

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

447 Kernel used to PSF-match the convolved image. 

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

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

450 PSF-matching kernel. 

451 """ 

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

453 

454 convolveTemplate = self.chooseConvolutionMethod(template, science) 

455 

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

457 convolveTemplate=convolveTemplate, 

458 runSourceDetection=False) 

459 

460 if self.config.doSubtractBackground: 

461 backgroundModel = kernelResult.backgroundModel 

462 else: 

463 backgroundModel = None 

464 if convolveTemplate: 

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

466 backgroundModel=backgroundModel) 

467 else: 

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

469 backgroundModel=backgroundModel) 

470 subtractResults.kernelSources = kernelResult.kernelSources 

471 

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

473 

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

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

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

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

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

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

480 self.metadata["differenceFootprintRatioMean"], 

481 self.metadata["differenceFootprintRatioStdev"]) 

482 

483 return subtractResults 

484 

485 def chooseConvolutionMethod(self, template, science): 

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

487 matching kernel. 

488 

489 Parameters 

490 ---------- 

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

492 Template exposure, warped to match the science exposure. 

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

494 Science exposure to subtract from the template. 

495 

496 Returns 

497 ------- 

498 convolveTemplate : `bool` 

499 Convolve the template to match the two images? 

500 

501 Raises 

502 ------ 

503 RuntimeError 

504 If an unsupported convolution mode is supplied. 

505 """ 

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

507 convolveTemplate = _shapeTest(template, 

508 science, 

509 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer, 

510 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid) 

511 if convolveTemplate: 

512 if self.sciencePsfSize < self.templatePsfSize: 

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

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

515 else: 

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

517 else: 

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

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

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

521 convolveTemplate = True 

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

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

524 convolveTemplate = False 

525 else: 

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

527 return convolveTemplate 

528 

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

530 """Construct the PSF-matching kernel. 

531 

532 Parameters 

533 ---------- 

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

535 Template exposure, warped to match the science exposure. 

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

537 Science exposure to subtract from the template. 

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

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

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

541 images around them. 

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

543 convolveTemplate : `bool`, optional 

544 Construct the matching kernel to convolve the template? 

545 runSourceDetection : `bool`, optional 

546 Run a minimal version of source detection to determine kernel 

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

548 from must be supplied. 

549 

550 Returns 

551 ------- 

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

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

554 Background model that was fit while solving for the 

555 PSF-matching kernel 

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

557 Kernel used to PSF-match the convolved image. 

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

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

560 PSF-matching kernel. 

561 """ 

562 

563 if convolveTemplate: 

564 reference = template 

565 target = science 

566 referenceFwhmPix = self.templatePsfSize 

567 targetFwhmPix = self.sciencePsfSize 

568 else: 

569 reference = science 

570 target = template 

571 referenceFwhmPix = self.sciencePsfSize 

572 targetFwhmPix = self.templatePsfSize 

573 try: 

574 if runSourceDetection: 

575 kernelSources = self.runKernelSourceDetection(template, science) 

576 else: 

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

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

579 preconvolved=False, 

580 templateFwhmPix=referenceFwhmPix, 

581 scienceFwhmPix=targetFwhmPix) 

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

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

584 # Raise NoWorkFound if template fraction is insufficient 

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

586 self.config.minTemplateFractionForExpectedSuccess, 

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

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

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

590 raise e 

591 

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

593 psfMatchingKernel=kernelResult.psfMatchingKernel, 

594 kernelSources=kernelSources) 

595 

596 def runKernelSourceDetection(self, template, science): 

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

598 to reject candidate sources. 

599 

600 Parameters 

601 ---------- 

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

603 Template exposure, warped to match the science exposure. 

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

605 Science exposure to subtract from the template. 

606 

607 Returns 

608 ------- 

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

610 Sources from the input catalog to use to construct the 

611 PSF-matching kernel. 

612 """ 

613 kernelSize = self.makeKernel.makeKernelBasisList( 

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

615 sigmaToFwhm = 2*np.log(2*np.sqrt(2)) 

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

617 candidateList=None, 

618 sigma=self.sciencePsfSize/sigmaToFwhm) 

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)