Coverage for tests/test_getTemplate.py: 91%
292 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-09 10:57 +0000
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/>.
22import collections
23import itertools
24import unittest
25import unittest.mock
27from astropy.stats import gaussian_sigma_to_fwhm
28import numpy as np
30import lsst.afw.detection
31import lsst.afw.geom
32import lsst.afw.image
33import lsst.afw.math
34import lsst.afw.table
35from lsst.daf.butler import DataCoordinate, DimensionUniverse
36import lsst.geom
37import lsst.ip.diffim
38import lsst.meas.algorithms
39import lsst.meas.base.tests
40import lsst.pex.exceptions
41import lsst.pipe.base as pipeBase
42import lsst.skymap
43import lsst.utils.tests
45from utils import generate_data_id
47# Change this to True, `setup display_ds9`, and open ds9 (or use another afw
48# display backend) to show the tract/patch layouts on the image.
49debug = False
50if debug: 50 ↛ 51line 50 didn't jump to line 51 because the condition on line 50 was never true
51 import lsst.afw.display
52 display = lsst.afw.display.Display()
53 display.frame = 1
56def _showTemplate(box, template):
57 """Show the corners of the template we made in this test."""
58 for point in box.getCorners():
59 display.dot("+", point.x, point.y, ctype="orange", size=40)
60 display.frame = 2
61 display.image(template, "warped template")
62 display.frame = 3
63 display.image(template.variance, "warped variance")
66class GetTemplateTaskTestCase(lsst.utils.tests.TestCase):
67 """Test that GetTemplateTask works on both one tract and multiple tract
68 input coadd exposures.
70 Makes a synthetic exposure large enough to fit four small tracts with 2x2
71 (300x300 pixel) patches each, extracts pixels for those patches by warping,
72 and tests GetTemplateTask's output against boxes that overlap various
73 combinations of one or multiple tracts.
74 """
75 def setUp(self):
76 self.scale = 0.2 # arcsec/pixel
77 self.skymap = self._makeSkymap()
78 self.patches = collections.defaultdict(list)
79 self.dataIds = collections.defaultdict(list)
80 self.exposure = self._makeExposure()
81 self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
83 if debug: 83 ↛ 84line 83 didn't jump to line 84 because the condition on line 83 was never true
84 display.image(self.exposure, "base exposure")
86 for tract_id in range(4):
87 tract = self.skymap.generateTract(tract_id)
88 self._makePatches(tract)
90 def _makeSkymap(self):
91 """Make a Skymap with 4 tracts with 4 patches each.
92 """
93 tractScale = 0.02 # degrees
94 # On-sky coordinates of the tract centers.
95 coords = [(0, 0),
96 (0, tractScale),
97 (tractScale, 0),
98 (tractScale, tractScale),
99 ]
100 config = lsst.skymap.DiscreteSkyMap.ConfigClass()
101 config.raList = [c[0] for c in coords]
102 config.decList = [c[1] for c in coords]
103 # Half the tract center step size, to keep the tract overlap small.
104 config.radiusList = [tractScale/2 for c in coords]
105 config.projection = "TAN"
106 config.pixelScale = self.scale
107 config.tractOverlap = 0.0005
108 config.tractBuilder = "legacy"
109 config.tractBuilder["legacy"].patchInnerDimensions = (300, 300)
110 config.tractBuilder["legacy"].patchBorder = 10
111 return lsst.skymap.DiscreteSkyMap(config=config)
113 def _makeExposure(self):
114 """Create a large image to break up into tracts and patches.
116 The image will have a source every 100 pixels in x and y, and a WCS
117 that results in the tracts all fitting in the image, with tract=0
118 in the lower left, tract=1 to the right, tract=2 above, and tract=3
119 to the upper right.
120 """
121 box = lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(800, 800))
122 # This WCS was constructed so that tract 0 mostly fills the lower left
123 # quadrant of the image, and the other tracts fill the rest; slight
124 # extra rotation as a check on the final warp layout, scaled by 5%
125 # from the patch pixel scale.
126 cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*self.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees)
127 wcs = lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150),
128 lsst.geom.SpherePoint(0, 0, lsst.geom.radians),
129 cd_matrix)
130 dataset = lsst.meas.base.tests.TestDataset(box, wcs=wcs)
131 for x, y in itertools.product(np.arange(0, 500, 100), np.arange(0, 500, 100)):
132 dataset.addSource(1e5, lsst.geom.Point2D(x, y))
133 exposure, _ = dataset.realize(2, dataset.makeMinimalSchema())
134 exposure.setFilter(lsst.afw.image.FilterLabel("a", "a_test"))
135 return exposure
137 def _makePatches(self, tract):
138 """Populate the patches and dataId dicts, keyed on tract id, with the
139 warps of the main exposure and minimal dataIds, respectively.
140 """
141 if debug: 141 ↛ 142line 141 didn't jump to line 142 because the condition on line 141 was never true
142 color = ['red', 'green', 'cyan', 'yellow'][tract.tract_id]
143 point = self.exposure.wcs.skyToPixel(tract.ctr_coord)
144 # Show the tract center, colored by tract id.
145 display.dot("x", point.x, point.y, ctype=color, size=30)
147 # Use 5th order to minimize artifacts on the templates.
148 config = lsst.afw.math.Warper.ConfigClass()
149 config.warpingKernelName = "lanczos5"
150 warper = lsst.afw.math.Warper.fromConfig(config)
151 for patchId in range(tract.num_patches.x*tract.num_patches.y):
152 patch = tract.getPatchInfo(patchId)
153 box = patch.getOuterBBox()
155 if debug: 155 ↛ 157line 155 didn't jump to line 157 because the condition on line 155 was never true
156 # Show the patch corners as patch ids, colored by tract id.
157 points = self.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x)
158 for x in box.getCorners()]))
159 for p in points:
160 display.dot(patchId, p.x, p.y, ctype=color)
162 # This is mostly taken from drp_tasks makePsfMatchedWarp, but
163 # ip_diffim cannot depend on drp_tasks.
164 xyTransform = lsst.afw.geom.makeWcsPairTransform(self.exposure.wcs, patch.wcs)
165 warpedPsf = lsst.meas.algorithms.WarpedPsf(self.exposure.psf, xyTransform)
166 warped = warper.warpExposure(patch.wcs, self.exposure, destBBox=box)
167 warped.setPsf(warpedPsf)
169 warped.getInfo().setCoaddInputs(
170 self._makeCoaddInputs([(self.exposure.wcs, self.exposure.getBBox())]))
171 dataRef = pipeBase.InMemoryDatasetHandle(
172 warped,
173 storageClass="ExposureF",
174 copy=True,
175 dataId=generate_data_id(
176 tract=tract,
177 patch=patch,
178 )
179 )
180 self.patches[tract.tract_id].append(dataRef)
181 dataCoordinate = DataCoordinate.standardize({"tract": tract.tract_id,
182 "patch": patchId,
183 "band": "a",
184 "skymap": "skymap"},
185 universe=DimensionUniverse())
186 self.dataIds[tract.tract_id].append(dataCoordinate)
188 def _checkMetadata(self, template, config, box, wcs, nPsfs):
189 """Check that the various metadata components were set correctly.
190 """
191 expectedBox = lsst.geom.Box2I(box)
192 expectedBox.grow(config.templateBorderSize)
193 self.assertEqual(template.getBBox(), expectedBox)
194 # WCS should match our exposure, not any of the coadd tracts.
195 for tract in self.patches:
196 self.assertNotEqual(template.wcs, self.patches[tract][0].get().wcs)
197 self.assertEqual(template.wcs, self.exposure.wcs)
198 # The template pixels are calibrated to nJy.
199 self.assertEqual(template.photoCalib, lsst.afw.image.PhotoCalib(1.0))
200 self.assertEqual(template.getXY0(), expectedBox.getMin())
201 self.assertEqual(template.filter.bandLabel, "a")
202 self.assertEqual(template.filter.physicalLabel, "a_test")
203 # The template PSF is a Gaussian with the width of the CoaddPsf of
204 # its inputs.
205 self.assertIsInstance(template.psf, lsst.afw.detection.GaussianPsf)
206 coaddPsf = self._makeCoaddPsf(template, config)
207 self.assertEqual(coaddPsf.getComponentCount(), nPsfs)
208 position = coaddPsf.getAveragePosition()
209 self.assertFloatsAlmostEqual(template.psf.getSigma(),
210 coaddPsf.computeShape(position).getTraceRadius())
211 self.assertTrue(template.getInfo().hasCoaddInputs())
212 self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs)
214 def _makeCoaddPsf(self, template, config):
215 """Return the CoaddPsf that ``run`` approximated with a Gaussian,
216 rebuilt from the coadd inputs recorded on the template.
217 """
218 task = lsst.ip.diffim.GetTemplateTask(config=config)
219 return task._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.wcs)
221 def _checkPixels(self, template, config, box):
222 """Check that the pixel values in the template are close to the
223 original image, calibrated to nJy.
224 """
225 # All pixels should have real values!
226 expectedBox = lsst.geom.Box2I(box)
227 expectedBox.grow(config.templateBorderSize)
228 calibration = self.exposure.photoCalib.getCalibrationMean()
229 expected = self.exposure.photoCalib.calibrateImage(self.exposure[expectedBox].maskedImage)
231 if debug: 231 ↛ 232line 231 didn't jump to line 232 because the condition on line 231 was never true
232 _showTemplate(expectedBox, template)
234 # Check that we fully filled the template from the patches.
235 self.assertTrue(np.all(np.isfinite(template.image.array)))
236 # Because of the scale changes, there will be some ringing in the
237 # difference between the template and the original image; pick
238 # tolerances large enough to account for that.
239 self.assertImagesAlmostEqual(template.image, expected.image, rtol=.1, atol=4*calibration)
240 # Variance plane ==4 in the original image (realize() takes a noise
241 # sigma). Warping sets the level from the pixel areas and the two
242 # warping kernels, which `_correctVariance` corrects up to the
243 # difference between the configured coaddWarpKernel and the lanczos5
244 # `_makePatches` really used. A per-pixel ripple that no scalar
245 # correction can remove remains on top of that, so check the level
246 # tightly and allow for the ripple around it.
247 variance = template.variance.array[np.isfinite(template.variance.array)]
248 median = np.median(variance)
249 self.assertFloatsAlmostEqual(median, np.median(expected.variance.array),
250 rtol=0.35, msg="variance level differs")
251 self.assertLess(np.percentile(variance, 99)/median, 1.6, msg="variance ripple too large")
252 self.assertGreater(np.percentile(variance, 1)/median, 0.6, msg="variance ripple too large")
253 # Not checking the mask, as warping changes the sizes of the masks.
255 def testRunOneTractInput(self):
256 """Test a bounding box that fully fits inside one tract, with only
257 that tract passed as input. This checks that the code handles a single
258 tract input correctly.
259 """
260 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
261 task = lsst.ip.diffim.GetTemplateTask()
262 # Restrict to tract 0, since the box fits in just that tract.
263 # Task modifies the input bbox, so pass a copy.
264 result = task.run(coaddExposureHandles={0: self.patches[0]},
265 bbox=lsst.geom.Box2I(box),
266 wcs=self.exposure.wcs,
267 dataIds={0: self.dataIds[0]},
268 physical_filter="a_test")
270 # All 4 patches from tract 0 are included in this template.
271 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 4)
272 self._checkPixels(result.template, task.config, box)
274 def testRunOneTractMultipleInputs(self):
275 """Test a bounding box that fully fits inside one tract but where
276 multiple tracts were passed in. This checks that patches that are
277 mostly NaN after warping are merged correctly in the output.
278 """
279 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
280 task = lsst.ip.diffim.GetTemplateTask()
281 # Task modifies the input bbox, so pass a copy.
282 result = task.run(coaddExposureHandles=self.patches,
283 bbox=lsst.geom.Box2I(box),
284 wcs=self.exposure.wcs,
285 dataIds=self.dataIds,
286 physical_filter="a_test")
288 # All 4 patches from two tracts are included in this template.
289 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 6)
290 self._checkPixels(result.template, task.config, box)
292 def testRunTwoTracts(self):
293 """Test a bounding box that crosses tract boundaries.
294 """
295 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))
296 task = lsst.ip.diffim.GetTemplateTask()
297 # Task modifies the input bbox, so pass a copy.
298 result = task.run(coaddExposureHandles=self.patches,
299 bbox=lsst.geom.Box2I(box),
300 wcs=self.exposure.wcs,
301 dataIds=self.dataIds,
302 physical_filter="a_test")
304 # All 4 patches from all 4 tracts are included in this template
305 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9)
306 self._checkPixels(result.template, task.config, box)
308 def testRunCalibratesToNanojansky(self):
309 """Test that the template is in nJy whatever the calibration of the
310 coadds, and that coadds already in nJy are not rescaled.
311 """
312 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
313 calibration = self.exposure.photoCalib.getCalibrationMean()
314 self.assertNotEqual(calibration, 1.0)
315 task = lsst.ip.diffim.GetTemplateTask()
316 templates = {}
317 for photoCalib in (self.exposure.photoCalib, lsst.afw.image.PhotoCalib(1.0)):
318 self._remakePatches(photoCalib)
319 # Task modifies the input bbox, so pass a copy.
320 result = task.run(coaddExposureHandles={0: self.patches[0]},
321 bbox=lsst.geom.Box2I(box),
322 wcs=self.exposure.wcs,
323 dataIds={0: self.dataIds[0]},
324 physical_filter="a_test")
325 self.assertEqual(result.template.photoCalib, lsst.afw.image.PhotoCalib(1.0))
326 templates[photoCalib.getCalibrationMean()] = result.template
328 self.assertFloatsAlmostEqual(templates[calibration].image.array,
329 calibration*templates[1.0].image.array, rtol=1e-6)
330 self.assertFloatsAlmostEqual(templates[calibration].variance.array,
331 calibration**2*templates[1.0].variance.array, rtol=1e-6)
333 def _remakePatches(self, photoCalib):
334 """Rebuild the coadd patches from the main exposure with a different
335 photometric calibration.
336 """
337 self.exposure.setPhotoCalib(photoCalib)
338 self.patches.clear()
339 self.dataIds.clear()
340 for tract_id in range(4):
341 self._makePatches(self.skymap.generateTract(tract_id))
343 def testGaussianPsfGridFallback(self):
344 """If the CoaddPsf cannot be evaluated at its average position, the
345 Gaussian width is averaged over a grid of positions instead.
346 """
347 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
348 task = lsst.ip.diffim.GetTemplateTask()
349 result = task.run(coaddExposureHandles={0: self.patches[0]},
350 bbox=lsst.geom.Box2I(box),
351 wcs=self.exposure.wcs,
352 dataIds={0: self.dataIds[0]},
353 physical_filter="a_test")
354 template = result.template
355 template.setPsf(self._makeCoaddPsf(template, task.config))
356 expected = lsst.ip.diffim.utils.evaluateMeanPsfFwhm(
357 template,
358 fwhmExposureBuffer=task.config.fwhmExposureBuffer,
359 fwhmExposureGrid=task.config.fwhmExposureGrid,
360 )/gaussian_sigma_to_fwhm
362 error = lsst.pex.exceptions.InvalidParameterError("No inputs at the average position.")
363 with unittest.mock.patch("lsst.ip.diffim.getTemplate.getPsfFwhm", side_effect=error):
364 with self.assertLogs(task.log.name, level="INFO") as cm:
365 psf = task._makeGaussianPsf(template)
367 self.assertIn("grid of points", "\n".join(cm.output))
368 self.assertIsInstance(psf, lsst.afw.detection.GaussianPsf)
369 self.assertFloatsAlmostEqual(psf.getSigma(), expected)
370 self.assertFloatsAlmostEqual(task.metadata["templatePsfSigma"], expected)
372 def testRunNoTemplate(self):
373 """A bounding box that doesn't overlap the patches will raise.
374 """
375 box = lsst.geom.Box2I(lsst.geom.Point2I(1200, 1200), lsst.geom.Point2I(1600, 1600))
376 task = lsst.ip.diffim.GetTemplateTask()
377 with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"):
378 task.run(coaddExposureHandles=self.patches,
379 bbox=lsst.geom.Box2I(box),
380 wcs=self.exposure.wcs,
381 dataIds=self.dataIds,
382 physical_filter="a_test")
384 def testMissingPatches(self):
385 """Test that a missing patch results in an appropriate mask.
387 This fixes the bug reported on DM-44997 (image and variance were NaN
388 but the mask was not set to NO_DATA for those pixels).
389 """
390 # tract=0, patch=1 is the lower-left corner, as displayed in DS9.
391 self.patches[0].pop(1)
392 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
393 task = lsst.ip.diffim.GetTemplateTask()
394 # Task modifies the input bbox, so pass a copy.
395 result = task.run(coaddExposureHandles=self.patches,
396 bbox=lsst.geom.Box2I(box),
397 wcs=self.exposure.wcs,
398 dataIds=self.dataIds,
399 physical_filter="a_test")
400 no_data = (result.template.mask.array & result.template.mask.getPlaneBitMask("NO_DATA")) != 0
401 self.assertTrue(np.isfinite(result.template.image.array).all())
402 self.assertTrue(np.isfinite(result.template.variance.array).all())
403 self.assertEqual(no_data.sum(), 20990)
405 @lsst.utils.tests.methodParameters(
406 box=[
407 lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)),
408 lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)),
409 ],
410 nInput=[8, 16],
411 )
412 def testNanInputs(self, box=None, nInput=None):
413 """Test that the template has finite values when some of the input
414 pixels have NaN as variance.
415 """
416 for tract, patchRefs in self.patches.items():
417 for patchRef in patchRefs:
418 patchCoadd = patchRef.get()
419 bbox = lsst.geom.Box2I()
420 bbox.include(lsst.geom.Point2I(patchCoadd.getBBox().getCenter()))
421 bbox.grow(3)
422 patchCoadd.variance[bbox].array *= np.nan
424 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))
425 task = lsst.ip.diffim.GetTemplateTask()
426 result = task.run(coaddExposureHandles=self.patches,
427 bbox=lsst.geom.Box2I(box),
428 wcs=self.exposure.wcs,
429 dataIds=self.dataIds,
430 physical_filter="a_test")
431 if debug: 431 ↛ 432line 431 didn't jump to line 432 because the condition on line 431 was never true
432 _showTemplate(box, result.template)
433 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9)
434 # We just check that the pixel values are all finite. We cannot check that pixel values
435 # in the template are closer to the original anymore.
436 self.assertTrue(np.isfinite(result.template.image.array).all())
438 def _runCorrection(self, doCorrectVariancePlateScale=False, doScaleVariance=False,
439 handles=None):
440 """Run the task on tract 0 with the named variance corrections.
442 Everything defaults to off, so each test turns on exactly what it is
443 exercising.
444 """
445 config = lsst.ip.diffim.GetTemplateTask.ConfigClass()
446 config.doCorrectVariancePlateScale = doCorrectVariancePlateScale
447 config.doScaleVariance = doScaleVariance
448 task = lsst.ip.diffim.GetTemplateTask(config=config)
449 result = task.run(coaddExposureHandles={0: handles or self.patches[0]},
450 bbox=lsst.geom.Box2I(self.varianceBox),
451 wcs=self.exposure.wcs,
452 dataIds={0: self.dataIds[0]},
453 physical_filter="a_test")
454 return task, result.template
456 def _makeScaledWcs(self, factor):
457 """Make a WCS like the base exposure's, but with its pixel scale
458 multiplied by ``factor``.
459 """
460 cdMatrix = lsst.afw.geom.makeCdMatrix(factor*1.05*self.scale*lsst.geom.arcseconds,
461 93*lsst.geom.degrees)
462 return lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150),
463 lsst.geom.SpherePoint(0, 0, lsst.geom.radians),
464 cdMatrix)
466 @staticmethod
467 def _makeCoaddInputs(records):
468 """Make a CoaddInputs holding the given input records.
470 Parameters
471 ----------
472 records : `list` [`tuple` [`lsst.afw.geom.SkyWcs` or `None`, \
473 `lsst.geom.Box2I`]]
474 The WCS and bbox to record for each input. A `None` WCS makes a
475 record that `_plateScaleFactor` has to skip. May be empty, to
476 simulate a coadd whose plate scale cannot be reconstructed.
477 """
478 ccdSchema = lsst.afw.table.ExposureTable.makeMinimalSchema()
479 weightKey = ccdSchema.addField("weight", type=float, doc="Coadd weight")
480 coaddInputs = lsst.afw.image.CoaddInputs(
481 lsst.afw.table.ExposureTable.makeMinimalSchema(), ccdSchema)
482 for wcs, bbox in records:
483 record = coaddInputs.ccds.addNew()
484 record.setWcs(wcs)
485 record.setBBox(bbox)
486 # Included because real coadds have it, though a single-record
487 # correction does not use it.
488 record.set(weightKey, 1.0)
489 return coaddInputs
491 def _patchHandles(self, tract, records):
492 """Return handles for a tract's patches, with their CoaddInputs
493 replaced by ``records``.
494 """
495 handles = []
496 for ref in self.patches[tract]:
497 coadd = ref.get()
498 coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records))
499 handles.append(pipeBase.InMemoryDatasetHandle(
500 coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId))
501 return handles
503 def testCorrectVariancePlateScale(self):
504 """The plate scale correction is the total pixel area change from the
505 images the coadds were built from to the science image.
506 """
507 _, off = self._runCorrection()
509 # The fixture's records are the science image itself, so there is no
510 # net change in pixel area: the same-instrument case.
511 task, on = self._runCorrection(doCorrectVariancePlateScale=True)
512 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.0, rtol=1e-6)
514 # Coarser original pixels than science pixels, the DECam-template
515 # case: the correction is the ratio of their areas.
516 for scaleFactor in (1.315, 0.5):
517 with self.subTest(scaleFactor=scaleFactor):
518 handles = self._patchHandles(
519 0, [(self._makeScaledWcs(scaleFactor), self.exposure.getBBox())])
520 task, on = self._runCorrection(doCorrectVariancePlateScale=True,
521 handles=handles)
522 factor = task.metadata["variancePlateScaleFactor"]
523 self.assertFloatsAlmostEqual(factor, scaleFactor**2, rtol=1e-6)
524 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor,
525 rtol=1e-5, ignoreNaNs=True)
527 def testCorrectVariancePlateScaleUsesOneRecord(self):
528 """Only the first usable coadd input is read.
530 The spread between records is just the local pixel scale, a few
531 tenths of a percent across a real focal plane, so one stands for all
532 of them.
533 """
534 handles = self._patchHandles(
535 0, [(None, self.exposure.getBBox()),
536 (self._makeScaledWcs(2.0), self.exposure.getBBox()),
537 (self.exposure.wcs, self.exposure.getBBox())])
538 task, _ = self._runCorrection(doCorrectVariancePlateScale=True, handles=handles)
540 # The first record with a WCS: the one without is skipped, and the
541 # last (which would give 1.0) is never reached.
542 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 4.0, rtol=1e-6)
544 def testCorrectVariancePlateScaleNeedsCoaddInputs(self):
545 """Without usable coadd inputs the plate scale cannot be
546 reconstructed, and the task must say so rather than silently applying
547 only part of the correction.
548 """
549 handles = self._patchHandles(0, [])
550 with self.assertRaisesRegex(RuntimeError, "doCorrectVariancePlateScale"):
551 self._runCorrection(doCorrectVariancePlateScale=True, handles=handles)
553 def _scaleInputVariance(self, tract, factor):
554 """Return fresh handles for one tract's patches, with their variance
555 planes multiplied by ``factor``.
557 Parameters
558 ----------
559 tract : `int`
560 Id of the tract whose patches should be copied.
561 factor : `float`
562 Factor to multiply the input variance planes by.
564 Returns
565 -------
566 handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`]
567 Handles to the modified patches.
568 """
569 handles = []
570 for handle in self.patches[tract]:
571 # ``copy=True`` on the original handles means this is a copy, so
572 # the patches shared with the other tests are left untouched.
573 patch = handle.get()
574 patch.variance.array *= factor
575 handles.append(pipeBase.InMemoryDatasetHandle(patch,
576 storageClass="ExposureF",
577 copy=True,
578 dataId=handle.dataId))
579 return handles
581 def testScaleVariance(self):
582 """Test that the template variance plane is rescaled to match the
583 empirical pixel noise, and that the factor used is recorded in the
584 task metadata.
585 """
586 scaleFactor = 1.345
587 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
589 def _configureAndRunTask(doScaleVariance, varianceScale=1.):
590 """Build a template from tract 0, optionally rescaling the input
591 variance planes by ``varianceScale`` first.
592 """
593 config = lsst.ip.diffim.GetTemplateTask.ConfigClass()
594 config.doScaleVariance = doScaleVariance
595 task = lsst.ip.diffim.GetTemplateTask(config=config)
596 # Task modifies the input bbox, so pass a copy.
597 result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)},
598 bbox=lsst.geom.Box2I(box),
599 wcs=self.exposure.wcs,
600 dataIds={0: self.dataIds[0]},
601 physical_filter="a_test")
602 return task, result.template
604 # With scaling disabled the subtask is never constructed, and nothing
605 # is recorded in the metadata.
606 taskOff, templateOff = _configureAndRunTask(False)
607 self.assertFalse(hasattr(taskOff, "scaleVariance"))
608 self.assertNotIn("scaleTemplateVarianceFactor", taskOff.metadata)
610 # Both warps -- lanczos5 in ``_makePatches`` and lanczos3 in the
611 # task -- correlate the noise. The variance plane tracks only the
612 # per-pixel diagonal, which the second warp leaves too low, so
613 # ``scaleVariance`` measures a factor well above 1 even though the
614 # input variance planes are correct.
615 #
616 taskOn, templateOn = _configureAndRunTask(True)
617 factor = taskOn.metadata["scaleTemplateVarianceFactor"]
618 # TODO DM-55879: this value is pinned on purpose. The lanczos warping
619 # kernels introduce small correlations that artificially suppress the
620 # image pixel stddev and inflate the variance scaling factor. This
621 # should be changed to 1.0 after DM-55879 is merged.
622 self.assertFloatsAlmostEqual(factor, 1.1465, atol=0.01,
623 msg="Measured template variance scaling changed; see the"
624 " comment above if the correlation correction landed.")
625 # The only difference from the unscaled template is the constant
626 # factor applied to the variance plane.
627 self.assertFloatsAlmostEqual(templateOn.variance.array,
628 templateOff.variance.array*factor, rtol=1e-5)
629 # Tolerance here is float32 round-off: repeated runs of the task are
630 # not bitwise identical.
631 self.assertImagesAlmostEqual(templateOn.image, templateOff.image, rtol=1e-5, atol=1e-5)
633 # If the input variance planes under-estimate the noise by a known
634 # factor, the measured factor grows by that amount and the same
635 # output variance plane is recovered.
636 taskLow, templateLow = _configureAndRunTask(True, varianceScale=1/scaleFactor)
637 self.assertFloatsAlmostEqual(taskLow.metadata["scaleTemplateVarianceFactor"],
638 factor*scaleFactor, rtol=1e-5)
639 self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5)
642def setup_module(module):
643 lsst.utils.tests.init()
646class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
647 pass
650if __name__ == "__main__": 650 ↛ 651line 650 didn't jump to line 651 because the condition on line 650 was never true
651 lsst.utils.tests.init()
652 unittest.main()