Coverage for tests/test_fits_projection.py: 94%

1# This file is part of dax_images_cutout. 

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/>. 

21 

22from __future__ import annotations 

23 

24import unittest 

25 

26import astropy.io.fits 

27import numpy as np 

28 

29from lsst.dax.images.cutout._fits_projection import projection_and_bbox_from_fits_header 

30from lsst.images import Box, SkyProjection 

31 

32try: 

33 from lsst.afw.geom import getImageXY0FromMetadata, makeSkyWcs 

34 from lsst.daf.base import PropertyList 

35 from lsst.geom import Box2I, Extent2I 

36 

37 HAVE_AFW = True 

38except ImportError: 

39 HAVE_AFW = False 

40 

41 

42def _make_header() -> astropy.io.fits.Header: 

43 """Make a primary gnomonic sky WCS plus an 'A' alternate WCS for XY0.""" 

44 header = astropy.io.fits.Header() 

45 cards = { 

46 "WCSAXES": 2, 

47 "CRPIX1": 10.0, 

48 "CRPIX2": 20.0, 

49 "CRVAL1": 12.0, 

50 "CRVAL2": 13.0, 

51 "CD1_1": -2.78e-05, 

52 "CD2_2": 2.78e-05, 

53 "CTYPE1": "RA---TAN", 

54 "CTYPE2": "DEC--TAN", 

55 "CRPIX1A": 1.0, 

56 "CRPIX2A": 1.0, 

57 "CRVAL1A": 100.0, 

58 "CRVAL2A": 200.0, 

59 "CD1_1A": 1.0, 

60 "CD2_2A": 1.0, 

61 "CTYPE1A": "LINEAR", 

62 "CTYPE2A": "LINEAR", 

63 } 

64 for key, value in cards.items(): 

65 header[key] = value 

66 return header 

67 

68 

69class FitsProjectionTestCase(unittest.TestCase): 

70 """Tests for the FITS-header projection helper (no afw required).""" 

71 

72 def setUp(self) -> None: 

73 self.header = _make_header() 

74 # (ny, nx) as returned by astropy ``hdu.shape``. 

75 self.shape = (64, 48) 

76 

77 def test_bbox_origin_and_size(self) -> None: 

78 _, bbox = projection_and_bbox_from_fits_header(self.header, self.shape) 

79 self.assertIsInstance(bbox, Box) 

80 # The "A" WCS has CRPIX*A = 1 and CRVAL*A = (100, 200), so grid (1, 1) 

81 # maps to the parent origin (100, 200); the size comes from ``shape``. 

82 self.assertEqual(bbox.x.start, 100) 

83 self.assertEqual(bbox.y.start, 200) 

84 self.assertEqual(bbox.x.size, self.shape[1]) 

85 self.assertEqual(bbox.y.size, self.shape[0]) 

86 

87 def test_reference_pixel_maps_to_crval(self) -> None: 

88 projection, bbox = projection_and_bbox_from_fits_header(self.header, self.shape) 

89 self.assertIsInstance(projection, SkyProjection) 

90 # The projection is in parent pixel coordinates; the primary reference 

91 # pixel CRPIX (1-based grid) sits at parent CRPIX - 1 + origin and must 

92 # map to CRVAL. 

93 ref_x = self.header["CRPIX1"] - 1 + bbox.x.start 

94 ref_y = self.header["CRPIX2"] - 1 + bbox.y.start 

95 sky = projection.pixel_to_sky(x=np.array([ref_x]), y=np.array([ref_y])) 

96 self.assertAlmostEqual(float(sky.ra.deg[0]), self.header["CRVAL1"], places=6) 

97 self.assertAlmostEqual(float(sky.dec.deg[0]), self.header["CRVAL2"], places=6) 

98 

99 

100@unittest.skipUnless(HAVE_AFW, "lsst.afw/lsst.geom not available") 

101class FitsProjectionAfwReferenceTestCase(unittest.TestCase): 

102 """Cross-check the helper against the legacy afw implementation.""" 

103 

104 def setUp(self) -> None: 

105 self.header = _make_header() 

106 self.shape = (64, 48) 

107 

108 def test_projection_matches_afw(self) -> None: 

109 projection, bbox = projection_and_bbox_from_fits_header(self.header, self.shape) 

110 pl = PropertyList() 

111 pl.update(self.header) 

112 wcs = makeSkyWcs(pl) 

113 xs = np.array([bbox.x.start + 5.0, bbox.x.start + 30.0]) 

114 ys = np.array([bbox.y.start + 7.0, bbox.y.start + 25.0]) 

115 sky = projection.pixel_to_sky(x=xs, y=ys) 

116 for i in range(len(xs)): 

117 reference = wcs.pixelToSky(float(xs[i]), float(ys[i])) 

118 self.assertAlmostEqual(float(sky.ra.deg[i]), reference.getRa().asDegrees(), places=6) 

119 self.assertAlmostEqual(float(sky.dec.deg[i]), reference.getDec().asDegrees(), places=6) 

120 

121 def test_bbox_matches_afw_reference(self) -> None: 

122 _, bbox = projection_and_bbox_from_fits_header(self.header, self.shape) 

123 # afw reference for the parent bbox (XY0 from 'A' WCS + dimensions). 

124 pl = PropertyList() 

125 pl.update(self.header) 

126 xy0 = getImageXY0FromMetadata(pl, "A", strip=False) 

127 dimensions = Extent2I(self.shape[1], self.shape[0]) 

128 self.assertEqual(bbox, Box.from_legacy(Box2I(xy0, dimensions))) 

129 

130 

131if __name__ == "__main__": 131 ↛ 132line 131 didn't jump to line 132 because the condition on line 131 was never true

132 unittest.main()