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