Coverage for tests/test_getTemplate.py: 95%

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, handles=None, instrument="SciCam"): 

441 """Run the task on tract 0, with a science image from ``instrument`` 

442 and variance scaling turned off. 

443 

444 The plate scale correction is applied only when ``instrument`` is not 

445 the instrument of the template coadds, which is ``"TestCam"`` unless 

446 ``handles`` says otherwise. 

447 """ 

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

449 config.doScaleVariance = False 

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 instrument=instrument) 

457 return task, result.template 

458 

459 def _makeScaledWcs(self, factor): 

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

461 multiplied by ``factor``. 

462 """ 

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

464 93*lsst.geom.degrees) 

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

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

467 cdMatrix) 

468 

469 @staticmethod 

470 def _makeCoaddInputs(records, instrument="TestCam"): 

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

472 

473 Parameters 

474 ---------- 

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

476 `lsst.geom.Box2I`]] 

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

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

479 simulate a coadd whose plate scale cannot be reconstructed. 

480 instrument : `str` or `None`, optional 

481 Instrument label to record for every input. `None` leaves the 

482 inputs without a visit info, and an empty string gives them a 

483 visit info with no label; either way the instrument is unknown. 

484 """ 

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

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

487 coaddInputs = lsst.afw.image.CoaddInputs( 

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

489 for wcs, bbox in records: 

490 record = coaddInputs.ccds.addNew() 

491 record.setWcs(wcs) 

492 record.setBBox(bbox) 

493 if instrument is not None: 

494 record.setVisitInfo(lsst.afw.image.VisitInfo(instrumentLabel=instrument)) 

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

496 # correction does not use it. 

497 record.set(weightKey, 1.0) 

498 return coaddInputs 

499 

500 def _patchHandles(self, tract, records, instrument="TestCam"): 

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

502 replaced by ``records`` from ``instrument``. 

503 """ 

504 handles = [] 

505 for ref in self.patches[tract]: 

506 coadd = ref.get() 

507 coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records, instrument=instrument)) 

508 handles.append(pipeBase.InMemoryDatasetHandle( 

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

510 return handles 

511 

512 def testCorrectVariancePlateScale(self): 

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

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

515 """ 

516 # The science image and the coadds are from the same instrument, so 

517 # the factor is 1 without reading any coadd inputs. 

518 task, off = self._runCorrection(instrument="TestCam") 

519 self.assertEqual(task.metadata["variancePlateScaleFactor"], 1.0) 

520 self.assertNotIn("templateVarianceCorrectionFactor", task.metadata) 

521 

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

523 # net change in pixel area even when the instruments differ. 

524 task, on = self._runCorrection() 

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

526 

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

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

529 for scaleFactor in (1.315, 0.5): 

530 with self.subTest(scaleFactor=scaleFactor): 

531 handles = self._patchHandles( 

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

533 task, on = self._runCorrection(handles=handles) 

534 factor = task.metadata["variancePlateScaleFactor"] 

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

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

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

538 

539 def testCorrectVariancePlateScaleUsesOneRecord(self): 

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

541 

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

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

544 of them. 

545 """ 

546 handles = self._patchHandles( 

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

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

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

550 task, _ = self._runCorrection(handles=handles) 

551 

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

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

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

555 

556 def testCorrectVariancePlateScaleSkipsDroppedCoadd(self): 

557 """The plate scale comes from the first coadd that is used, not from 

558 one dropped for its instrument. 

559 """ 

560 bbox = self.exposure.getBBox() 

561 handles = (self._patchHandles(0, [(self._makeScaledWcs(2.0), bbox)], instrument="OtherCam")[:1] 

562 + self._patchHandles(0, [(self._makeScaledWcs(1.315), bbox)])[1:]) 

563 task, _ = self._runCorrection(handles=handles) 

564 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6) 

565 

566 def testCorrectVariancePlateScaleNextTract(self): 

567 """If the first coadd kept in a tract cannot give the plate scale, the 

568 first coadd kept in the next tract is used. 

569 """ 

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

571 handles = {0: self._patchHandles(0, []), 

572 1: self._patchHandles(1, [(self._makeScaledWcs(1.315), self.exposure.getBBox())])} 

573 task.run(coaddExposureHandles=handles, 

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

575 wcs=self.exposure.wcs, 

576 dataIds={0: self.dataIds[0], 1: self.dataIds[1]}, 

577 physical_filter="a_test", 

578 instrument="SciCam") 

579 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6) 

580 

581 def testCorrectVariancePlateScaleNeedsCoaddInputs(self): 

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

583 reconstructed, so the template cannot be built. 

584 """ 

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

586 with self.assertRaisesRegex(RuntimeError, "cannot be corrected for the change in pixel area"): 

587 self._runCorrection(handles=handles) 

588 

589 def _runTract0(self, handles): 

590 """Build a template from the given tract 0 patch handles. 

591 """ 

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

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

594 return task.run(coaddExposureHandles={0: handles}, 

595 bbox=box, 

596 wcs=self.exposure.wcs, 

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

598 physical_filter="a_test") 

599 

600 def testRunRecordsInstrument(self): 

601 """The instrument of each coadd is recorded with the template's coadd 

602 inputs. 

603 """ 

604 result = self._runTract0(self.patches[0]) 

605 ccds = result.template.getInfo().getCoaddInputs().ccds 

606 self.assertEqual(len(ccds), 4) 

607 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) 

608 

609 def testRunUnknownInstrument(self): 

610 """A coadd whose inputs record no instrument is dropped when the 

611 other coadds record one, and kept, with the template's instrument 

612 ``"UNKNOWN"``, when none of them records one. 

613 """ 

614 records = [(self.exposure.wcs, self.exposure.getBBox())] 

615 # Inputs with no visit info, and with a visit info but no label. 

616 for unknown in (None, ""): 

617 with self.subTest(label=unknown): 

618 handles = (self._patchHandles(0, records, instrument=unknown)[:1] 

619 + self._patchHandles(0, records)[1:]) 

620 ccds, dropped = self._runLogged({0: handles}) 

621 self.assertEqual(len(ccds), 3) 

622 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) 

623 self.assertEqual(len(dropped), 1) 

624 self.assertIn("made from UNKNOWN", dropped[0]) 

625 

626 # No coadd records an instrument. 

627 result = self._runTract0(self._patchHandles(0, records, instrument=unknown)) 

628 ccds = result.template.getInfo().getCoaddInputs().ccds 

629 self.assertEqual(len(ccds), 4) 

630 self.assertEqual({record["instrument"] for record in ccds}, {"UNKNOWN"}) 

631 

632 def _runLogged(self, handles, dataIds=None, box=None, instrument=None): 

633 """Build a template, capturing the messages about dropped coadds. 

634 

635 Parameters 

636 ---------- 

637 handles : `dict` [`int`, `list`] 

638 Patch handles, keyed by tract. 

639 dataIds : `dict` [`int`, `list`], optional 

640 Data ids matching ``handles``; defaults to the fixture's. 

641 box : `lsst.geom.Box2I`, optional 

642 Box to build the template on; defaults to one in tract 0. 

643 instrument : `str`, optional 

644 Instrument of the science image. 

645 

646 Returns 

647 ------- 

648 ccds : `lsst.afw.table.ExposureCatalog` 

649 The coadd inputs recorded on the template. 

650 dropped : `list` [`str`] 

651 The info messages about dropped coadds. 

652 """ 

653 if dataIds is None: 

654 dataIds = {tract: self.dataIds[tract] for tract in handles} 

655 if box is None: 

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

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

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

659 result = task.run(coaddExposureHandles=handles,