Coverage for tests / test_tableUtils.py: 16%

1# LSST Data Management System 

2# Copyright 2016 LSST Corporation. 

3# 

4# This product includes software developed by the 

5# LSST Project (http://www.lsst.org/). 

6# 

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

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

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

10# (at your option) any later version. 

11# 

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

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

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

15# GNU General Public License for more details. 

16# 

17# You should have received a copy of the LSST License Statement and 

18# the GNU General Public License along with this program. If not, 

19# see <http://www.lsstcorp.org/LegalNotices/>. 

20# 

21# The classes in this test are a little non-standard to reduce code 

22# duplication and support automated unittest discovery. 

23# A base class includes all the code that implements the testing and 

24# itself inherits from unittest.TestCase. unittest automated discovery 

25# will scan all classes that inherit from unittest.TestCase and invoke 

26# any test methods found. To prevent this base class from being executed 

27# the test methods are placed in a different class that does not inherit 

28# from unittest.TestCase. The actual test classes then inherit from 

29# both the testing class and the implementation class allowing test 

30# discovery to only run tests found in the subclasses. 

31 

32import math 

33import unittest 

34 

35import numpy as np 

36 

37import lsst.utils.tests 

38import lsst.geom 

39import lsst.afw.geom as afwGeom 

40import lsst.afw.table as afwTable 

41 

42 

43class UpdateTestCase(lsst.utils.tests.TestCase): 

44 """A test case for the lsst.afw.table.updateRefCentroids and updateSourceCoords 

45 """ 

46 

47 def setUp(self): 

48 self.crval = lsst.geom.SpherePoint(44.0, 45.0, lsst.geom.degrees) 

49 self.crpix = lsst.geom.Point2D(15000, 4000) 

50 

51 arcsecPerPixel = 1/3600.0 

52 cdMatrix = afwGeom.makeCdMatrix(arcsecPerPixel * lsst.geom.arcseconds) 

53 self.wcs = afwGeom.makeSkyWcs(crval=self.crval, crpix=self.crpix, cdMatrix=cdMatrix) 

54 

55 refSchema = afwTable.SimpleTable.makeMinimalSchema() 

56 self.refCentroidKey = afwTable.Point2DKey.addFields( 

57 refSchema, "centroid", "centroid", "pixels") 

58 self.refCoordKey = afwTable.CoordKey(refSchema["coord"]) 

59 self.refHasCentroidKey = refSchema.addField("hasCentroid", type="Flag") 

60 self.refCat = afwTable.SimpleCatalog(refSchema) 

61 

62 # an alias is required to make src.getCentroid() work; 

63 # simply defining a field named "slot_Centroid" doesn't suffice 

64 srcSchema = afwTable.SourceTable.makeMinimalSchema() 

65 self.srcCentroidKey = afwTable.Point2DKey.addFields(srcSchema, "base_SdssCentroid", 

66 "centroid", "pixels") 

67 self.srcCentroidErrKey = afwTable.CovarianceMatrix2fKey.addFields(srcSchema, "base_SdssCentroid", 

68 ["x", "y"], "pixels") 

69 self.srcCoordErrKey = afwTable.CoordKey.addErrorFields(srcSchema) 

70 

71 srcAliases = srcSchema.getAliasMap() 

72 srcAliases.set("slot_Centroid", "base_SdssCentroid") 

73 self.srcCoordKey = afwTable.CoordKey(srcSchema["coord"]) 

74 self.sourceCat = afwTable.SourceCatalog(srcSchema) 

75 

76 def tearDown(self): 

77 del self.wcs 

78 del self.refCat 

79 del self.sourceCat 

80 

81 def testNull(self): 

82 """Check that an empty list causes no problems for either function""" 

83 afwTable.updateRefCentroids(self.wcs, []) 

84 afwTable.updateSourceCoords(self.wcs, []) 

85 

86 def testRefCenter(self): 

87 """Check that a ref obj at the center is handled as expected""" 

88 refObj = self.refCat.addNew() 

89 refObj.set(self.refCoordKey, self.crval) 

90 

91 # initial centroid should be nan and hasCentroid False 

92 nanRefCentroid = self.refCat[0].get(self.refCentroidKey) 

93 for val in nanRefCentroid: 

94 self.assertTrue(math.isnan(val)) 

95 self.assertFalse(self.refCat[0].get(self.refHasCentroidKey)) 

96 

97 # computed centroid should be crpix and hasCentroid True 

98 afwTable.updateRefCentroids(self.wcs, self.refCat) 

99 refCentroid = self.refCat[0].get(self.refCentroidKey) 

100 self.assertPairsAlmostEqual(refCentroid, self.crpix) 

101 self.assertTrue(self.refCat[0].get(self.refHasCentroidKey)) 

102 

103 # coord should not be changed 

104 self.assertEqual(self.refCat[0].get(self.refCoordKey), self.crval) 

105 

106 def testSourceCenter(self): 

107 """Check that a source at the center is handled as expected""" 

108 src = self.sourceCat.addNew() 

109 src.set(self.srcCentroidKey, self.crpix) 

110 

111 # initial coord should be nan; as a sanity-check 

112 nanSourceCoord = self.sourceCat[0].get(self.srcCoordKey) 

113 for val in nanSourceCoord: 

114 self.assertTrue(math.isnan(val)) 

115 

116 # compute coord should be crval 

117 afwTable.updateSourceCoords(self.wcs, self.sourceCat) 

118 srcCoord = self.sourceCat[0].get(self.srcCoordKey) 

119 self.assertPairsAlmostEqual(srcCoord, self.crval) 

120 

121 # centroid should not be changed; also make sure that getCentroid words 

122 self.assertEqual(self.sourceCat[0].getCentroid(), self.crpix) 

123 

124 def testLists(self): 

125 """Check updating lists of reference objects and sources""" 

126 # arbitrary but reasonable values that are intentionally different than 

127 # testCatalogs 

128 maxPix = 1000 

129 numPoints = 10 

130 self.setCatalogs(maxPix=maxPix, numPoints=numPoints) 

131 

132 # update the catalogs as lists 

133 afwTable.updateSourceCoords(self.wcs, [s for s in self.sourceCat]) 

134 afwTable.updateRefCentroids(self.wcs, [r for r in self.refCat]) 

135 

136 self.checkCatalogs() 

137 

138 def testCatalogs(self): 

139 """Check updating catalogs of reference objects and sources""" 

140 # arbitrary but reasonable values that are intentionally different than 

141 # testLists 

142 maxPix = 2000 

143 numPoints = 9 

144 self.setCatalogs(maxPix=maxPix, numPoints=numPoints) 

145 

146 # update the catalogs 

147 afwTable.updateSourceCoords(self.wcs, self.sourceCat) 

148 afwTable.updateRefCentroids(self.wcs, self.refCat) 

149 

150 # check that centroids and coords match 

151 self.checkCatalogs() 

152 

153 def testCoordErrors(self): 

154 """Check that updateSourceCoords has correctly propagated the centroid 

155 errors. 

156 """ 

157 maxPix = 2000 

158 numPoints = 10 

159 

160 self.setCatalogs(maxPix=maxPix, numPoints=numPoints) 

161 scale = (1.0 * lsst.geom.arcseconds).asDegrees() 

162 # update the catalogs 

163 afwTable.updateSourceCoords(self.wcs, self.sourceCat) 

164 factor = math.pi/180. 

165 factor_sq = factor**2 

166 for src in self.sourceCat: 

167 center = src.get(self.srcCentroidKey) 

168 skyCenter = self.wcs.pixelToSky(center) 

169 localGnomonicWcs = lsst.afw.geom.makeSkyWcs( 

170 center, skyCenter, np.diag((scale, scale))) 

171 measurementToLocalGnomonic = self.wcs.getTransform().then( 

172 localGnomonicWcs.getTransform().inverted() 

173 ) 

174 localMatrix = measurementToLocalGnomonic.getJacobian(center) 

175 radMatrix = np.radians(localMatrix / 3600) 

176 

177 centroidErr = src.get(self.srcCentroidErrKey) 

178 coordErr = radMatrix.dot(centroidErr.dot(radMatrix.T)) 

179 catCoordErr = src.get(self.srcCoordErrKey) 

180 np.testing.assert_allclose(coordErr, catCoordErr, rtol=1e-10, atol=1e-10) 

181 

182 (ra, dec), (ra_err, dec_err, ra_dec_cov) = afwTable.convertCentroid( 

183 self.wcs, 

184 center.getX(), center.getY(), 

185 math.sqrt(centroidErr[0, 0]), math.sqrt(centroidErr[1, 1]), centroidErr[0, 1], 

186 ) 

187 np.testing.assert_allclose( 

188 (ra, dec), 

189 (skyCenter[0].asDegrees(), skyCenter[1].asDegrees()), 

190 rtol=1e-10, atol=1e-10, 

191 ) 

192 np.testing.assert_allclose( 

193 ((ra_err*factor)**2, (dec_err*factor)**2, ra_dec_cov*factor_sq), 

194 (catCoordErr[0, 0], catCoordErr[1, 1], catCoordErr[0, 1]), 

195 rtol=1e-10, atol=1e-10, 

196 ) 

197 

198 def checkCatalogs(self, maxPixDiff=1e-5, maxSkyDiff=0.001*lsst.geom.arcseconds): 

199 """Check that the source and reference object catalogs have equal centroids and coords""" 

200 self.assertEqual(len(self.sourceCat), len(self.refCat)) 

201 

202 for src, refObj in zip(self.sourceCat, self.refCat): 

203 self.assertTrue(refObj.get(self.refHasCentroidKey)) 

204 srcCentroid = src.get(self.srcCentroidKey) 

205 refCentroid = refObj.get(self.refCentroidKey) 

206 self.assertPairsAlmostEqual( 

207 srcCentroid, refCentroid, maxDiff=maxPixDiff) 

208 

209 srcCoord = src.get(self.srcCoordKey) 

210 refCoord = refObj.get(self.refCoordKey) 

211 self.assertSpherePointsAlmostEqual( 

212 srcCoord, refCoord, maxSep=maxSkyDiff) 

213 

214 def setCatalogs(self, maxPix, numPoints): 

215 """Set the source centroids and reference object coords 

216 

217 Set self.sourceCat centroids to a square grid of points 

218 and set self.refCat coords to the corresponding sky positions 

219 

220 The catalogs must be empty to start 

221 

222 @param[in] maxPix maximum pixel position; used for both x and y; 

223 the min is the negative of maxPix 

224 @param[in] numPoints number of points in x or y; total points = numPoints*numPoints 

225 """ 

226 if len(self.sourceCat) != 0: 

227 raise RuntimeError("self.sourceCat must be empty") 

228 if len(self.refCat) != 0: 

229 raise RuntimeError("self.refCat must be empty") 

230 

231 for i in np.linspace(-maxPix, maxPix, numPoints): 

232 for j in np.linspace(-maxPix, maxPix, numPoints): 

233 centroid = lsst.geom.Point2D(i, j) 

234 centroidErr = np.array([[0.05, 0], [0, 0.05]]) 

235 src = self.sourceCat.addNew() 

236 src.set(self.srcCentroidKey, centroid) 

237 src.set(self.srcCentroidErrKey, centroidErr) 

238 

239 refObj = self.refCat.addNew() 

240 coord = self.wcs.pixelToSky(centroid) 

241 refObj.set(self.refCoordKey, coord) 

242 

243 

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

245 pass 

246 

247 

248def setup_module(module): 

249 lsst.utils.tests.init() 

250 

251 

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

253 lsst.utils.tests.init() 

254 unittest.main()