Coverage for tests/test_ringsSkyMap.py: 98%

1import unittest 

2import math 

3import numpy as np 

4 

5import lsst.utils.tests 

6import lsst.geom 

7 

8from lsst.skymap.ringsSkyMap import RingsSkyMap 

9from helper import skyMapTestCase 

10 

11 

12class RingsTestCase(skyMapTestCase.SkyMapTestCase): 

13 

14 def setUp(self): 

15 config = RingsSkyMap.ConfigClass() 

16 config.numRings = 3 

17 self.setAttributes( 

18 SkyMapClass=RingsSkyMap, 

19 name="rings", 

20 config=config, 

21 numTracts=26, 

22 neighborAngularSeparation=None, # no uniform tract separation 

23 numNeighbors=None, # ignored because neighborAngularSeparation=None 

24 ) 

25 

26 def testPoles(self): 

27 """Test that findAllTracts behaves at the poles 

28 

29 Testing fix to DM-10686. 

30 """ 

31 skymap = self.getSkyMap() 

32 for ra in (0, 123, 321, 359.9): 

33 tracts = skymap.findAllTracts(lsst.geom.SpherePoint(ra, 90, lsst.geom.degrees)) 

34 self.assertListEqual(tracts, [skymap[len(skymap) - 1]]) 

35 tracts = skymap.findAllTracts(lsst.geom.SpherePoint(ra, -90, lsst.geom.degrees)) 

36 self.assertListEqual(tracts, [skymap[0]]) 

37 

38 def testSha1Compare(self): 

39 """Test that RingsSkyMap's extra state is included in its hash.""" 

40 defaultSkyMap = self.getSkyMap() 

41 for numRings in (4, 5): 

42 config = self.getConfig() 

43 config.numRings = numRings 

44 skyMap = self.getSkyMap(config=config) 

45 self.assertNotEqual(skyMap, defaultSkyMap) 

46 for raStart in (60.0, 75.0): 

47 config = self.getConfig() 

48 config.raStart = raStart 

49 skyMap = self.getSkyMap(config=config) 

50 self.assertNotEqual(skyMap, defaultSkyMap) 

51 

52 def testCorners(self): 

53 """Test that corners of a tract can be found in the tract""" 

54 skymap = self.getSkyMap() 

55 for tract in skymap: 

56 vertices = tract.getVertexList() 

57 for coord in vertices: 

58 self.assertIn(tract.getId(), [tt.getId() for tt in skymap.findAllTracts(coord)]) 

59 

60 def testInnerSkyRegion(self): 

61 """Test that the inner_sky_region covers the inner sky region.""" 

62 skymap = self.getSkyMap() 

63 for tract in skymap: 

64 innerSkyRegion = tract.inner_sky_region 

65 bbox = tract.bbox 

66 wcs = tract.wcs 

67 x = np.linspace(bbox.getMinX(), bbox.getMaxX(), 100) 

68 y = np.linspace(bbox.getMinY(), bbox.getMaxY(), 100) 

69 xx, yy = np.meshgrid(x, y) 

70 

71 ra, dec = wcs.pixelToSkyArray(xx.ravel(), yy.ravel()) 

72 np.testing.assert_array_equal( 

73 innerSkyRegion.contains(ra, dec), 

74 skymap.findTractIdArray(ra, dec) == tract.getId(), 

75 ) 

76 

77 

78class NonzeroRaStartRingsTestCase(RingsTestCase): 

79 """Test that setting raStart != 0 works""" 

80 def getConfig(self): 

81 config = super().getConfig() 

82 config.raStart = 234 

83 return config 

84 

85 

86class HscRingsTestCase(lsst.utils.tests.TestCase): 

87 def getConfig(self): 

88 """Return a configuration matching that used for the HSC SSP""" 

89 config = RingsSkyMap.ConfigClass() 

90 config.numRings = 120 

91 config.projection = "TAN" 

92 config.tractOverlap = 1.0/60 # Overlap between tracts (degrees) 

93 config.pixelScale = 0.168 

94 return config 

95 

96 def setUp(self): 

97 self.skymap = RingsSkyMap(self.getConfig()) 

98 

99 def tearDown(self): 

100 del self.skymap 

101 

102 def testDm7770(self): 

103 """Test that DM-7770 has been fixed 

104 

105 These operations previously caused: 

106 lsst::pex::exceptions::RuntimeError: 'Error: wcslib 

107 returned a status code of 9 at sky 30.18, -3.8 deg: 

108 One or more of the world coordinates were invalid' 

109 

110 We are only testing function, and not the actual results. 

111 """ 

112 coordList = [lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees) for 

113 ra, dec in [(30.18, -3.8), (31.3, -3.8), (31.3, -2.7), (30.18, -2.7)]] 

114 for coord in coordList: 

115 self.skymap.findAllTracts(coord) 

116 self.skymap.findTractPatchList(coordList) 

117 

118 def testDm14809(self): 

119 """Test that DM-14809 has been fixed""" 

120 skyMapTestCase.checkDm14809(self, self.skymap) 

121 

122 # Check that the first tract in the last ring exists 

123 coord = self.getFirstTractLastRingCoord() 

124 tract = self.skymap.findTract(coord) 

125 self.assertTrue(tract.contains(coord)) 

126 

127 tractId = self.skymap.findTractIdArray(coord.getLongitude().asRadians(), 

128 coord.getLatitude().asRadians(), 

129 degrees=False) 

130 self.assertEqual(tractId[0], tract.getId()) 

131 

132 def testWraparound(self): 

133 """Check wrapping at RA=0 

134 

135 How-to-reproduce of a bug identified by Sogo Mineo. 

136 """ 

137 tractId = 9712 

138 deviation = 10 / 3600.0 # 10 arcsec 

139 tract = self.skymap[tractId] 

140 center = tract.getCtrCoord() 

141 centerRa = center.getRa().asDegrees() 

142 centerDec = center.getDec().asDegrees() 

143 for devRa in [-deviation, deviation]: 

144 coord = lsst.geom.SpherePoint(centerRa + devRa, centerDec, lsst.geom.degrees) 

145 foundTractId = self.skymap.findTract(coord).getId() 

146 self.assertEqual(tractId, foundTractId) 

147 

148 foundTractArrayId = self.skymap.findTractIdArray([centerRa + devRa], 

149 [centerDec], 

150 degrees=True) 

151 self.assertEqual(tractId, foundTractArrayId[0]) 

152 

153 def testFindTractIdArray(self): 

154 """Test findTractIdArray. 

155 

156 Test an array of positions to ensure that ``findTract`` and 

157 ``findTractIdArray`` give the same answers. 

158 """ 

159 np.random.seed(12345) 

160 

161 ras = np.random.uniform(low=0.0, high=360.0, size=1000) 

162 decs = np.random.uniform(low=-90.0, high=90.0, size=1000) 

163 

164 coords = [lsst.geom.SpherePoint(ra*lsst.geom.degrees, dec*lsst.geom.degrees) 

165 for ra, dec in zip(ras, decs)] 

166 

167 indexes = [self.skymap.findTract(coord).getId() for coord in coords] 

168 indexes2 = self.skymap.findTractIdArray(ras, decs, degrees=True) 

169 

170 np.testing.assert_array_equal(indexes2, indexes) 

171 

172 def testFindTractIdPatchIdArray(self): 

173 """Test findTractIdPatchIdArray.""" 

174 np.random.seed(12345) 

175 

176 ras = np.random.uniform(low=0.0, high=360.0, size=1000) 

177 decs = np.random.uniform(low=-90.0, high=90.0, size=1000) 

178 

179 coords = [lsst.geom.SpherePoint(ra*lsst.geom.degrees, dec*lsst.geom.degrees) 

180 for ra, dec in zip(ras, decs)] 

181 

182 tractIds = [self.skymap.findTract(coord).getId() for coord in coords] 

183 patchIds = np.zeros(len(tractIds), dtype=np.int32) 

184 for i, tractId in enumerate(tractIds): 

185 tract = self.skymap[tractId] 

186 patch = tract.findPatch(coords[i]) 

187 patchIds[i] = patch.sequential_index 

188 

189 tractIds2, patchIds2 = self.skymap.findTractIdPatchIdArray(ras, decs, degrees=True) 

190 

191 np.testing.assert_array_equal(tractIds2, tractIds) 

192 np.testing.assert_array_equal(patchIds2, patchIds) 

193 

194 def getFirstTractLastRingCoord(self): 

195 """Return the coordinates of the first tract in the last ring 

196 

197 This tract is missing in version=0, but this is fixed in version=1.in 

198 """ 

199 ringNum = self.skymap.config.numRings - 1 

200 ringSize = math.pi/(self.skymap.config.numRings + 1) 

201 firstRingStart = ringSize*0.5 - 0.5*math.pi 

202 dec = ringNum*ringSize + firstRingStart 

203 return lsst.geom.SpherePoint(self.skymap.config.raStart*lsst.geom.degrees, 

204 dec*lsst.geom.radians) 

205 

206 

207class Version0HscRingsTestCase(HscRingsTestCase): 

208 """Testing that the version=0 RingsSkyMap works in the expected way""" 

209 def setUp(self): 

210 self.skymap = RingsSkyMap(self.getConfig(), version=0) 

211 

212 def testDm14809(self): 

213 """Test that DM-14809 has been partially fixed 

214 

215 The observed behaviour was: 

216 

217 skyMap.findTract(skyMap[9712].getCtrCoord()).getId() != 9712 

218 

219 and 

220 

221 skyMap[1].getCtrCoord() == skyMap[11].getCtrCoord() 

222 

223 Specifically for version=0, we fixed the ``findTract`` behaviour but 

224 left the tract duplication (tract 11 duplicates tract 1) and the 

225 missing tract (the first tract in the last ring) so that the tract 

226 numbering would remain unchanged. 

227 """ 

228 # Check that the tract found for central coordinate of a tract is that 

229 # tract. 

230 expect = [tract.getId() for tract in self.skymap] 

231 expect[self.skymap._ringNums[0] + 1] = 1 # Due to the bug 

232 got = [self.skymap.findTract(tract.getCtrCoord()).getId() for tract in self.skymap] 

233 self.assertEqual(got, expect) 

234 

235 # Check that the tract central coordinates are unique 

236 # Round to integer arcminutes so differences are relatively immune to 

237 # small numerical inaccuracies. 

238 centers = set([(int(coord.getRa().asArcminutes()), int(coord.getDec().asArcminutes())) for 

239 coord in (tract.getCtrCoord() for tract in self.skymap)]) 

240 self.assertEqual(len(centers), len(self.skymap) - 1) # One tract is duplicated 

241 self.assertEqual(self.skymap[1].getCtrCoord(), 

242 self.skymap[self.skymap._ringNums[0] + 1].getCtrCoord()) # This is the duplicate 

243 

244 # Check that some particular tracts we know and love haven't moved 

245 degrees = lsst.geom.degrees 

246 # 9712 is at RA=0, and was identified as problematic in DM-14809 

247 self.assertEqual(self.skymap[9712].getCtrCoord(), 

248 lsst.geom.SpherePoint(0.0*degrees, 0.7438016528925696*degrees)) 

249 # The Cosmos field 

250 self.assertEqual(self.skymap[9813].getCtrCoord(), 

251 lsst.geom.SpherePoint(150.2479338842975*degrees, 2.2314049586776834*degrees)) 

252 

253 # Check that the first tract in the last ring does NOT exist 

254 # (due to the bug). 

255 coord = self.getFirstTractLastRingCoord() 

256 tract = self.skymap.findTract(coord) 

257 self.assertFalse(tract.contains(coord)) 

258 

259 

260class MemoryTester(lsst.utils.tests.MemoryTestCase): 

261 pass 

262 

263 

264def setup_module(module): 

265 lsst.utils.tests.init() 

266 

267 

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

269 lsst.utils.tests.init() 

270 unittest.main()