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.matcher.retarget(measAstrom.MatchOptimisticBTask) 

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

89 sourceSelection = tempSolver.sourceSelector.run(sourceCat) 

90 

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

92 

93 if doPlot: 

94 import matplotlib.pyplot as plt 

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

96 ss in distortedCat] 

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

98 

99 def plot(catalog, symbol): 

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

101 

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

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

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

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

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

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

108 plt.show() 

109 

110 sourceCat = distortedCat 

111 

112 matchRes = self.matchOptimisticB.matchObjectsToSources( 

113 refCat=refCat, 

114 sourceCat=sourceCat, 

115 wcs=self.distortedWcs, 

116 sourceFluxField='slot_ApFlux_instFlux', 

117 refFluxField="r_flux", 

118 ) 

119 matches = matchRes.matches 

120 if doPlot: 

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

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

123 

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

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

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

127 maxDistErr = 0*lsst.geom.radians 

128 for refObj, source, distRad in matches: 

129 sourceCoord = source.get(srcCoordKey) 

130 refCoord = refObj.get(refCoordKey) 

131 predDist = sourceCoord.separation(refCoord) 

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

133 maxDistErr = max(distErr, maxDistErr) 

134 

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

136 refCentroid = refObj.get(refCentroidKey) 

137 sourceCentroid = source.getCentroid() 

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

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

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

141 

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

143 

144 def computePosRefCatalog(self, sourceCat): 

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

146 """ 

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

148 refCat = afwTable.SimpleCatalog(minimalPosRefSchema) 

149 for source in sourceCat: 

150 refObj = refCat.addNew() 

151 refObj.setCoord(source.getCoord()) 

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

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

154 refObj.set("hasCentroid", True) 

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

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

157 refObj.setId(source.getId()) 

158 return refCat 

159 

160 def loadSourceCatalog(self, filename): 

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

162 """ 

163 sourceCat = afwTable.SourceCatalog.readFits(filename) 

164 aliasMap = sourceCat.schema.getAliasMap() 

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

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

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

168 

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

170 

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

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

173 for src in sourceCat: 

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

175 src.set(centroidKey, adjCentroid) 

176 src.set(instFluxKey, 1000) 

177 src.set(instFluxErrKey, 1) 

178 

179 # Set catalog coord 

180 for src in sourceCat: 

181 src.updateCoord(self.wcs) 

182 return sourceCat 

183 

184 def testArgumentErrors(self): 

185 """Test argument sanity checking in matchOptimisticB 

186 """ 

187 matchControl = matchOptimisticB.MatchOptimisticBControl() 

188 

189 sourceCat = self.loadSourceCatalog(self.filename) 

190 emptySourceCat = afwTable.SourceCatalog(sourceCat.schema) 

191 

192 refCat = self.computePosRefCatalog(sourceCat) 

193 emptyRefCat = afwTable.SimpleCatalog(refCat.schema) 

194 

195 with self.assertRaises(pexExcept.InvalidParameterError): 

196 matchOptimisticB.matchOptimisticB( 

197 emptyRefCat, 

198 sourceCat, 

199 matchControl, 

200 self.wcs, 

201 0, 

202 ) 

203 with self.assertRaises(pexExcept.InvalidParameterError): 

204 matchOptimisticB.matchOptimisticB( 

205 refCat, 

206 emptySourceCat, 

207 matchControl, 

208 self.wcs, 

209 0, 

210 ) 

211 with self.assertRaises(pexExcept.InvalidParameterError): 

212 matchOptimisticB.matchOptimisticB( 

213 refCat, 

214 sourceCat, 

215 matchControl, 

216 self.wcs, 

217 len(refCat), 

218 ) 

219 with self.assertRaises(pexExcept.InvalidParameterError): 

220 matchOptimisticB.matchOptimisticB( 

221 refCat, 

222 sourceCat, 

223 matchControl, 

224 self.wcs, 

225 -1, 

226 ) 

227 

228 def testConfigPickle(self): 

229 """Test that we can pickle the Config 

230 

231 This is required for use in singleFrameDriver. 

232 See DM-18314. 

233 """ 

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

235 self.assertEqual(config, self.config) 

236 

237 

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

239 pass 

240 

241 

242def setup_module(module): 

243 lsst.utils.tests.init() 

244 

245 

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

247 

248 lsst.utils.tests.init() 

249 unittest.main()