Coverage for python / lsst / sphgeom / _healpixPixelization.py: 23%

1# This file is part of sphgeom. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (http://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 software is dual licensed under the GNU General Public License and also 

10# under a 3-clause BSD license. Recipients may choose which of these licenses 

11# to use; please see the files gpl-3.0.txt and/or bsd_license.txt, 

12# respectively. If you choose the GPL option then the following text applies 

13# (but note that there is still no warranty even if you opt for BSD instead): 

14# 

15# This program is free software: you can redistribute it and/or modify 

16# it under the terms of the GNU General Public License as published by 

17# the Free Software Foundation, either version 3 of the License, or 

18# (at your option) any later version. 

19# 

20# This program is distributed in the hope that it will be useful, 

21# but WITHOUT ANY WARRANTY; without even the implied warranty of 

22# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

23# GNU General Public License for more details. 

24# 

25# You should have received a copy of the GNU General Public License 

26# along with this program. If not, see <http://www.gnu.org/licenses/>. 

27__all__ = ["HealpixPixelization"] 

28 

29import hpgeom as hpg 

30import numpy as np 

31 

32from ._sphgeom import Box, Circle, ConvexPolygon, Ellipse, LonLat, RangeSet, Region, UnitVector3d 

33from .pixelization_abc import PixelizationABC 

34 

35 

36class HealpixPixelization(PixelizationABC): 

37 """HEALPix pixelization class. 

38 

39 Parameters 

40 ---------- 

41 level : `int` 

42 Pixelization level. HEALPix nside = 2**level 

43 """ 

44 

45 MAX_LEVEL = 17 

46 

47 def __init__(self, level: int): 

48 if level < 0: 

49 raise ValueError("HealPix level must be >= 0.") 

50 

51 self._level = level 

52 self._nside = 2**level 

53 

54 self._npix = hpg.nside_to_npixel(self._nside) 

55 

56 # Values used to do pixel/region intersections 

57 self._bit_shift = 8 

58 self._nside_highres = self._nside * (2 ** (self._bit_shift // 2)) 

59 

60 @property 

61 def nside(self): 

62 return self._nside 

63 

64 def getLevel(self): 

65 return self._level 

66 

67 level = property(getLevel) 

68 

69 def universe(self) -> RangeSet: 

70 return RangeSet(0, self._npix) 

71 

72 def pixel(self, i) -> Region: 

73 # This is arbitrarily returning 4 points on a side 

74 # to approximate the pixel shape. 

75 varr = hpg.angle_to_vector(*hpg.boundaries(self._nside, i, step=4)) 

76 return ConvexPolygon([UnitVector3d(*c) for c in varr]) 

77 

78 def index(self, v: UnitVector3d) -> int: 

79 return hpg.vector_to_pixel(self._nside, v.x(), v.y(), v.z()) 

80 

81 def toString(self, i: int) -> str: 

82 return str(i) 

83 

84 def envelope(self, region: Region, maxRanges: int = 0): 

85 if maxRanges > 0: 

86 # If this is important, the rangeset can be consolidated. 

87 raise NotImplementedError("HealpixPixelization: maxRanges not implemented") 

88 pixels_highres = self._interior_pixels_from_region(self._nside_highres, region) 

89 

90 # Dilate the high resolution pixels by one to ensure that the full 

91 # region is completely covered at high resolution. 

92 neighbors = hpg.neighbors(self._nside_highres, pixels_highres).ravel() 

93 # Shift back to the original resolution and uniquify 

94 pixels = np.unique(np.right_shift(neighbors[neighbors >= 0], self._bit_shift)) 

95 

96 return RangeSet(pixels) 

97 

98 def interior(self, region: Region, maxRanges: int = 0): 

99 if maxRanges > 0: 

100 # If this is important, the rangeset can be consolidated. 

101 raise NotImplementedError("HealpixPixelization: maxRanges not implemented") 

102 pixels = self._interior_pixels_from_region(self._nside, region) 

103 

104 # Check that the corners of the pixels are entirely enclosed in 

105 # the region 

106 

107 # Returns arrays [npixels, ncorners], where ncorners is 4. 

108 corners_lon, corners_lat = hpg.boundaries(self._nside, pixels, step=1, degrees=False) 

109 

110 corners_int = region.contains(corners_lon.ravel(), corners_lat.ravel()).reshape((len(pixels), 4)) 

111 interior = np.sum(corners_int, axis=1) == 4 

112 pixels = pixels[interior] 

113 

114 return RangeSet(pixels) 

115 

116 def _interior_pixels_from_region(self, nside: int, region: Region): 

117 """Get interior pixels from a region. 

118 

119 Parameters 

120 ---------- 

121 nside : `int` 

122 Healpix nside to retrieve interior pixels. 

123 region : `lsst.sphgeom.Region` 

124 Sphgeom region to find interior pixels. 

125 

126 Returns 

127 ------- 

128 pixels : `np.ndarray` 

129 Array of pixels at resolution nside, nest ordering. 

130 """ 

131 if isinstance(region, Circle): 

132 center = LonLat(region.getCenter()) 

133 radius_rad = region.getOpeningAngle().asRadians() 

134 if radius_rad > 0: 

135 pixels = hpg.query_circle( 

136 nside, 

137 center.getLon().asRadians(), 

138 center.getLat().asRadians(), 

139 radius_rad, 

140 degrees=False, 

141 inclusive=True, 

142 ) 

143 elif radius_rad == 0: 

144 pixels = hpg.angle_to_pixel( 

145 nside, 

146 center.getLon().asRadians(), 

147 center.getLat().asRadians(), 

148 degrees=False, 

149 ) 

150 else: 

151 raise ValueError("Invalid circle radius.") 

152 elif isinstance(region, ConvexPolygon): 

153 vertices = np.array([[v.x(), v.y(), v.z()] for v in region.getVertices()]) 

154 retry = False 

155 try: 

156 pixels = hpg.query_polygon_vec(nside, vertices, inclusive=True) 

157 except RuntimeError: 

158 # This happens if the polygon is determined to have a 

159 # degenerate corner (which can happen with rounding). 

160 # In this case, redo with a bounding circle. 

161 retry = True 

162 

163 if retry: 

164 pixels = self._interior_pixels_from_region(nside, region.getBoundingCircle()) 

165 elif isinstance(region, Box): 

166 pixels = hpg.query_box( 

167 nside, 

168 region.getLon().getA().asRadians(), 

169 region.getLon().getB().asRadians(), 

170 region.getLat().getA().asRadians(), 

171 region.getLat().getB().asRadians(), 

172 degrees=False, 

173 inclusive=True, 

174 ) 

175 elif isinstance(region, Ellipse): 

176 # hpgeom supports query_ellipse given center, alpha, beta, 

177 # and orientation. However, until we figure out how to get 

178 # the orientation out of the Ellipse region, we will use the 

179 # bounding circle as was done with healpy. 

180 _circle = region.getBoundingCircle() 

181 center = LonLat(_circle.getCenter()) 

182 pixels = hpg.query_circle( 

183 nside, 

184 center.getLon().asRadians(), 

185 center.getLat().asRadians(), 

186 _circle.getOpeningAngle().asRadians(), 

187 degrees=False, 

188 inclusive=True, 

189 ) 

190 else: 

191 raise ValueError("Invalid region.") 

192 

193 return pixels 

194 

195 def __eq__(self, other): 

196 if isinstance(other, HealpixPixelization): 

197 return self._level == other._level 

198 

199 def __repr__(self): 

200 return f"HealpixPixelization({self._level})" 

201 

202 def __reduce__(self): 

203 return (self.__class__, (self._level,))