31def makeHpxWcs(hips_order, pixel, shift_order=9):
32 """
33 Make a SkyWcs object with HEALPix grid projection (HPX).
34
35 The SkyWcs generated by this function is suitable to be used with a
36 Hierarchical Progressive Survey (HiPS) FITS image as described in
37 https://www.ivoa.net/documents/HiPS/20170519/REC-HIPS-1.0-20170519.pdf
38
39 A HiPS image covers one HEALPix cell, with the HEALPix nside equal to
40 2**hips_order. Each cell is 'shift_order' orders deeper than the HEALPix
41 cell, with 2**shift_order x 2**shift_order sub-pixels on a side, which
42 defines the target resolution of the HiPS image. The IVOA recommends
43 shift_order=9, for 2**9=512 pixels on a side. See Notes below to
44 convert from hips_order to image resolution.
45
46 Parameters
47 ----------
48 hips_order : `int`
49 HiPS order, such that HEALPix nside=2**hips_order.
50 Must be a positive integer.
51 pixel : `int`
52 Pixel number in the nest ordering scheme.
53 shift_order : `int`, optional
54 Shift order for subpixels, such that there are 2**shift_order
55 sub-pixels on a side of the HiPS cell.
56 Must be a positive integer.
57
58 Returns
59 -------
60 wcs : `lsst.geom.SkyWcs`
61
62 Raises
63 ------
64 `ValueError`: Raise if hips_order is <=0, or if shift_order is <=0, or
65 if pixel number is out of range for the given hips_order
66 (0 <= pixel < 12*nside*nside).
67
68 Notes
69 -----
70 Table 5 from
71 https://www.ivoa.net/documents/HiPS/20170519/REC-HIPS-1.0-20170519.pdf
72 shows the relationship between hips_order, number of tiles (full
73 sky coverage), cell size, and sub-pixel size/image resolution (with
74 the default shift_order=9):
75
76 +------------+-----------------+--------------+------------------+
77 | hips_order | Number of Tiles | Cell Size | Image Resolution |
78 +============+=================+==============+==================+
79 | 0 | 12 | 58.63 deg | 6.871 arcmin |
80 | 1 | 48 | 29.32 deg | 3.435 arcmin |
81 | 2 | 192 | 14.66 deg | 1.718 arcmin |
82 | 3 | 768 | 7.329 deg | 51.53 arcsec |
83 | 4 | 3072 | 3.665 deg | 25.77 arcsec |
84 | 5 | 12288 | 1.832 deg | 12.88 arcsec |
85 | 6 | 49152 | 54.97 arcmin | 6.442 arcsec |
86 | 7 | 196608 | 27.48 arcmin | 3.221 arcsec |
87 | 8 | 786432 | 13.74 arcmin | 1.61 arcsec |
88 | 9 | 3145728 | 6.871 arcmin | 805.2mas |
89 | 10 | 12582912 | 3.435 arcmin | 402.6mas |
90 | 11 | 50331648 | 1.718 arcmin | 201.3mas |
91 | 12 | 201326592 | 51.53 arcsec | 100.6mas |
92 | 13 | 805306368 | 25.77 arcsec | 50.32mas |
93 +------------+-----------------+--------------+------------------+
94 """
95 if shift_order <= 0:
96 raise ValueError(f"shift_order {shift_order} must be positive.")
97 hips_tilepix = 2**shift_order
98
99 if hips_order <= 0:
100 raise ValueError(f"order {hips_order} must be positive.")
101 nside_cell = 2**hips_order
102
103 if pixel < 0 or pixel >= 12*nside_cell*nside_cell:
104 raise ValueError(f"pixel value {pixel} out of range.")
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122 nside_pix = nside_cell*hips_tilepix
123
124 cos45 = np.sqrt(2.0)/2.0
125
126 scale = 90.0/nside_pix/np.sqrt(2.0)
127 cos45_scale = cos45*scale
128
129
130
131
132 cent_ra_proj, cent_dec_proj = _hpx_projected_center(hips_order, pixel)
133
134 md = PropertySet()
135 md['CD1_1'] = -cos45_scale
136 md['CD1_2'] = -cos45_scale
137 md['CD2_1'] = cos45_scale
138 md['CD2_2'] = -cos45_scale
139 md['CTYPE1'] = 'RA---HPX'
140 md['CTYPE2'] = 'DEC--HPX'
141 md['CRVAL1'] = 0.0
142 md['CRVAL2'] = 0.0
143 md['PV2_1'] = 4
144 md['PV2_2'] = 3
145 md['CRPIX1'] = ((hips_tilepix + 1)/2.0) - 0.5*(-cent_ra_proj/cos45_scale + cent_dec_proj/cos45_scale)
146 md['CRPIX2'] = ((hips_tilepix + 1)/2.0) - 0.5*(-cent_ra_proj/cos45_scale - cent_dec_proj/cos45_scale)
147
148 return makeSkyWcs(md)
149
150
