Coverage for tests/test_getTemplate.py: 96%
684 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-09 17:53 +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
27import astropy.units as u
28from astropy.stats import gaussian_sigma_to_fwhm
29import astshim
30import numpy as np
32import lsst.afw.cameraGeom.testUtils
33import lsst.afw.detection
34import lsst.afw.geom
35import lsst.afw.image
36import lsst.afw.math
37import lsst.afw.table
38from lsst.daf.base import DateTime
39from lsst.daf.butler import DataCoordinate, DatasetRef, DatasetType, DimensionUniverse
40import lsst.geom
41import lsst.images
42import lsst.images.psfs
43from lsst.images.tests import compare_masked_image_to_legacy
44import lsst.ip.diffim
45import lsst.meas.algorithms
46import lsst.meas.base.tests
47import lsst.pex.exceptions
48import lsst.pipe.base as pipeBase
49import lsst.pipe.base.testUtils
50import lsst.skymap
51import lsst.utils.tests
53from utils import generate_data_id, makeTestExposureRecord
55# Change this to True, `setup display_ds9`, and open ds9 (or use another afw
56# display backend) to show the tract/patch layouts on the image.
57debug = False
58if debug: 58 ↛ 59line 58 didn't jump to line 59 because the condition on line 58 was never true
59 import lsst.afw.display
60 display = lsst.afw.display.Display()
61 display.frame = 1
63# Sum of squared warping kernel values along one axis, with the kernel
64# normalized to unit sum at each sub-pixel phase, averaged over phase.
65KERNEL_SUM_SQ_1D = {
66 "bilinear": 2/3,
67 "lanczos3": 0.8933699084,
68 "lanczos4": 0.9141249535,
69 "lanczos5": 0.9342245513,
70}
73def _showTemplate(box, template):
74 """Show the corners of the template we made in this test."""
75 for point in box.getCorners():
76 display.dot("+", point.x, point.y, ctype="orange", size=40)
77 display.frame = 2
78 display.image(template, "warped template")
79 display.frame = 3
80 display.image(template.variance, "warped variance")
83class GetTemplateTaskTestCase(lsst.utils.tests.TestCase):
84 """Test that GetTemplateTask works on both one tract and multiple tract
85 input coadd exposures.
87 Makes a synthetic exposure large enough to fit four small tracts with 2x2
88 (300x300 pixel) patches each, extracts pixels for those patches by warping,
89 and tests GetTemplateTask's output against boxes that overlap various
90 combinations of one or multiple tracts.
91 """
92 scale = 0.2 # arcsec/pixel
93 # Use 5th order to minimize artifacts on the templates.
94 coaddWarpKernel = "lanczos5"
96 @classmethod
97 def setUpClass(cls):
98 # Warping the patches is the slowest part of every test, so it is done
99 # once. Tests must not modify the exposure; the patch handles return
100 # copies, and setUp gives each test its own lists of them.
101 cls.skymap = cls._makeSkymap()
102 cls.patches = collections.defaultdict(list)
103 cls.dataIds = collections.defaultdict(list)
104 cls.exposure = cls._makeExposure()
106 if debug: 106 ↛ 107line 106 didn't jump to line 107 because the condition on line 106 was never true
107 display.image(cls.exposure, "base exposure")
109 for tract_id in range(4):
110 tract = cls.skymap.generateTract(tract_id)
111 cls._makePatches(tract)
113 @classmethod
114 def tearDownClass(cls):
115 del cls.skymap, cls.patches, cls.dataIds, cls.exposure
117 def setUp(self):
118 cls = type(self)
119 self.patches = collections.defaultdict(list, {k: list(v) for k, v in cls.patches.items()})
120 self.dataIds = collections.defaultdict(list, {k: list(v) for k, v in cls.dataIds.items()})
121 self.coaddRefs = {}
122 self.varianceBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
124 @classmethod
125 def _makeSkymap(cls):
126 """Make a Skymap with 4 tracts with 4 patches each.
127 """
128 tractScale = 0.02 # degrees
129 # On-sky coordinates of the tract centers.
130 coords = [(0, 0),
131 (0, tractScale),
132 (tractScale, 0),
133 (tractScale, tractScale),
134 ]
135 config = lsst.skymap.DiscreteSkyMap.ConfigClass()
136 config.raList = [c[0] for c in coords]
137 config.decList = [c[1] for c in coords]
138 # Half the tract center step size, to keep the tract overlap small.
139 config.radiusList = [tractScale/2 for c in coords]
140 config.projection = "TAN"
141 config.pixelScale = cls.scale
142 config.tractOverlap = 0.0005
143 config.tractBuilder = "legacy"
144 config.tractBuilder["legacy"].patchInnerDimensions = (300, 300)
145 config.tractBuilder["legacy"].patchBorder = 10
146 return lsst.skymap.DiscreteSkyMap(config=config)
148 @classmethod
149 def _makeExposure(cls):
150 """Create a large image to break up into tracts and patches.
152 The image will have a source every 100 pixels in x and y, and a WCS
153 that results in the tracts all fitting in the image, with tract=0
154 in the lower left, tract=1 to the right, tract=2 above, and tract=3
155 to the upper right.
156 """
157 box = lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(800, 800))
158 # This WCS was constructed so that tract 0 mostly fills the lower left
159 # quadrant of the image, and the other tracts fill the rest; slight
160 # extra rotation as a check on the final warp layout, scaled by 5%
161 # from the patch pixel scale.
162 cd_matrix = lsst.afw.geom.makeCdMatrix(1.05*cls.scale*lsst.geom.arcseconds, 93*lsst.geom.degrees)
163 wcs = lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150),
164 lsst.geom.SpherePoint(0, 0, lsst.geom.radians),
165 cd_matrix)
166 dataset = lsst.meas.base.tests.TestDataset(box, wcs=wcs)
167 for x, y in itertools.product(np.arange(0, 500, 100), np.arange(0, 500, 100)):
168 dataset.addSource(1e5, lsst.geom.Point2D(x, y))
169 exposure, _ = dataset.realize(2, dataset.makeMinimalSchema())
170 exposure.setFilter(lsst.afw.image.FilterLabel("a", "a_test"))
171 return exposure
173 @classmethod
174 def _makePatches(cls, tract):
175 """Populate the patches and dataId dicts, keyed on tract id, with the
176 warps of the main exposure and minimal dataIds, respectively.
177 """
178 if debug: 178 ↛ 179line 178 didn't jump to line 179 because the condition on line 178 was never true
179 color = ['red', 'green', 'cyan', 'yellow'][tract.tract_id]
180 point = cls.exposure.wcs.skyToPixel(tract.ctr_coord)
181 # Show the tract center, colored by tract id.
182 display.dot("x", point.x, point.y, ctype=color, size=30)
184 config = lsst.afw.math.Warper.ConfigClass()
185 config.warpingKernelName = cls.coaddWarpKernel
186 warper = lsst.afw.math.Warper.fromConfig(config)
187 for patchId in range(tract.num_patches.x*tract.num_patches.y):
188 patch = tract.getPatchInfo(patchId)
189 box = patch.getOuterBBox()
191 if debug: 191 ↛ 193line 191 didn't jump to line 193 because the condition on line 191 was never true
192 # Show the patch corners as patch ids, colored by tract id.
193 points = cls.exposure.wcs.skyToPixel(patch.wcs.pixelToSky([lsst.geom.Point2D(x)
194 for x in box.getCorners()]))
195 for p in points:
196 display.dot(patchId, p.x, p.y, ctype=color)
198 # This is mostly taken from drp_tasks makePsfMatchedWarp, but
199 # ip_diffim cannot depend on drp_tasks.
200 xyTransform = lsst.afw.geom.makeWcsPairTransform(cls.exposure.wcs, patch.wcs)
201 warpedPsf = lsst.meas.algorithms.WarpedPsf(cls.exposure.psf, xyTransform)
202 warped = warper.warpExposure(patch.wcs, cls.exposure, destBBox=box)
203 warped.setPsf(warpedPsf)
205 warped.getInfo().setCoaddInputs(
206 cls._makeCoaddInputs([(cls.exposure.wcs, cls.exposure.getBBox())]))
207 dataRef = pipeBase.InMemoryDatasetHandle(
208 warped,
209 storageClass="ExposureF",
210 copy=True,
211 dataId=generate_data_id(
212 tract=tract,
213 patch=patch,
214 )
215 )
216 cls.patches[tract.tract_id].append(dataRef)
217 dataCoordinate = DataCoordinate.standardize({"tract": tract.tract_id,
218 "patch": patchId,
219 "band": "a",
220 "skymap": "skymap"},
221 universe=DimensionUniverse())
222 cls.dataIds[tract.tract_id].append(dataCoordinate)
224 def _checkMetadata(self, template, config, box, wcs, nPsfs):
225 """Check that the various metadata components were set correctly.
226 """
227 expectedBox = lsst.geom.Box2I(box)
228 expectedBox.grow(config.templateBorderSize)
229 self.assertEqual(template.getBBox(), expectedBox)
230 # WCS should match our exposure, not any of the coadd tracts.
231 for tract in self.patches:
232 self.assertNotEqual(template.wcs, self.patches[tract][0].get().wcs)
233 self.assertEqual(template.wcs, self.exposure.wcs)
234 # The template pixels are calibrated to nJy.
235 self.assertEqual(template.photoCalib, lsst.afw.image.PhotoCalib(1.0))
236 self.assertEqual(template.getXY0(), expectedBox.getMin())
237 self.assertEqual(template.filter.bandLabel, "a")
238 self.assertEqual(template.filter.physicalLabel, "a_test")
239 # The template PSF is a Gaussian with the width of the CoaddPsf of
240 # its inputs.
241 self.assertIsInstance(template.psf, lsst.afw.detection.GaussianPsf)
242 coaddPsf = self._makeCoaddPsf(template, config)
243 self.assertEqual(coaddPsf.getComponentCount(), nPsfs)
244 position = coaddPsf.getAveragePosition()
245 self.assertFloatsAlmostEqual(template.psf.getSigma(),
246 coaddPsf.computeShape(position).getTraceRadius())
247 self.assertTrue(template.getInfo().hasCoaddInputs())
248 self.assertEqual(len(template.getInfo().getCoaddInputs().ccds), nPsfs)
250 def _makeConfig(self):
251 """Make a GetTemplateTask config whose coaddWarpKernel is the kernel
252 the test patches were made with.
253 """
254 config = lsst.ip.diffim.GetTemplateTask.ConfigClass()
255 config.coaddWarpKernel = self.coaddWarpKernel
256 return config
258 def _makeCoaddPsf(self, template, config):
259 """Return the CoaddPsf that ``run`` approximated with a Gaussian,
260 rebuilt from the coadd inputs recorded on the template.
261 """
262 task = lsst.ip.diffim.GetTemplateTask(config=config)
263 return task._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.wcs)
265 def _checkPixels(self, template, config, box):
266 """Check that the pixel values in the template are close to the
267 original image, calibrated to nJy.
268 """
269 # All pixels should have real values!
270 expectedBox = lsst.geom.Box2I(box)
271 expectedBox.grow(config.templateBorderSize)
272 calibration = self.exposure.photoCalib.getCalibrationMean()
273 expected = self.exposure.photoCalib.calibrateImage(self.exposure[expectedBox].maskedImage)
275 if debug: 275 ↛ 276line 275 didn't jump to line 276 because the condition on line 275 was never true
276 _showTemplate(expectedBox, template)
278 # Check that we fully filled the template from the patches.
279 self.assertTrue(np.all(np.isfinite(template.image.array)))
280 # Because of the scale changes, there will be some ringing in the
281 # difference between the template and the original image; pick
282 # tolerances large enough to account for that.
283 self.assertImagesAlmostEqual(template.image, expected.image, rtol=.1, atol=4*calibration)
284 # Variance plane ==4 in the original image (realize() takes a noise
285 # sigma). Warping sets the level from the pixel areas and the two
286 # warping kernels, which `_correctVariance` corrects. A per-pixel
287 # ripple that no scalar correction can remove remains on top of that,
288 # so check the level tightly and allow for the ripple around it.
289 variance = template.variance.array[np.isfinite(template.variance.array)]
290 median = np.median(variance)
291 self.assertFloatsAlmostEqual(median, np.median(expected.variance.array),
292 rtol=0.02, msg="variance level differs")
293 self.assertLess(np.percentile(variance, 99)/median, 1.6, msg="variance ripple too large")
294 self.assertGreater(np.percentile(variance, 1)/median, 0.6, msg="variance ripple too large")
295 # Not checking the mask, as warping changes the sizes of the masks.
297 def testRunOneTractInput(self):
298 """Test a bounding box that fully fits inside one tract, with only
299 that tract passed as input. This checks that the code handles a single
300 tract input correctly.
301 """
302 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
303 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
304 # Restrict to tract 0, since the box fits in just that tract.
305 result = task.run(coaddExposureHandles={0: self.patches[0]},
306 bbox=box,
307 wcs=self.exposure.wcs,
308 dataIds={0: self.dataIds[0]},
309 physical_filter="a_test")
311 # All 4 patches from tract 0 are included in this template.
312 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 4)
313 self._checkPixels(result.template, task.config, box)
315 def testRunOneTractMultipleInputs(self):
316 """Test a bounding box that fully fits inside one tract but where
317 multiple tracts were passed in. This checks that patches that are
318 mostly NaN after warping are merged correctly in the output.
319 """
320 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
321 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
322 result = task.run(coaddExposureHandles=self.patches,
323 bbox=box,
324 wcs=self.exposure.wcs,
325 dataIds=self.dataIds,
326 physical_filter="a_test")
328 # All 4 patches from two tracts are included in this template.
329 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 6)
330 self._checkPixels(result.template, task.config, box)
332 def testRunTwoTracts(self):
333 """Test a bounding box that crosses tract boundaries.
334 """
335 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))
336 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
337 result = task.run(coaddExposureHandles=self.patches,
338 bbox=box,
339 wcs=self.exposure.wcs,
340 dataIds=self.dataIds,
341 physical_filter="a_test")
343 # All 4 patches from all 4 tracts are included in this template
344 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9)
345 self._checkPixels(result.template, task.config, box)
347 def testRunCalibratesToNanojansky(self):
348 """Test that the template is in nJy whatever the calibration of the
349 coadds, and that coadds already in nJy are not rescaled.
350 """
351 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
352 calibration = self.exposure.photoCalib.getCalibrationMean()
353 self.assertNotEqual(calibration, 1.0)
354 # The same coadd pixels, labelled as already in nJy.
355 uncalibrated = []
356 for handle in self.patches[0]:
357 coadd = handle.get()
358 coadd.setPhotoCalib(lsst.afw.image.PhotoCalib(1.0))
359 uncalibrated.append(pipeBase.InMemoryDatasetHandle(coadd, storageClass="ExposureF", copy=True,
360 dataId=handle.dataId))
361 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
362 templates = {}
363 for key, handles in ((calibration, self.patches[0]), (1.0, uncalibrated)):
364 result = task.run(coaddExposureHandles={0: handles},
365 bbox=box,
366 wcs=self.exposure.wcs,
367 dataIds={0: self.dataIds[0]},
368 physical_filter="a_test")
369 self.assertEqual(result.template.photoCalib, lsst.afw.image.PhotoCalib(1.0))
370 templates[key] = result.template
372 self.assertFloatsAlmostEqual(templates[calibration].image.array,
373 calibration*templates[1.0].image.array, rtol=1e-6)
374 self.assertFloatsAlmostEqual(templates[calibration].variance.array,
375 calibration**2*templates[1.0].variance.array, rtol=1e-6)
377 def testGaussianPsfGridFallback(self):
378 """If the CoaddPsf cannot be evaluated at its average position, the
379 Gaussian width is averaged over a grid of positions instead.
380 """
381 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
382 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
383 result = task.run(coaddExposureHandles={0: self.patches[0]},
384 bbox=box,
385 wcs=self.exposure.wcs,
386 dataIds={0: self.dataIds[0]},
387 physical_filter="a_test")
388 template = result.template
389 template.setPsf(self._makeCoaddPsf(template, task.config))
390 expected = lsst.ip.diffim.utils.evaluateMeanPsfFwhm(
391 template,
392 fwhmExposureBuffer=task.config.fwhmExposureBuffer,
393 fwhmExposureGrid=task.config.fwhmExposureGrid,
394 )/gaussian_sigma_to_fwhm
396 error = lsst.pex.exceptions.InvalidParameterError("No inputs at the average position.")
397 with unittest.mock.patch("lsst.ip.diffim.getTemplate.getPsfFwhm", side_effect=error):
398 with self.assertLogs(task.log.name, level="INFO") as cm:
399 psf = task._makeGaussianPsf(template)
401 self.assertIn("grid of points", "\n".join(cm.output))
402 self.assertIsInstance(psf, lsst.afw.detection.GaussianPsf)
403 self.assertFloatsAlmostEqual(psf.getSigma(), expected)
404 self.assertFloatsAlmostEqual(task.metadata["templatePsfSigma"], expected)
406 def testRunNoTemplate(self):
407 """A bounding box that doesn't overlap the patches will raise.
408 """
409 box = lsst.geom.Box2I(lsst.geom.Point2I(1200, 1200), lsst.geom.Point2I(1600, 1600))
410 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
411 with self.assertRaisesRegex(lsst.pipe.base.NoWorkFound, "No patches found"):
412 task.run(coaddExposureHandles=self.patches,
413 bbox=box,
414 wcs=self.exposure.wcs,
415 dataIds=self.dataIds,
416 physical_filter="a_test")
418 def testMissingPatches(self):
419 """Test that a missing patch results in an appropriate mask.
421 This fixes the bug reported on DM-44997 (image and variance were NaN
422 but the mask was not set to NO_DATA for those pixels).
423 """
424 # tract=0, patch=1 is the lower-left corner, as displayed in DS9.
425 self.patches[0].pop(1)
426 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
427 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
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 no_data = (result.template.mask.array & result.template.mask.getPlaneBitMask("NO_DATA")) != 0
434 self.assertTrue(np.isfinite(result.template.image.array).all())
435 self.assertTrue(np.isfinite(result.template.variance.array).all())
436 self.assertEqual(no_data.sum(), 20990)
438 @lsst.utils.tests.methodParameters(
439 box=[
440 lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)),
441 lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)),
442 ],
443 nInput=[8, 16],
444 )
445 def testNanInputs(self, box=None, nInput=None):
446 """Test that the template has finite values when some of the input
447 pixels have NaN as variance.
448 """
449 for tract, patchRefs in self.patches.items():
450 for patchRef in patchRefs:
451 patchCoadd = patchRef.get()
452 bbox = lsst.geom.Box2I()
453 bbox.include(lsst.geom.Point2I(patchCoadd.getBBox().getCenter()))
454 bbox.grow(3)
455 patchCoadd.variance[bbox].array *= np.nan
457 box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))
458 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
459 result = task.run(coaddExposureHandles=self.patches,
460 bbox=box,
461 wcs=self.exposure.wcs,
462 dataIds=self.dataIds,
463 physical_filter="a_test")
464 if debug: 464 ↛ 465line 464 didn't jump to line 465 because the condition on line 464 was never true
465 _showTemplate(box, result.template)
466 self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 9)
467 # We just check that the pixel values are all finite. We cannot check that pixel values
468 # in the template are closer to the original anymore.
469 self.assertTrue(np.isfinite(result.template.image.array).all())
471 def _runCorrection(self, doCorrectVarianceWarpCorrelation=False, doScaleVariance=False,
472 handles=None, instrument="SciCam", coaddWarpKernel=None,
473 templateWarpKernel="lanczos3"):
474 """Run the task on tract 0, with a science image from ``instrument``.
476 The plate scale correction is applied only when ``instrument`` is not
477 the instrument of the template coadds, which is ``"TestCam"`` unless
478 ``handles`` says otherwise. ``coaddWarpKernel`` defaults to the kernel
479 the patches were made with, so the reconstructed factors are exact.
480 """
481 config = self._makeConfig()
482 config.doCorrectVarianceWarpCorrelation = doCorrectVarianceWarpCorrelation
483 config.doScaleVariance = doScaleVariance
484 if coaddWarpKernel is not None:
485 config.coaddWarpKernel = coaddWarpKernel
486 config.warp.warpingKernelName = templateWarpKernel
487 task = lsst.ip.diffim.GetTemplateTask(config=config)
488 result = task.run(coaddExposureHandles={0: handles or self.patches[0]},
489 bbox=lsst.geom.Box2I(self.varianceBox),
490 wcs=self.exposure.wcs,
491 dataIds={0: self.dataIds[0]},
492 physical_filter="a_test",
493 instrument=instrument)
494 return task, result.template
496 def _makeScaledWcs(self, factor):
497 """Make a WCS like the base exposure's, but with its pixel scale
498 multiplied by ``factor``.
499 """
500 cdMatrix = lsst.afw.geom.makeCdMatrix(factor*1.05*self.scale*lsst.geom.arcseconds,
501 93*lsst.geom.degrees)
502 return lsst.afw.geom.makeSkyWcs(lsst.geom.Point2D(120, 150),
503 lsst.geom.SpherePoint(0, 0, lsst.geom.radians),
504 cdMatrix)
506 @staticmethod
507 def _makeCoaddInputs(records, instrument="TestCam", epochs=None):
508 """Make a CoaddInputs holding the given input records.
510 Parameters
511 ----------
512 records : `list` [`tuple` [`lsst.afw.geom.SkyWcs` or `None`, \
513 `lsst.geom.Box2I`]]
514 The WCS and bbox to record for each input. A `None` WCS makes a
515 record that `_plateScaleFactor` has to skip. May be empty, to
516 simulate a coadd whose plate scale cannot be reconstructed.
517 instrument : `str` or `None`, optional
518 Instrument label to record for every input. `None` leaves the
519 inputs without a visit info, and an empty string gives them a
520 visit info with no label; either way the instrument is unknown.
521 epochs : `list` [`tuple` [`int`, `float`, `float`]], optional
522 Visit, mid-exposure MJD (TAI) and coadd weight of each input. By
523 default every input is visit 0 with weight 1 and no observation
524 time. The time is not recorded if ``instrument`` is `None`.
525 """
526 ccdSchema = lsst.afw.table.ExposureTable.makeMinimalSchema()
527 weightKey = ccdSchema.addField("weight", type=float, doc="Coadd weight")
528 visitKey = ccdSchema.addField("visit", type=np.int64, doc="Visit of the input")
529 ccdSchema.addField("ccd", type=np.int32, doc="Detector of the input")
530 coaddInputs = lsst.afw.image.CoaddInputs(
531 lsst.afw.table.ExposureTable.makeMinimalSchema(), ccdSchema)
532 if epochs is None:
533 epochs = [(0, None, 1.0)]*len(records)
534 for (wcs, bbox), (visit, mjd, weight) in zip(records, epochs, strict=True):
535 record = coaddInputs.ccds.addNew()
536 record.setWcs(wcs)
537 record.setBBox(bbox)
538 record.set(visitKey, visit)
539 record.set(weightKey, weight)
540 if instrument is not None:
541 date = DateTime(mjd, DateTime.MJD, DateTime.TAI) if mjd is not None else DateTime()
542 record.setVisitInfo(lsst.afw.image.VisitInfo(instrumentLabel=instrument, date=date))
543 return coaddInputs
545 def _patchHandles(self, tract, records, instrument="TestCam", epochs=None):
546 """Return handles for a tract's patches, with their CoaddInputs
547 replaced by ``records`` from ``instrument``, observed at ``epochs``.
548 """
549 handles = []
550 for ref in self.patches[tract]:
551 coadd = ref.get()
552 coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records, instrument=instrument,
553 epochs=epochs))
554 handles.append(pipeBase.InMemoryDatasetHandle(
555 coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId))
556 return handles
558 def testCorrectVarianceWarpCorrelation(self):
559 """The correlation correction covers both warps: the one that built
560 the coadds and the one this task performs.
561 """
562 _, off = self._runCorrection(instrument="TestCam")
563 for coaddKernel, templateKernel in (("lanczos3", "lanczos3"),
564 ("lanczos5", "lanczos3"),
565 ("lanczos5", "lanczos5"),
566 ("bilinear", "lanczos3")):
567 with self.subTest(coaddWarpKernel=coaddKernel, templateWarpKernel=templateKernel):
568 task, on = self._runCorrection(doCorrectVarianceWarpCorrelation=True,
569 instrument="TestCam",
570 coaddWarpKernel=coaddKernel,
571 templateWarpKernel=templateKernel)
572 factor = task.metadata["varianceWarpCorrelationFactor"]
573 # A kernel spreads power, so this always makes the plane
574 # larger.
575 self.assertGreater(factor, 1.0)
576 # The phase average converges slowly for bilinear, because
577 # its sum of squares has a kink where the phase wraps.
578 rtol = 1e-4 if "bilinear" in (coaddKernel, templateKernel) else 1e-5
579 self.assertFloatsAlmostEqual(
580 factor,
581 1/(KERNEL_SUM_SQ_1D[coaddKernel]*KERNEL_SUM_SQ_1D[templateKernel])**2,
582 rtol=rtol)
583 self.assertEqual(task.metadata["variancePlateScaleFactor"], 1.0)
584 if templateKernel == "lanczos3":
585 # A pure scalar, with no per-pixel structure of its own;
586 # only comparable when `off` used the same template warp.
587 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor,
588 rtol=1e-5, ignoreNaNs=True)
590 def testCorrectVariancePlateScale(self):
591 """The plate scale correction is the total pixel area change from the
592 images the coadds were built from to the science image.
593 """
594 # The science image and the coadds are from the same instrument, so
595 # the factor is 1 without reading any coadd inputs.
596 task, off = self._runCorrection(instrument="TestCam")
597 self.assertEqual(task.metadata["variancePlateScaleFactor"], 1.0)
598 self.assertNotIn("templateVarianceCorrectionFactor", task.metadata)
600 # The fixture's records are the science image itself, so there is no
601 # net change in pixel area even when the instruments differ.
602 task, on = self._runCorrection()
603 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.0, rtol=1e-6)
604 self.assertNotIn("varianceWarpCorrelationFactor", task.metadata)
606 # Coarser original pixels than science pixels, the DECam-template
607 # case: the correction is the ratio of their areas.
608 for scaleFactor in (1.315, 0.5):
609 with self.subTest(scaleFactor=scaleFactor):
610 handles = self._patchHandles(
611 0, [(self._makeScaledWcs(scaleFactor), self.exposure.getBBox())])
612 task, on = self._runCorrection(handles=handles)
613 factor = task.metadata["variancePlateScaleFactor"]
614 self.assertFloatsAlmostEqual(factor, scaleFactor**2, rtol=1e-6)
615 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor,
616 rtol=1e-5, ignoreNaNs=True)
618 def testVarianceCorrectionsCompose(self):
619 """The two corrections are independent and multiply together.
621 The patches are warped from the exposure and the template is warped
622 back to its pixel grid, so the two warps cancel geometrically, leaving
623 only the two kernels' correlation to undo. The substituted coadd
624 inputs claim coarser original pixels than the exposure really has, so
625 with both corrections applied the variance must come back to its
626 original value times the claimed change in pixel area.
627 """
628 plateScale = 1.315
629 handles = self._patchHandles(0, [(self._makeScaledWcs(plateScale), self.exposure.getBBox())])
630 # A science image from the coadds' instrument applies no correction
631 # at all; one from another instrument applies both.
632 _, off = self._runCorrection(handles=handles, instrument="TestCam")
633 task, on = self._runCorrection(doCorrectVarianceWarpCorrelation=True, handles=handles)
635 # Both together is exactly the product of the two factors.
636 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], plateScale**2, rtol=1e-6)
637 scale = (task.metadata["varianceWarpCorrelationFactor"]
638 * task.metadata["variancePlateScaleFactor"])
639 self.assertFloatsAlmostEqual(task.metadata["templateVarianceCorrectionFactor"],
640 scale, rtol=1e-12)
641 self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*scale,
642 rtol=1e-5, ignoreNaNs=True)
644 # Only the variance plane is touched. Nothing about the correction
645 # feeds back into how the patches are combined.
646 self.assertImagesEqual(on.image, off.image)
647 self.assertMasksEqual(on.mask, off.mask)
649 # realize() takes a noise sigma, so the original variance is 4, not 2,
650 # before calibrating to nJy.
651 calibration = self.exposure.photoCalib.getCalibrationMean()
652 original = calibration**2*np.median(self.exposure.variance.array)
653 self.assertFloatsAlmostEqual(np.median(on.variance.array), original*plateScale**2, rtol=0.02)
655 # Which is a real correction: uncorrected, the plane is too small by
656 # the two warping kernels' phase-averaged sums of squares.
657 self.assertFloatsAlmostEqual(
658 np.median(off.variance.array)/original,
659 (KERNEL_SUM_SQ_1D[self.coaddWarpKernel]*KERNEL_SUM_SQ_1D["lanczos3"])**2, rtol=0.02)
661 def testCorrectVariancePlateScaleUsesOneRecord(self):
662 """Only the first usable coadd input is read.
664 The spread between records is just the local pixel scale, a few
665 tenths of a percent across a real focal plane, so one stands for all
666 of them.
667 """
668 handles = self._patchHandles(
669 0, [(None, self.exposure.getBBox()),
670 (self._makeScaledWcs(2.0), self.exposure.getBBox()),
671 (self.exposure.wcs, self.exposure.getBBox())])
672 task, _ = self._runCorrection(handles=handles)
674 # The first record with a WCS: the one without is skipped, and the
675 # last (which would give 1.0) is never reached.
676 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 4.0, rtol=1e-6)
678 def testCorrectVariancePlateScaleSkipsDroppedCoadd(self):
679 """The plate scale comes from the first coadd that is used, not from
680 one dropped for its instrument.
681 """
682 bbox = self.exposure.getBBox()
683 handles = (self._patchHandles(0, [(self._makeScaledWcs(2.0), bbox)], instrument="OtherCam")[:1]
684 + self._patchHandles(0, [(self._makeScaledWcs(1.315), bbox)])[1:])
685 task, _ = self._runCorrection(handles=handles)
686 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6)
688 def testCorrectVariancePlateScaleNextTract(self):
689 """If the first coadd kept in a tract cannot give the plate scale, the
690 first coadd kept in the next tract is used.
691 """
692 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
693 handles = {0: self._patchHandles(0, []),
694 1: self._patchHandles(1, [(self._makeScaledWcs(1.315), self.exposure.getBBox())])}
695 task.run(coaddExposureHandles=handles,
696 bbox=lsst.geom.Box2I(self.varianceBox),
697 wcs=self.exposure.wcs,
698 dataIds={0: self.dataIds[0], 1: self.dataIds[1]},
699 physical_filter="a_test",
700 instrument="SciCam")
701 self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6)
703 def testCorrectVariancePlateScaleNeedsCoaddInputs(self):
704 """Without usable coadd inputs the plate scale cannot be
705 reconstructed, so the template cannot be built.
706 """
707 handles = self._patchHandles(0, [])
708 with self.assertRaisesRegex(RuntimeError, "cannot be corrected for the change in pixel area"):
709 self._runCorrection(handles=handles)
711 def testRunReturnsTemplateInputs(self):
712 """`run` returns the inputs of each coadd that overlap the science
713 image, with the part of the science image each one covers. Inputs
714 with no weight or no observation time are left out, as are inputs
715 that cover only the template border.
716 """
717 wcs = self.exposure.wcs
718 whole = (wcs, self.exposure.getBBox())
719 right = (wcs, lsst.geom.Box2I(lsst.geom.Point2I(90, -200), lsst.geom.Point2I(800, 800)))
720 outside = (wcs, lsst.geom.Box2I(lsst.geom.Point2I(600, 600), lsst.geom.Point2I(800, 800)))
721 border = (wcs, lsst.geom.Box2I(lsst.geom.Point2I(185, 0), lsst.geom.Point2I(195, 180)))
722 handles = (self._patchHandles(0, [whole, right, outside, whole, whole, border],
723 epochs=[(1, 60000.25, 1.0), (2, 60010.5, 3.0), (3, 60100.0, 1.0),
724 (4, 60200.0, 0.0), (5, None, 1.0), (6, 60300.0, 1.0)])[:1]
725 + self._patchHandles(0, [whole], epochs=[(1, 60000.25, 2.0)])[1:])
726 result = self._runTract0(handles, outputImageType="future")
728 templateInputs = result.templateInputs
729 self.assertEqual(set(templateInputs), {(0, patch) for patch in range(4)})
730 first = templateInputs[0, 0]
731 self.assertEqual([(i.visit, i.detector, i.weight) for i in first], [(1, 0, 1.0), (2, 0, 3.0)])
732 self.assertFloatsAlmostEqual(np.array([i.mid_mjd_tai for i in first]), np.array([60000.25, 60010.5]),
733 rtol=0, atol=1e-6)
734 # The box that _runTract0 requests.
735 scienceBox = lsst.geom.Box2D(lsst.geom.Point2D(-0.5, -0.5), lsst.geom.Point2D(180.5, 180.5))
736 self.assertFloatsAlmostEqual(first[0].bounds.area, scienceBox.getArea(), rtol=1e-6)
737 # Warping with lanczos5 left the 4.5 pixels along the input's edges
738 # out of the coadd.
739 rightArea = (scienceBox.getMaxX() - 94.0)*scienceBox.getHeight()
740 self.assertFloatsAlmostEqual(first[1].bounds.area, rightArea, rtol=1e-6)
741 self.assertFloatsAlmostEqual(first[1].bounds.x_vertices.min(), 94.0, rtol=0, atol=1e-6)
742 # Each coadd lists its own record of an input shared with another.
743 for patch in range(1, 4):
744 self.assertEqual([(i.visit, i.weight) for i in templateInputs[0, patch]], [(1, 2.0)])
746 def testRunTemplateInputsOnlyForFuture(self):
747 """`run` records no template inputs for the legacy output type.
748 """
749 records = [(self.exposure.wcs, self.exposure.getBBox())]
750 result = self._runTract0(self._patchHandles(0, records, epochs=[(1, 60000.25, 1.0)]))
751 self.assertEqual(result.templateInputs, {})
753 def testRunTemplateInputsSkipDroppedCoadd(self):
754 """A coadd dropped for its instrument has no template inputs.
755 """
756 records = [(self.exposure.wcs, self.exposure.getBBox())]
757 handles = (self._patchHandles(0, records, epochs=[(1, 60000.25, 1.0)])[:3]
758 + self._patchHandles(0, records, instrument="OtherCam",
759 epochs=[(2, 60100.0, 1.0)])[3:])
760 result = self._runTract0(handles, outputImageType="future")
761 self.assertEqual(set(result.templateInputs), {(0, 0), (0, 1), (0, 2)})
763 def testRunTemplateInputsLeaveOutWarpBorder(self):
764 """An input's footprint leaves out the border along its edges that
765 the coadd's warping kernel left out.
766 """
767 right = (self.exposure.wcs, lsst.geom.Box2I(lsst.geom.Point2I(90, -200), lsst.geom.Point2I(800, 800)))
768 for kernel, border in (("bilinear", 0.5), ("lanczos3", 2.5)):
769 with self.subTest(coaddWarpKernel=kernel):
770 result = self._runTract0(self._patchHandles(0, [right], epochs=[(1, 60000.25, 1.0)]),
771 outputImageType="future", coaddWarpKernel=kernel)
772 bounds = result.templateInputs[0, 0][0].bounds
773 self.assertFloatsAlmostEqual(bounds.x_vertices.min(), 89.5 + border, rtol=0, atol=1e-6)
775 def testRunTemplateInputsCoverScienceEdges(self):
776 """An input's footprint reaches the edges of a distorted science image
777 wherever the input covers them.
778 """
779 # Science pixel (x, y) lands on input pixel (x, y + a*(x - xc)**2),
780 # which distorts the top and bottom edges of the science image by about
781 # 2 pixels between its corners.
782 xc = 90.0
783 left = (self.exposure.wcs, lsst.geom.Box2I(lsst.geom.Point2I(-200, -200), lsst.geom.Point2I(89, 800)))
784 # The centers of the edge pixels that the input's warp covers; warping
785 # with lanczos5 left the 4.5 pixels along the input's edges out.
786 x = np.arange(0.0, 85.0)
787 y = np.arange(0.0, 181.0)
788 edgeX = np.concatenate([x, x, np.zeros_like(y)])
789 edgeY = np.concatenate([np.zeros_like(x), np.full_like(x, 180.0), y])
790 for a in (-2.5e-4, 2.5e-4):
791 with self.subTest(a=a):
792 forward = np.array([[1.0, 1, 1, 0], [1.0, 2, 0, 1], [a, 2, 2, 0], [-2*a*xc, 2, 1, 0],
793 [a*xc**2, 2, 0, 0]])
794 inverse = np.array([[1.0, 1, 1, 0], [1.0, 2, 0, 1], [-a, 2, 2, 0], [2*a*xc, 2, 1, 0],
795 [-a*xc**2, 2, 0, 0]])
796 distortion = lsst.afw.geom.TransformPoint2ToPoint2(astshim.PolyMap(forward, inverse))
797 wcs = lsst.afw.geom.makeModifiedWcs(distortion, self.exposure.wcs, False)
798 result = self._runTract0(self._patchHandles(0, [left], epochs=[(1, 60000.25, 1.0)]),
799 outputImageType="future", wcs=wcs)
800 bounds = result.templateInputs[0, 0][0].bounds
801 self.assertTrue(np.all(bounds.contains(x=edgeX, y=edgeY)))
803 def _runTract0(self, handles, outputImageType="legacy", wcs=None, coaddWarpKernel=None):
804 """Build a template from the given tract 0 patch handles, with a task
805 configured for ``outputImageType`` and ``coaddWarpKernel``, for a
806 science image with ``wcs``, by default that of the test exposure.
807 """
808 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
809 config = self._makeConfig()
810 config.output_image_type = outputImageType
811 if coaddWarpKernel is not None:
812 config.coaddWarpKernel = coaddWarpKernel
813 task = lsst.ip.diffim.GetTemplateTask(config=config)
814 return task.run(coaddExposureHandles={0: handles},
815 bbox=box,
816 wcs=self.exposure.wcs if wcs is None else wcs,
817 dataIds={0: self.dataIds[0]},
818 physical_filter="a_test")
820 def testRunRecordsInstrument(self):
821 """The instrument of each coadd is recorded with the template's coadd
822 inputs.
823 """
824 result = self._runTract0(self.patches[0])
825 ccds = result.template.getInfo().getCoaddInputs().ccds
826 self.assertEqual(len(ccds), 4)
827 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"})
829 def testRunUnknownInstrument(self):
830 """A coadd whose inputs record no instrument is dropped when the
831 other coadds record one, and kept, with the template's instrument
832 ``"UNKNOWN"``, when none of them records one.
833 """
834 records = [(self.exposure.wcs, self.exposure.getBBox())]
835 # Inputs with no visit info, and with a visit info but no label.
836 for unknown in (None, ""):
837 with self.subTest(label=unknown):
838 handles = (self._patchHandles(0, records, instrument=unknown)[:1]
839 + self._patchHandles(0, records)[1:])
840 ccds, dropped = self._runLogged({0: handles})
841 self.assertEqual(len(ccds), 3)
842 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"})
843 self.assertEqual(len(dropped), 1)
844 self.assertIn("made from UNKNOWN", dropped[0])
846 # No coadd records an instrument.
847 result = self._runTract0(self._patchHandles(0, records, instrument=unknown))
848 ccds = result.template.getInfo().getCoaddInputs().ccds
849 self.assertEqual(len(ccds), 4)
850 self.assertEqual({record["instrument"] for record in ccds}, {"UNKNOWN"})
852 def _runLogged(self, handles, dataIds=None, box=None, instrument=None):
853 """Build a template, capturing the messages about dropped coadds.
855 Parameters
856 ----------
857 handles : `dict` [`int`, `list`]
858 Patch handles, keyed by tract.
859 dataIds : `dict` [`int`, `list`], optional
860 Data ids matching ``handles``; defaults to the fixture's.
861 box : `lsst.geom.Box2I`, optional
862 Box to build the template on; defaults to one in tract 0.
863 instrument : `str`, optional
864 Instrument of the science image.
866 Returns
867 -------
868 ccds : `lsst.afw.table.ExposureCatalog`
869 The coadd inputs recorded on the template.
870 dropped : `list` [`str`]
871 The info messages about dropped coadds.
872 """
873 if dataIds is None:
874 dataIds = {tract: self.dataIds[tract] for tract in handles}
875 if box is None:
876 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
877 task = lsst.ip.diffim.GetTemplateTask(config=self._makeConfig())
878 with self.assertLogs(task.log.name, level="INFO") as cm:
879 result = task.run(coaddExposureHandles=handles,
880 bbox=box,
881 wcs=self.exposure.wcs,
882 dataIds=dataIds,
883 physical_filter="a_test",
884 instrument=instrument)
885 dropped = [record for record in cm.records
886 if record.getMessage().startswith("Dropping template coadd")]
887 self.assertEqual({record.levelname for record in dropped}, {"INFO"} if dropped else set())
888 return result.template.getInfo().getCoaddInputs().ccds, [record.getMessage() for record in dropped]
890 def testRunDropsMinorityInstrument(self):
891 """Coadds made from an instrument other than the one most coadds were
892 made from are dropped.
893 """
894 records = [(self.exposure.wcs, self.exposure.getBBox())]
895 handles = (self._patchHandles(0, records)[:3]
896 + self._patchHandles(0, records, instrument="OtherCam")[3:])
898 ccds, dropped = self._runLogged({0: handles})
899 self.assertEqual({record["patch"] for record in ccds}, {0, 1, 2})
900 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"})
901 self.assertEqual(len(dropped), 1)
902 self.assertIn(str(self.dataIds[0][3]), dropped[0])
903 self.assertIn("made from OtherCam", dropped[0])
905 def testRunPrefersScienceInstrument(self):
906 """Coadds made from the science image's instrument are used even when
907 they are in the minority; otherwise the majority instrument is used.
908 """
909 records = [(self.exposure.wcs, self.exposure.getBBox())]
910 handles = (self._patchHandles(0, records, instrument="OtherCam")[:3]
911 + self._patchHandles(0, records)[3:])
913 ccds, dropped = self._runLogged({0: handles}, instrument="TestCam")
914 self.assertEqual({record["patch"] for record in ccds}, {3})
915 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"})
916 self.assertEqual(len(dropped), 3)
917 self.assertTrue(all("made from OtherCam" in message for message in dropped))
919 # No coadd was made from the science image's instrument.
920 ccds, dropped = self._runLogged({0: handles}, instrument="SciCam")
921 self.assertEqual({record["patch"] for record in ccds}, {0, 1, 2})
922 self.assertEqual({record["instrument"] for record in ccds}, {"OtherCam"})
923 self.assertEqual(len(dropped), 1)
925 def testRunInstrumentTie(self):
926 """A tie goes to the instrument whose name sorts first.
927 """
928 records = [(self.exposure.wcs, self.exposure.getBBox())]
929 handles = (self._patchHandles(0, records)[:2]
930 + self._patchHandles(0, records, instrument="OtherCam")[2:])
932 ccds, dropped = self._runLogged({0: handles})
933 self.assertEqual({record["patch"] for record in ccds}, {2, 3})
934 self.assertEqual({record["instrument"] for record in ccds}, {"OtherCam"})
935 self.assertEqual(len(dropped), 2)
937 def testRunDropsNonOverlappingCoadd(self):
938 """A coadd from another instrument is dropped even if it does not
939 overlap the template.
940 """
941 records = [(self.exposure.wcs, self.exposure.getBBox())]
942 # Tract 3 does not overlap a box in the lower left of tract 0.
943 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(50, 50))
944 handles = {0: self.patches[0], 3: self._patchHandles(3, records, instrument="OtherCam")[:1]}
945 dataIds = {0: self.dataIds[0], 3: self.dataIds[3][:1]}
947 ccds, dropped = self._runLogged(handles, dataIds=dataIds, box=box)
948 self.assertEqual({record["tract"] for record in ccds}, {0})
949 self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"})
950 self.assertEqual(len(dropped), 1)
951 self.assertIn(str(self.dataIds[3][0]), dropped[0])
953 def _scaleInputVariance(self, tract, factor):
954 """Return fresh handles for one tract's patches, with their variance
955 planes multiplied by ``factor``.
957 Parameters
958 ----------
959 tract : `int`
960 Id of the tract whose patches should be copied.
961 factor : `float`
962 Factor to multiply the input variance planes by.
964 Returns
965 -------
966 handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`]
967 Handles to the modified patches.
968 """
969 handles = []
970 for handle in self.patches[tract]:
971 # ``copy=True`` on the original handles means this is a copy, so
972 # the patches shared with the other tests are left untouched.
973 patch = handle.get()
974 patch.variance.array *= factor
975 handles.append(pipeBase.InMemoryDatasetHandle(patch,
976 storageClass="ExposureF",
977 copy=True,
978 dataId=handle.dataId))
979 return handles
981 def testScaleVariance(self):
982 """Test that the template variance plane is rescaled to match the
983 empirical pixel noise, and that the factor used is recorded in the
984 task metadata.
985 """
986 scaleFactor = 1.345
987 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
989 def _configureAndRunTask(doScaleVariance, varianceScale=1.):
990 """Build a template from tract 0, optionally rescaling the input
991 variance planes by ``varianceScale`` first.
992 """
993 config = self._makeConfig()
994 config.doScaleVariance = doScaleVariance
995 task = lsst.ip.diffim.GetTemplateTask(config=config)
996 result = task.run(coaddExposureHandles={0: self._scaleInputVariance(0, varianceScale)},
997 bbox=box,
998 wcs=self.exposure.wcs,
999 dataIds={0: self.dataIds[0]},
1000 physical_filter="a_test")
1001 return task, result.template
1003 # With scaling disabled the subtask is never constructed, and nothing
1004 # is recorded in the metadata.
1005 taskOff, templateOff = _configureAndRunTask(doScaleVariance=False)
1006 self.assertFalse(hasattr(taskOff, "scaleVariance"))
1007 self.assertNotIn("scaleTemplateVarianceFactor", taskOff.metadata)
1009 # The image is correlated after warping, so ``doScaleVariance`` will
1010 # apply a correction. The value here is empirical from a warped image
1011 # with white noise and no sources.
1012 taskOn, templateOn = _configureAndRunTask(doScaleVariance=True)
1013 factor = taskOn.metadata["scaleTemplateVarianceFactor"]
1014 self.assertFloatsAlmostEqual(factor, 1.145, atol=0.01,
1015 msg="Measured template variance scaling changed; check the"
1016 " warping kernels used by _makePatches and the task,"
1017 " and the order of the variance corrections.")
1018 # The only difference from the unscaled template is the constant
1019 # factor applied to the variance plane.
1020 self.assertFloatsAlmostEqual(templateOn.variance.array,
1021 templateOff.variance.array*factor, rtol=1e-5)
1022 # Tolerance here is float32 round-off: repeated runs of the task are
1023 # not bitwise identical.
1024 self.assertImagesAlmostEqual(templateOn.image, templateOff.image, rtol=1e-5, atol=1e-5)
1026 # If the input variance planes under-estimate the noise by a known
1027 # factor, the measured factor grows by that amount and the same
1028 # output variance plane is recovered.
1029 taskLow, templateLow = _configureAndRunTask(doScaleVariance=True, varianceScale=1/scaleFactor)
1030 self.assertFloatsAlmostEqual(taskLow.metadata["scaleTemplateVarianceFactor"],
1031 factor*scaleFactor, rtol=1e-5)
1032 self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5)
1034 def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True, handles=None):
1035 """Build a template from tract 0 with a task configured for the
1036 future output type, without converting it.
1038 The task is kept as ``self.futureTask``.
1040 Parameters
1041 ----------
1042 raiseOnUndefinedMaskMap : `bool`, optional
1043 Value of the task config field of the same name.
1044 handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`], optional
1045 Handles to the tract 0 patches to use instead of the fixture's.
1047 Returns
1048 -------
1049 result : `lsst.pipe.base.Struct`
1050 The legacy output struct.
1051 box : `lsst.geom.Box2I`
1052 The bounding box the template was requested on, before the task
1053 grew it by the template border.
1054 """
1055 config = self._makeConfig()
1056 config.output_image_type = "future"
1057 config.raiseOnUndefinedMaskMap = raiseOnUndefinedMaskMap
1058 task = lsst.ip.diffim.GetTemplateTask(config=config)
1059 self.visit = 9876
1060 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180))
1061 result = task.run(coaddExposureHandles={0: self.patches[0] if handles is None else handles},
1062 bbox=box,
1063 wcs=self.exposure.wcs,
1064 dataIds={0: self.dataIds[0]},
1065 physical_filter="a_test",
1066 visit=self.visit)
1067 # run() always produces legacy types, whatever output_image_type is.
1068 self.assertIsInstance(result.template, lsst.afw.image.ExposureF)
1069 self.futureTask = task
1070 return result, box
1072 def _runFuture(self, handles=None):
1073 """Build a template and convert it, the way ``runQuantum`` does.
1075 Parameters
1076 ----------
1077 handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`], optional
1078 Handles to the tract 0 patches to use instead of the fixture's.
1080 Returns
1081 -------
1082 result : `lsst.pipe.base.Struct`
1083 The converted output struct.
1084 box : `lsst.geom.Box2I`
1085 The bounding box the template was requested on, before the task
1086 grew it by the template border.
1087 """
1088 # lsst.images needs per-amplifier raw geometry and a field angle
1089 # transform, which the trivial test detectors do not have.
1090 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0]
1091 result, box = self._runLegacyForFuture(handles=handles)
1092 result.template = self.futureTask.convert_outputs_to_future(
1093 result.template, self._coaddRefs(0), detector=detector,
1094 exposureRecord=self._exposureRecord(), templateInputs=result.templateInputs)
1095 return result, box, detector
1097 def _exposureRecord(self):
1098 """Return the exposure record that ``runQuantum`` rebuilds from the
1099 observation metadata of the science image.
1100 """
1101 return makeTestExposureRecord(DimensionUniverse(), visit=self.visit)
1103 def _coaddRefs(self, tract):
1104 """Return butler references for the coadds of one tract, like the
1105 ones ``runQuantum`` passes on from its ``coaddExposures`` input.
1107 Parameters
1108 ----------
1109 tract : `int`
1110 Id of the tract whose coadds to make references for.
1112 Returns
1113 -------
1114 refs : `list` [`lsst.daf.butler.DatasetRef`]
1115 One reference per patch of that tract. The same references are
1116 returned on every call, because each new one would get a new
1117 dataset id.
1118 """
1119 if tract not in self.coaddRefs:
1120 datasetType = DatasetType("template_coadd", ("tract", "patch", "band", "skymap"),
1121 "ExposureF", universe=DimensionUniverse())
1122 self.coaddRefs[tract] = [DatasetRef(datasetType, dataId, run="test_run")
1123 for dataId in self.dataIds[tract]]
1124 return self.coaddRefs[tract]
1126 def testConvertOutputsToFuture(self):
1127 """Test that the output template is converted to a DifferenceImage
1128 that keeps the pixels and the detector it was built on.
1129 """
1130 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0]
1131 result, box = self._runLegacyForFuture()
1132 legacy = result.template.clone()
1133 result.template = self.futureTask.convert_outputs_to_future(
1134 result.template, self._coaddRefs(0), detector=detector,
1135 exposureRecord=self._exposureRecord())
1137 self.assertIsInstance(result.template, lsst.images.DifferenceImage)
1138 compare_masked_image_to_legacy(result.template, legacy.maskedImage,
1139 plane_map=lsst.images.get_legacy_template_mask_planes())
1140 self.assertEqual(result.template.unit, u.nJy)
1141 self.assertIsNotNone(result.template.detector)
1142 self.assertEqual(result.template.detector.name, detector.getName())
1143 # The template is grown by the border, so it is larger than both the
1144 # requested box and the detector.
1145 expectedBox = lsst.geom.Box2I(box)
1146 expectedBox.grow(lsst.ip.diffim.GetTemplateTask.ConfigClass().templateBorderSize)
1147 self.assertEqual(result.template.bbox.to_legacy(), expectedBox)
1148 self.assertEqual(result.template.obs_info.visit_id, self.visit)
1150 def testConvertOutputsToFutureCoaddMaskPlanes(self):
1151 """Test that the coadd mask planes a warped template carries are
1152 converted, and come back when the image is converted to legacy.
1154 The template is converted with the template plane map, so the coadd
1155 planes are named there; the source injection planes are optional, and
1156 are added because this template has pixels set in them.
1157 """
1158 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0]
1159 planes = ("CLIPPED", "REJECTED", "INEXACT_PSF", "SENSOR_EDGE", "HIGH_VARIANCE",
1160 "INJECTED", "INJECTED_CORE")
1161 result, _ = self._runLegacyForFuture()
1162 mask = result.template.mask
1163 for n, plane in enumerate(planes):
1164 mask.addMaskPlane(plane)
1165 mask.array[0, n] |= mask.getPlaneBitMask(plane)
1166 legacyTemplate = result.template.clone()
1168 result.template = self.futureTask.convert_outputs_to_future(
1169 result.template, self._coaddRefs(0), detector=detector,
1170 exposureRecord=self._exposureRecord())
1172 template = result.template
1173 self.assertIsInstance(template, lsst.images.DifferenceImage)
1174 # Name the optional planes that have pixels set in the map, so that
1175 # dropping one fails the comparison instead of skipping it. The others
1176 # are left out because they may be registered by unrelated tests.
1177 optional = lsst.images.get_legacy_optional_mask_planes()
1178 planeMap = lsst.images.get_legacy_template_mask_planes()
1179 planeMap.update({plane: optional[plane] for plane in planes if plane in optional})
1180 compare_masked_image_to_legacy(template, legacyTemplate.maskedImage, plane_map=planeMap)
1181 # The butler converts a DifferenceImage to an ExposureF with the
1182 # difference image plane map, which does not name the coadd planes;
1183 # they keep their own names instead of being dropped.
1184 legacy = template.to_legacy()
1185 for n, plane in enumerate(planes):
1186 with self.subTest(plane=plane):
1187 self.assertIn(plane, legacy.mask.getMaskPlaneDict())
1188 self.assertEqual(
1189 np.count_nonzero(legacy.mask.array & legacy.mask.getPlaneBitMask(plane)), 1)
1191 def _addUnmappedMaskPlane(self, template):
1192 """Set one pixel in the template's BRIGHT_OBJECT plane, which the
1193 template plane map does not include, and one in its CLIPPED plane,
1194 which it does.
1195 """
1196 mask = template.mask
1197 for n, plane in enumerate(("BRIGHT_OBJECT", "CLIPPED")):
1198 mask.addMaskPlane(plane)
1199 mask.array[0, n] |= mask.getPlaneBitMask(plane)
1201 def testConvertOutputsToFutureUnmappedMaskPlaneRaises(self):
1202 """With raiseOnUndefinedMaskMap=True, a template with pixels set in
1203 an unmapped mask plane cannot be converted.
1204 """
1205 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0]
1206 result, _ = self._runLegacyForFuture()
1207 self._addUnmappedMaskPlane(result.template)
1208 with self.assertRaisesRegex(RuntimeError, "BRIGHT_OBJECT"):
1209 self.futureTask.convert_outputs_to_future(
1210 result.template, self._coaddRefs(0), detector=detector,
1211 exposureRecord=self._exposureRecord())
1213 def testConvertOutputsToFutureUnmappedMaskPlaneWarns(self):
1214 """With raiseOnUndefinedMaskMap=False, an unmapped mask plane is
1215 dropped with a warning and the mapped planes are kept.
1216 """
1217 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0]
1218 result, _ = self._runLegacyForFuture(raiseOnUndefinedMaskMap=False)
1219 self._addUnmappedMaskPlane(result.template)
1220 with self.assertLogs(self.futureTask.log.name, level="WARNING") as cm:
1221 template = self.futureTask.convert_outputs_to_future(
1222 result.template, self._coaddRefs(0), detector=detector,
1223 exposureRecord=self._exposureRecord())
1224 self.assertIn("BRIGHT_OBJECT", "\n".join(cm.output))
1225 self.assertIsInstance(template, lsst.images.DifferenceImage)
1226 self.assertNotIn("BRIGHT_OBJECT", template.mask.schema.names)
1227 self.assertEqual(np.count_nonzero(template.mask.get("CLIPPED")), 1)
1229 def testConvertOutputsToFutureLosesProvenance(self):
1230 """The coadd inputs that `run` attaches are not carried by
1231 `lsst.images.DifferenceImage`.
1232 """
1233 legacy, _ = self._runLegacyForFuture()
1234 self.assertTrue(legacy.template.getInfo().hasCoaddInputs())
1235 result, _, _ = self._runFuture()
1236 self.assertFalse(result.template.to_legacy().getInfo().hasCoaddInputs())
1238 def testConvertOutputsToFutureTemplates(self):
1239 """The coadds that went into the template are recorded on the
1240 converted image.
1241 """
1242 result, _, _ = self._runFuture()
1243 refs = {(ref.dataId["tract"], ref.dataId["patch"]): ref for ref in self._coaddRefs(0)}
1245 templates = result.template.templates
1246 self.assertGreater(len(templates), 0)
1247 for template in templates:
1248 with self.subTest(tract=template.tract, patch=template.patch):
1249 ref = refs[template.tract, template.patch]
1250 self.assertEqual(template.skymap, "skymap")
1251 self.assertEqual(template.dataset_id, ref.id)
1252 self.assertEqual(template.dataset_run, ref.run)
1253 self.assertEqual(template.instrument, "TestCam")
1254 self.assertFalse(template.psf_shape_flag)
1255 self.assertGreater(template.psf_shape_xx, 0)
1256 self.assertGreater(template.bounds.area, 0)
1258 def testConvertOutputsToFutureInstrument(self):
1259 """The template instrument is the one recorded by the coadds that
1260 know it, or ``"UNKNOWN"`` if none of them do.
1261 """
1262 records = [(self.exposure.wcs, self.exposure.getBBox())]
1263 handles = (self._patchHandles(0, records, instrument=None)[:1]
1264 + self._patchHandles(0, records)[1:])
1265 result, _, _ = self._runFuture(handles=handles)
1266 self.assertEqual({template.instrument for template in result.template.templates}, {"TestCam"})
1268 result, _, _ = self._runFuture(handles=self._patchHandles(0, records, instrument=None))
1269 self.assertEqual({template.instrument for template in result.template.templates}, {"UNKNOWN"})
1271 def testConvertOutputsToFutureTemplateInputs(self):
1272 """Each coadd recorded on the converted template lists its inputs.
1273 """
1274 records = [(self.exposure.wcs, self.exposure.getBBox())]
1275 handles = self._patchHandles(0, records, epochs=[(1, 60000.25, 1.0)])
1276 result, _, _ = self._runFuture(handles=handles)
1277 templates = result.template.templates
1278 self.assertGreater(len(templates), 0)
1279 for template in templates:
1280 with self.subTest(tract=template.tract, patch=template.patch):
1281 self.assertEqual([(i.visit, i.detector, i.weight) for i in template.inputs], [(1, 0, 1.0)])
1282 self.assertFloatsAlmostEqual(template.inputs[0].mid_mjd_tai, 60000.25, rtol=0, atol=1e-6)
1284 def testConvertOutputsToFuturePsf(self):
1285 """The Gaussian PSF set by `run` is kept, defined over the detector.
1286 """
1287 result, _ = self._runLegacyForFuture()
1288 detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0]
1289 legacyPsf = result.template.getPsf()
1291 result.template = self.futureTask.convert_outputs_to_future(
1292 result.template, self._coaddRefs(0), detector=detector,
1293 exposureRecord=self._exposureRecord())
1295 psf = result.template.psf
1296 self.assertIsInstance(psf, lsst.images.psfs.GaussianPointSpreadFunction)
1297 self.assertFloatsAlmostEqual(psf.sigma, legacyPsf.getSigma())
1298 self.assertEqual(psf.bounds.to_legacy(), detector.getBBox())
1299 self.assertEqual(psf.kernel_bbox.to_legacy().getDimensions(), legacyPsf.getDimensions())
1301 def testConvertOutputsToFuturePsfToLegacy(self):
1302 """The Gaussian survives conversion back to an Exposure, which is how
1303 a task that has not been converted reads the template.
1304 """
1305 result, _, _ = self._runFuture()
1307 psf = result.template.to_legacy().getPsf()
1308 self.assertIsInstance(psf, lsst.afw.detection.GaussianPsf)
1309 self.assertFloatsAlmostEqual(psf.getSigma(), result.template.psf.sigma)
1312class WarpCorrelationFactorTestCase(lsst.utils.tests.TestCase):
1313 """Test the warping correlation correction on pure noise, with no other
1314 variance correction.
1315 """
1317 def setUp(self):
1318 self.size = 1024
1319 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(self.size, self.size))
1321 def _makeWcs(self, angle, offset):
1322 """Make a WCS with a fixed plate scale, rotated by ``angle`` degrees
1323 and with its reference pixel moved by ``offset`` pixels.
1324 """
1325 center = lsst.geom.Point2D(self.size/2 + offset[0], self.size/2 + offset[1])
1326 cdMatrix = lsst.afw.geom.makeCdMatrix(0.2*lsst.geom.arcseconds, angle*lsst.geom.degrees)
1327 return lsst.afw.geom.makeSkyWcs(center, lsst.geom.SpherePoint(45, 30, lsst.geom.degrees),
1328 cdMatrix)
1330 def testCorrelationFactorRestoresNoisePower(self):
1331 """Warp pure noise to a rotated, offset frame with the coadd kernel,
1332 then to a third frame with the task's own warper.
1334 Warping leaves the variance of a flux summed over a region (the noise
1335 power) unchanged, but lowers the per-pixel scatter and the variance
1336 plane. The variance plane is lowered by a larger factor because the
1337 input to the second warp is already correlated. The correction must
1338 bring the variance plane back to the noise power of the original image.
1339 """
1340 sigma = 2.0
1341 rng = np.random.default_rng(42)
1342 original = lsst.afw.image.ExposureF(self.bbox, self._makeWcs(0, (0, 0)))
1343 original.image.array[:] = rng.normal(0, sigma, original.image.array.shape)
1344 original.variance.array[:] = sigma**2
1346 config = lsst.ip.diffim.GetTemplateTask.ConfigClass()
1347 config.doScaleVariance = False
1348 config.coaddWarpKernel = "lanczos3"
1349 task = lsst.ip.diffim.GetTemplateTask(config=config)
1350 coadd = lsst.afw.math.Warper(config.coaddWarpKernel).warpExposure(
1351 self._makeWcs(37, (0.31, 0.72)), original, destBBox=self.bbox)
1352 template = task.warper.warpExposure(self._makeWcs(-21, (0.58, 0.23)), coadd, destBBox=self.bbox)
1354 # The central quarter is covered in every frame.
1355 inner = lsst.geom.Box2I(lsst.geom.Point2I(self.size//4, self.size//4),
1356 lsst.geom.Extent2I(self.size//2, self.size//2))
1357 originalStd = original.image[inner].array.std()
1358 image = template.image[inner].array.astype(float)
1359 pixelStd = image.std()
1360 # Each warp moves flux between neighboring pixels but doesn't create
1361 # or destroy any, so the noise in a flux summed over a region is the
1362 # same before and after. Summing the image over blocks and measuring
1363 # how much those sums scatter gives that quantity directly, without
1364 # depending on the variance plane. Every placement of the block is
1365 # used, from a summed-area table, to reduce the noise.
1366 blockSize = 8
1367 table = np.pad(image.cumsum(axis=0).cumsum(axis=1), ((1, 0), (1, 0)))
1368 sums = (table[blockSize:, blockSize:] - table[:-blockSize, blockSize:]
1369 - table[blockSize:, :-blockSize] + table[:-blockSize, :-blockSize])
1370 blockStd = sums.std()/blockSize
1371 planeStdBefore = np.sqrt(template.variance[inner].array.mean())
1372 task._correctVariance(template, plateScaleFactor=1.0)
1373 planeStdAfter = np.sqrt(template.variance[inner].array.mean())
1374 factor = task.metadata["varianceWarpCorrelationFactor"]
1376 # afw propagates the variance as if each warp's input pixels were
1377 # independent, so the plane drops by exactly the factor's reciprocal.
1378 self.assertFloatsAlmostEqual(planeStdBefore, sigma/np.sqrt(factor), rtol=1e-5)
1379 self.assertFloatsAlmostEqual(planeStdAfter, sigma, rtol=1e-5)
1381 # The warps conserve the noise power, which the corrected plane
1382 # matches.
1383 self.assertFloatsAlmostEqual(blockStd, originalStd, rtol=0.03)
1384 self.assertFloatsAlmostEqual(blockStd, planeStdAfter, rtol=0.03)
1386 # The per-pixel scatter drops, though by less than the uncorrected
1387 # plane, because the second warp receives correlated input. The
1388 # corrected plane deliberately exceeds it.
1389 self.assertLess(pixelStd, 0.9*originalStd)
1390 self.assertGreater(pixelStd, planeStdBefore)
1391 self.assertLess(pixelStd, planeStdAfter)
1394class GetTemplateConnectionsTestCase(lsst.utils.tests.TestCase):
1395 """Test the connections that ``output_image_type`` switches on.
1397 These only need the config classes, not a built template.
1398 """
1400 def testOutputImageType(self):
1401 """Test that the future output type changes the template storage
1402 class, and leaves the input connections alone.
1403 """
1404 Connections = lsst.ip.diffim.GetTemplateTask.ConfigClass.ConnectionsClass
1405 config = lsst.ip.diffim.GetTemplateTask.ConfigClass()
1407 connections = Connections(config=config)
1408 self.assertEqual(connections.template.storageClass, "ExposureF")
1409 self.assertNotIn("obs_info", connections.inputs)
1411 config.output_image_type = "future"
1412 connections = Connections(config=config)
1413 self.assertEqual(connections.template.storageClass, "DifferenceImage")
1414 # The dataset name is unchanged; only its storage class differs.
1415 self.assertEqual(connections.template.name, "goodSeeingDiff_templateExp")
1416 # The detector the conversion needs comes from the science image.
1417 self.assertEqual(connections.detector.name, "calexp.detector")
1418 self.assertEqual(connections.detector.storageClass, "Detector")
1419 # The observation metadata is only read in future mode, where it is
1420 # the one input the conversion adds.
1421 self.assertEqual(connections.obs_info.name, "calexp.obs_info")
1422 self.assertEqual(connections.obs_info.storageClass, "ObservationInfo")
1424 def testLintConnections(self):
1425 """Check that the connections are self-consistent in both modes.
1426 """
1427 for task in (lsst.ip.diffim.GetTemplateTask, lsst.ip.diffim.GetDcrTemplateTask):
1428 with self.subTest(task=task.__name__):
1429 lsst.pipe.base.testUtils.lintConnections(task.ConfigClass.ConnectionsClass)
1432def setup_module(module):
1433 lsst.utils.tests.init()
1436class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
1437 pass
1440if __name__ == "__main__": 1440 ↛ 1441line 1440 didn't jump to line 1441 because the condition on line 1440 was never true
1441 lsst.utils.tests.init()
1442 unittest.main()