Coverage for tests / test_warpExposure.py: 16%
336 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 01:50 -0700
1# This file is part of afw.
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/>.
22"""Test warpExposure
23"""
24import os
25import unittest
27import numpy as np
29import lsst.utils
30import lsst.utils.logging
31import lsst.utils.tests
32import lsst.daf.base as dafBase
33from lsst.afw.coord import Observatory, Weather
34import lsst.geom
35import lsst.afw.geom as afwGeom
36import lsst.afw.image as afwImage
37import lsst.afw.math as afwMath
38import lsst.pex.exceptions as pexExcept
39import lsst.afw.display as afwDisplay
40from lsst.log import Log
42# Change the level to Log.DEBUG to see debug messages
43Log.getLogger("lsst.afw.image.Mask").setLevel(Log.INFO)
45# Change integer argument to adjust trace logging.
46lsst.utils.logging.trace_set_at("lsst.afw.math.warp", -1)
48afwDisplay.setDefaultMaskTransparency(75)
50display = False
51# set True to save afw-warped images as FITS files
52SAVE_FITS_FILES = False
53# set True to save failed afw-warped images as FITS files even if
54# SAVE_FITS_FILES is False
55SAVE_FAILED_FITS_FILES = True
57try:
58 afwdataDir = lsst.utils.getPackageDir("afwdata")
59except LookupError:
60 afwdataDir = None
61else:
62 dataDir = os.path.join(afwdataDir, "data")
64 originalExposureName = "medexp.fits"
65 originalExposurePath = os.path.join(dataDir, originalExposureName)
66 subExposureName = "medsub.fits"
67 subExposurePath = os.path.join(dataDir, originalExposureName)
68 originalFullExposureName = os.path.join(
69 "CFHT", "D4", "cal-53535-i-797722_1.fits")
70 originalFullExposurePath = os.path.join(dataDir, originalFullExposureName)
73def makeVisitInfo():
74 """Return a non-NaN visitInfo."""
75 return afwImage.VisitInfo(exposureTime=10.01,
76 darkTime=11.02,
77 date=dafBase.DateTime(65321.1, dafBase.DateTime.MJD, dafBase.DateTime.TAI),
78 ut1=12345.1,
79 era=45.1*lsst.geom.degrees,
80 boresightRaDec=lsst.geom.SpherePoint(23.1, 73.2, lsst.geom.degrees),
81 boresightAzAlt=lsst.geom.SpherePoint(134.5, 33.3, lsst.geom.degrees),
82 boresightAirmass=1.73,
83 boresightRotAngle=73.2*lsst.geom.degrees,
84 rotType=afwImage.RotType.SKY,
85 observatory=Observatory(11.1*lsst.geom.degrees, 22.2*lsst.geom.degrees, 0.333),
86 weather=Weather(1.1, 2.2, 34.5),
87 )
90class WarpExposureTestCase(lsst.utils.tests.TestCase):
91 """Test case for warpExposure
92 """
94 def setUp(self):
95 np.random.seed(0)
97 @unittest.skipIf(afwdataDir is None, "afwdata not setup")
98 def testNullWarpExposure(self, interpLength=10):
99 """Test that warpExposure maps an image onto itself.
101 Note:
102 - NO_DATA and off-CCD pixels must be ignored
103 - bad mask pixels get smeared out so we have to excluded all bad mask pixels
104 from the output image when comparing masks.
105 """
106 originalExposure = afwImage.ExposureF(originalExposurePath)
107 originalExposure.getInfo().setId(10313423)
108 originalExposure.getInfo().setVisitInfo(makeVisitInfo())
109 originalFilterLabel = afwImage.FilterLabel(band="i")
110 originalPhotoCalib = afwImage.PhotoCalib(1.0e5, 1.0e3)
111 originalExposure.setFilter(originalFilterLabel)
112 originalExposure.setPhotoCalib(originalPhotoCalib)
113 afwWarpedExposure = afwImage.ExposureF(
114 originalExposure.getBBox(),
115 originalExposure.getWcs())
116 warpingControl = afwMath.WarpingControl(
117 "lanczos4", "", 0, interpLength)
118 afwMath.warpExposure(
119 afwWarpedExposure, originalExposure, warpingControl)
120 if SAVE_FITS_FILES:
121 afwWarpedExposure.writeFits("afwWarpedExposureNull.fits")
123 self.assertEqual(afwWarpedExposure.getFilter().bandLabel,
124 originalFilterLabel.bandLabel)
125 self.assertEqual(afwWarpedExposure.getPhotoCalib(), originalPhotoCalib)
126 self.assertEqual(afwWarpedExposure.getInfo().getVisitInfo(),
127 originalExposure.getInfo().getVisitInfo())
129 afwWarpedMaskedImage = afwWarpedExposure.getMaskedImage()
130 afwWarpedMask = afwWarpedMaskedImage.getMask()
131 noDataBitMask = afwWarpedMask.getPlaneBitMask("NO_DATA")
133 # compare all non-DATA pixels of image and variance, but relax specs a bit
134 # because of minor noise introduced by bad pixels
135 noDataMaskArr = afwWarpedMaskedImage.mask.array & noDataBitMask
136 msg = "afw null-warped MaskedImage (all pixels, relaxed tolerance)"
137 self.assertMaskedImagesAlmostEqual(afwWarpedMaskedImage, originalExposure.getMaskedImage(),
138 doMask=False, skipMask=noDataMaskArr, atol=1e-5, msg=msg)
140 # compare good pixels (mask=0) of image, mask and variance using full
141 # tolerance
142 msg = "afw null-warped MaskedImage (good pixels, max tolerance)"
143 self.assertMaskedImagesAlmostEqual(afwWarpedMaskedImage, originalExposure.getMaskedImage(),
144 skipMask=afwWarpedMask, msg=msg)
146 @unittest.skipIf(afwdataDir is None, "afwdata not setup")
147 def testNullWarpImage(self, interpLength=10):
148 """Test that warpImage maps an image onto itself.
149 """
150 originalExposure = afwImage.ExposureF(originalExposurePath)
151 afwWarpedExposure = afwImage.ExposureF(originalExposurePath)
152 originalImage = originalExposure.getMaskedImage().getImage()
153 afwWarpedImage = afwWarpedExposure.getMaskedImage().getImage()
154 originalWcs = originalExposure.getWcs()
155 afwWarpedWcs = afwWarpedExposure.getWcs()
156 warpingControl = afwMath.WarpingControl(
157 "lanczos4", "", 0, interpLength)
158 afwMath.warpImage(afwWarpedImage, afwWarpedWcs,
159 originalImage, originalWcs, warpingControl)
160 if SAVE_FITS_FILES:
161 afwWarpedImage.writeFits("afwWarpedImageNull.fits")
162 afwWarpedImageArr = afwWarpedImage.getArray()
163 noDataMaskArr = np.isnan(afwWarpedImageArr)
164 # relax specs a bit because of minor noise introduced by bad pixels
165 msg = "afw null-warped Image"
166 self.assertImagesAlmostEqual(originalImage, afwWarpedImage, skipMask=noDataMaskArr,
167 atol=1e-5, msg=msg)
169 @unittest.skipIf(afwdataDir is None, "afwdata not setup")
170 def testNullWcs(self, interpLength=10):
171 """Cannot warp from or into an exposure without a Wcs.
172 """
173 exposureWithWcs = afwImage.ExposureF(originalExposurePath)
174 mi = exposureWithWcs.getMaskedImage()
175 exposureWithoutWcs = afwImage.ExposureF(mi.getDimensions())
176 warpingControl = afwMath.WarpingControl(
177 "bilinear", "", 0, interpLength)
179 with self.assertRaises(pexExcept.InvalidParameterError):
180 afwMath.warpExposure(exposureWithWcs, exposureWithoutWcs, warpingControl)
182 with self.assertRaises(pexExcept.InvalidParameterError):
183 afwMath.warpExposure(exposureWithoutWcs, exposureWithWcs, warpingControl)
185 def testWarpIntoSelf(self, interpLength=10):
186 """Cannot warp in-place
187 """
188 wcs = afwGeom.makeSkyWcs(
189 crpix=lsst.geom.Point2D(0, 0),
190 crval=lsst.geom.SpherePoint(359, 0, lsst.geom.degrees),
191 cdMatrix=afwGeom.makeCdMatrix(1.0e-8*lsst.geom.degrees),
192 )
193 exposure = afwImage.ExposureF(lsst.geom.Extent2I(100, 100), wcs)
194 maskedImage = exposure.getMaskedImage()
195 warpingControl = afwMath.WarpingControl(
196 "bilinear", "", 0, interpLength)
198 with self.assertRaises(pexExcept.InvalidParameterError):
199 afwMath.warpExposure(exposure, exposure, warpingControl)
201 with self.assertRaises(pexExcept.InvalidParameterError):
202 afwMath.warpImage(maskedImage, wcs, maskedImage, wcs, warpingControl)
204 with self.assertRaises(pexExcept.InvalidParameterError):
205 afwMath.warpImage(maskedImage.getImage(), wcs, maskedImage.getImage(), wcs, warpingControl)
207 def testWarpingControl(self):
208 """Test the basic mechanics of WarpingControl
209 """
210 for interpLength in (0, 1, 52):
211 wc = afwMath.WarpingControl("lanczos3", "", 0, interpLength)
212 self.assertFalse(wc.hasMaskWarpingKernel())
213 self.assertEqual(wc.getInterpLength(), interpLength)
214 for newInterpLength in (3, 7, 9):
215 wc.setInterpLength(newInterpLength)
216 self.assertEqual(wc.getInterpLength(), newInterpLength)
218 for cacheSize in (0, 100):
219 wc = afwMath.WarpingControl("lanczos3", "bilinear", cacheSize)
220 self.assertTrue(wc.hasMaskWarpingKernel())
221 self.assertEqual(wc.getCacheSize(), cacheSize)
222 self.assertEqual(wc.getWarpingKernel().getCacheSize(), cacheSize)
223 self.assertEqual(
224 wc.getMaskWarpingKernel().getCacheSize(), cacheSize)
225 for newCacheSize in (1, 50):
226 wc.setCacheSize(newCacheSize)
227 self.assertEqual(wc.getCacheSize(), newCacheSize)
228 self.assertEqual(
229 wc.getWarpingKernel().getCacheSize(), newCacheSize)
230 self.assertEqual(
231 wc.getMaskWarpingKernel().getCacheSize(), newCacheSize)
233 wc = afwMath.WarpingControl("lanczos4", "nearest", 64, 7, 42)
234 self.assertTrue(wc.isPersistable())
235 with lsst.utils.tests.getTempFilePath(".fits", expectOutput=True) as tempFile:
236 wc.writeFits(tempFile)
237 wc2 = afwMath.WarpingControl.readFits(tempFile)
238 self.assertEqual(wc.getCacheSize(), wc2.getCacheSize())
239 self.assertEqual(wc.getInterpLength(), wc2.getInterpLength())
240 self.assertEqual(wc.getWarpingKernel().getBBox(), wc2.getWarpingKernel().getBBox())
241 self.assertEqual(wc.getWarpingKernel().getKernelParameters(),
242 wc2.getWarpingKernel().getKernelParameters())
243 self.assertEqual(wc.hasMaskWarpingKernel(), wc2.hasMaskWarpingKernel())
244 self.assertEqual(wc.getMaskWarpingKernel().getBBox(), wc2.getMaskWarpingKernel().getBBox())
245 self.assertEqual(wc.getMaskWarpingKernel().getKernelParameters(),
246 wc2.getMaskWarpingKernel().getKernelParameters())
247 self.assertEqual(wc.getGrowFullMask(), wc2.getGrowFullMask())
249 def testWarpingControlError(self):
250 """Test error handling of WarpingControl
251 """
252 # error: mask kernel smaller than main kernel
253 for kernelName, maskKernelName in (
254 ("bilinear", "lanczos3"),
255 ("bilinear", "lanczos4"),
256 ("lanczos3", "lanczos4"),
257 ):
258 with self.assertRaises(pexExcept.InvalidParameterError):
259 afwMath.WarpingControl(kernelName, maskKernelName)
261 # error: new mask kernel larger than main kernel
262 warpingControl = afwMath.WarpingControl("bilinear")
263 for maskKernelName in ("lanczos3", "lanczos4"):
264 with self.assertRaises(pexExcept.InvalidParameterError):
265 warpingControl.setMaskWarpingKernelName(maskKernelName)
267 # error: new kernel smaller than mask kernel
268 warpingControl = afwMath.WarpingControl("lanczos4", "lanczos4")
269 for kernelName in ("bilinear", "lanczos3"):
270 with self.assertRaises(pexExcept.InvalidParameterError):
271 warpingControl.setWarpingKernelName(kernelName)
273 # OK: main kernel at least as big as mask kernel
274 for kernelName, maskKernelName in (
275 ("bilinear", "bilinear"),
276 ("lanczos3", "lanczos3"),
277 ("lanczos3", "bilinear"),
278 ("lanczos4", "lanczos3"),
279 ):
280 # this should not raise any exception
281 afwMath.WarpingControl(kernelName, maskKernelName)
283 # invalid kernel names
284 for kernelName, maskKernelName in (
285 ("badname", ""),
286 ("lanczos", ""), # no digit after lanczos
287 ("lanczos3", "badname"),
288 ("lanczos3", "lanczos"),
289 ):
290 with self.assertRaises(pexExcept.InvalidParameterError):
291 afwMath.WarpingControl(kernelName, maskKernelName)
293 def testWarpMask(self):
294 """Test that warping the mask plane with a different kernel does the right thing
295 """
296 for kernelName, maskKernelName in (
297 ("bilinear", "bilinear"),
298 ("lanczos3", "lanczos3"),
299 ("lanczos3", "bilinear"),
300 ("lanczos4", "lanczos3"),
301 ):
302 for growFullMask in (0, 1, 3, 0xFFFF):
303 self.verifyMaskWarp(
304 kernelName=kernelName,
305 maskKernelName=maskKernelName,
306 growFullMask=growFullMask,
307 )
309 def testMatchSwarpBilinearImage(self):
310 """Test that warpExposure matches swarp using a bilinear warping kernel
311 """
312 self.compareToSwarp("bilinear", useWarpExposure=False, atol=0.15)
314 def testMatchSwarpBilinearExposure(self):
315 """Test that warpExposure matches swarp using a bilinear warping kernel
316 """
317 self.compareToSwarp("bilinear", useWarpExposure=True,
318 useSubregion=False, useDeepCopy=True)
320 def testMatchSwarpLanczos2Image(self):
321 """Test that warpExposure matches swarp using a lanczos2 warping kernel
322 """
323 self.compareToSwarp("lanczos2", useWarpExposure=False)
325 def testMatchSwarpLanczos2Exposure(self):
326 """Test that warpExposure matches swarp using a lanczos2 warping kernel.
327 """
328 self.compareToSwarp("lanczos2", useWarpExposure=True)
330 def testMatchSwarpLanczos2SubExposure(self):
331 """Test that warpExposure matches swarp using a lanczos2 warping kernel with a subexposure
332 """
333 for useDeepCopy in (False, True):
334 self.compareToSwarp("lanczos2", useWarpExposure=True,
335 useSubregion=True, useDeepCopy=useDeepCopy)
337 def testMatchSwarpLanczos3Image(self):
338 """Test that warpExposure matches swarp using a lanczos2 warping kernel
339 """
340 self.compareToSwarp("lanczos3", useWarpExposure=False)
342 def testMatchSwarpLanczos3(self):
343 """Test that warpExposure matches swarp using a lanczos4 warping kernel.
344 """
345 self.compareToSwarp("lanczos3", useWarpExposure=True)
347 def testMatchSwarpLanczos4Image(self):
348 """Test that warpExposure matches swarp using a lanczos2 warping kernel
349 """
350 self.compareToSwarp("lanczos4", useWarpExposure=False)
352 def testMatchSwarpLanczos4(self):
353 """Test that warpExposure matches swarp using a lanczos4 warping kernel.
354 """
355 self.compareToSwarp("lanczos4", useWarpExposure=True)
357 def testMatchSwarpNearestExposure(self):
358 """Test that warpExposure matches swarp using a nearest neighbor warping kernel
359 """
360 self.compareToSwarp("nearest", useWarpExposure=True, atol=60)
362 @unittest.skipIf(afwdataDir is None, "afwdata not setup")
363 def testTransformBasedWarp(self):
364 """Test warping using TransformPoint2ToPoint2
365 """
366 for interpLength in (0, 1, 2, 4):
367 kernelName = "lanczos3"
368 rtol = 4e-5
369 atol = 1e-2
370 warpingControl = afwMath.WarpingControl(
371 warpingKernelName=kernelName,
372 interpLength=interpLength,
373 )
375 originalExposure = afwImage.ExposureF(originalExposurePath)
376 originalMetadata = afwImage.DecoratedImageF(originalExposurePath).getMetadata()
377 originalSkyWcs = afwGeom.makeSkyWcs(originalMetadata)
379 swarpedImageName = f"medswarp1{kernelName}.fits"
380 swarpedImagePath = os.path.join(dataDir, swarpedImageName)
381 swarpedDecoratedImage = afwImage.DecoratedImageF(swarpedImagePath)
382 swarpedImage = swarpedDecoratedImage.getImage()
384 swarpedMetadata = swarpedDecoratedImage.getMetadata()
385 warpedSkyWcs = afwGeom.makeSkyWcs(swarpedMetadata)
387 # original image is source, warped image is destination
388 srcToDest = afwGeom.makeWcsPairTransform(originalSkyWcs, warpedSkyWcs)
390 afwWarpedMaskedImage = afwImage.MaskedImageF(swarpedImage.getDimensions())
391 originalMaskedImage = originalExposure.getMaskedImage()
393 numGoodPix = afwMath.warpImage(afwWarpedMaskedImage, originalMaskedImage,
394 srcToDest, warpingControl)
395 self.assertGreater(numGoodPix, 50)
397 afwWarpedImage = afwWarpedMaskedImage.getImage()
398 afwWarpedImageArr = afwWarpedImage.getArray()
399 noDataMaskArr = np.isnan(afwWarpedImageArr)
400 self.assertImagesAlmostEqual(afwWarpedImage, swarpedImage,
401 skipMask=noDataMaskArr, rtol=rtol, atol=atol)
403 def testTicket2441(self):
404 """Test ticket 2441: warpExposure sometimes mishandles zero-extent dest exposures"""
405 fromWcs = afwGeom.makeSkyWcs(
406 crpix=lsst.geom.Point2D(0, 0),
407 crval=lsst.geom.SpherePoint(359, 0, lsst.geom.degrees),
408 cdMatrix=afwGeom.makeCdMatrix(scale=1.0e-8*lsst.geom.degrees),
409 )
410 fromExp = afwImage.ExposureF(afwImage.MaskedImageF(10, 10), fromWcs)
412 toWcs = afwGeom.makeSkyWcs(
413 crpix=lsst.geom.Point2D(410000, 11441),
414 crval=lsst.geom.SpherePoint(45, 0, lsst.geom.degrees),
415 cdMatrix=afwGeom.makeCdMatrix(scale=0.00011*lsst.geom.degrees, flipX=True),
416 projection="CEA",
417 )
418 toExp = afwImage.ExposureF(afwImage.MaskedImageF(0, 0), toWcs)
420 warpControl = afwMath.WarpingControl("lanczos3")
421 # if a bug described in ticket #2441 is present, this will raise an
422 # exception:
423 numGoodPix = afwMath.warpExposure(toExp, fromExp, warpControl)
424 self.assertEqual(numGoodPix, 0)
426 def testTicketDM4063(self):
427 """Test that a uint16 array can be cast to a bool array, to avoid DM-4063
428 """
429 a = np.array([0, 1, 0, 23], dtype=np.uint16)
430 b = np.array([True, True, False, False], dtype=bool)
431 acast = np.array(a != 0, dtype=bool)
432 orArr = acast | b
433 desOrArr = np.array([True, True, False, True], dtype=bool)
434 # Note: assertEqual(bool arr, bool arr) fails with:
435 # ValueError: The truth value of an array with more than one element is
436 # ambiguous
437 try:
438 self.assertTrue(np.all(orArr == desOrArr))
439 except Exception as e:
440 print(f"Failed: {orArr!r} != {desOrArr!r}: {e}")
441 raise
443 def testSmallSrc(self):
444 """Verify that a source image that is too small will not raise an exception
446 This tests another bug that was fixed in ticket #2441
447 """
448 fromWcs = afwGeom.makeSkyWcs(
449 crpix=lsst.geom.Point2D(0, 0),
450 crval=lsst.geom.SpherePoint(359, 0, lsst.geom.degrees),
451 cdMatrix=afwGeom.makeCdMatrix(scale=1.0e-8*lsst.geom.degrees),
452 )
453 fromExp = afwImage.ExposureF(afwImage.MaskedImageF(1, 1), fromWcs)
455 toWcs = afwGeom.makeSkyWcs(
456 crpix=lsst.geom.Point2D(0, 0),
457 crval=lsst.geom.SpherePoint(358, 0, lsst.geom.degrees),
458 cdMatrix=afwGeom.makeCdMatrix(scale=1.1e-8*lsst.geom.degrees),
459 )
460 toExp = afwImage.ExposureF(afwImage.MaskedImageF(10, 10), toWcs)
462 warpControl = afwMath.WarpingControl("lanczos3")
463 # if a bug described in ticket #2441 is present, this will raise an
464 # exception:
465 numGoodPix = afwMath.warpExposure(toExp, fromExp, warpControl)
466 self.assertEqual(numGoodPix, 0)
467 self.assertTrue(np.all(np.isnan(toExp.image.array)))
468 self.assertTrue(np.all(np.isinf(toExp.variance.array)))
469 noDataBitMask = afwImage.Mask.getPlaneBitMask("NO_DATA")
470 self.assertTrue(np.all(toExp.mask.array == noDataBitMask))
472 def verifyMaskWarp(self, kernelName, maskKernelName, growFullMask, interpLength=10, cacheSize=100000,
473 rtol=4e-05, atol=1e-2):
474 """Verify that using a separate mask warping kernel produces the correct results
476 Inputs:
477 - kernelName: name of warping kernel in the form used by afwImage.makeKernel
478 - maskKernelName: name of mask warping kernel in the form used by afwImage.makeKernel
479 - interpLength: interpLength argument for lsst.afw.math.WarpingControl
480 - cacheSize: cacheSize argument for lsst.afw.math.WarpingControl;
481 0 disables the cache
482 10000 gives some speed improvement but less accurate results (atol must be increased)
483 100000 gives better accuracy but no speed improvement in this test
484 - rtol: relative tolerance as used by np.allclose
485 - atol: absolute tolerance as used by np.allclose
486 """
487 srcWcs = afwGeom.makeSkyWcs(
488 crpix=lsst.geom.Point2D(10, 11),
489 crval=lsst.geom.SpherePoint(41.7, 32.9, lsst.geom.degrees),
490 cdMatrix=afwGeom.makeCdMatrix(scale=0.2*lsst.geom.degrees),
491 )
492 destWcs = afwGeom.makeSkyWcs(
493 crpix=lsst.geom.Point2D(9, 10),
494 crval=lsst.geom.SpherePoint(41.65, 32.95, lsst.geom.degrees),
495 cdMatrix=afwGeom.makeCdMatrix(scale=0.17*lsst.geom.degrees),
496 )
498 srcMaskedImage = afwImage.MaskedImageF(100, 101)
499 srcExposure = afwImage.ExposureF(srcMaskedImage, srcWcs)
501 shape = srcMaskedImage.image.array.shape
502 srcMaskedImage.image.array[:] = np.random.normal(10000, 1000, size=shape)
503 srcMaskedImage.variance.array[:] = np.random.normal(9000, 900, size=shape)
504 srcMaskedImage.mask.array[:] = np.reshape(
505 np.arange(0, shape[0] * shape[1], 1, dtype=np.uint16), shape)
507 warpControl = afwMath.WarpingControl(
508 kernelName,
509 maskKernelName,
510 cacheSize,
511 interpLength,
512 growFullMask
513 )
514 destMaskedImage = afwImage.MaskedImageF(110, 121)
515 destExposure = afwImage.ExposureF(destMaskedImage, destWcs)
516 afwMath.warpExposure(destExposure, srcExposure, warpControl)
518 # now compute with two separate mask planes
519 warpControl.setGrowFullMask(0)
520 narrowMaskedImage = afwImage.MaskedImageF(110, 121)
521 narrowExposure = afwImage.ExposureF(narrowMaskedImage, destWcs)
522 afwMath.warpExposure(narrowExposure, srcExposure, warpControl)
524 warpControl.setMaskWarpingKernelName("")
525 broadMaskedImage = afwImage.MaskedImageF(110, 121)
526 broadExposure = afwImage.ExposureF(broadMaskedImage, destWcs)
527 afwMath.warpExposure(broadExposure, srcExposure, warpControl)
529 if (kernelName != maskKernelName) and (growFullMask != 0xFFFF):
530 # we expect the mask planes to differ
531 if np.all(narrowExposure.mask.array == broadExposure.mask.array):
532 self.fail("No difference between broad and narrow mask")
534 predMask = (broadExposure.mask.array & growFullMask) | (
535 narrowExposure.mask.array & ~growFullMask).astype(np.uint16)
536 predArraySet = (broadExposure.image.array, predMask, broadExposure.variance.array)
537 predExposure = afwImage.makeMaskedImageFromArrays(*predArraySet)
539 msg = f"Separate mask warping failed; warpingKernel={kernelName}; maskWarpingKernel={maskKernelName}"
540 self.assertMaskedImagesAlmostEqual(destExposure.getMaskedImage(), predExposure,
541 doImage=True, doMask=True, doVariance=True,
542 rtol=rtol, atol=atol, msg=msg)
544 @unittest.skipIf(afwdataDir is None, "afwdata not setup")
545 def compareToSwarp(self, kernelName,
546 useWarpExposure=True, useSubregion=False, useDeepCopy=False,
547 interpLength=10, cacheSize=100000,
548 rtol=4e-05, atol=1e-2):
549 """Compare warpExposure to swarp for given warping kernel.
551 Note that swarp only warps the image plane, so only test that plane.
553 Inputs:
554 - kernelName: name of kernel in the form used by afwImage.makeKernel
555 - useWarpExposure: if True, call warpExposure to warp an ExposureF,
556 else call warpImage to warp an ImageF and also call the Transform version
557 - useSubregion: if True then the original source exposure (from which the usual
558 test exposure was extracted) is read and the correct subregion extracted
559 - useDeepCopy: if True then the copy of the subimage is a deep copy,
560 else it is a shallow copy; ignored if useSubregion is False
561 - interpLength: interpLength argument for lsst.afw.math.WarpingControl
562 - cacheSize: cacheSize argument for lsst.afw.math.WarpingControl;
563 0 disables the cache
564 10000 gives some speed improvement but less accurate results (atol must be increased)
565 100000 gives better accuracy but no speed improvement in this test
566 - rtol: relative tolerance as used by np.allclose
567 - atol: absolute tolerance as used by np.allclose
568 """
569 warpingControl = afwMath.WarpingControl(
570 kernelName,
571 "", # there is no point to a separate mask kernel since we aren't testing the mask plane
572 cacheSize,
573 interpLength,
574 )
575 if useSubregion:
576 originalFullExposure = afwImage.ExposureF(originalExposurePath)
577 # "medsub" is a subregion of med starting at 0-indexed pixel (40, 150) of size 145 x 200
578 bbox = lsst.geom.Box2I(lsst.geom.Point2I(40, 150),
579 lsst.geom.Extent2I(145, 200))
580 originalExposure = afwImage.ExposureF(
581 originalFullExposure, bbox, afwImage.LOCAL, useDeepCopy)
582 swarpedImageName = f"medsubswarp1{kernelName}.fits"
583 else:
584 originalExposure = afwImage.ExposureF(originalExposurePath)
585 swarpedImageName = f"medswarp1{kernelName}.fits"
587 swarpedImagePath = os.path.join(dataDir, swarpedImageName)
588 swarpedDecoratedImage = afwImage.DecoratedImageF(swarpedImagePath)
589 swarpedImage = swarpedDecoratedImage.getImage()
590 swarpedMetadata = swarpedDecoratedImage.getMetadata()
591 warpedWcs = afwGeom.makeSkyWcs(swarpedMetadata)
593 if useWarpExposure:
594 # path for saved afw-warped image
595 afwWarpedImagePath = f"afwWarpedExposure1{kernelName}.fits"
597 afwWarpedMaskedImage = afwImage.MaskedImageF(
598 swarpedImage.getDimensions())
599 afwWarpedExposure = afwImage.ExposureF(
600 afwWarpedMaskedImage, warpedWcs)
601 afwMath.warpExposure(
602 afwWarpedExposure, originalExposure, warpingControl)
603 afwWarpedMask = afwWarpedMaskedImage.getMask()
604 if SAVE_FITS_FILES:
605 afwWarpedExposure.writeFits(afwWarpedImagePath)
606 if display:
607 afwDisplay.Display(frame=1).mtv(afwWarpedExposure, title="Warped")
609 swarpedMaskedImage = afwImage.MaskedImageF(swarpedImage)
611 if display:
612 afwDisplay.Display(frame=2).mtv(swarpedMaskedImage, title="SWarped")
614 msg = f"afw and swarp {kernelName}-warped differ (ignoring bad pixels)"
615 try:
616 self.assertMaskedImagesAlmostEqual(afwWarpedMaskedImage, swarpedMaskedImage,
617 doImage=True, doMask=False, doVariance=False,
618 skipMask=afwWarpedMask, rtol=rtol, atol=atol, msg=msg)
619 except Exception:
620 if SAVE_FAILED_FITS_FILES:
621 afwWarpedExposure.writeFits(afwWarpedImagePath)
622 print(f"Saved failed afw-warped exposure as: {afwWarpedImagePath}")
623 raise
624 else:
625 # path for saved afw-warped image
626 afwWarpedImagePath = f"afwWarpedImage1{kernelName}.fits"
627 afwWarpedImage2Path = f"afwWarpedImage1{kernelName}_xyTransform.fits"
629 afwWarpedImage = afwImage.ImageF(swarpedImage.getDimensions())
630 originalImage = originalExposure.getMaskedImage().getImage()
631 originalWcs = originalExposure.getWcs()
632 afwMath.warpImage(afwWarpedImage, warpedWcs, originalImage,
633 originalWcs, warpingControl)
634 if display:
635 afwDisplay.Display(frame=1).mtv(afwWarpedImage, title="Warped")
636 afwDisplay.Display(frame=2).mtv(swarpedImage, title="SWarped")
637 diff = swarpedImage.Factory(swarpedImage, True)
638 diff -= afwWarpedImage
639 afwDisplay.Display(frame=3).mtv(diff, title="swarp - afw")
640 if SAVE_FITS_FILES:
641 afwWarpedImage.writeFits(afwWarpedImagePath)
643 afwWarpedImageArr = afwWarpedImage.getArray()
644 noDataMaskArr = np.isnan(afwWarpedImageArr)
645 msg = f"afw and swarp {kernelName}-warped images do not match (ignoring NaN pixels)"
646 try:
647 self.assertImagesAlmostEqual(afwWarpedImage, swarpedImage,
648 skipMask=noDataMaskArr, rtol=rtol, atol=atol, msg=msg)
649 except Exception:
650 if SAVE_FAILED_FITS_FILES:
651 # save the image anyway
652 afwWarpedImage.writeFits(afwWarpedImagePath)
653 print(f"Saved failed afw-warped image as: {afwWarpedImagePath}")
654 raise
656 afwWarpedImage2 = afwImage.ImageF(swarpedImage.getDimensions())
657 srcToDest = afwGeom.makeWcsPairTransform(originalWcs, warpedWcs)
658 afwMath.warpImage(afwWarpedImage2, originalImage,
659 srcToDest, warpingControl)
660 msg = f"afw transform-based and WCS-based {kernelName}-warped images do not match"
661 try:
662 self.assertImagesAlmostEqual(afwWarpedImage2, afwWarpedImage,
663 rtol=rtol, atol=atol, msg=msg)
664 except Exception:
665 if SAVE_FAILED_FITS_FILES:
666 # save the image anyway
667 afwWarpedImage.writeFits(afwWarpedImage2)
668 print(f"Saved failed afw-warped image as: {afwWarpedImage2Path}")
669 raise
672class MemoryTester(lsst.utils.tests.MemoryTestCase):
673 pass
676def setup_module(module):
677 lsst.utils.tests.init()
680if __name__ == "__main__": 680 ↛ 681line 680 didn't jump to line 681 because the condition on line 680 was never true
681 lsst.utils.tests.init()
682 unittest.main()