Coverage for tests / test_skyWcs.py: 12%
566 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 01:18 -0700
1import itertools
2import os
3import sys
4import unittest
6import astropy.io.fits
7import astropy.coordinates
8import astropy.wcs
9import astshim as ast
10import numpy as np
11from numpy.testing import assert_allclose
13from lsst.daf.base import PropertyList
14import lsst.geom
15import lsst.afw.cameraGeom as cameraGeom
16from lsst.afw.geom import (
17 TransformPoint2ToPoint2, TransformPoint2ToSpherePoint, makeRadialTransform,
18 SkyWcs, makeSkyWcs, makeCdMatrix, makeWcsPairTransform,
19 makeFlippedWcs, makeModifiedWcs, makeTanSipWcs,
20 getIntermediateWorldCoordsToSky, getPixelToIntermediateWorldCoords,
21 stripWcsMetadata
22)
23from lsst.afw.geom import getCdMatrixFromMetadata, getSipMatrixFromMetadata, makeSimpleWcsMetadata
24from lsst.afw.geom.testUtils import makeSipIwcToPixel, makeSipPixelToIwc
25from lsst.afw.fits import makeLimitedFitsHeader
26from lsst.afw.image import ExposureF
29def addActualPixelsFrame(skyWcs, actualPixelsToPixels):
30 """Add an "ACTUAL_PIXELS" frame to a SkyWcs and return the result
32 Parameters
33 ----------
34 skyWcs : `lsst.afw.geom.SkyWcs`
35 The WCS to which you wish to add an ACTUAL_PIXELS frame
36 actualPixelsToPixels : `lsst.afw.geom.TransformPoint2ToPoint2`
37 The transform from ACTUAL_PIXELS to PIXELS
38 """
39 actualPixelsToPixelsMap = actualPixelsToPixels.getMapping()
40 actualPixelsFrame = ast.Frame(2, "Domain=ACTUAL_PIXELS")
41 frameDict = skyWcs.getFrameDict()
42 frameDict.addFrame("PIXELS", actualPixelsToPixelsMap.inverted(), actualPixelsFrame)
43 frameDict.setBase("ACTUAL_PIXELS")
44 frameDict.setCurrent("SKY")
45 return SkyWcs(frameDict)
48class SkyWcsBaseTestCase(lsst.utils.tests.TestCase):
49 def checkPersistence(self, skyWcs, bbox):
50 """Check persistence of a SkyWcs
51 """
52 className = "SkyWcs"
54 # check writeString and readString
55 skyWcsStr = skyWcs.writeString()
56 serialVersion, serialClassName, serialRest = skyWcsStr.split(" ", 2)
57 self.assertEqual(int(serialVersion), 1)
58 self.assertEqual(serialClassName, className)
59 badStr1 = " ".join(["2", serialClassName, serialRest])
60 with self.assertRaises(lsst.pex.exceptions.TypeError):
61 skyWcs.readString(badStr1)
62 badClassName = "x" + serialClassName
63 badStr2 = " ".join(["1", badClassName, serialRest])
64 with self.assertRaises(lsst.pex.exceptions.TypeError):
65 skyWcs.readString(badStr2)
66 skyWcsFromStr1 = skyWcs.readString(skyWcsStr)
67 self.assertEqual(skyWcs, skyWcsFromStr1)
68 self.assertEqual(type(skyWcs), type(skyWcsFromStr1))
69 self.assertEqual(skyWcs.getFrameDict(), skyWcsFromStr1.getFrameDict())
71 pixelPoints = [
72 lsst.geom.Point2D(0, 0),
73 lsst.geom.Point2D(1000, 0),
74 lsst.geom.Point2D(0, 1000),
75 lsst.geom.Point2D(-50, -50),
76 ]
77 skyPoints = skyWcs.pixelToSky(pixelPoints)
78 pixelPoints2 = skyWcs.skyToPixel(skyPoints)
79 assert_allclose(pixelPoints, pixelPoints2, atol=1e-7)
81 # check that WCS is properly saved as part of an exposure FITS file
82 exposure = ExposureF(100, 100, skyWcs)
83 with lsst.utils.tests.getTempFilePath(".fits") as outFile:
84 exposure.writeFits(outFile)
85 exposureRoundTrip = ExposureF(outFile)
86 wcsFromExposure = exposureRoundTrip.getWcs()
87 self.assertWcsAlmostEqualOverBBox(skyWcs, wcsFromExposure, bbox, maxDiffPix=0,
88 maxDiffSky=0*lsst.geom.radians)
90 def checkFrameDictConstructor(self, skyWcs, bbox):
91 """Check that the FrameDict constructor works
92 """
93 frameDict = skyWcs.getFrameDict()
94 wcsFromFrameDict = SkyWcs(frameDict)
95 self.assertWcsAlmostEqualOverBBox(skyWcs, wcsFromFrameDict, bbox, maxDiffPix=0,
96 maxDiffSky=0*lsst.geom.radians)
98 self.checkPersistence(wcsFromFrameDict, bbox)
100 # check that it is impossible to build a SkyWcs if a required frame is missing
101 for domain in ("PIXELS", "IWC", "SKY"):
102 badFrameDict = skyWcs.getFrameDict()
103 badFrameDict.removeFrame(domain)
104 with self.assertRaises(lsst.pex.exceptions.TypeError):
105 SkyWcs(badFrameDict)
107 def checkMakeFlippedWcs(self, skyWcs, skyAtol=1e-7*lsst.geom.arcseconds):
108 """Check makeFlippedWcs on the provided WCS
109 """
110 # make an arbitrary bbox, but one that includes zero in one axis
111 # and does not include zero in the other axis
112 # the center of the bbox is used as the center of flipping
113 # and the corners of the bbox are the input positions that are tested
114 bbox = lsst.geom.Box2D(lsst.geom.Point2D(-100, 1000), lsst.geom.Extent2D(2000, 1501))
115 # dict of (isRight, isTop): position
116 minPos = bbox.getMin()
117 maxPos = bbox.getMax()
118 center = bbox.getCenter()
119 cornerDict = {
120 (False, False): minPos,
121 (False, True): lsst.geom.Point2D(minPos[0], maxPos[1]),
122 (True, False): lsst.geom.Point2D(maxPos[0], minPos[1]),
123 (True, True): maxPos,
124 }
125 for flipLR, flipTB in itertools.product((False, True), (False, True)):
126 flippedWcs = makeFlippedWcs(wcs=skyWcs, flipLR=flipLR, flipTB=flipTB, center=center)
127 # the center is unchanged
128 self.assertSpherePointsAlmostEqual(skyWcs.pixelToSky(center),
129 flippedWcs.pixelToSky(center), maxSep=skyAtol)
131 for isR, isT in itertools.product((False, True), (False, True)):
132 origPos = cornerDict[(isR, isT)]
133 flippedPos = cornerDict[(isR ^ flipLR, isT ^ flipTB)]
134 self.assertSpherePointsAlmostEqual(skyWcs.pixelToSky(origPos),
135 flippedWcs.pixelToSky(flippedPos), maxSep=skyAtol)
137 def assertSkyWcsAstropyWcsAlmostEqual(self, skyWcs, astropyWcs, bbox,
138 pixAtol=1e-4, skyAtol=1e-4*lsst.geom.arcseconds,
139 checkRoundTrip=True):
140 """Assert that a SkyWcs and the corresponding astropy.wcs.WCS agree over a specified bounding box
141 """
142 bbox = lsst.geom.Box2D(bbox)
143 center = bbox.getCenter()
144 xArr = bbox.getMinX(), center[0], bbox.getMaxX()
145 yArr = bbox.getMinY(), center[1], bbox.getMaxY()
146 pixPosList = [lsst.geom.Point2D(x, y) for x, y in itertools.product(xArr, yArr)]
148 # pixelToSky
149 skyPosList = skyWcs.pixelToSky(pixPosList)
150 astropySkyPosList = self.astropyPixelsToSky(astropyWcs=astropyWcs, pixPosList=pixPosList)
151 self.assertSpherePointListsAlmostEqual(skyPosList, astropySkyPosList, maxSep=skyAtol)
153 if not checkRoundTrip:
154 return
156 # astropy round trip
157 astropyPixPosRoundTrip = self.astropySkyToPixels(astropyWcs=astropyWcs, skyPosList=astropySkyPosList)
158 self.assertPairListsAlmostEqual(pixPosList, astropyPixPosRoundTrip, maxDiff=pixAtol)
160 # SkyWcs round trip
161 pixPosListRoundTrip = skyWcs.skyToPixel(skyPosList)
162 self.assertPairListsAlmostEqual(pixPosList, pixPosListRoundTrip, maxDiff=pixAtol)
164 # skyToPixel astropy vs SkyWcs
165 astropyPixPosList2 = self.astropySkyToPixels(astropyWcs=astropyWcs, skyPosList=skyPosList)
166 self.assertPairListsAlmostEqual(pixPosListRoundTrip, astropyPixPosList2, maxDiff=pixAtol)
168 def astropyPixelsToSky(self, astropyWcs, pixPosList):
169 """Use an astropy wcs to convert pixels to sky
171 @param[in] astropyWcs a celestial astropy.wcs.WCS with 2 axes in RA, Dec order
172 @param[in] pixPosList 0-based pixel positions as lsst.geom.Point2D or similar pairs
173 @returns sky coordinates as a list of lsst.geom.SpherePoint
175 Converts the output to ICRS
176 """
177 xarr = [p[0] for p in pixPosList]
178 yarr = [p[1] for p in pixPosList]
179 skyCoordList = astropy.wcs.utils.pixel_to_skycoord(xp=xarr,
180 yp=yarr,
181 wcs=astropyWcs,
182 origin=0,
183 mode="all")
184 icrsList = [sc.transform_to("icrs") for sc in skyCoordList]
185 return [lsst.geom.SpherePoint(sc.ra.deg, sc.dec.deg, lsst.geom.degrees) for sc in icrsList]
187 def astropySkyToPixels(self, astropyWcs, skyPosList):
188 """Use an astropy wcs to convert pixels to sky
190 @param[in] astropyWcs a celestial astropy.wcs.WCS with 2 axes in RA, Dec order
191 @param[in] skyPosList ICRS sky coordinates as a list of lsst.geom.SpherePoint
192 @returns a list of lsst.geom.Point2D, 0-based pixel positions
194 Converts the input from ICRS to the coordinate system of the wcs
195 """
196 skyCoordList = [astropy.coordinates.SkyCoord(c[0].asDegrees(),
197 c[1].asDegrees(),
198 frame="icrs",
199 unit="deg") for c in skyPosList]
200 xyArr = [astropy.wcs.utils.skycoord_to_pixel(coords=sc,
201 wcs=astropyWcs,
202 origin=0,
203 mode="all") for sc in skyCoordList]
204 # float is needed to avoid truncation to int
205 return [lsst.geom.Point2D(float(x), float(y)) for x, y in xyArr]
208class SimpleSkyWcsTestCase(SkyWcsBaseTestCase):
209 """Test the simple FITS version of makeSkyWcs
210 """
212 def setUp(self):
213 self.crpix = lsst.geom.Point2D(100, 100)
214 self.crvalList = [
215 lsst.geom.SpherePoint(0, 45, lsst.geom.degrees),
216 lsst.geom.SpherePoint(0.00001, 45, lsst.geom.degrees),
217 lsst.geom.SpherePoint(359.99999, 45, lsst.geom.degrees),
218 lsst.geom.SpherePoint(30, 89.99999, lsst.geom.degrees),
219 lsst.geom.SpherePoint(30, -89.99999, lsst.geom.degrees),
220 ]
221 self.orientationList = [
222 0 * lsst.geom.degrees,
223 0.00001 * lsst.geom.degrees,
224 -0.00001 * lsst.geom.degrees,
225 -45 * lsst.geom.degrees,
226 90 * lsst.geom.degrees,
227 ]
228 self.scale = 1.0 * lsst.geom.arcseconds
229 self.tinyPixels = 1.0e-10
230 self.tinyAngle = 1.0e-10 * lsst.geom.radians
231 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(-1000, -1000),
232 lsst.geom.Extent2I(2000, 2000)) # arbitrary but reasonable
234 def checkTanWcs(self, crval, orientation, flipX):
235 """Construct a pure TAN SkyWcs and check that it operates as specified
237 Parameters
238 ----------
239 crval : `lsst.geom.SpherePoint`
240 Desired reference sky position.
241 Must not be at either pole.
242 orientation : `lsst.geom.Angle`
243 Position angle of pixel +Y, measured from N through E.
244 At 0 degrees, +Y is along N and +X is along E/W if flipX false/true
245 At 90 degrees, +Y is along E and +X is along S/N if flipX false/true
246 flipX : `bool`
247 Flip x axis? See `orientation` for details.
249 Returns
250 -------
251 wcs : `lsst.afw.geom.SkyWcs`
252 The generated pure TAN SkyWcs
253 """
254 cdMatrix = makeCdMatrix(scale=self.scale, orientation=orientation, flipX=flipX)
255 wcs = makeSkyWcs(crpix=self.crpix, crval=crval, cdMatrix=cdMatrix)
256 self.checkPersistence(wcs, bbox=self.bbox)
257 self.checkMakeFlippedWcs(wcs)
259 self.assertTrue(wcs.isFits)
260 self.assertEqual(wcs.isFlipped, bool(flipX))
262 xoffAng = 0*lsst.geom.degrees if flipX else 180*lsst.geom.degrees
264 pixelList = [
265 lsst.geom.Point2D(self.crpix[0], self.crpix[1]),
266 lsst.geom.Point2D(self.crpix[0] + 1, self.crpix[1]),
267 lsst.geom.Point2D(self.crpix[0], self.crpix[1] + 1),
268 ]
269 skyList = wcs.pixelToSky(pixelList)
271 # check pixels to sky
272 predSkyList = [
273 crval,
274 crval.offset(xoffAng - orientation, self.scale),
275 crval.offset(90*lsst.geom.degrees - orientation, self.scale),
276 ]
277 self.assertSpherePointListsAlmostEqual(predSkyList, skyList)
278 self.assertSpherePointListsAlmostEqual(predSkyList, wcs.pixelToSky(pixelList))
279 for pixel, predSky in zip(pixelList, predSkyList):
280 self.assertSpherePointsAlmostEqual(predSky, wcs.pixelToSky(pixel))
281 self.assertSpherePointsAlmostEqual(predSky, wcs.pixelToSky(pixel[0], pixel[1]))
283 # check sky to pixels
284 self.assertPairListsAlmostEqual(pixelList, wcs.skyToPixel(skyList))
285 self.assertPairListsAlmostEqual(pixelList, wcs.skyToPixel(skyList))
286 for pixel, sky in zip(pixelList, skyList):
287 self.assertPairsAlmostEqual(pixel, wcs.skyToPixel(sky))
288 # self.assertPairsAlmostEqual(pixel, wcs.skyToPixel(sky[0], sky[1]))
290 # check CRVAL round trip
291 self.assertSpherePointsAlmostEqual(wcs.getSkyOrigin(), crval,
292 maxSep=self.tinyAngle)
294 crpix = wcs.getPixelOrigin()
295 self.assertPairsAlmostEqual(crpix, self.crpix, maxDiff=self.tinyPixels)
297 self.assertFloatsAlmostEqual(wcs.getCdMatrix(), cdMatrix, atol=1e-15, rtol=1e-11)
299 pixelScale = wcs.getPixelScale()
300 self.assertAnglesAlmostEqual(self.scale, pixelScale, maxDiff=self.tinyAngle)
302 pixelScale = wcs.getPixelScale(self.crpix)
303 self.assertAnglesAlmostEqual(self.scale, pixelScale, maxDiff=self.tinyAngle)
305 # check that getFitsMetadata can operate at high precision
306 # and has axis order RA, Dec
307 fitsMetadata = wcs.getFitsMetadata(True)
308 self.assertEqual(fitsMetadata.getScalar("CTYPE1")[0:4], "RA--")
309 self.assertEqual(fitsMetadata.getScalar("CTYPE2")[0:4], "DEC-")
311 # Compute a WCS with the pixel origin shifted by an arbitrary amount
312 # The resulting sky origin should not change
313 offset = lsst.geom.Extent2D(500, -322) # arbitrary
314 shiftedWcs = wcs.copyAtShiftedPixelOrigin(offset)
315 self.assertTrue(shiftedWcs.isFits)
316 predShiftedPixelOrigin = self.crpix + offset
317 self.assertPairsAlmostEqual(shiftedWcs.getPixelOrigin(), predShiftedPixelOrigin,
318 maxDiff=self.tinyPixels)
319 self.assertSpherePointsAlmostEqual(shiftedWcs.getSkyOrigin(), crval, maxSep=self.tinyAngle)
321 shiftedPixelList = [p + offset for p in pixelList]
322 shiftedSkyList = shiftedWcs.pixelToSky(shiftedPixelList)
323 self.assertSpherePointListsAlmostEqual(skyList, shiftedSkyList, maxSep=self.tinyAngle)
325 # Check that the shifted WCS can be round tripped as FITS metadata
326 shiftedMetadata = shiftedWcs.getFitsMetadata(precise=True)
327 shiftedWcsCopy = makeSkyWcs(shiftedMetadata)
328 shiftedBBox = lsst.geom.Box2D(predShiftedPixelOrigin,
329 predShiftedPixelOrigin + lsst.geom.Extent2I(2000, 2000))
330 self.assertWcsAlmostEqualOverBBox(shiftedWcs, shiftedWcsCopy, shiftedBBox)
332 wcsCopy = SkyWcs.readString(wcs.writeString())
333 self.assertTrue(wcsCopy.isFits)
335 return wcs
337 def checkNonFitsWcs(self, wcs):
338 """Check SkyWcs.getFitsMetadata for a WCS that cannot be represented as a FITS-WCS
339 """
340 # the modified WCS should not be representable as pure FITS-WCS
341 self.assertFalse(wcs.isFits)
342 with self.assertRaises(RuntimeError):
343 wcs.getFitsMetadata(True)
344 with self.assertRaises(RuntimeError):
345 wcs.getFitsMetadata(False)
347 # When a WCS is not valid FITS itself, we can explicitly install a FITS
348 # approximation. The 'wcsApprox' below is just arbitrary, not really
349 # an approximation in this test, but since SkyWcs shouldn't be
350 # responsible for the accuracy of the approximation that's fine.
351 wcsApprox = makeSkyWcs(
352 lsst.geom.Point2D(5, 4),
353 lsst.geom.SpherePoint(20.0, 30.0, lsst.geom.degrees),
354 np.identity(2),
355 )
356 wcsWithApproximation = wcs.copyWithFitsApproximation(wcsApprox)
357 self.assertTrue(wcsWithApproximation.hasFitsApproximation())
358 self.assertEqual(wcsWithApproximation.getFitsApproximation(), wcsApprox)
359 # When we save a SkyWcs with a FITS approximation it round-trips.
360 with lsst.utils.tests.getTempFilePath(".fits") as filePath:
361 wcsWithApproximation.writeFits(filePath)
362 wcsFromFits = SkyWcs.readFits(filePath)
363 self.assertTrue(wcsFromFits.hasFitsApproximation())
364 self.assertEqual(wcsFromFits.getFitsApproximation(), wcsApprox)
365 # When we save an Exposure with a non-FITS SkyWcs that has a FITS
366 # approximation, the approximation appears in the headers.
367 exposure = ExposureF(lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(3, 2)))
368 exposure.setWcs(wcsWithApproximation)
369 with lsst.utils.tests.getTempFilePath(".fits") as filePath:
370 exposure.writeFits(filePath)
371 exposureFromFits = ExposureF(filePath)
372 self.assertEqual(exposureFromFits.wcs, wcsWithApproximation)
373 self.assertEqual(exposureFromFits.wcs.getFitsApproximation(), wcsApprox)
374 with astropy.io.fits.open(filePath) as fits:
375 self.assertEqual(fits[1].header["CRVAL1"], 20.0)
376 self.assertEqual(fits[1].header["CRVAL2"], 30.0)
378 def testTanWcs(self):
379 """Check a variety of TanWcs, with crval not at a pole.
380 """
381 for crval, orientation, flipX in itertools.product(self.crvalList,
382 self.orientationList,
383 (False, True)):
384 self.checkTanWcs(crval=crval,
385 orientation=orientation,
386 flipX=flipX,
387 )
389 def testTanWcsFromFrameDict(self):
390 """Test making a TAN WCS from a FrameDict
391 """
392 cdMatrix = makeCdMatrix(scale=self.scale)
393 skyWcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
394 self.checkFrameDictConstructor(skyWcs, bbox=self.bbox)
396 def testGetFrameDict(self):
397 """Test that getFrameDict returns a deep copy
398 """
399 cdMatrix = makeCdMatrix(scale=self.scale)
400 skyWcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
401 for domain in ("PIXELS", "IWC", "SKY"):
402 frameDict = skyWcs.getFrameDict()
403 frameDict.removeFrame(domain)
404 self.assertFalse(frameDict.hasDomain(domain))
405 self.assertTrue(skyWcs.getFrameDict().hasDomain(domain))
407 def testMakeModifiedWcsNoActualPixels(self):
408 """Test makeModifiedWcs on a SkyWcs that has no ACTUAL_PIXELS frame
409 """
410 cdMatrix = makeCdMatrix(scale=self.scale)
411 originalWcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
412 originalFrameDict = originalWcs.getFrameDict()
414 # make an arbitrary but reasonable transform to insert using makeModifiedWcs
415 pixelTransform = makeRadialTransform([0.0, 1.0, 0.0, 0.0011])
416 # the result of the insertion should be as follows
417 desiredPixelsToSky = pixelTransform.then(originalWcs.getTransform())
419 pixPointList = ( # arbitrary but reasonable
420 lsst.geom.Point2D(0.0, 0.0),
421 lsst.geom.Point2D(1000.0, 0.0),
422 lsst.geom.Point2D(0.0, 2000.0),
423 lsst.geom.Point2D(-1111.0, -2222.0),
424 )
425 for modifyActualPixels in (False, True):
426 modifiedWcs = makeModifiedWcs(pixelTransform=pixelTransform,
427 wcs=originalWcs,
428 modifyActualPixels=modifyActualPixels)
429 modifiedFrameDict = modifiedWcs.getFrameDict()
430 skyList = modifiedWcs.pixelToSky(pixPointList)
432 # compare pixels to sky
433 desiredSkyList = desiredPixelsToSky.applyForward(pixPointList)
434 self.assertSpherePointListsAlmostEqual(skyList, desiredSkyList)
436 # compare pixels to IWC
437 pixelsToIwc = TransformPoint2ToPoint2(modifiedFrameDict.getMapping("PIXELS", "IWC"))
438 desiredPixelsToIwc = TransformPoint2ToPoint2(
439 pixelTransform.getMapping().then(originalFrameDict.getMapping("PIXELS", "IWC")))
440 self.assertPairListsAlmostEqual(pixelsToIwc.applyForward(pixPointList),
441 desiredPixelsToIwc.applyForward(pixPointList))
443 self.checkNonFitsWcs(modifiedWcs)
445 def testMakeModifiedWcsWithActualPixels(self):
446 """Test makeModifiedWcs on a SkyWcs that has an ACTUAL_PIXELS frame
447 """
448 cdMatrix = makeCdMatrix(scale=self.scale)
449 baseWcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
450 # model actual pixels to pixels as an arbitrary zoom factor;
451 # this is not realistic, but is fine for a unit test
452 actualPixelsToPixels = TransformPoint2ToPoint2(ast.ZoomMap(2, 0.72))
453 originalWcs = addActualPixelsFrame(baseWcs, actualPixelsToPixels)
454 originalFrameDict = originalWcs.getFrameDict()
456 # make an arbitrary but reasonable transform to insert using makeModifiedWcs
457 pixelTransform = makeRadialTransform([0.0, 1.0, 0.0, 0.0011]) # arbitrary but reasonable
459 pixPointList = ( # arbitrary but reasonable
460 lsst.geom.Point2D(0.0, 0.0),
461 lsst.geom.Point2D(1000.0, 0.0),
462 lsst.geom.Point2D(0.0, 2000.0),
463 lsst.geom.Point2D(-1111.0, -2222.0),
464 )
465 for modifyActualPixels in (True, False):
466 modifiedWcs = makeModifiedWcs(pixelTransform=pixelTransform,
467 wcs=originalWcs,
468 modifyActualPixels=modifyActualPixels)
469 modifiedFrameDict = modifiedWcs.getFrameDict()
470 self.assertEqual(modifiedFrameDict.getFrame(modifiedFrameDict.BASE).domain, "ACTUAL_PIXELS")
471 modifiedActualPixelsToPixels = \
472 TransformPoint2ToPoint2(modifiedFrameDict.getMapping("ACTUAL_PIXELS", "PIXELS"))
473 modifiedPixelsToIwc = TransformPoint2ToPoint2(modifiedFrameDict.getMapping("PIXELS", "IWC"))
475 # compare pixels to sky
476 skyList = modifiedWcs.pixelToSky(pixPointList)
477 if modifyActualPixels:
478 desiredPixelsToSky = pixelTransform.then(originalWcs.getTransform())
479 else:
480 originalPixelsToSky = \
481 TransformPoint2ToSpherePoint(originalFrameDict.getMapping("PIXELS", "SKY"))
482 desiredPixelsToSky = actualPixelsToPixels.then(pixelTransform).then(originalPixelsToSky)
483 desiredSkyList = desiredPixelsToSky.applyForward(pixPointList)
484 self.assertSpherePointListsAlmostEqual(skyList, desiredSkyList)
486 # compare ACTUAL_PIXELS to PIXELS and PIXELS to IWC
487 if modifyActualPixels:
488 # check that ACTUAL_PIXELS to PIXELS has been modified as expected
489 desiredActualPixelsToPixels = pixelTransform.then(actualPixelsToPixels)
490 self.assertPairListsAlmostEqual(modifiedActualPixelsToPixels.applyForward(pixPointList),
491 desiredActualPixelsToPixels.applyForward(pixPointList))
493 # check that PIXELS to IWC is unchanged
494 originalPixelsToIwc = TransformPoint2ToPoint2(originalFrameDict.getMapping("PIXELS", "IWC"))
495 self.assertPairListsAlmostEqual(modifiedPixelsToIwc.applyForward(pixPointList),
496 originalPixelsToIwc.applyForward(pixPointList))
498 else:
499 # check that ACTUAL_PIXELS to PIXELS is unchanged
500 self.assertPairListsAlmostEqual(actualPixelsToPixels.applyForward(pixPointList),
501 actualPixelsToPixels.applyForward(pixPointList))
503 # check that PIXELS to IWC has been modified as expected
504 desiredPixelsToIwc = TransformPoint2ToPoint2(
505 pixelTransform.getMapping().then(originalFrameDict.getMapping("PIXELS", "IWC")))
506 self.assertPairListsAlmostEqual(modifiedPixelsToIwc.applyForward(pixPointList),
507 desiredPixelsToIwc.applyForward(pixPointList))
509 self.checkNonFitsWcs(modifiedWcs)
511 def testMakeSkyWcsFromPixelsToFieldAngle(self):
512 """Test makeSkyWcs from a pixelsToFieldAngle transform
513 """
514 pixelSizeMm = 25e-3
515 # place the detector in several positions at several orientations
516 # use fewer CRVAL and orientations to speed up the test
517 for fpPosition, yaw, addOpticalDistortion, crval, pixelOrientation, \
518 flipX, projection in itertools.product(
519 (lsst.geom.Point3D(0, 0, 0), lsst.geom.Point3D(-100, 500, 1.5)),
520 (0*lsst.geom.degrees, 71*lsst.geom.degrees), (False, True),
521 self.crvalList[0:2], self.orientationList[0:2], (False, True), ("TAN", "STG")):
522 with self.subTest(fpPosition=repr(fpPosition), yaw=repr(yaw),
523 addOpticalDistortion=repr(addOpticalDistortion),
524 crval=repr(crval), orientation=repr(pixelOrientation)):
525 pixelsToFocalPlane = cameraGeom.Orientation(
526 fpPosition=fpPosition,
527 yaw=yaw,
528 ).makePixelFpTransform(lsst.geom.Extent2D(pixelSizeMm, pixelSizeMm))
529 # Compute crpix before adding optical distortion,
530 # since it is not affected by such distortion
531 crpix = pixelsToFocalPlane.applyInverse(lsst.geom.Point2D(0, 0))
532 radiansPerMm = self.scale.asRadians() / pixelSizeMm
533 focalPlaneToFieldAngle = lsst.afw.geom.makeTransform(
534 lsst.geom.AffineTransform(lsst.geom.LinearTransform.makeScaling(radiansPerMm)))
535 pixelsToFieldAngle = pixelsToFocalPlane.then(focalPlaneToFieldAngle)
537 cdMatrix = makeCdMatrix(scale=self.scale, orientation=pixelOrientation, flipX=flipX)
538 wcs1 = makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix, projection=projection)
540 if addOpticalDistortion:
541 # Model optical distortion as a pixel transform,
542 # so it can be added to the WCS created from crpix,
543 # cdMatrix, etc. using makeModifiedWcs
544 pixelTransform = makeRadialTransform([0.0, 1.0, 0.0, 0.0011])
545 pixelsToFieldAngle = pixelTransform.then(pixelsToFieldAngle)
546 wcs1 = makeModifiedWcs(pixelTransform=pixelTransform, wcs=wcs1, modifyActualPixels=False)
548 # orientation is with respect to detector x, y
549 # but this flavor of makeSkyWcs needs it with respect to focal plane x, y
550 focalPlaneOrientation = pixelOrientation + (yaw if flipX else -yaw)
551 wcs2 = makeSkyWcs(pixelsToFieldAngle=pixelsToFieldAngle,
552 orientation=focalPlaneOrientation,
553 flipX=flipX,
554 boresight=crval,
555 projection=projection)
556 self.assertWcsAlmostEqualOverBBox(wcs1, wcs2, self.bbox)
558 @unittest.skipIf(sys.version_info[0] < 3, "astropy.wcs rejects the header on py2")
559 def testAgainstAstropyWcs(self):
560 bbox = lsst.geom.Box2D(lsst.geom.Point2D(-1000, -1000), lsst.geom.Extent2D(2000, 2000))
561 for crval, orientation, flipX, projection in itertools.product(self.crvalList,
562 self.orientationList,
563 (False, True),
564 ("TAN", "STG", "CEA", "AIT")):
565 cdMatrix = makeCdMatrix(scale=self.scale, orientation=orientation, flipX=flipX)
566 metadata = makeSimpleWcsMetadata(crpix=self.crpix, crval=crval, cdMatrix=cdMatrix,
567 projection=projection)
568 header = makeLimitedFitsHeader(metadata)
569 astropyWcs = astropy.wcs.WCS(header)
570 skyWcs = makeSkyWcs(crpix=self.crpix, crval=crval, cdMatrix=cdMatrix, projection=projection)
571 # Most projections only seem to agree to within 1e-4 in the round trip test
572 self.assertSkyWcsAstropyWcsAlmostEqual(skyWcs=skyWcs, astropyWcs=astropyWcs, bbox=bbox)
574 def testPixelToSkyArray(self):
575 """Test the numpy-array version of pixelToSky
576 """
577 cdMatrix = makeCdMatrix(scale=self.scale)
578 wcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
580 xPoints = np.array([0.0, 1000.0, 0.0, -1111.0])
581 yPoints = np.array([0.0, 0.0, 2000.0, -2222.0])
583 pixPointList = [lsst.geom.Point2D(x, y) for x, y in zip(xPoints, yPoints)]
585 spherePoints = wcs.pixelToSky(pixPointList)
587 ra, dec = wcs.pixelToSkyArray(xPoints, yPoints, degrees=False)
588 for r, d, spherePoint in zip(ra, dec, spherePoints):
589 self.assertAlmostEqual(r, spherePoint.getRa().asRadians())
590 self.assertAlmostEqual(d, spherePoint.getDec().asRadians())
592 ra, dec = wcs.pixelToSkyArray(xPoints, yPoints, degrees=True)
593 for r, d, spherePoint in zip(ra, dec, spherePoints):
594 self.assertAlmostEqual(r, spherePoint.getRa().asDegrees())
595 self.assertAlmostEqual(d, spherePoint.getDec().asDegrees())
597 def testSkyToPixelArray(self):
598 """Test the numpy-array version of skyToPixel
599 """
600 cdMatrix = makeCdMatrix(scale=self.scale)
601 wcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
603 raPoints = np.array([3.92646679e-02, 3.59646622e+02,
604 3.96489283e-02, 4.70419353e-01])
605 decPoints = np.array([44.9722155, 44.97167735,
606 45.52775599, 44.3540619])
608 spherePointList = [lsst.geom.SpherePoint(ra*lsst.geom.degrees,
609 dec*lsst.geom.degrees)
610 for ra, dec in zip(raPoints, decPoints)]
612 pixPoints = wcs.skyToPixel(spherePointList)
614 x, y = wcs.skyToPixelArray(np.deg2rad(raPoints), np.deg2rad(decPoints))
615 for x0, y0, pixPoint in zip(x, y, pixPoints):
616 self.assertAlmostEqual(x0, pixPoint.getX())
617 self.assertAlmostEqual(y0, pixPoint.getY())
619 x, y = wcs.skyToPixelArray(raPoints, decPoints, degrees=True)
620 for x0, y0, pixPoint in zip(x, y, pixPoints):
621 self.assertAlmostEqual(x0, pixPoint.getX())
622 self.assertAlmostEqual(y0, pixPoint.getY())
624 def testStr(self):
625 """Test that we can get something coherent when printing a SkyWcs.
626 """
627 cdMatrix = makeCdMatrix(scale=self.scale)
628 skyWcs = makeSkyWcs(crpix=self.crpix, crval=self.crvalList[0], cdMatrix=cdMatrix)
629 self.assertIn(f"Sky Origin: {self.crvalList[0]}", str(skyWcs))
630 self.assertIn(f"Pixel Origin: {self.crpix}", str(skyWcs))
631 self.assertIn("Pixel Scale: ", str(skyWcs))
634class MetadataWcsTestCase(SkyWcsBaseTestCase):
635 """Test metadata constructor of SkyWcs
636 """
638 def setUp(self):
639 metadata = PropertyList()
640 for name, value in (
641 ("RADESYS", "ICRS"),
642 ("EQUINOX", 2000.),
643 ("CRVAL1", 215.604025685476),
644 ("CRVAL2", 53.1595451514076),
645 ("CRPIX1", 1109.99981456774),
646 ("CRPIX2", 560.018167811613),
647 ("CTYPE1", "RA---TAN"),
648 ("CTYPE2", "DEC--TAN"),
649 ("CUNIT1", "deg"),
650 ("CUNIT2", "deg"),
651 ("CD1_1", 5.10808596133527E-05),
652 ("CD1_2", 1.85579539217196E-07),
653 ("CD2_2", -5.10281493481982E-05),
654 ("CD2_1", -1.85579539217196E-07),
655 ):
656 metadata.set(name, value)
657 self.metadata = metadata
659 def tearDown(self):
660 del self.metadata
662 def checkWcs(self, skyWcs):
663 pixelOrigin = skyWcs.getPixelOrigin()
664 skyOrigin = skyWcs.getSkyOrigin()
665 for i in range(2):
666 # subtract 1 from FITS CRPIX to get LSST convention
667 self.assertAlmostEqual(pixelOrigin[i], self.metadata.getScalar(f"CRPIX{i+1}") - 1)
668 self.assertAnglesAlmostEqual(skyOrigin[i],
669 self.metadata.getScalar(f"CRVAL{i+1}")*lsst.geom.degrees)
670 cdMatrix = skyWcs.getCdMatrix()
671 for i, j in itertools.product(range(2), range(2)):
672 self.assertAlmostEqual(cdMatrix[i, j], self.metadata.getScalar(f"CD{i+1}_{j+1}"))
674 self.assertTrue(skyWcs.isFits)
676 skyWcsCopy = SkyWcs.readString(skyWcs.writeString())
677 self.assertTrue(skyWcsCopy.isFits)
678 self.checkMakeFlippedWcs(skyWcs)
680 @unittest.skipIf(sys.version_info[0] < 3, "astropy.wcs rejects the header on py2")
681 def testAgainstAstropyWcs(self):
682 skyWcs = makeSkyWcs(self.metadata, strip=False)
683 header = makeLimitedFitsHeader(self.metadata)
684 astropyWcs = astropy.wcs.WCS(header)
685 bbox = lsst.geom.Box2D(lsst.geom.Point2D(-1000, -1000), lsst.geom.Extent2D(3000, 3000))
686 self.assertSkyWcsAstropyWcsAlmostEqual(skyWcs=skyWcs, astropyWcs=astropyWcs, bbox=bbox)
688 def testLinearizeMethods(self):
689 skyWcs = makeSkyWcs(self.metadata)
690 # use a sky position near, but not at, the WCS origin
691 sky00 = skyWcs.getSkyOrigin().offset(45 * lsst.geom.degrees, 1.2 * lsst.geom.degrees)
692 pix00 = skyWcs.skyToPixel(sky00)
693 for skyUnit in (lsst.geom.degrees, lsst.geom.radians):
694 linPixToSky1 = skyWcs.linearizePixelToSky(sky00, skyUnit) # should match inverse of linSkyToPix1
695 linPixToSky2 = skyWcs.linearizePixelToSky(pix00, skyUnit) # should match inverse of linSkyToPix1
696 linSkyToPix1 = skyWcs.linearizeSkyToPixel(sky00, skyUnit)
697 linSkyToPix2 = skyWcs.linearizeSkyToPixel(pix00, skyUnit) # should match linSkyToPix1
699 for pixel in (pix00, pix00 + lsst.geom.Extent2D(1000, -1230)):
700 linSky = linPixToSky1(pixel)
701 self.assertPairsAlmostEqual(linPixToSky2(pixel), linSky)
702 self.assertPairsAlmostEqual(linSkyToPix1(linSky), pixel)
703 self.assertPairsAlmostEqual(linSkyToPix2(linSky), pixel)
705 sky00Doubles = sky00.getPosition(skyUnit)
706 pix00gApprox = linSkyToPix1(sky00Doubles)
707 self.assertPairsAlmostEqual(pix00gApprox, pix00)
708 self.assertAlmostEqual(pix00.getX(), pix00gApprox.getX())
709 self.assertAlmostEqual(pix00.getY(), pix00gApprox.getY())
710 pixelScale = skyWcs.getPixelScale(pix00)
711 pixelArea = pixelScale.asAngularUnits(skyUnit)**2
712 predictedPixelArea = 1 / linSkyToPix1.getLinear().computeDeterminant()
713 self.assertAlmostEqual(pixelArea, predictedPixelArea)
715 def testBasics(self):
716 skyWcs = makeSkyWcs(self.metadata, strip=False)
717 self.assertEqual(len(self.metadata.names(False)), 14)
718 self.checkWcs(skyWcs)
719 makeSkyWcs(self.metadata, strip=True)
720 self.assertEqual(len(self.metadata.names(False)), 0)
722 def testBasicsStrip(self):
723 stripWcsMetadata(self.metadata)
724 self.assertEqual(len(self.metadata.names(False)), 0)
725 # The metadata should be unchanged if we attempt to strip it again
726 metadataCopy = self.metadata.deepCopy()
727 stripWcsMetadata(self.metadata)
728 for key in self.metadata.keys():
729 self.assertEqual(self.metadata[key], metadataCopy[key])
731 def testNormalizationFk5(self):
732 """Test that readLsstSkyWcs correctly normalizes FK5 1975 to ICRS
733 """
734 equinox = 1975.0
735 metadata = self.metadata
737 metadata.set("RADESYS", "FK5")
738 metadata.set("EQUINOX", equinox)
739 crpix = lsst.geom.Point2D(metadata.getScalar("CRPIX1") - 1, metadata.getScalar("CRPIX2") - 1)
740 # record the original CRVAL before reading and stripping metadata
741 crvalFk5Deg = (metadata.getScalar("CRVAL1"), metadata.getScalar("CRVAL2"))
743 # create the wcs and retrieve crval
744 skyWcs = makeSkyWcs(metadata)
745 crval = skyWcs.getSkyOrigin()
747 # compare to computed crval
748 computedCrval = skyWcs.pixelToSky(crpix)
749 self.assertSpherePointsAlmostEqual(crval, computedCrval)
751 # get predicted crval by converting with astropy
752 crvalFk5 = astropy.coordinates.SkyCoord(crvalFk5Deg[0], crvalFk5Deg[1], frame="fk5",
753 equinox=f"J{equinox}", unit="deg")
754 predictedCrvalIcrs = crvalFk5.icrs
755 predictedCrval = lsst.geom.SpherePoint(predictedCrvalIcrs.ra.radian, predictedCrvalIcrs.dec.radian,
756 lsst.geom.radians)
757 self.assertSpherePointsAlmostEqual(crval, predictedCrval, maxSep=0.002*lsst.geom.arcseconds)
759 def testNormalizationDecRa(self):
760 """Test that a Dec, RA WCS is normalized to RA, Dec
761 """
762 crpix = lsst.geom.Point2D(self.metadata.getScalar("CRPIX1") - 1,
763 self.metadata.getScalar("CRPIX2") - 1)
765 # swap RA, Decaxes in metadata
766 crvalIn = lsst.geom.SpherePoint(self.metadata.getScalar("CRVAL1"),
767 self.metadata.getScalar("CRVAL2"), lsst.geom.degrees)
768 self.metadata.set("CRVAL1", crvalIn[1].asDegrees())
769 self.metadata.set("CRVAL2", crvalIn[0].asDegrees())
770 self.metadata.set("CTYPE1", "DEC--TAN")
771 self.metadata.set("CTYPE2", "RA---TAN")
773 # create the wcs
774 skyWcs = makeSkyWcs(self.metadata)
776 # compare pixel origin to input crval
777 crval = skyWcs.getSkyOrigin()
778 self.assertSpherePointsAlmostEqual(crval, crvalIn)
780 # compare to computed crval
781 computedCrval = skyWcs.pixelToSky(crpix)
782 self.assertSpherePointsAlmostEqual(crval, computedCrval)
784 def testReadDESHeader(self):
785 """Verify that we can read a DES header"""
786 self.metadata.set("RADESYS", "ICRS ") # note trailing white space
787 self.metadata.set("CTYPE1", "RA---TPV")
788 self.metadata.set("CTYPE2", "DEC--TPV")
790 skyWcs = makeSkyWcs(self.metadata, strip=False)
791 self.checkWcs(skyWcs)
793 def testCD_PC(self):
794 """Test that we can read a FITS file with both CD and PC keys (like early Suprimecam files)"""
795 md = PropertyList()
796 for k, v in (
797 ("EQUINOX", 2000.0),
798 ("RADESYS", "ICRS"),
799 ("CRPIX1", 5353.0),
800 ("CRPIX2", -35.0),
801 ("CD1_1", 0.0),
802 ("CD1_2", -5.611E-05),
803 ("CD2_1", -5.611E-05),
804 ("CD2_2", -0.0),
805 ("CRVAL1", 4.5789875),
806 ("CRVAL2", 16.30004444),
807 ("CUNIT1", "deg"),
808 ("CUNIT2", "deg"),
809 ("CTYPE1", "RA---TAN"),
810 ("CTYPE2", "DEC--TAN"),
811 ("CDELT1", -5.611E-05),
812 ("CDELT2", 5.611E-05),
813 ):
814 md.set(k, v)
816 wcs = makeSkyWcs(md, strip=False)
818 pixPos = lsst.geom.Point2D(1000, 2000)
819 pred_skyPos = lsst.geom.SpherePoint(4.459815023498577, 16.544199850984768, lsst.geom.degrees)
821 skyPos = wcs.pixelToSky(pixPos)
822 self.assertSpherePointsAlmostEqual(skyPos, pred_skyPos)
824 for badPC in (False, True):
825 for k, v in (
826 ("PC001001", 0.0),
827 ("PC001002", -1.0 if badPC else 1.0),
828 ("PC002001", 1.0 if badPC else -1.0),
829 ("PC002002", 0.0),
830 ):
831 md.set(k, v)
833 # Check Greisen and Calabretta A&A 395 1061 (2002), Eq. 3
834 if not badPC:
835 for i in (1, 2,):
836 for j in (1, 2,):
837 self.assertEqual(md.getScalar(f"CD{i}_{j}"),
838 md.getScalar(f"CDELT{i}")*md.getScalar(f"PC00{i}00{j}"))
840 wcs2 = makeSkyWcs(md, strip=False)
841 skyPos2 = wcs2.pixelToSky(pixPos)
842 self.assertSpherePointsAlmostEqual(skyPos2, pred_skyPos)
844 def testNoEpoch(self):
845 """Ensure we're not writing epoch headers (DATE-OBS, MJD-OBS)"""
846 self.metadata.set("EQUINOX", 2000.0) # Triggers AST writing DATE-OBS, MJD-OBS
847 skyWcs = makeSkyWcs(self.metadata)
848 header = skyWcs.getFitsMetadata()
849 self.assertFalse(header.exists("DATE-OBS"))
850 self.assertFalse(header.exists("MJD-OBS"))
852 def testCdMatrix(self):
853 """Ensure we're writing CD matrix elements even if they're zero"""
854 self.metadata.remove("CD1_2")
855 self.metadata.remove("CD2_1")
856 skyWcs = makeSkyWcs(self.metadata, strip=False)
857 header = skyWcs.getFitsMetadata()
858 for keyword in ("CD1_1", "CD2_2"):
859 # There's some mild rounding going on
860 self.assertFloatsAlmostEqual(header.get(keyword), self.metadata.get(keyword), atol=1.0e-16)
861 for keyword in ("CD1_2", "CD2_1"):
862 self.assertTrue(header.exists(keyword))
863 self.assertEqual(header.get(keyword), 0.0)
866class TestTanSipTestCase(SkyWcsBaseTestCase):
868 def setUp(self):
869 metadata = PropertyList()
870 # the following was fit using CreateWcsWithSip from meas_astrom
871 # and is valid over this bbox: (minimum=(0, 0), maximum=(3030, 3030))
872 # This same metadata was used to create testdata/oldTanSipwWs.fits
873 for name, value in (
874 ("RADESYS", "ICRS"),
875 ("CTYPE1", "RA---TAN-SIP"),
876 ("CTYPE2", "DEC--TAN-SIP"),
877 ("CRPIX1", 1531.1824767147),
878 ("CRPIX2", 1531.1824767147),
879 ("CRVAL1", 43.035511801383),
880 ("CRVAL2", 44.305697682784),
881 ("CUNIT1", "deg"),
882 ("CUNIT2", "deg"),
883 ("CD1_1", 0.00027493991598151),
884 ("CD1_2", -3.2758487104158e-06),
885 ("CD2_1", 3.2301310675830e-06),
886 ("CD2_2", 0.00027493937506632),
887 ("A_ORDER", 5),
888 ("A_0_2", -1.7769487466972e-09),
889 ("A_0_3", 5.3745894718340e-13),
890 ("A_0_4", -7.2921116596880e-17),
891 ("A_0_5", 8.6947236956136e-21),
892 ("A_1_1", 5.4246387438098e-08),
893 ("A_1_2", -1.5689083084641e-12),
894 ("A_1_3", 1.2424130500997e-16),
895 ("A_1_4", 3.9982572658006e-20),
896 ("A_2_0", 4.9268299826160e-08),
897 ("A_2_1", 1.6365657558495e-12),
898 ("A_2_2", 1.1976983061953e-16),
899 ("A_2_3", -1.7262037266467e-19),
900 ("A_3_0", -5.9235031179999e-13),
901 ("A_3_1", -3.4444326387310e-16),
902 ("A_3_2", 1.4377441160800e-19),
903 ("A_4_0", 1.8736407845095e-16),
904 ("A_4_1", 2.9213314172884e-20),
905 ("A_5_0", -5.3601346091084e-20),
906 ("B_ORDER", 5),
907 ("B_0_2", 4.9268299822979e-08),
908 ("B_0_3", -5.9235032026906e-13),
909 ("B_0_4", 1.8736407776035e-16),
910 ("B_0_5", -5.3601341373220e-20),
911 ("B_1_1", 5.4246387435453e-08),
912 ("B_1_2", 1.6365657531115e-12),
913 ("B_1_3", -3.4444326228808e-16),
914 ("B_1_4", 2.9213312399941e-20),
915 ("B_2_0", -1.7769487494962e-09),
916 ("B_2_1", -1.5689082999319e-12),
917 ("B_2_2", 1.1976983393279e-16),
918 ("B_2_3", 1.4377441169892e-19),
919 ("B_3_0", 5.3745894237186e-13),
920 ("B_3_1", 1.2424130479929e-16),
921 ("B_3_2", -1.7262036838229e-19),
922 ("B_4_0", -7.2921117326608e-17),
923 ("B_4_1", 3.9982566975450e-20),
924 ("B_5_0", 8.6947240592408e-21),
925 ("AP_ORDER", 6),
926 ("AP_0_0", -5.4343024221207e-11),
927 ("AP_0_1", 5.5722265946666e-12),
928 ("AP_0_2", 1.7769484042400e-09),
929 ("AP_0_3", -5.3773609554820e-13),
930 ("AP_0_4", 7.3035278852156e-17),
931 ("AP_0_5", -8.7151153799062e-21),
932 ("AP_0_6", 3.2535945427624e-27),
933 ("AP_1_0", -3.8944805432871e-12),
934 ("AP_1_1", -5.4246388067582e-08),
935 ("AP_1_2", 1.5741716194971e-12),
936 ("AP_1_3", -1.2447067748187e-16),
937 ("AP_1_4", -3.9960260822306e-20),
938 ("AP_1_5", 1.1297941471380e-26),
939 ("AP_2_0", -4.9268299293185e-08),
940 ("AP_2_1", -1.6256111849359e-12),
941 ("AP_2_2", -1.1973373130440e-16),
942 ("AP_2_3", 1.7266948205700e-19),
943 ("AP_2_4", -3.7059606160753e-26),
944 ("AP_3_0", 5.9710911995811e-13),
945 ("AP_3_1", 3.4464427650041e-16),
946 ("AP_3_2", -1.4381853884204e-19),
947 ("AP_3_3", -7.6527426974322e-27),
948 ("AP_4_0", -1.8748435698960e-16),
949 ("AP_4_1", -2.9267280226373e-20),
950 ("AP_4_2", 4.8004317051259e-26),
951 ("AP_5_0", 5.3657330221120e-20),
952 ("AP_5_1", -1.6904065766661e-27),
953 ("AP_6_0", -1.9484495120493e-26),
954 ("BP_ORDER", 6),
955 ("BP_0_0", -5.4291220607725e-11),
956 ("BP_0_1", -3.8944871307931e-12),
957 ("BP_0_2", -4.9268299290361e-08),
958 ("BP_0_3", 5.9710912831833e-13),
959 ("BP_0_4", -1.8748435594265e-16),
960 ("BP_0_5", 5.3657325543368e-20),
961 ("BP_0_6", -1.9484577299247e-26),
962 ("BP_1_0", 5.5722051513577e-12),
963 ("BP_1_1", -5.4246388065000e-08),
964 ("BP_1_2", -1.6256111821465e-12),
965 ("BP_1_3", 3.4464427499767e-16),
966 ("BP_1_4", -2.9267278448109e-20),
967 ("BP_1_5", -1.6904244067295e-27),
968 ("BP_2_0", 1.7769484069376e-09),
969 ("BP_2_1", 1.5741716110182e-12),
970 ("BP_2_2", -1.1973373446176e-16),
971 ("BP_2_3", -1.4381853893526e-19),
972 ("BP_2_4", 4.8004294492911e-26),
973 ("BP_3_0", -5.3773609074713e-13),
974 ("BP_3_1", -1.2447067726801e-16),
975 ("BP_3_2", 1.7266947774875e-19),
976 ("BP_3_3", -7.6527556667042e-27),
977 ("BP_4_0", 7.3035279660505e-17),
978 ("BP_4_1", -3.9960255158200e-20),
979 ("BP_4_2", -3.7059659675039e-26),
980 ("BP_5_0", -8.7151157361284e-21),
981 ("BP_5_1", 1.1297944388060e-26),
982 ("BP_6_0", 3.2535788867488e-27),
983 ):
984 metadata.set(name, value)
985 self.metadata = metadata
986 self.bbox = lsst.geom.Box2D(lsst.geom.Point2D(-1000, -1000), lsst.geom.Extent2D(3000, 3000))
988 def testTanSipFromFrameDict(self):
989 """Test making a TAN-SIP WCS from a FrameDict
990 """
991 skyWcs = makeSkyWcs(self.metadata, strip=False)
992 self.checkFrameDictConstructor(skyWcs, bbox=self.bbox)
994 def testFitsMetadata(self):
995 """Test that getFitsMetadata works for TAN-SIP
996 """
997 skyWcs = makeSkyWcs(self.metadata, strip=False)
998 self.assertTrue(skyWcs.isFits)
999 fitsMetadata = skyWcs.getFitsMetadata(precise=True)
1000 skyWcsCopy = makeSkyWcs(fitsMetadata)
1001 self.assertWcsAlmostEqualOverBBox(skyWcs, skyWcsCopy, self.bbox)
1002 self.checkPersistence(skyWcs, bbox=self.bbox)
1004 def testGetIntermediateWorldCoordsToSky(self):
1005 """Test getIntermediateWorldCoordsToSky and getPixelToIntermediateWorldCoords
1006 """
1007 crpix = lsst.geom.Extent2D(self.metadata.getScalar("CRPIX1") - 1,
1008 self.metadata.getScalar("CRPIX2") - 1)
1009 skyWcs = makeSkyWcs(self.metadata, strip=False)
1010 for simplify in (False, True):
1011 pixelToIwc = getPixelToIntermediateWorldCoords(skyWcs, simplify)
1012 iwcToSky = getIntermediateWorldCoordsToSky(skyWcs, simplify)
1013 self.assertTrue(isinstance(pixelToIwc, TransformPoint2ToPoint2))
1014 self.assertTrue(isinstance(iwcToSky, TransformPoint2ToSpherePoint))
1015 if simplify:
1016 self.assertTrue(pixelToIwc.getMapping().isSimple)
1017 self.assertTrue(iwcToSky.getMapping().isSimple)
1018 # else the mapping may have already been simplified inside the WCS,
1019 # so don't assert isSimple is false
1021 # check that the chained transforms produce the same results as the WCS
1022 # in the forward and inverse direction
1023 pixPosList = []
1024 for dx in (0, 1000):
1025 for dy in (0, 1000):
1026 pixPosList.append(lsst.geom.Point2D(dx, dy) + crpix)
1027 iwcPosList = pixelToIwc.applyForward(pixPosList)
1028 skyPosList = iwcToSky.applyForward(iwcPosList)
1029 self.assertSpherePointListsAlmostEqual(skyPosList, skyWcs.pixelToSky(pixPosList))
1030 self.assertPairListsAlmostEqual(pixelToIwc.applyInverse(iwcToSky.applyInverse(skyPosList)),
1031 skyWcs.skyToPixel(skyPosList))
1033 self.assertPairListsAlmostEqual(iwcPosList, iwcToSky.applyInverse(skyPosList))
1034 self.assertPairListsAlmostEqual(pixPosList, pixelToIwc.applyInverse(iwcPosList))
1036 # compare extracted pixelToIwc to a version of pixelToIwc computed directly from the metadata
1037 ourPixelToIwc = makeSipPixelToIwc(self.metadata)
1038 self.assertPairListsAlmostEqual(pixelToIwc.applyForward(pixPosList),
1039 ourPixelToIwc.applyForward(pixPosList))
1041 # compare extracted iwcToPixel to a version of iwcToPixel computed directly from the metadata
1042 ourIwcToPixel = makeSipIwcToPixel(self.metadata)
1043 self.assertPairListsAlmostEqual(pixelToIwc.applyInverse(iwcPosList),
1044 ourIwcToPixel.applyForward(iwcPosList))
1046 @unittest.skipIf(sys.version_info[0] < 3, "astropy.wcs rejects the header on py2")
1047 def testAgainstAstropyWcs(self):
1048 skyWcs = makeSkyWcs(self.metadata, strip=False)
1049 header = makeLimitedFitsHeader(self.metadata)
1050 astropyWcs = astropy.wcs.WCS(header)
1051 self.assertSkyWcsAstropyWcsAlmostEqual(skyWcs=skyWcs, astropyWcs=astropyWcs, bbox=self.bbox)
1053 def testMakeTanSipWcs(self):
1054 referenceWcs = makeSkyWcs(self.metadata, strip=False)
1056 crpix = lsst.geom.Point2D(self.metadata.getScalar("CRPIX1") - 1,
1057 self.metadata.getScalar("CRPIX2") - 1)
1058 crval = lsst.geom.SpherePoint(self.metadata.getScalar("CRVAL1"),
1059 self.metadata.getScalar("CRVAL2"), lsst.geom.degrees)
1060 cdMatrix = getCdMatrixFromMetadata(self.metadata)
1061 sipA = getSipMatrixFromMetadata(self.metadata, "A")
1062 sipB = getSipMatrixFromMetadata(self.metadata, "B")
1063 sipAp = getSipMatrixFromMetadata(self.metadata, "AP")
1064 sipBp = getSipMatrixFromMetadata(self.metadata, "BP")
1065 skyWcs1 = makeTanSipWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix, sipA=sipA, sipB=sipB)
1066 skyWcs2 = makeTanSipWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix, sipA=sipA, sipB=sipB,
1067 sipAp=sipAp, sipBp=sipBp)
1069 self.assertWcsAlmostEqualOverBBox(referenceWcs, skyWcs1, self.bbox)
1070 self.assertWcsAlmostEqualOverBBox(referenceWcs, skyWcs2, self.bbox)
1071 self.checkMakeFlippedWcs(skyWcs1)
1072 self.checkMakeFlippedWcs(skyWcs2)
1074 def testReadWriteFits(self):
1075 wcsFromMetadata = makeSkyWcs(self.metadata)
1076 with lsst.utils.tests.getTempFilePath(".fits") as filePath:
1077 wcsFromMetadata.writeFits(filePath)
1078 wcsFromFits = SkyWcs.readFits(filePath)
1080 self.assertWcsAlmostEqualOverBBox(wcsFromFits, wcsFromMetadata, self.bbox, maxDiffPix=0,
1081 maxDiffSky=0*lsst.geom.radians)
1083 def testReadOldTanSipFits(self):
1084 """Test reading a FITS file containing data for an lsst::afw::image::TanWcs
1086 That file was made using the same metadata as this test
1087 """
1088 dataDir = os.path.join(os.path.split(__file__)[0], "data")
1089 filePath = os.path.join(dataDir, "oldTanSipWcs.fits")
1090 wcsFromFits = SkyWcs.readFits(filePath)
1092 wcsFromMetadata = makeSkyWcs(self.metadata)
1094 bbox = lsst.geom.Box2D(lsst.geom.Point2D(-1000, -1000), lsst.geom.Extent2D(3000, 3000))
1095 self.assertWcsAlmostEqualOverBBox(wcsFromFits, wcsFromMetadata, bbox)
1097 def testReadOldTanFits(self):
1098 """Test reading a FITS file containing data for an lsst::afw::image::TanWcs
1100 That file was made using the same metadata follows
1101 (like self.metadata without the distortion)
1102 """
1103 tanMetadata = PropertyList()
1104 # the following was fit using CreateWcsWithSip from meas_astrom
1105 # and is valid over this bbox: (minimum=(0, 0), maximum=(3030, 3030))
1106 # This same metadata was used to create testdata/oldTanSipwWs.fits
1107 for name, value in (
1108 ("RADESYS", "ICRS"),
1109 ("CTYPE1", "RA---TAN"),
1110 ("CTYPE2", "DEC--TAN"),
1111 ("CRPIX1", 1531.1824767147),
1112 ("CRPIX2", 1531.1824767147),
1113 ("CRVAL1", 43.035511801383),
1114 ("CRVAL2", 44.305697682784),
1115 ("CUNIT1", "deg"),
1116 ("CUNIT2", "deg"),
1117 ("CD1_1", 0.00027493991598151),
1118 ("CD1_2", -3.2758487104158e-06),
1119 ("CD2_1", 3.2301310675830e-06),
1120 ("CD2_2", 0.00027493937506632),
1121 ):
1122 tanMetadata.set(name, value)
1124 dataDir = os.path.join(os.path.split(__file__)[0], "data")
1125 filePath = os.path.join(dataDir, "oldTanWcs.fits")
1126 wcsFromFits = SkyWcs.readFits(filePath)
1128 wcsFromMetadata = makeSkyWcs(tanMetadata)
1130 bbox = lsst.geom.Box2D(lsst.geom.Point2D(-1000, -1000), lsst.geom.Extent2D(3000, 3000))
1131 self.assertWcsAlmostEqualOverBBox(wcsFromFits, wcsFromMetadata, bbox)
1134class WcsPairTransformTestCase(SkyWcsBaseTestCase):
1135 """Test functionality of makeWcsPairTransform.
1136 """
1137 def setUp(self):
1138 SkyWcsBaseTestCase.setUp(self)
1139 crpix = lsst.geom.Point2D(100, 100)
1140 crvalList = [
1141 lsst.geom.SpherePoint(0, 45, lsst.geom.degrees),
1142 lsst.geom.SpherePoint(0.00001, 45, lsst.geom.degrees),
1143 lsst.geom.SpherePoint(359.99999, 45, lsst.geom.degrees),
1144 lsst.geom.SpherePoint(30, 89.99999, lsst.geom.degrees),
1145 ]
1146 orientationList = [
1147 0 * lsst.geom.degrees,
1148 0.00001 * lsst.geom.degrees,
1149 -0.00001 * lsst.geom.degrees,
1150 -45 * lsst.geom.degrees,
1151 90 * lsst.geom.degrees,
1152 ]
1153 scale = 1.0 * lsst.geom.arcseconds
1155 self.wcsList = []
1156 for crval in crvalList:
1157 for orientation in orientationList:
1158 cd = makeCdMatrix(scale=scale, orientation=orientation)
1159 self.wcsList.append(makeSkyWcs(
1160 crpix=crpix,
1161 crval=crval,
1162 cdMatrix=cd))
1163 self.pixelPoints = [lsst.geom.Point2D(x, y) for x, y in
1164 itertools.product((0.0, -2.0, 42.5, 1042.3),
1165 (27.6, -0.1, 0.0, 196.0))]
1167 def testGenericWcs(self):
1168 """Test that input and output points represent the same sky position.
1170 Would prefer a black-box test, but don't have the numbers for it.
1171 """
1172 inPoints = self.pixelPoints
1173 for wcs1 in self.wcsList:
1174 for wcs2 in self.wcsList:
1175 transform = makeWcsPairTransform(wcs1, wcs2)
1176 outPoints = transform.applyForward(inPoints)
1177 inPointsRoundTrip = transform.applyInverse(outPoints)
1178 self.assertPairListsAlmostEqual(inPoints, inPointsRoundTrip)
1179 self.assertSpherePointListsAlmostEqual(wcs1.pixelToSky(inPoints),
1180 wcs2.pixelToSky(outPoints))
1182 def testSameWcs(self):
1183 """Confirm that pairing two identical Wcs gives an identity transform.
1184 """
1185 for wcs in self.wcsList:
1186 transform = makeWcsPairTransform(wcs, wcs)
1187 # check that the transform has been simplified
1188 self.assertTrue(transform.getMapping().isSimple)
1189 # check the transform
1190 outPoints1 = transform.applyForward(self.pixelPoints)
1191 outPoints2 = transform.applyInverse(outPoints1)
1192 self.assertPairListsAlmostEqual(self.pixelPoints, outPoints1)
1193 self.assertPairListsAlmostEqual(outPoints1, outPoints2)
1196class TestMemory(lsst.utils.tests.MemoryTestCase):
1197 pass
1200def setup_module(module):
1201 lsst.utils.tests.init()
1204if __name__ == "__main__": 1204 ↛ 1205line 1204 didn't jump to line 1205 because the condition on line 1204 was never true
1205 lsst.utils.tests.init()
1206 unittest.main()