Coverage for python/lsst/dax/images/cutout/_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 

24__all__ = ("projection_and_bbox_from_fits_header",) 

25 

26from collections.abc import Sequence 

27 

28import astropy.io.fits 

29import astropy.units as u 

30import numpy as np 

31 

32from lsst.images import Box, GeneralFrame, SkyProjection 

33from lsst.images._transforms import _ast 

34 

35# Identifier AST gives the alternate ("A") FITS WCS that encodes the parent 

36# pixel origin (XY0) for LSST image data. 

37_PARENT_WCS_IDENT = "A" 

38 

39# Pixel coordinate frame for the projection. Only the unit is significant. 

40_PIXEL_FRAME = GeneralFrame(unit=u.pix) 

41 

42 

43def projection_and_bbox_from_fits_header( 

44 header: astropy.io.fits.Header, shape: Sequence[int] 

45) -> tuple[SkyProjection, Box]: 

46 """Build a sky projection and parent bounding box from a FITS header. 

47 

48 Parameters 

49 ---------- 

50 header : `astropy.io.fits.Header` 

51 FITS header carrying a primary celestial WCS and an ``"A"`` alternate 

52 WCS that encodes the parent pixel origin. 

53 shape : `~collections.abc.Sequence` [ `int` ] 

54 Array shape ``(ny, nx)`` of the image the header describes, as returned 

55 by ``astropy.io.fits.ImageHDU.shape``. 

56 

57 Returns 

58 ------- 

59 projection : `lsst.images.SkyProjection` 

60 Sky projection for the primary WCS, in ICRS. 

61 bbox : `lsst.images.Box` 

62 Parent bounding box: the origin comes from the ``"A"`` WCS and the 

63 size from ``shape``. 

64 

65 Notes 

66 ----- 

67 All AST handling goes through the wrapper ``lsst.images`` exposes in 

68 ``lsst.images._transforms._ast``, which presents the astshim interface 

69 regardless of whether astshim or starlink-pyast backs it. This function is 

70 deliberately confined to that wrapper plus the public ``lsst.images`` APIs 

71 so it can move into ``lsst.images`` later (e.g. as 

72 ``SkyProjection.from_fits_header``). 

73 """ 

74 frame_set = _ast.FitsChan(_ast.StringStream(header.tostring())).read() 

75 parent_index = _frame_index_by_ident(frame_set, _PARENT_WCS_IDENT) 

76 

77 # The base frame is the raw 1-based pixel grid; the "A" frame is the parent 

78 # pixel system that carries the XY0 origin. Read XY0 as the parent 

79 # coordinate of grid (1, 1) before retargeting, then build the projection 

80 # on the parent frame so its pixels match the parent bounding box (and the 

81 # legacy ``makeSkyWcs`` convention). 

82 grid_to_parent = frame_set.getMapping(frame_set.base, parent_index) 

83 origin = grid_to_parent.applyForward(np.array([[1.0], [1.0]])).ravel() 

84 x0 = int(round(float(origin[0]))) 

85 y0 = int(round(float(origin[1]))) 

86 

87 frame_set.base = parent_index 

88 projection = SkyProjection.from_ast_frame_set(frame_set, _PIXEL_FRAME) 

89 

90 ny, nx = int(shape[0]), int(shape[1]) 

91 bbox = Box.factory[y0 : y0 + ny, x0 : x0 + nx] 

92 return projection, bbox 

93 

94 

95def _frame_index_by_ident(frame_set: _ast.FrameSet, ident: str) -> int: 

96 """Return the 1-based index of the frame whose identifier is ``ident``.""" 

97 for index in range(1, frame_set.nFrame + 1): 97 ↛ 100line 97 didn't jump to line 100 because the loop on line 97 didn't complete

98 if frame_set.getFrame(index).ident.strip() == ident: 

99 return index 

100 raise LookupError(f"FITS header has no AST frame with identifier {ident!r}.")