Coverage for tests/test_getTemplate.py: 94%

1# This file is part of ip_diffim. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import collections 

23import itertools 

24import unittest 

25import unittest.mock 

26 

27import astropy.units as u 

28from astropy.stats import gaussian_sigma_to_fwhm 

29import numpy as np 

30 

31import lsst.afw.cameraGeom.testUtils 

32import lsst.afw.detection 

33import lsst.afw.geom 

34import lsst.afw.image 

35import lsst.afw.math 

36import lsst.afw.table 

37from lsst.daf.butler import DataCoordinate, DatasetRef, DatasetType, DimensionUniverse 

38import lsst.geom 

39import lsst.images 

40import lsst.images.psfs 

41from lsst.images.tests import compare_masked_image_to_legacy 

42import lsst.ip.diffim 

43import lsst.meas.algorithms 

44import lsst.meas.base.tests 

45import lsst.pex.exceptions 

46import lsst.pipe.base as pipeBase 

47import lsst.pipe.base.testUtils 

48import lsst.skymap 

49import lsst.utils.tests 

50 

51from utils import generate_data_id, makeTestExposureRecord 

52 

53# Change this to True, `setup display_ds9`, and open ds9 (or use another afw 

54# display backend) to show the tract/patch layouts on the image. 

55debug = False 

56if debug: 56 ↛ 57line 56 didn't jump to line 57 because the condition on line 56 was never true

57 import lsst.afw.display 

58 display = lsst.afw.display.Display() 

59 display.frame = 1 

60 

61 

62def _showTemplate(box, template): 

63 """Show the corners of the template we made in this test.""" 

64 for point in box.getCorners(): 

65 display.dot("+", point.x, point.y, ctype="orange", size=40) 

66 display.frame = 2 

67 display.image(template, "warped template") 

68 display.frame = 3 

69 display.image(template.variance, "warped variance") 

70 

71 

72class GetTemplateTaskTestCase(lsst.utils.tests.TestCase): 

73 """Test that GetTemplateTask works on both one tract and multiple tract 

74 input coadd exposures. 

75 

76 Makes a synthetic exposure large enough to fit four small tracts with 2x2 

77 (300x300 pixel) patches each, extracts pixels for those patches by warping, 

78 and tests GetTemplateTask's output against boxes that overlap various 

79 combinations of one or multiple tracts. 

80 """ 

81 def setUp(self): 

82 self.scale = 0.2 # arcsec/pixel 

83 self.skymap = self._makeSkymap() 

84 self.patches = collections.defaultdict(list) 

85 self.dataIds = collections.defaultdict(list) 

86 self.coaddRefs = {} 

87 self.exposure = self._makeExposure() 

88 self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

89 

90 if debug: 90 ↛ 91line 90 didn't jump to line 91 because the condition on line 90 was never true

91 display.image(self.exposure, "base exposure") 

92 

93 for tract_id in range(4): 

94 tract = self.skymap.generateTract(tract_id) 

95 self._makePatches(tract) 

96 

97 def _makeSkymap(self): 

98 """Make a Skymap with 4 tracts with 4 patches each. 

99 """ 

100 tractScale = 0.02 # degrees 

101 # On-sky coordinates of the tract centers. 

102 coords = [(0, 0), 

103 (0, tractScale), 

104 (tractScale, 0), 

105 (tractScale, tractScale), 

106 ] 

107 config = lsst.skymap.DiscreteSkyMap.ConfigClass() 

108 config.raList = [c[0] for c in coords] 

109 config.decList = [c[1] for c in coords] 

110 # Half the tract center step size, to keep the tract overlap small. 

111 config.radiusList = [tractScale/2 for c in coords] 

112 config.projection = "TAN" 

113 config.pixelScale = self.scale 

114 config.tractOverlap = 0.0005 

115 config.tractBuilder = "legacy" 

116 config.tractBuilder["legacy"].patchInnerDimensions = (300, 300) 

117 config.tractBuilder["legacy"].patchBorder = 10 

118 return lsst.skymap.DiscreteSkyMap(config=config) 

119 

120 def _makeExposure(self): 

121 """Create a large image to break up into tracts and patches. 

122 

123 The image will have a source every 100 pixels in x and y, and a WCS 

124 that results in the tracts all fitting in the image, with tract=0 

125 in the lower left, tract=1 to the right, tract=2 above, and tract=3 

126 to the upper right. 

127 """ 

128 box = lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(800, 800)) 

129 # This WCS was constructed so that tract 0 mostly fills the lower left 

130 # quadrant of the image, and the other tracts fill the rest; slight 

131 # extra rotation as a check on the final warp layout, scaled by 5% 

132 # from the patch pixel scale. 

133 cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*self.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees) 

134 wcs = lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150), 

135 lsst.geom.SpherePoint(0, 0, lsst.geom.radians), 

136 cd_matrix) 

137 dataset = lsst.meas.base.tests.TestDataset(box, wcs=wcs) 

138 for x, y in itertools.product(np.arange(0, 500, 100), np.arange(0, 500, 100)): 

139 dataset.addSource(1e5, lsst.geom.Point2D(x, y)) 

140 exposure, _ = dataset.realize(2, dataset.makeMinimalSchema()) 

141 exposure.setFilter(lsst.afw.image.FilterLabel("a", "a_test")) 

142 return exposure 

143 

144 def _makePatches(self, tract): 

145 """Populate the patches and dataId dicts, keyed on tract id, with the 

146 warps of the main exposure and minimal dataIds, respectively. 

147 """ 

148 if debug: 148 ↛ 149line 148 didn't jump to line 149 because the condition on line 148 was never true

149 color = ['red', 'green', 'cyan', 'yellow'][tract.tract_id] 

150 point = self.exposure.wcs.skyToPixel(tract.ctr_coord) 

151 # Show the tract center, colored by tract id. 

152 display.dot("x", point.x, point.y, ctype=color, size=30) 

153 

154 # Use 5th order to minimize artifacts on the templates. 

155 config = lsst.afw.math.Warper.ConfigClass() 

156 config.warpingKernelName = "lanczos5" 

157 warper = lsst.afw.math.Warper.fromConfig(config) 

158 for patchId in range(tract.num_patches.x*tract.num_patches.y): 

159 patch = tract.getPatchInfo(patchId) 

160 box = patch.getOuterBBox() 

161 

162 if debug: 162 ↛ 164line 162 didn't jump to line 164 because the condition on line 162 was never true

163 # Show the patch corners as patch ids, colored by tract id. 

164 points = self.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x) 

165 for x in box.getCorners()])) 

166 for p in points: 

167 display.dot(patchId, p.x, p.y, ctype=color) 

168 

169 # This is mostly taken from drp_tasks makePsfMatchedWarp, but 

170 # ip_diffim cannot depend on drp_tasks. 

171 xyTransform = lsst.afw.geom.makeWcsPairTransform(self.exposure.wcs, patch.wcs) 

172 warpedPsf = lsst.meas.algorithms.WarpedPsf(self.exposure.psf, xyTransform) 

173 warped = warper.warpExposure(patch.wcs, self.exposure, destBBox=box) 

174 warped.setPsf(warpedPsf) 

175 

176 warped.getInfo().setCoaddInputs( 

177 self._makeCoaddInputs([(self.exposure.wcs, self.exposure.getBBox())])) 

178 dataRef = pipeBase.InMemoryDatasetHandle( 

179 warped, 

180 storageClass="ExposureF", 

181 copy=True, 

182 dataId=generate_data_id( 

183 tract=tract, 

184 patch=patch, 

185 ) 

186 ) 

187 self.patches[tract.tract_id].append(dataRef) 

188 dataCoordinate = DataCoordinate.standardize({"tract": tract.tract_id, 

189 "patch": patchId, 

190 "band": "a", 

191 "skymap": "skymap"}, 

192 universe=DimensionUniverse()) 

193 self.dataIds[tract.tract_id].append(dataCoordinate) 

194 

195 def _checkMetadata(self, template, config, box, wcs, nPsfs): 

196 """Check that the various metadata components were set correctly. 

197 """ 

198 expectedBox = lsst.geom.Box2I(box) 

199 expectedBox.grow(config.templateBorderSize) 

200 self.assertEqual(template.getBBox(), expectedBox) 

201 # WCS should match our exposure, not any of the coadd tracts. 

202 for tract in self.patches: 

203 self.assertNotEqual(template.wcs, self.patches[tract][0].get().wcs) 

204 self.assertEqual(template.wcs, self.exposure.wcs) 

205 # The template pixels are calibrated to nJy. 

206 self.assertEqual(template.photoCalib, lsst.afw.image.PhotoCalib(1.0)) 

207 self.assertEqual(template.getXY0(), expectedBox.getMin()) 

208 self.assertEqual(template.filter.bandLabel, "a") 

209 self.assertEqual(template.filter.physicalLabel, "a_test") 

210 # The template PSF is a Gaussian with the width of the CoaddPsf of 

211 # its inputs. 

212 self.assertIsInstance(template.psf, lsst.afw.detection.GaussianPsf) 

213 coaddPsf = self._makeCoaddPsf(template, config) 

214 self.assertEqual(coaddPsf.getComponentCount(), nPsfs) 

215 position = coaddPsf.getAveragePosition() 

216 self.assertFloatsAlmostEqual(template.psf.getSigma(), 

217 coaddPsf.computeShape(position).getTraceRadius()) 

218 self.assertTrue(template.getInfo().hasCoaddInputs()) 

219 self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs) 

220 

221 def _makeCoaddPsf(self, template, config): 

222 """Return the CoaddPsf that ``run`` approximated with a Gaussian, 

223 rebuilt from the coadd inputs recorded on the template. 

224 """ 

225 task = lsst.ip.diffim.GetTemplateTask(config=config) 

226 return task._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.wcs) 

227 

228 def _checkPixels(self, template, config, box): 

229 """Check that the pixel values in the template are close to the 

230 original image, calibrated to nJy. 

231 """ 

232 # All pixels should have real values! 

233 expectedBox = lsst.geom.Box2I(box) 

234 expectedBox.grow(config.templateBorderSize) 

235 calibration = self.exposure.photoCalib.getCalibrationMean() 

236 expected = self.exposure.photoCalib.calibrateImage(self.exposure[expectedBox].maskedImage) 

237 

238 if debug: 238 ↛ 239line 238 didn't jump to line 239 because the condition on line 238 was never true

239 _showTemplate(expectedBox, template) 

240 

241 # Check that we fully filled the template from the patches. 

242 self.assertTrue(np.all(np.isfinite(template.image.array))) 

243 # Because of the scale changes, there will be some ringing in the 

244 # difference between the template and the original image; pick 

245 # tolerances large enough to account for that. 

246 self.assertImagesAlmostEqual(template.image, expected.image, rtol=.1, atol=4*calibration) 

247 # Variance plane ==4 in the original image (realize() takes a noise 

248 # sigma). Warping sets the level from the pixel areas and the two 

249 # warping kernels, which `_correctVariance` corrects up to the 

250 # difference between the configured coaddWarpKernel and the lanczos5 

251 # `_makePatches` really used. A per-pixel ripple that no scalar 

252 # correction can remove remains on top of that, so check the level 

253 # tightly and allow for the ripple around it. 

254 variance = template.variance.array[np.isfinite(template.variance.array)] 

255 median = np.median(variance) 

256 self.assertFloatsAlmostEqual(median, np.median(expected.variance.array), 

257 rtol=0.35, msg="variance level differs") 

258 self.assertLess(np.percentile(variance, 99)/median, 1.6, msg="variance ripple too large") 

259 self.assertGreater(np.percentile(variance, 1)/median, 0.6, msg="variance ripple too large") 

260 # Not checking the mask, as warping changes the sizes of the masks. 

261 

262 def testRunOneTractInput(self): 

263 """Test a bounding box that fully fits inside one tract, with only 

264 that tract passed as input. This checks that the code handles a single 

265 tract input correctly. 

266 """ 

267 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

268 task = lsst.ip.diffim.GetTemplateTask() 

269 # Restrict to tract 0, since the box fits in just that tract. 

270 result = task.run(coaddExposureHandles={0: self.patches[0]}, 

271 bbox=box, 

272 wcs=self.exposure.wcs, 

273 dataIds={0: self.dataIds[0]}, 

274 physical_filter="a_test") 

275 

276 # All 4 patches from tract 0 are included in this template. 

277 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 4) 

278 self._checkPixels(result.template, task.config, box) 

279 

280 def testRunOneTractMultipleInputs(self): 

281 """Test a bounding box that fully fits inside one tract but where 

282 multiple tracts were passed in. This checks that patches that are 

283 mostly NaN after warping are merged correctly in the output. 

284 """ 

285 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

286 task = lsst.ip.diffim.GetTemplateTask() 

287 result = task.run(coaddExposureHandles=self.patches, 

288 bbox=box, 

289 wcs=self.exposure.wcs, 

290 dataIds=self.dataIds, 

291 physical_filter="a_test") 

292 

293 # All 4 patches from two tracts are included in this template. 

294 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 6) 

295 self._checkPixels(result.template, task.config, box) 

296 

297 def testRunTwoTracts(self): 

298 """Test a bounding box that crosses tract boundaries. 

299 """ 

300 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) 

301 task = lsst.ip.diffim.GetTemplateTask() 

302 result = task.run(coaddExposureHandles=self.patches, 

303 bbox=box, 

304 wcs=self.exposure.wcs, 

305 dataIds=self.dataIds, 

306 physical_filter="a_test") 

307 

308 # All 4 patches from all 4 tracts are included in this template 

309 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9) 

310 self._checkPixels(result.template, task.config, box) 

311 

312 def testRunCalibratesToNanojansky(self): 

313 """Test that the template is in nJy whatever the calibration of the 

314 coadds, and that coadds already in nJy are not rescaled. 

315 """ 

316 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

317 calibration = self.exposure.photoCalib.getCalibrationMean() 

318 self.assertNotEqual(calibration, 1.0) 

319 task = lsst.ip.diffim.GetTemplateTask() 

320 templates = {} 

321 for photoCalib in (self.exposure.photoCalib, lsst.afw.image.PhotoCalib(1.0)): 

322 self._remakePatches(photoCalib) 

323 result = task.run(coaddExposureHandles={0: self.patches[0]}, 

324 bbox=box, 

325 wcs=self.exposure.wcs, 

326 dataIds={0: self.dataIds[0]}, 

327 physical_filter="a_test") 

328 self.assertEqual(result.template.photoCalib, lsst.afw.image.PhotoCalib(1.0)) 

329 templates[photoCalib.getCalibrationMean()] = result.template 

330 

331 self.assertFloatsAlmostEqual(templates[calibration].image.array, 

332 calibration*templates[1.0].image.array, rtol=1e-6) 

333 self.assertFloatsAlmostEqual(templates[calibration].variance.array, 

334 calibration**2*templates[1.0].variance.array, rtol=1e-6) 

335 

336 def _remakePatches(self, photoCalib): 

337 """Rebuild the coadd patches from the main exposure with a different 

338 photometric calibration. 

339 """ 

340 self.exposure.setPhotoCalib(photoCalib) 

341 self.patches.clear() 

342 self.dataIds.clear() 

343 for tract_id in range(4): 

344 self._makePatches(self.skymap.generateTract(tract_id)) 

345 

346 def testGaussianPsfGridFallback(self): 

347 """If the CoaddPsf cannot be evaluated at its average position, the 

348 Gaussian width is averaged over a grid of positions instead. 

349 """ 

350 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

351 task = lsst.ip.diffim.GetTemplateTask() 

352 result = task.run(coaddExposureHandles={0: self.patches[0]}, 

353 bbox=box, 

354 wcs=self.exposure.wcs, 

355 dataIds={0: self.dataIds[0]}, 

356 physical_filter="a_test") 

357 template = result.template 

358 template.setPsf(self._makeCoaddPsf(template, task.config)) 

359 expected = lsst.ip.diffim.utils.evaluateMeanPsfFwhm( 

360 template, 

361 fwhmExposureBuffer=task.config.fwhmExposureBuffer, 

362 fwhmExposureGrid=task.config.fwhmExposureGrid, 

363 )/gaussian_sigma_to_fwhm 

364 

365 error = lsst.pex.exceptions.InvalidParameterError("No inputs at the average position.") 

366 with unittest.mock.patch("lsst.ip.diffim.getTemplate.getPsfFwhm", side_effect=error): 

367 with self.assertLogs(task.log.name, level="INFO") as cm: 

368 psf = task._makeGaussianPsf(template) 

369 

370 self.assertIn("grid of points", "\n".join(cm.output)) 

371 self.assertIsInstance(psf, lsst.afw.detection.GaussianPsf) 

372 self.assertFloatsAlmostEqual(psf.getSigma(), expected) 

373 self.assertFloatsAlmostEqual(task.metadata["templatePsfSigma"], expected) 

374 

375 def testRunNoTemplate(self): 

376 """A bounding box that doesn't overlap the patches will raise. 

377 """ 

378 box = lsst.geom.Box2I(lsst.geom.Point2I(1200, 1200), lsst.geom.Point2I(1600, 1600)) 

379 task = lsst.ip.diffim.GetTemplateTask() 

380 with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"): 

381 task.run(coaddExposureHandles=self.patches, 

382 bbox=box, 

383 wcs=self.exposure.wcs, 

384 dataIds=self.dataIds, 

385 physical_filter="a_test") 

386 

387 def testMissingPatches(self): 

388 """Test that a missing patch results in an appropriate mask. 

389 

390 This fixes the bug reported on DM-44997 (image and variance were NaN 

391 but the mask was not set to NO_DATA for those pixels). 

392 """ 

393 # tract=0, patch=1 is the lower-left corner, as displayed in DS9. 

394 self.patches[0].pop(1) 

395 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

396 task = lsst.ip.diffim.GetTemplateTask() 

397 result = task.run(coaddExposureHandles=self.patches, 

398 bbox=box, 

399 wcs=self.exposure.wcs, 

400 dataIds=self.dataIds, 

401 physical_filter="a_test") 

402 no_data = (result.template.mask.array & result.template.mask.getPlaneBitMask("NO_DATA")) != 0 

403 self.assertTrue(np.isfinite(result.template.image.array).all()) 

404 self.assertTrue(np.isfinite(result.template.variance.array).all()) 

405 self.assertEqual(no_data.sum(), 20990) 

406 

407 @lsst.utils.tests.methodParameters( 

408 box=[ 

409 lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)), 

410 lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)), 

411 ], 

412 nInput=[8, 16], 

413 ) 

414 def testNanInputs(self, box=None, nInput=None): 

415 """Test that the template has finite values when some of the input 

416 pixels have NaN as variance. 

417 """ 

418 for tract, patchRefs in self.patches.items(): 

419 for patchRef in patchRefs: 

420 patchCoadd = patchRef.get() 

421 bbox = lsst.geom.Box2I() 

422 bbox.include(lsst.geom.Point2I(patchCoadd.getBBox().getCenter())) 

423 bbox.grow(3) 

424 patchCoadd.variance[bbox].array *= np.nan 

425 

426 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)) 

427 task = lsst.ip.diffim.GetTemplateTask() 

428 result = task.run(coaddExposureHandles=self.patches, 

429 bbox=box, 

430 wcs=self.exposure.wcs, 

431 dataIds=self.dataIds, 

432 physical_filter="a_test") 

433 if debug: 433 ↛ 434line 433 didn't jump to line 434 because the condition on line 433 was never true

434 _showTemplate(box, result.template) 

435 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9) 

436 # We just check that the pixel values are all finite. We cannot check that pixel values 

437 # in the template are closer to the original anymore. 

438 self.assertTrue(np.isfinite(result.template.image.array).all()) 

439 

440 def _runCorrection(self, doCorrectVariancePlateScale=False, doScaleVariance=False, 

441 handles=None): 

442 """Run the task on tract 0 with the named variance corrections. 

443 

444 Everything defaults to off, so each test turns on exactly what it is 

445 exercising. 

446 """ 

447 config = lsst.ip.diffim.GetTemplateTask.ConfigClass() 

448 config.doCorrectVariancePlateScale = doCorrectVariancePlateScale 

449 config.doScaleVariance = doScaleVariance 

450 task = lsst.ip.diffim.GetTemplateTask(config=config) 

451 result = task.run(coaddExposureHandles={0: handles or self.patches[0]}, 

452 bbox=lsst.geom.Box2I(self.varianceBox), 

453 wcs=self.exposure.wcs, 

454 dataIds={0: self.dataIds[0]}, 

455 physical_filter="a_test") 

456 return task, result.template 

457 

458 def _makeScaledWcs(self, factor): 

459 """Make a WCS like the base exposure's, but with its pixel scale 

460 multiplied by ``factor``. 

461 """ 

462 cdMatrix = lsst.afw.geom.makeCdMatrix(factor*1.05*self.scale*lsst.geom.arcseconds, 

463 93*lsst.geom.degrees) 

464 return lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150), 

465 lsst.geom.SpherePoint(0, 0, lsst.geom.radians), 

466 cdMatrix) 

467 

468 @staticmethod 

469 def _makeCoaddInputs(records): 

470 """Make a CoaddInputs holding the given input records. 

471 

472 Parameters 

473 ---------- 

474 records : `list` [`tuple` [`lsst.afw.geom.SkyWcs` or `None`, \ 

475 `lsst.geom.Box2I`]] 

476 The WCS and bbox to record for each input. A `None` WCS makes a 

477 record that `_plateScaleFactor` has to skip. May be empty, to 

478 simulate a coadd whose plate scale cannot be reconstructed. 

479 """ 

480 ccdSchema = lsst.afw.table.ExposureTable.makeMinimalSchema() 

481 weightKey = ccdSchema.addField("weight", type=float, doc="Coadd weight") 

482 coaddInputs = lsst.afw.image.CoaddInputs( 

483 lsst.afw.table.ExposureTable.makeMinimalSchema(), ccdSchema) 

484 for wcs, bbox in records: 

485 record = coaddInputs.ccds.addNew() 

486 record.setWcs(wcs) 

487 record.setBBox(bbox) 

488 # Included because real coadds have it, though a single-record 

489 # correction does not use it. 

490 record.set(weightKey, 1.0) 

491 return coaddInputs 

492 

493 def _patchHandles(self, tract, records): 

494 """Return handles for a tract's patches, with their CoaddInputs 

495 replaced by ``records``. 

496 """ 

497 handles = [] 

498 for ref in self.patches[tract]: 

499 coadd = ref.get() 

500 coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records)) 

501 handles.append(pipeBase.InMemoryDatasetHandle( 

502 coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId)) 

503 return handles 

504 

505 def testCorrectVariancePlateScale(self): 

506 """The plate scale correction is the total pixel area change from the 

507 images the coadds were built from to the science image. 

508 """ 

509 _, off = self._runCorrection() 

510 

511 # The fixture's records are the science image itself, so there is no 

512 # net change in pixel area: the same-instrument case. 

513 task, on = self._runCorrection(doCorrectVariancePlateScale=True) 

514 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.0, rtol=1e-6) 

515 

516 # Coarser original pixels than science pixels, the DECam-template 

517 # case: the correction is the ratio of their areas. 

518 for scaleFactor in (1.315, 0.5): 

519 with self.subTest(scaleFactor=scaleFactor): 

520 handles = self._patchHandles( 

521 0, [(self._makeScaledWcs(scaleFactor), self.exposure.getBBox())]) 

522 task, on = self._runCorrection(doCorrectVariancePlateScale=True, 

523 handles=handles) 

524 factor = task.metadata["variancePlateScaleFactor"] 

525 self.assertFloatsAlmostEqual(factor, scaleFactor**2, rtol=1e-6) 

526 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor, 

527 rtol=1e-5, ignoreNaNs=True) 

528 

529 def testCorrectVariancePlateScaleUsesOneRecord(self): 

530 """Only the first usable coadd input is read. 

531 

532 The spread between records is just the local pixel scale, a few 

533 tenths of a percent across a real focal plane, so one stands for all 

534 of them. 

535 """ 

536 handles = self._patchHandles( 

537 0, [(None, self.exposure.getBBox()), 

538 (self._makeScaledWcs(2.0), self.exposure.getBBox()), 

539 (self.exposure.wcs, self.exposure.getBBox())]) 

540 task, _ = self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) 

541 

542 # The first record with a WCS: the one without is skipped, and the 

543 # last (which would give 1.0) is never reached. 

544 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 4.0, rtol=1e-6) 

545 

546 def testCorrectVariancePlateScaleNeedsCoaddInputs(self): 

547 """Without usable coadd inputs the plate scale cannot be 

548 reconstructed, and the task must say so rather than silently applying 

549 only part of the correction. 

550 """ 

551 handles = self._patchHandles(0, []) 

552 with self.assertRaisesRegex(RuntimeError, "doCorrectVariancePlateScale"): 

553 self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) 

554 

555 def _scaleInputVariance(self, tract, factor): 

556 """Return fresh handles for one tract's patches, with their variance 

557 planes multiplied by ``factor``. 

558 

559 Parameters 

560 ---------- 

561 tract : `int` 

562 Id of the tract whose patches should be copied. 

563 factor : `float` 

564 Factor to multiply the input variance planes by. 

565 

566 Returns 

567 ------- 

568 handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`] 

569 Handles to the modified patches. 

570 """ 

571 handles = [] 

572 for handle in self.patches[tract]: 

573 # ``copy=True`` on the original handles means this is a copy, so 

574 # the patches shared with the other tests are left untouched. 

575 patch = handle.get() 

576 patch.variance.array *= factor 

577 handles.append(pipeBase.InMemoryDatasetHandle(patch, 

578 storageClass="ExposureF", 

579 copy=True, 

580 dataId=handle.dataId)) 

581 return handles 

582 

583 def testScaleVariance(self): 

584 """Test that the template variance plane is rescaled to match the 

585 empirical pixel noise, and that the factor used is recorded in the 

586 task metadata. 

587 """ 

588 scaleFactor = 1.345 

589 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

590 

591 def _configureAndRunTask(doScaleVariance, varianceScale=1.): 

592 """Build a template from tract 0, optionally rescaling the input 

593 variance planes by ``varianceScale`` first. 

594 """ 

595 config = lsst.ip.diffim.GetTemplateTask.ConfigClass() 

596 config.doScaleVariance = doScaleVariance 

597 task = lsst.ip.diffim.GetTemplateTask(config=config) 

598 result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)}, 

599 bbox=box, 

600 wcs=self.exposure.wcs, 

601 dataIds={0: self.dataIds[0]}, 

602 physical_filter="a_test") 

603 return task, result.template 

604 

605 # With scaling disabled the subtask is never constructed, and nothing 

606 # is recorded in the metadata. 

607 taskOff, templateOff = _configureAndRunTask(False) 

608 self.assertFalse(hasattr(taskOff, "scaleVariance")) 

609 self.assertNotIn("scaleTemplateVarianceFactor", taskOff.metadata) 

610 

611 # Both warps -- lanczos5 in ``_makePatches`` and lanczos3 in the 

612 # task -- correlate the noise. The variance plane tracks only the 

613 # per-pixel diagonal, which the second warp leaves too low, so 

614 # ``scaleVariance`` measures a factor well above 1 even though the 

615 # input variance planes are correct. 

616 # 

617 taskOn, templateOn = _configureAndRunTask(True) 

618 factor = taskOn.metadata["scaleTemplateVarianceFactor"] 

619 # TODO DM-55879: this value is pinned on purpose. The lanczos warping 

620 # kernels introduce small correlations that artificially suppress the 

621 # image pixel stddev and inflate the variance scaling factor. This 

622 # should be changed to 1.0 after DM-55879 is merged. 

623 self.assertFloatsAlmostEqual(factor, 1.1465, atol=0.01, 

624 msg="Measured template variance scaling changed; see the" 

625 " comment above if the correlation correction landed.") 

626 # The only difference from the unscaled template is the constant 

627 # factor applied to the variance plane. 

628 self.assertFloatsAlmostEqual(templateOn.variance.array, 

629 templateOff.variance.array*factor, rtol=1e-5) 

630 # Tolerance here is float32 round-off: repeated runs of the task are 

631 # not bitwise identical. 

632 self.assertImagesAlmostEqual(templateOn.image, templateOff.image, rtol=1e-5, atol=1e-5) 

633 

634 # If the input variance planes under-estimate the noise by a known 

635 # factor, the measured factor grows by that amount and the same 

636 # output variance plane is recovered. 

637 taskLow, templateLow = _configureAndRunTask(True, varianceScale=1/scaleFactor) 

638 self.assertFloatsAlmostEqual(taskLow.metadata["scaleTemplateVarianceFactor"], 

639 factor*scaleFactor, rtol=1e-5) 

640 self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) 

641 

642 def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): 

643 """Build a template from tract 0 with a task configured for the 

644 future output type, without converting it. 

645 

646 The task is kept as ``self.futureTask``. 

647 

648 Parameters 

649 ---------- 

650 raiseOnUndefinedMaskMap : `bool`, optional 

651 Value of the task config field of the same name. 

652 

653 Returns 

654 ------- 

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

656 The legacy output struct. 

657 box : `lsst.geom.Box2I` 

658 The bounding box the template was requested on, before the task 

659 grew it by the template border. 

660 """ 

661 config = lsst.ip.diffim.GetTemplateTask.ConfigClass() 

662 config.output_image_type = "future" 

663 config.raiseOnUndefinedMaskMap = raiseOnUndefinedMaskMap 

664 task = lsst.ip.diffim.GetTemplateTask(config=config) 

665 self.visit = 9876 

666 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) 

667 result = task.run(coaddExposureHandles={0: self.patches[0]}, 

668 bbox=box, 

669 wcs=self.exposure.wcs, 

670 dataIds={0: self.dataIds[0]}, 

671 physical_filter="a_test", 

672 visit=self.visit) 

673 # run() always produces legacy types, whatever output_image_type is. 

674 self.assertIsInstance(result.template, lsst.afw.image.ExposureF) 

675 self.futureTask = task 

676 return result, box 

677 

678 def _runFuture(self): 

679 """Build a template and convert it, the way ``runQuantum`` does. 

680 

681 Returns 

682 ------- 

683 result : `lsst.pipe.base.Struct` 

684 The converted output struct. 

685 box : `lsst.geom.Box2I` 

686 The bounding box the template was requested on, before the task 

687 grew it by the template border. 

688 """ 

689 # lsst.images needs per-amplifier raw geometry and a field angle 

690 # transform, which the trivial test detectors do not have. 

691 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] 

692 result, box = self._runLegacyForFuture() 

693 result.template = self.futureTask.convert_outputs_to_future(