Coverage for tests / test_matchOptimisticB.py: 20%

1# 

2# LSST Data Management System 

3# Copyright 2008, 2009, 2010 LSST Corporation. 

4# 

5# This product includes software developed by the 

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

7# 

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

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

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

11# (at your option) any later version. 

12# 

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

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

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

16# GNU General Public License for more details. 

17# 

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

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

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

21# 

22 

23import math 

24import os 

25import unittest 

26import pickle 

27 

28import lsst.geom 

29import lsst.afw.geom as afwGeom 

30import lsst.afw.table as afwTable 

31import lsst.utils.tests 

32import lsst.pex.exceptions as pexExcept 

33from lsst.meas.algorithms import convertReferenceCatalog 

34import lsst.meas.astrom.sip.genDistortedImage as distort 

35import lsst.meas.astrom as measAstrom 

36import lsst.meas.astrom.matchOptimisticB as matchOptimisticB 

37 

38 

39class TestMatchOptimisticB(unittest.TestCase): 

40 

41 def setUp(self): 

42 self.config = measAstrom.MatchOptimisticBTask.ConfigClass() 

43 self.matchOptimisticB = measAstrom.MatchOptimisticBTask(config=self.config) 

44 self.wcs = afwGeom.makeSkyWcs(crpix=lsst.geom.Point2D(791.4, 559.7), 

45 crval=lsst.geom.SpherePoint(36.930640, -4.939560, lsst.geom.degrees), 

46 cdMatrix=afwGeom.makeCdMatrix(scale=5.17e-5*lsst.geom.degrees)) 

47 self.distortedWcs = self.wcs 

48 

49 self.filename = os.path.join(os.path.dirname(__file__), "cat.xy.fits") 

50 self.tolArcsec = .4 

51 self.tolPixel = .1 

52 

53 def testLinearXDistort(self): 

54 self.singleTestInstance(self.filename, distort.linearXDistort) 

55 

56 def testLinearYDistort(self): 

57 self.singleTestInstance(self.filename, distort.linearYDistort) 

58 

59 def testQuadraticDistort(self): 

60 self.singleTestInstance(self.filename, distort.quadraticDistort) 

61 

62 def testLargeDistortion(self): 

63 # This transform is about as extreme as I can get: 

64 # using 0.0005 in the last value appears to produce numerical issues. 

65 # It produces a maximum deviation of 459 pixels, which should be sufficient. 

66 pixelsToTanPixels = afwGeom.makeRadialTransform([0.0, 1.1, 0.0004]) 

67 self.distortedWcs = afwGeom.makeModifiedWcs(pixelTransform=pixelsToTanPixels, 

68 wcs=self.wcs, 

69 modifyActualPixels=False) 

70 

71 def applyDistortion(src): 

72 out = src.table.copyRecord(src) 

73 out.set(out.table.getCentroidSlot().getMeasKey(), 

74 pixelsToTanPixels.applyInverse(src.getCentroid())) 

75 return out 

76 

77 self.singleTestInstance(self.filename, applyDistortion) 

78 

79 def singleTestInstance(self, filename, distortFunc, doPlot=False): 

80 sourceCat = self.loadSourceCatalog(self.filename) 

81 refCat = self.computePosRefCatalog(sourceCat) 

82 

83 # Apply source selector to sourceCat, using the astrometry config defaults 

84 tempConfig = measAstrom.AstrometryTask.ConfigClass() 

85 # This field isn't in the old test catalog. 

86 tempConfig.sourceSelector["science"].flags.bad.remove("base_PixelFlags_flag_nodata") 

87 tempConfig.sourceSelector["science"].doCentroidErrorLimit = False 

88 tempConfig.matcher.retarget(measAstrom.MatchOptimisticBTask) 

89 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None) 

90 sourceSelection = tempSolver.sourceSelector.run(sourceCat) 

91 

92 distortedCat = distort.distortList(sourceSelection.sourceCat, distortFunc) 

93 

94 if doPlot: 

95 import matplotlib.pyplot as plt 

96 undistorted = [self.wcs.skyToPixel(self.distortedWcs.pixelToSky(ss.getCentroid())) for 

97 ss in distortedCat] 

98 refs = [self.wcs.skyToPixel(ss.getCoord()) for ss in refCat] 

99 

100 def plot(catalog, symbol): 

101 plt.plot([ss.getX() for ss in catalog], [ss.getY() for ss in catalog], symbol) 

102 

103 # plot(sourceCat, 'k+') # Original positions: black + 

104 plot(distortedCat, 'b+') # Distorted positions: blue + 

105 plot(undistorted, 'g+') # Undistorted positions: green + 

106 plot(refs, 'rx') # Reference catalog: red x 

107 # The green + should overlap with the red x, because that's how matchOptimisticB does it. 

108 # The black + happens to overlap with those also, but that's beside the point. 

109 plt.show() 

110 

111 sourceCat = distortedCat 

112 

113 matchRes = self.matchOptimisticB.matchObjectsToSources( 

114 refCat=refCat, 

115 sourceCat=sourceCat, 

116 wcs=self.distortedWcs, 

117 sourceFluxField='slot_ApFlux_instFlux', 

118 refFluxField="r_flux", 

119 ) 

120 matches = matchRes.matches 

121 if doPlot: 

122 measAstrom.plotAstrometry(matches=matches, refCat=refCat, sourceCat=sourceCat) 

123 self.assertEqual(len(matches), 183) 

124 

125 refCoordKey = afwTable.CoordKey(refCat.schema["coord"]) 

126 srcCoordKey = afwTable.CoordKey(sourceCat.schema["coord"]) 

127 refCentroidKey = afwTable.Point2DKey(refCat.getSchema()["centroid"]) 

128 maxDistErr = 0*lsst.geom.radians 

129 for refObj, source, distRad in matches: 

130 sourceCoord = source.get(srcCoordKey) 

131 refCoord = refObj.get(refCoordKey) 

132 predDist = sourceCoord.separation(refCoord) 

133 distErr = abs(predDist - distRad*lsst.geom.radians) 

134 maxDistErr = max(distErr, maxDistErr) 

135 

136 if refObj.getId() != source.getId(): 

137 refCentroid = refObj.get(refCentroidKey) 

138 sourceCentroid = source.getCentroid() 

139 radius = math.hypot(*(refCentroid - sourceCentroid)) 

140 self.fail("ID mismatch: %s at %s != %s at %s; error = %0.1f pix" % 

141 (refObj.getId(), refCentroid, source.getId(), sourceCentroid, radius)) 

142 

143 self.assertLess(maxDistErr.asArcseconds(), 1e-7) 

144 

145 def computePosRefCatalog(self, sourceCat): 

146 """Generate a position reference catalog from a source catalog 

147 """ 

148 minimalPosRefSchema = convertReferenceCatalog._makeSchema(filterNameList=["r"], addCentroid=True) 

149 refCat = afwTable.SimpleCatalog(minimalPosRefSchema) 

150 for source in sourceCat: 

151 refObj = refCat.addNew() 

152 refObj.setCoord(source.getCoord()) 

153 refObj.set("centroid_x", source.getX()) 

154 refObj.set("centroid_y", source.getY()) 

155 refObj.set("hasCentroid", True) 

156 refObj.set("r_flux", source.get("slot_ApFlux_instFlux")) 

157 refObj.set("r_fluxErr", source.get("slot_ApFlux_instFluxErr")) 

158 refObj.setId(source.getId()) 

159 return refCat 

160 

161 def loadSourceCatalog(self, filename): 

162 """Load a list of xy points from a file, set coord, and return a SourceSet of points 

163 """ 

164 sourceCat = afwTable.SourceCatalog.readFits(filename) 

165 aliasMap = sourceCat.schema.getAliasMap() 

166 aliasMap.set("slot_ApFlux", "base_PsfFlux") 

167 instFluxKey = sourceCat.schema["slot_ApFlux_instFlux"].asKey() 

168 instFluxErrKey = sourceCat.schema["slot_ApFlux_instFluxErr"].asKey() 

169 

170 # print("schema=", sourceCat.schema) 

171 

172 # Source x,y positions are ~ (500,1500) x (500,1500) 

173 centroidKey = sourceCat.table.getCentroidSlot().getMeasKey() 

174 for src in sourceCat: 

175 adjCentroid = src.get(centroidKey) - lsst.geom.Extent2D(500, 500) 

176 src.set(centroidKey, adjCentroid) 

177 src.set(instFluxKey, 1000) 

178 src.set(instFluxErrKey, 1) 

179 

180 # Set catalog coord 

181 for src in sourceCat: 

182 src.updateCoord(self.wcs) 

183 return sourceCat 

184 

185 def testArgumentErrors(self): 

186 """Test argument sanity checking in matchOptimisticB 

187 """ 

188 matchControl = matchOptimisticB.MatchOptimisticBControl() 

189 

190 sourceCat = self.loadSourceCatalog(self.filename) 

191 emptySourceCat = afwTable.SourceCatalog(sourceCat.schema) 

192 

193 refCat = self.computePosRefCatalog(sourceCat) 

194 emptyRefCat = afwTable.SimpleCatalog(refCat.schema) 

195 

196 with self.assertRaises(pexExcept.InvalidParameterError): 

197 matchOptimisticB.matchOptimisticB( 

198 emptyRefCat, 

199 sourceCat, 

200 matchControl, 

201 self.wcs, 

202 0, 

203 ) 

204 with self.assertRaises(pexExcept.InvalidParameterError): 

205 matchOptimisticB.matchOptimisticB( 

206 refCat, 

207 emptySourceCat, 

208 matchControl, 

209 self.wcs, 

210 0, 

211 ) 

212 with self.assertRaises(pexExcept.InvalidParameterError): 

213 matchOptimisticB.matchOptimisticB( 

214 refCat, 

215 sourceCat, 

216 matchControl, 

217 self.wcs, 

218 len(refCat), 

219 ) 

220 with self.assertRaises(pexExcept.InvalidParameterError): 

221 matchOptimisticB.matchOptimisticB( 

222 refCat, 

223 sourceCat, 

224 matchControl, 

225 self.wcs, 

226 -1, 

227 ) 

228 

229 def testConfigPickle(self): 

230 """Test that we can pickle the Config 

231 

232 This is required for use in singleFrameDriver. 

233 See DM-18314. 

234 """ 

235 config = pickle.loads(pickle.dumps(self.config)) 

236 self.assertEqual(config, self.config) 

237 

238 

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

240 pass 

241 

242 

243def setup_module(module): 

244 lsst.utils.tests.init() 

245 

246 

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

248 

249 lsst.utils.tests.init() 

250 unittest.main()