Coverage for tests / test_HealpixPixelization.py: 19%

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 

28import pickle 

29import unittest 

30 

31import hpgeom as hpg 

32import numpy as np 

33 

34try: 

35 import yaml 

36except ImportError: 

37 yaml = None 

38 

39from lsst.sphgeom import Angle, Box, Circle, ConvexPolygon, Ellipse, HealpixPixelization, LonLat, UnitVector3d 

40 

41 

42class HealpixPixelizationTestCase(unittest.TestCase): 

43 """Test HEALPix pixelization.""" 

44 

45 def test_construction(self): 

46 """Test construction of a HealpixPixelization.""" 

47 with self.assertRaises(ValueError): 

48 HealpixPixelization(-1) 

49 h1 = HealpixPixelization(5) 

50 self.assertEqual(h1.level, 5) 

51 self.assertEqual(h1.getLevel(), 5) 

52 self.assertEqual(h1.nside, 32) 

53 h2 = HealpixPixelization(6) 

54 h3 = HealpixPixelization(h2.level) 

55 self.assertNotEqual(h1, h2) 

56 self.assertEqual(h2, h3) 

57 

58 def test_indexing(self): 

59 """Test indexing of HealpixPixelization.""" 

60 h = HealpixPixelization(5) 

61 vec = UnitVector3d(1, 1, 1) 

62 lonlat = LonLat(vec) 

63 pix = hpg.angle_to_pixel(h.nside, lonlat.getLon().asDegrees(), lonlat.getLat().asDegrees()) 

64 self.assertEqual(h.index(UnitVector3d(1, 1, 1)), pix) 

65 

66 def test_pixel(self): 

67 """Test pixel polygon of HealpixPixelization.""" 

68 h = HealpixPixelization(5) 

69 pix_poly = h.pixel(10) 

70 self.assertIsInstance(pix_poly, ConvexPolygon) 

71 

72 def test_envelope(self): 

73 """Test envelope method of HealpixPixelization.""" 

74 # Make the hardest intersection: a region that _just_ 

75 # touches a healpix pixel. 

76 h = HealpixPixelization(5) 

77 pix = hpg.angle_to_pixel(h.nside, 50.0, 20.0) 

78 

79 corners_ra, corners_dec = hpg.boundaries(h.nside, pix, step=1) 

80 

81 # Take the southernmost corner... 

82 smost = np.argmin(corners_dec) 

83 

84 # Choose challenging comparison box corners: 

85 ra_range = np.array([corners_ra[smost] - 0.5, corners_ra[smost] + 0.5]) 

86 dec_range = np.array([corners_dec[smost] - 0.5, corners_dec[smost] + 1e-8]) 

87 

88 # Test the box region 

89 box = Box( 

90 point1=LonLat.fromDegrees(ra_range[0], dec_range[0]), 

91 point2=LonLat.fromDegrees(ra_range[1], dec_range[1]), 

92 ) 

93 # These pixels have been checked to completely overlap the region 

94 self._check_envelope(h, box, [98, 99, 104, 105]) 

95 

96 # Try a polygon region: 

97 poly = ConvexPolygon( 

98 [ 

99 UnitVector3d(LonLat.fromDegrees(ra_range[0], dec_range[0])), 

100 UnitVector3d(LonLat.fromDegrees(ra_range[1], dec_range[0])), 

101 UnitVector3d(LonLat.fromDegrees(ra_range[1], dec_range[1])), 

102 UnitVector3d(LonLat.fromDegrees(ra_range[0], dec_range[1])), 

103 UnitVector3d(LonLat.fromDegrees(ra_range[0], dec_range[0])), 

104 ] 

105 ) 

106 self._check_envelope(h, poly, [98, 99, 104, 105]) 

107 

108 # Try a circle region 

109 circle = Circle( 

110 center=UnitVector3d( 

111 LonLat.fromDegrees((ra_range[0] + ra_range[1]) / 2.0, (dec_range[0] + dec_range[1]) / 2.0) 

112 ), 

113 angle=Angle.fromDegrees((dec_range[1] - dec_range[0]) / 2.0), 

114 ) 

115 self._check_envelope(h, circle, [98, 99, 104, 105]) 

116 

117 # Try a circle region with zero-radius (representing a point) 

118 pt_circle = Circle( 

119 center=UnitVector3d( 

120 LonLat.fromDegrees((ra_range[0] + ra_range[1]) / 2.0, (dec_range[0] + dec_range[1]) / 2.0) 

121 ), 

122 ) 

123 self._check_envelope(h, pt_circle, [98]) 

124 

125 # Try an ellipse region 

126 circle = Ellipse( 

127 center=UnitVector3d( 

128 LonLat.fromDegrees((ra_range[0] + ra_range[1]) / 2.0, (dec_range[0] + dec_range[1]) / 2.0) 

129 ), 

130 alpha=Angle.fromDegrees(1.5 * (dec_range[1] - dec_range[0]) / 2.0), 

131 beta=Angle.fromDegrees(0.5 * (dec_range[1] - dec_range[0]) / 2.0), 

132 orientation=Angle.fromDegrees(45.0), 

133 ) 

134 self._check_envelope(h, circle, [98, 99, 104, 105]) 

135 

136 def test_envelope_missing_neighbors(self): 

137 """Test envelope, with a pixel with missing neighbors. 

138 

139 Testing DM-47043. 

140 """ 

141 h = HealpixPixelization(0) 

142 h_highres = HealpixPixelization(4) 

143 

144 self.assertEqual(h._nside_highres, h_highres._nside) 

145 

146 # Find a pixel with a missing neighbor at high res. 

147 n = hpg.neighbors(h_highres._nside, np.arange(hpg.nside_to_npixel(h_highres._nside))) 

148 ind = np.argmin(n) 

149 self.assertEqual(n.ravel()[ind], -1) 

150 

151 # This is the pixel with a bad (-1) neighbor. 

152 badpix = ind // 8 

153 ra, dec = hpg.pixel_to_angle(h_highres._nside, badpix) 

154 

155 box = Box( 

156 point1=LonLat.fromDegrees(ra - 0.1, dec - 0.1), 

157 point2=LonLat.fromDegrees(ra + 0.1, dec + 0.1), 

158 ) 

159 

160 self._check_envelope(h, box, [0, 1, 5]) 

161 

162 def test_envelope_small(self): 

163 """Test envelope, DM-54274.""" 

164 h = HealpixPixelization(3) 

165 

166 # Test the box region 

167 box = Box( 

168 point1=LonLat.fromDegrees(0.6, 0.6), 

169 point2=LonLat.fromDegrees(0.61, 0.61), 

170 ) 

171 self._check_envelope(h, box, [282, 304]) 

172 

173 # Try a polygon region: 

174 poly = ConvexPolygon( 

175 [ 

176 UnitVector3d(-0.9879246204440798, -0.08691162558237045, 0.12826267445773423), 

177 UnitVector3d(-0.9885269066487583, -0.08474209761828523, 0.1250333224492182), 

178 UnitVector3d(-0.9885167542515381, -0.08795173158697882, 0.12287847441621431), 

179 UnitVector3d(-0.9879142847819917, -0.09012140316514382, 0.1261090742779096), 

180 ], 

181 ) 

182 

183 self._check_envelope(h, poly, [433]) 

184 

185 # Try a circle region 

186 circle = Circle( 

187 center=UnitVector3d(LonLat.fromDegrees(0.6, 0.6)), 

188 angle=Angle.fromDegrees(0.01), 

189 ) 

190 self._check_envelope(h, circle, [282, 304]) 

191 

192 # Try an ellipse region 

193 circle = Ellipse( 

194 center=UnitVector3d(LonLat.fromDegrees(0.6, 0.6)), 

195 alpha=Angle.fromDegrees(0.02), 

196 beta=Angle.fromDegrees(0.005), 

197 orientation=Angle.fromDegrees(45.0), 

198 ) 

199 self._check_envelope(h, circle, [282, 304]) 

200 

201 def _check_envelope(self, pixelization, region, check_pixels): 

202 """Check the envelope from a region. 

203 

204 Parameters 

205 ---------- 

206 pixelization : `lsst.sphgeom.HealpixPixelization` 

207 region : `lsst.sphgeom.Region` 

208 check_pixels : `list` [`int`] 

209 """ 

210 pixel_range = pixelization.envelope(region) 

211 

212 pixels = [] 

213 for r in pixel_range.ranges(): 

214 pixels.extend(range(r[0], r[1])) 

215 

216 self.assertEqual(pixels, check_pixels) 

217 

218 def test_interior(self): 

219 """Test interior method of HealpixPixelization.""" 

220 h = HealpixPixelization(5) 

221 pix = hpg.angle_to_pixel(h.nside, 50.0, 20.0) 

222 

223 corners_ra, corners_dec = hpg.boundaries(h.nside, pix, step=1) 

224 

225 ra_range = np.array([corners_ra.min() - 1.0, corners_ra.max() + 1.0]) 

226 dec_range = np.array([corners_dec.min() - 1.0, corners_dec.max() + 1.0]) 

227 

228 # Test the box region 

229 box = Box( 

230 point1=LonLat.fromDegrees(ra_range[0], dec_range[0]), 

231 point2=LonLat.fromDegrees(ra_range[1], dec_range[1]), 

232 ) 

233 # These pixels have been checked to completely overlap the region 

234 self._check_interior(h, box, [pix]) 

235 

236 # Try a polygon region: 

237 poly = ConvexPolygon( 

238 [ 

239 UnitVector3d(LonLat.fromDegrees(ra_range[0], dec_range[0])), 

240 UnitVector3d(LonLat.fromDegrees(ra_range[1], dec_range[0])), 

241 UnitVector3d(LonLat.fromDegrees(ra_range[1], dec_range[1])), 

242 UnitVector3d(LonLat.fromDegrees(ra_range[0], dec_range[1])), 

243 UnitVector3d(LonLat.fromDegrees(ra_range[0], dec_range[0])), 

244 ] 

245 ) 

246 self._check_interior(h, poly, [pix]) 

247 

248 # Try a circle region 

249 circle = Circle( 

250 center=UnitVector3d( 

251 LonLat.fromDegrees((ra_range[0] + ra_range[1]) / 2.0, (dec_range[0] + dec_range[1]) / 2.0) 

252 ), 

253 angle=Angle.fromDegrees(2.5), 

254 ) 

255 self._check_interior(h, circle, [pix]) 

256 

257 # Try an ellipse region 

258 ellipse = Ellipse( 

259 center=UnitVector3d( 

260 LonLat.fromDegrees((ra_range[0] + ra_range[1]) / 2.0, (dec_range[0] + dec_range[1]) / 2.0) 

261 ), 

262 alpha=Angle.fromDegrees(1.5), 

263 beta=Angle.fromDegrees(2.5), 

264 orientation=Angle.fromDegrees(45.0), 

265 ) 

266 self._check_interior(h, ellipse, [pix]) 

267 

268 def test_interior_nearly_degenerate_polygon(self): 

269 """Check interior pixels with nearly degenerate polygon. 

270 

271 See DM-53933. 

272 """ 

273 poly = ConvexPolygon( 

274 [ 

275 UnitVector3d(-0.016144196513206758, -0.9235039286924489, -0.3832490816799893), 

276 UnitVector3d(-0.01988138785117867, -0.9228560783513854, -0.38463149775728545), 

277 UnitVector3d(-0.022348863683228713, -0.9209177038528237, -0.3891158066983548), 

278 UnitVector3d(-0.024841110338799998, -0.918545770328335, -0.3945333788782152), 

279 UnitVector3d(-0.026637426968850634, -0.9164294757774608, -0.39930873194901123), 

280 UnitVector3d(-0.02509225617228684, -0.9149335621264763, -0.4028237276709774), 

281 UnitVector3d(-0.0221015078570923, -0.911998223315482, -0.4095982959191202), 

282 UnitVector3d(-0.020553220684181986, -0.9104571830599119, -0.4130923418972051), 

283 UnitVector3d(-0.017371565494301272, -0.9072522623419655, -0.4202279871541906), 

284 UnitVector3d(-0.014347349160349979, -0.9041551644174409, -0.4269632211670476), 

285 UnitVector3d(-0.014303105070771129, -0.9041094310633312, -0.42706153871271785), 

286 UnitVector3d(-0.012802279210872009, -0.9025514012895127, -0.4303917630221838), 

287 UnitVector3d(-0.007862275539177377, -0.9015265103248066, -0.43265244227315125), 

288 UnitVector3d(-0.0019484762556141095, -0.900471339489277, -0.4349109911219402), 

289 UnitVector3d(0.0032282853511799436, -0.8996943229589627, -0.436508537613075), 

290 UnitVector3d(0.006969847869770022, -0.9003377930402403, -0.4351359323752772), 

291 UnitVector3d(0.022099717156107892, -0.9027944227880381, -0.42950417074160424), 

292 UnitVector3d(0.033589177231523826, -0.9045049275054283, -0.4251383342999184), 

293 UnitVector3d(0.03732521827078148, -0.9050320420971739, -0.4237025263772424), 

294 UnitVector3d(0.039795546158565454, -0.9070052460628258, -0.41923477684998206), 

295 UnitVector3d(0.04228779372661228, -0.9093772751891867, -0.41381724694752126), 

296 UnitVector3d(0.04409395578698435, -0.9114581951867685, -0.4090228373696683), 

297 UnitVector3d(0.04255005806475871, -0.9130743444969162, -0.4055671756691021), 

298 UnitVector3d(0.04109145924730725, -0.91458405813035, -0.4023027374885082), 

299 UnitVector3d(0.03791490987091337, -0.9178181221381234, -0.3951864803916365), 

300 UnitVector3d(0.03335712714052637, -0.9223387473238281, -0.38493965404208763), 

301 UnitVector3d(0.03180762744041508, -0.9238417869862389, -0.3814506881035668), 

302 UnitVector3d(0.030248610691484885, -0.9253347742243898, -0.3779425580195124), 

303 UnitVector3d(0.02530895801451361, -0.9263964395908428, -0.3756981412751856), 

304 UnitVector3d(0.01939515964342646, -0.9274517060282444, -0.37343963470379626), 

305 UnitVector3d(0.014215536919345694, -0.9281937989444676, -0.37182548340738003), 

306 UnitVector3d(0.010475126504663994, -0.9276710642902235, -0.3732514811803906), 

307 ], 

308 ) 

309 

310 h = HealpixPixelization(6) 

311 self._check_interior(h, poly, [28890, 28891, 28912, 28913, 28914]) 

312 

313 def _check_interior(self, pixelization, region, check_pixels): 

314 """Check the interior from a region. 

315 

316 Parameters 

317 ---------- 

318 pixelization : `lsst.sphgeom.HealpixPixelization` 

319 region : `lsst.sphgeom.Region` 

320 check_pixels : `list` [`int`] 

321 """ 

322 pixel_range = pixelization.interior(region) 

323 

324 pixels = [] 

325 for r in pixel_range.ranges(): 

326 pixels.extend(range(r[0], r[1])) 

327 

328 self.assertEqual(pixels, check_pixels) 

329 

330 def test_index_to_string(self): 

331 """Test converting index to string of HealpixPixelization.""" 

332 h = HealpixPixelization(5) 

333 self.assertEqual(h.toString(0), str(0)) 

334 self.assertEqual(h.toString(100), str(100)) 

335 

336 def test_string(self): 

337 """Test string representation of HealpixPixelization.""" 

338 h = HealpixPixelization(5) 

339 self.assertEqual(str(h), "HealpixPixelization(5)") 

340 self.assertEqual(str(h), repr(h)) 

341 self.assertEqual(h, eval(repr(h), {"HealpixPixelization": HealpixPixelization})) 

342 

343 def test_pickle(self): 

344 """Test pickling of HealpixPixelization.""" 

345 a = HealpixPixelization(5) 

346 b = pickle.loads(pickle.dumps(a)) 

347 self.assertEqual(a, b) 

348 

349 @unittest.skipIf(not yaml, "YAML module can not be imported") 

350 def test_yaml(self): 

351 """Test yaml representation of HealpixPixelization.""" 

352 a = HealpixPixelization(5) 

353 b = yaml.safe_load(yaml.dump(a)) 

354 self.assertEqual(a, b) 

355 

356 

357if __name__ == "__main__": 

358 unittest.main()