Coverage for tests/test_coordErrorConvention.py: 94%

1# This file is part of meas_base. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://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 program is free software: you can redistribute it and/or modify 

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

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

12# (at your option) any later version. 

13# 

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

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

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

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import unittest 

23 

24import numpy as np 

25 

26import lsst.geom 

27import lsst.afw.geom as afwGeom 

28import lsst.afw.table as afwTable 

29import lsst.utils.tests 

30 

31 

32class CoordErrorConventionTestCase(lsst.utils.tests.TestCase): 

33 """Verify that ``coord_raErr`` / ``coord_decErr`` / ``coord_ra_dec_Cov`` 

34 carry the tangent-plane uncertainty propagated through a local gnomonic 

35 Jacobian. 

36 """ 

37 

38 def setUp(self): 

39 # Reference position at high declination so cos(Dec) = 0.5 

40 self.referenceRa = 0.0 * lsst.geom.degrees 

41 self.referenceDec = 60.0 * lsst.geom.degrees 

42 self.cosDec = np.cos(self.referenceDec.asRadians()) 

43 

44 self.pixelScale = 0.2 * lsst.geom.arcseconds 

45 self.pixelScaleRadians = self.pixelScale.asRadians() 

46 

47 crpix = lsst.geom.Point2D(0.0, 0.0) 

48 crval = lsst.geom.SpherePoint(self.referenceRa, self.referenceDec) 

49 cdMatrix = afwGeom.makeCdMatrix(scale=self.pixelScale, 

50 orientation=0.0 * lsst.geom.radians, 

51 flipX=False) 

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

53 

54 # Minimal SourceCatalog with a centroid (slot + uncertainty fields) 

55 # and the global coord_raErr / coord_decErr / coord_ra_dec_Cov fields 

56 # that SourceRecord.updateCoord writes into. 

57 schema = afwTable.SourceTable.makeMinimalSchema() 

58 schema.addField("test_x", type="D", units="pixel", 

59 doc="centroid x") 

60 schema.addField("test_y", type="D", units="pixel", 

61 doc="centroid y") 

62 schema.addField("test_xErr", type="F", units="pixel", 

63 doc="1-sigma uncertainty on x") 

64 schema.addField("test_yErr", type="F", units="pixel", 

65 doc="1-sigma uncertainty on y") 

66 schema.addField("test_x_y_Cov", type="F", units="pixel^2", 

67 doc="x/y covariance") 

68 afwTable.CoordKey.addErrorFields(schema) 

69 schema.getAliasMap().set("slot_Centroid", "test") 

70 

71 self.catalog = afwTable.SourceCatalog(schema) 

72 

73 def _addSource(self, sigmaPix, xyCov=0.0): 

74 """Add one source at the WCS reference pixel with covariance 

75 ``[[sigmaPix**2, xyCov], [xyCov, sigmaPix**2]]``. 

76 """ 

77 source = self.catalog.addNew() 

78 source["test_x"] = 0.0 

79 source["test_y"] = 0.0 

80 source["test_xErr"] = sigmaPix 

81 source["test_yErr"] = sigmaPix 

82 source["test_x_y_Cov"] = xyCov 

83 return source 

84 

85 def test_skyPositionMatchesWcsReference(self): 

86 """Check the simple case where the recorded sky position is the WCS 

87 reference. 

88 """ 

89 source = self._addSource(sigmaPix=1.0) 

90 source.updateCoord(self.wcs) 

91 self.assertAlmostEqual(source["coord_ra"], self.referenceRa.asRadians(), 

92 places=10) 

93 self.assertAlmostEqual(source["coord_dec"], self.referenceDec.asRadians(), 

94 places=10) 

95 

96 def test_coordErrIsTangentPlaneSigma(self): 

97 """Check that the coordinate errors equal the pixel-scale sigma. 

98 """ 

99 sigmaPix = 1.0 

100 source = self._addSource(sigmaPix) 

101 source.updateCoord(self.wcs) 

102 expectedTangentErr = sigmaPix * self.pixelScaleRadians 

103 self.assertFloatsAlmostEqual(source["coord_decErr"], expectedTangentErr, rtol=1e-3) 

104 

105 # Note that there is no cos(Dec) factor. 

106 self.assertFloatsAlmostEqual(source["coord_raErr"], expectedTangentErr, rtol=1e-3) 

107 

108 def test_raDecCovIsTangentPlaneCov(self): 

109 """Check the calculations of the covariance matrix. 

110 

111 Note that ``coord_ra_dec_Cov`` is Cov(RA*cos(Dec), Dec). 

112 """ 

113 sigmaPix = 1.0 

114 

115 # Diagonal pixel covariance -> diagonal sky covariance. 

116 source = self._addSource(sigmaPix, xyCov=0.0) 

117 source.updateCoord(self.wcs) 

118 self.assertFloatsAlmostEqual(source["coord_ra_dec_Cov"], 0.0, 

119 atol=1e-20) 

120 

121 # Non-diagonal pixel covariance. 

122 rhoPix = 0.5 

123 source = self._addSource(sigmaPix, xyCov=rhoPix * sigmaPix**2) 

124 source.updateCoord(self.wcs) 

125 expectedCov = self.pixelScaleRadians**2 * (rhoPix * sigmaPix**2) 

126 if not self.wcs.isFlipped: 126 ↛ 128line 126 didn't jump to line 128 because the condition on line 126 was always true

127 expectedCov *= -1 

128 self.assertFloatsAlmostEqual(source["coord_ra_dec_Cov"], expectedCov, rtol=1e-3) 

129 

130 

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

132 pass 

133 

134 

135def setup_module(module): 

136 lsst.utils.tests.init() 

137 

138 

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

140 lsst.utils.tests.init() 

141 unittest.main()