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# Sum of squared warping kernel values along one axis, with the kernel 

62# normalized to unit sum at each sub-pixel phase, averaged over phase. 

63KERNEL_SUM_SQ_1D = { 

64 "bilinear": 2/3, 

65 "lanczos3": 0.8933699084, 

66 "lanczos4": 0.9141249535, 

67 "lanczos5": 0.9342245513, 

68} 

69 

70 

71def _showTemplate(box, template): 

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

73 for point in box.getCorners(): 

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

75 display.frame = 2 

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

77 display.frame = 3 

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

79 

80 

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

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

83 input coadd exposures. 

84 

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

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

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

88 combinations of one or multiple tracts. 

89 """ 

90 scale = 0.2 # arcsec/pixel 

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

92 coaddWarpKernel = "lanczos5" 

93 

94 @classmethod 

95 def setUpClass(cls): 

96 # Warping the patches is the slowest part of every test, so it is done 

97 # once. Tests must not modify the exposure; the patch handles return 

98 # copies, and setUp gives each test its own lists of them. 

99 cls.skymap = cls._makeSkymap() 

100 cls.patches = collections.defaultdict(list) 

101 cls.dataIds = collections.defaultdict(list) 

102 cls.exposure = cls._makeExposure() 

103 

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

105 display.image(cls.exposure, "base exposure") 

106 

107 for tract_id in range(4): 

108 tract = cls.skymap.generateTract(tract_id) 

109 cls._makePatches(tract) 

110 

111 @classmethod 

112 def tearDownClass(cls): 

113 del cls.skymap, cls.patches, cls.dataIds, cls.exposure 

114 

115 def setUp(self): 

116 cls = type(self) 

117 self.patches = collections.defaultdict(list, {k: list(v) for k, v in cls.patches.items()}) 

118 self.dataIds = collections.defaultdict(list, {k: list(v) for k, v in cls.dataIds.items()}) 

119 self.coaddRefs = {} 

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

121 

122 @classmethod 

123 def _makeSkymap(cls): 

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

125 """ 

126 tractScale = 0.02 # degrees 

127 # On-sky coordinates of the tract centers. 

128 coords = [(0, 0), 

129 (0, tractScale), 

130 (tractScale, 0), 

131 (tractScale, tractScale), 

132 ] 

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

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

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

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

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

138 config.projection = "TAN" 

139 config.pixelScale = cls.scale 

140 config.tractOverlap = 0.0005 

141 config.tractBuilder = "legacy" 

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

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

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

145 

146 @classmethod 

147 def _makeExposure(cls): 

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

149 

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

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

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

153 to the upper right. 

154 """ 

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

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

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

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

159 # from the patch pixel scale. 

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

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

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

163 cd_matrix) 

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

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

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

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

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

169 return exposure 

170 

171 @classmethod 

172 def _makePatches(cls, tract): 

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

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

175 """ 

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

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

178 point = cls.exposure.wcs.skyToPixel(tract.ctr_coord) 

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

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

181 

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

183 config.warpingKernelName = cls.coaddWarpKernel 

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

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

186 patch = tract.getPatchInfo(patchId) 

187 box = patch.getOuterBBox() 

188 

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

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

191 points = cls.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x) 

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

193 for p in points: 

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

195 

196 # This is mostly taken from drp_tasks makePsfMatchedWarp, but 

197 # ip_diffim cannot depend on drp_tasks. 

198 xyTransform = lsst.afw.geom.makeWcsPairTransform(cls.exposure.wcs, patch.wcs) 

199 warpedPsf = lsst.meas.algorithms.WarpedPsf(cls.exposure.psf, xyTransform) 

200 warped = warper.warpExposure(patch.wcs, cls.exposure, destBBox=box) 

201 warped.setPsf(warpedPsf) 

202 

203 warped.getInfo().setCoaddInputs( 

204 cls._makeCoaddInputs([(cls.exposure.wcs, cls.exposure.getBBox())])) 

205 dataRef = pipeBase.InMemoryDatasetHandle( 

206 warped, 

207 storageClass="ExposureF", 

208 copy=True, 

209 dataId=generate_data_id( 

210 tract=tract, 

211 patch=patch, 

212 ) 

213 ) 

214 cls.patches[tract.tract_id].append(dataRef) 

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

216 "patch": patchId, 

217 "band": "a", 

218 "skymap": "skymap"}, 

219 universe=DimensionUniverse()) 

220 cls.dataIds[tract.tract_id].append(dataCoordinate) 

221 

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

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

224 """ 

225 expectedBox = lsst.geom.Box2I(box) 

226 expectedBox.grow(config.templateBorderSize) 

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

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

229 for tract in self.patches: 

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

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

232 # The template pixels are calibrated to nJy. 

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

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

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

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

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

238 # its inputs. 

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

240 coaddPsf = self._makeCoaddPsf(template, config) 

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

242 position = coaddPsf.getAveragePosition() 

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

244 coaddPsf.computeShape(position).getTraceRadius()) 

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

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

247 

248 def _makeConfig(self): 

249 """Make a GetTemplateTask config whose coaddWarpKernel is the kernel 

250 the test patches were made with. 

251 """ 

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

253 config.coaddWarpKernel = self.coaddWarpKernel 

254 return config 

255 

256 def _makeCoaddPsf(self, template, config): 

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

258 rebuilt from the coadd inputs recorded on the template. 

259 """ 

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

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

262 

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

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

265 original image, calibrated to nJy. 

266 """ 

267 # All pixels should have real values! 

268 expectedBox = lsst.geom.Box2I(box) 

269 expectedBox.grow(config.templateBorderSize) 

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

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

272 

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

274 _showTemplate(expectedBox, template) 

275 

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

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

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

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

280 # tolerances large enough to account for that. 

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

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

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

284 # warping kernels, which `_correctVariance` corrects. A per-pixel 

285 # ripple that no scalar correction can remove remains on top of that, 

286 # so check the level tightly and allow for the ripple around it. 

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

288 median = np.median(variance) 

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

290 rtol=0.02, msg="variance level differs") 

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

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

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

294 

295 def testRunOneTractInput(self): 

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

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

298 tract input correctly. 

299 """ 

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

301 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

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

304 bbox=box, 

305 wcs=self.exposure.wcs, 

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

307 physical_filter="a_test") 

308 

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

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

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

312 

313 def testRunOneTractMultipleInputs(self): 

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

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

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

317 """ 

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

319 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

321 bbox=box, 

322 wcs=self.exposure.wcs, 

323 dataIds=self.dataIds, 

324 physical_filter="a_test") 

325 

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

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

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

329 

330 def testRunTwoTracts(self): 

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

332 """ 

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

334 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

336 bbox=box, 

337 wcs=self.exposure.wcs, 

338 dataIds=self.dataIds, 

339 physical_filter="a_test") 

340 

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

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

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

344 

345 def testRunCalibratesToNanojansky(self): 

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

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

348 """ 

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

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

351 self.assertNotEqual(calibration, 1.0) 

352 # The same coadd pixels, labelled as already in nJy. 

353 uncalibrated = [] 

354 for handle in self.patches[0]: 

355 coadd = handle.get() 

356 coadd.setPhotoCalib(lsst.afw.image.PhotoCalib(1.0)) 

357 uncalibrated.append(pipeBase.InMemoryDatasetHandle(coadd, storageClass="ExposureF", copy=True, 

358 dataId=handle.dataId)) 

359 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

360 templates = {} 

361 for key, handles in ((calibration, self.patches[0]), (1.0, uncalibrated)): 

362 result = task.run(coaddExposureHandles={0: handles}, 

363 bbox=box, 

364 wcs=self.exposure.wcs, 

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

366 physical_filter="a_test") 

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

368 templates[key] = result.template 

369 

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

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

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

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

374 

375 def testGaussianPsfGridFallback(self): 

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

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

378 """ 

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

380 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

382 bbox=box, 

383 wcs=self.exposure.wcs, 

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

385 physical_filter="a_test") 

386 template = result.template 

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

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

389 template, 

390 fwhmExposureBuffer=task.config.fwhmExposureBuffer, 

391 fwhmExposureGrid=task.config.fwhmExposureGrid, 

392 )/gaussian_sigma_to_fwhm 

393 

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

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

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

397 psf = task._makeGaussianPsf(template) 

398 

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

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

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

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

403 

404 def testRunNoTemplate(self): 

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

406 """ 

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

408 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

410 task.run(coaddExposureHandles=self.patches, 

411 bbox=box, 

412 wcs=self.exposure.wcs, 

413 dataIds=self.dataIds, 

414 physical_filter="a_test") 

415 

416 def testMissingPatches(self): 

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

418 

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

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

421 """ 

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

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

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

425 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

427 bbox=box, 

428 wcs=self.exposure.wcs, 

429 dataIds=self.dataIds, 

430 physical_filter="a_test") 

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

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

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

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

435 

436 @lsst.utils.tests.methodParameters( 

437 box=[ 

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

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

440 ], 

441 nInput=[8, 16], 

442 ) 

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

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

445 pixels have NaN as variance. 

446 """ 

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

448 for patchRef in patchRefs: 

449 patchCoadd = patchRef.get() 

450 bbox = lsst.geom.Box2I() 

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

452 bbox.grow(3) 

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

454 

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

456 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig()) 

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

458 bbox=box, 

459 wcs=self.exposure.wcs, 

460 dataIds=self.dataIds, 

461 physical_filter="a_test") 

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

463 _showTemplate(box, result.template) 

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

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

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

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

468 

469 def _runCorrection(self, doCorrectVarianceWarpCorrelation=False, doScaleVariance=False, 

470 handles=None, instrument="SciCam", coaddWarpKernel=None, 

471 templateWarpKernel="lanczos3"): 

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

473 

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

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

476 ``handles`` says otherwise. ``coaddWarpKernel`` defaults to the kernel 

477 the patches were made with, so the reconstructed factors are exact. 

478 """ 

479 config = self._makeConfig() 

480 config.doCorrectVarianceWarpCorrelation = doCorrectVarianceWarpCorrelation 

481 config.doScaleVariance = doScaleVariance 

482 if coaddWarpKernel is not None: 

483 config.coaddWarpKernel = coaddWarpKernel 

484 config.warp.warpingKernelName = templateWarpKernel 

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

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

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

488 wcs=self.exposure.wcs, 

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

490 physical_filter="a_test", 

491 instrument=instrument) 

492 return task, result.template 

493 

494 def _makeScaledWcs(self, factor): 

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

496 multiplied by ``factor``. 

497 """ 

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

499 93*lsst.geom.degrees) 

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

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

502 cdMatrix) 

503 

504 @staticmethod 

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

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

507 

508 Parameters 

509 ---------- 

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

511 `lsst.geom.Box2I`]] 

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

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

514 simulate a coadd whose plate scale cannot be reconstructed. 

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

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

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

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

519 """ 

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

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

522 coaddInputs = lsst.afw.image.CoaddInputs( 

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

524 for wcs, bbox in records: 

525 record = coaddInputs.ccds.addNew() 

526 record.setWcs(wcs) 

527 record.setBBox(bbox) 

528 if instrument is not None: 

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

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

531 # correction does not use it. 

532 record.set(weightKey, 1.0) 

533 return coaddInputs 

534 

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

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

537 replaced by ``records`` from ``instrument``. 

538 """ 

539 handles = [] 

540 for ref in self.patches[tract]: 

541 coadd = ref.get() 

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

543 handles.append(pipeBase.InMemoryDatasetHandle( 

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

545 return handles 

546 

547 def testCorrectVarianceWarpCorrelation(self): 

548 """The correlation correction covers both warps: the one that built 

549 the coadds and the one this task performs. 

550 """ 

551 _, off = self._runCorrection(instrument="TestCam") 

552 for coaddKernel, templateKernel in (("lanczos3", "lanczos3"), 

553 ("lanczos5", "lanczos3"), 

554 ("lanczos5", "lanczos5"), 

555 ("bilinear", "lanczos3")): 

556 with self.subTest(coaddWarpKernel=coaddKernel, templateWarpKernel=templateKernel): 

557 task, on = self._runCorrection(doCorrectVarianceWarpCorrelation=True, 

558 instrument="TestCam", 

559 coaddWarpKernel=coaddKernel, 

560 templateWarpKernel=templateKernel) 

561 factor = task.metadata["varianceWarpCorrelationFactor"] 

562 # A kernel spreads power, so this always makes the plane 

563 # larger. 

564 self.assertGreater(factor, 1.0) 

565 # The phase average converges slowly for bilinear, because 

566 # its sum of squares has a kink where the phase wraps. 

567 rtol = 1e-4 if "bilinear" in (coaddKernel, templateKernel) else 1e-5 

568 self.assertFloatsAlmostEqual( 

569 factor, 

570 1/(KERNEL_SUM_SQ_1D[coaddKernel]*KERNEL_SUM_SQ_1D[templateKernel])**2, 

571 rtol=rtol) 

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

573 if templateKernel == "lanczos3": 

574 # A pure scalar, with no per-pixel structure of its own; 

575 # only comparable when `off` used the same template warp. 

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

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

578 

579 def testCorrectVariancePlateScale(self): 

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

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

582 """ 

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

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

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

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

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

588 

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

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

591 task, on = self._runCorrection() 

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

593 self.assertNotIn("varianceWarpCorrelationFactor", task.metadata) 

594 

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

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

597 for scaleFactor in (1.315, 0.5): 

598 with self.subTest(scaleFactor=scaleFactor): 

599 handles = self._patchHandles( 

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

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

602 factor = task.metadata["variancePlateScaleFactor"] 

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

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

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

606 

607 def testVarianceCorrectionsCompose(self): 

608 """The two corrections are independent and multiply together. 

609 

610 The patches are warped from the exposure and the template is warped 

611 back to its pixel grid, so the two warps cancel geometrically, leaving 

612 only the two kernels' correlation to undo. The substituted coadd 

613 inputs claim coarser original pixels than the exposure really has, so 

614 with both corrections applied the variance must come back to its 

615 original value times the claimed change in pixel area. 

616 """ 

617 plateScale = 1.315 

618 handles = self._patchHandles(0, [(self._makeScaledWcs(plateScale), self.exposure.getBBox())]) 

619 # A science image from the coadds' instrument applies no correction 

620 # at all; one from another instrument applies both. 

621 _, off = self._runCorrection(handles=handles, instrument="TestCam") 

622 task, on = self._runCorrection(doCorrectVarianceWarpCorrelation=True, handles=handles) 

623 

624 # Both together is exactly the product of the two factors. 

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

626 scale = (task.metadata["varianceWarpCorrelationFactor"] 

627 * task.metadata["variancePlateScaleFactor"]) 

628 self.assertFloatsAlmostEqual(task.metadata["templateVarianceCorrectionFactor"], 

629 scale, rtol=1e-12) 

630 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*scale, 

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

632 

633 # Only the variance plane is touched. Nothing about the correction 

634 # feeds back into how the patches are combined. 

635 self.assertImagesEqual(on.image, off.image) 

636 self.assertMasksEqual(on.mask, off.mask) 

637 

638 # realize() takes a noise sigma, so the original variance is 4, not 2, 

639 # before calibrating to nJy. 

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

641 original = calibration**2*np.median(self.exposure.variance.array) 

642 self.assertFloatsAlmostEqual(np.median(on.variance.array), original*plateScale**2, rtol=0.02) 

643 

644 # Which is a real correction: uncorrected, the plane is too small by 

645 # the two warping kernels' phase-averaged sums of squares. 

646 self.assertFloatsAlmostEqual( 

647 np.median(off.variance.array)/original, 

648 (KERNEL_SUM_SQ_1D[self.coaddWarpKernel]*KERNEL_SUM_SQ_1D["lanczos3"])**2, rtol=0.02) 

649 

650 def testCorrectVariancePlateScaleUsesOneRecord(self): 

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

652 

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

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

655 of them. 

656 """ 

657 handles = self._patchHandles( 

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

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

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

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

662 

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

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

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

666 

667 def testCorrectVariancePlateScaleSkipsDroppedCoadd(self): 

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