Coverage for tests / test_fit_stellar_motion.py: 22%

1# This file is part of drp_tasks. 

2# 

3# LSST Data Management System 

4# This product includes software developed by the 

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

6# See COPYRIGHT file at the top of the source tree. 

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 <https://www.lsstcorp.org/LegalNotices/>. 

21# 

22 

23import unittest 

24 

25import astropy.time 

26import astropy.units as u 

27import numpy as np 

28import numpy.testing as npt 

29from astro_metadata_translator import makeObservationInfo 

30from astropy.coordinates import Angle, EarthLocation, SkyCoord 

31 

32import lsst.afw.table as afwTable 

33from lsst.drp.tasks.fit_stellar_motion import FitStellarMotionTask 

34from lsst.obs.base import MakeRawVisitInfoViaObsInfo 

35from lsst.pipe.base import InMemoryDatasetHandle 

36 

37 

38class FitStellarMotionTestCase(unittest.TestCase): 

39 def setUp(self): 

40 

41 visits = [2025051500122, 2025051500133, 2025051500144, 2025052300192, 2025052400145] 

42 self.time0 = astropy.time.Time(2025.5, format="jyear") 

43 timeDeltas = np.arange(-2, 3) * astropy.time.TimeDelta(300 * u.day) 

44 times = self.time0 + timeDeltas 

45 self.visitSummaryHandles = { 

46 visit: self._makeVisitSummaryTableHandle(time) for visit, time in zip(visits, times) 

47 } 

48 

49 nObjects = 10 

50 self.ras = np.linspace(150, 151, nObjects) 

51 self.decs = np.linspace(2, 3, nObjects) 

52 self.pmRAs = np.linspace(-5, 5, nObjects) 

53 self.pmDecs = np.linspace(-5, -3, nObjects) 

54 self.parallaxes = np.zeros(nObjects) 

55 

56 self.starCatalogHandle, self.visitCatalogHandles, self.starSources = self._makeCatalogs( 

57 visits, nObjects, timeDeltas 

58 ) 

59 

60 self.task = FitStellarMotionTask() 

61 

62 self.visitStars, self.visitInfo = self.task._load_sources( 

63 self.starSources, self.visitSummaryHandles, self.visitCatalogHandles 

64 ) 

65 

66 self.refCat = self._make_refCat(nObjects) 

67 

68 def _makeVisitSummaryTableHandle(self, time): 

69 """Make some arbitrary visit summary tables.""" 

70 schema = afwTable.ExposureTable.makeMinimalSchema() 

71 visitSummaryTable = afwTable.ExposureCatalog(schema) 

72 

73 record = visitSummaryTable.addNew() 

74 lsstLat = -30.244639 * u.degree 

75 lsstLon = -70.749417 * u.degree 

76 lsstAlt = 2663.0 * u.m 

77 loc = EarthLocation(lat=lsstLat, lon=lsstLon, height=lsstAlt) 

78 

79 obsInfo = makeObservationInfo( 

80 location=loc, 

81 datetime_begin=time - 15 * u.second, 

82 datetime_end=time + 15 * u.second, 

83 boresight_rotation_angle=Angle(0.0 * u.degree), 

84 boresight_rotation_coord="sky", 

85 observation_type="science", 

86 ) 

87 

88 visitInfo = MakeRawVisitInfoViaObsInfo.observationInfo2visitInfo(obsInfo) 

89 record.setVisitInfo(visitInfo) 

90 handle = InMemoryDatasetHandle(visitSummaryTable) 

91 return handle 

92 

93 def _makeCatalogs(self, visits, nObjects, timeDeltas): 

94 """Make catalogs to match the isolated_star, isolated_star_association, 

95 and source catalogs.""" 

96 nVisits = len(visits) 

97 

98 objectCoords = SkyCoord( 

99 self.ras * u.degree, 

100 self.decs * u.degree, 

101 pm_ra_cosdec=self.pmRAs * np.cos((self.decs * u.degree).to(u.radian)) * u.mas / u.yr, 

102 pm_dec=self.pmDecs * u.mas / u.yr, 

103 obstime=self.time0, 

104 ) 

105 

106 objIndices = np.arange(nObjects) 

107 

108 starCatalog = astropy.table.Table({"isolated_star_id": objIndices, "ra": self.ras, "dec": self.decs}) 

109 starCatalogHandle = InMemoryDatasetHandle(starCatalog, storageClass="ArrowAstropy") 

110 

111 sourceIds = [] 

112 visitCatalogHandles = {} 

113 for v, visit in enumerate(visits): 

114 visitCoords = objectCoords.apply_space_motion(dt=timeDeltas[v]) 

115 

116 # Make some arbitrary sourceIds 

117 visitSourceIds = visit * 100 + np.arange(nObjects) 

118 sourceIds.extend(visitSourceIds) 

119 catalog = { 

120 "sourceId": visitSourceIds, 

121 "ra": visitCoords.ra.degree, 

122 "dec": visitCoords.dec.degree, 

123 "raErr": np.ones(nObjects) * 1e-6, 

124 "decErr": np.ones(nObjects) * 1e-6, 

125 "ra_dec_Cov": np.ones(nObjects) * 1e-14, 

126 } 

127 catalog = astropy.table.Table(catalog) 

128 visitCatalogHandles[visit] = InMemoryDatasetHandle( 

129 catalog, 

130 storageClass="ArrowAstropy", 

131 parameters={ 

132 "columns": [ 

133 "sourceId", 

134 "ra", 

135 "dec", 

136 "raErr", 

137 "decErr", 

138 "ra_dec_Cov", 

139 ] 

140 }, 

141 ) 

142 starSourceDict = { 

143 "visit": np.repeat(visits, nObjects), 

144 "obj_index": np.tile(objIndices, nVisits), 

145 "sourceId": sourceIds, 

146 } 

147 starSources = astropy.table.Table(starSourceDict) 

148 starSources.add_index("sourceId") 

149 

150 return starCatalogHandle, visitCatalogHandles, starSources 

151 

152 def _make_refCat(self, nObjects): 

153 """Make a reference catalog.""" 

154 referenceIds = np.arange(nObjects) 

155 covariance = np.zeros((nObjects, 5, 5)) 

156 covariance[:, 0, 0] = 1 

157 covariance[:, 1, 1] = 1 

158 covariance[:, 2, 2] = 0.1 

159 covariance[:, 3, 3] = 0.1 

160 covariance[:, 4, 4] = 0.1 

161 refCat = { 

162 "ra": self.ras * u.degree, 

163 "dec": self.decs * u.degree, 

164 "raPM": self.pmRAs * u.mas / u.yr, 

165 "decPM": self.pmDecs * u.mas / u.yr, 

166 "parallax": self.parallaxes * u.mas, 

167 "covariance": covariance, 

168 "id": referenceIds, 

169 } 

170 return astropy.table.Table(refCat) 

171 

172 def test_fit_objects_without_reference(self): 

173 """Turn off task.config.includeReferenceCatalog to test fit without a 

174 reference catalog. 

175 """ 

176 

177 self.task.config.includeReferenceCatalog = False 

178 outCat, predictedRADec = self.task._fit_objects( 

179 self.visitStars, 

180 self.starCatalogHandle, 

181 self.starSources, 

182 self.visitInfo, 

183 self.time0, 

184 ) 

185 

186 npt.assert_allclose(self.pmRAs, outCat["raPM"], rtol=2e-3) 

187 npt.assert_allclose(self.pmDecs, outCat["decPM"], rtol=1e-3) 

188 npt.assert_allclose(self.parallaxes, outCat["parallax"], atol=2e-3) 

189 

190 npt.assert_allclose(predictedRADec["ra"], self.visitStars["ra"]) 

191 npt.assert_allclose(predictedRADec["dec"], self.visitStars["dec"]) 

192 

193 def test_fit_objects_with_reference(self): 

194 """Test fit with a reference catalog.""" 

195 

196 outCat, predictedRADec = self.task._fit_objects( 

197 self.visitStars, 

198 self.starCatalogHandle, 

199 self.starSources, 

200 self.visitInfo, 

201 self.time0, 

202 refCatalog=self.refCat, 

203 ) 

204 

205 npt.assert_allclose(self.pmRAs, outCat["raPM"], rtol=1e-4) 

206 npt.assert_allclose(self.pmDecs, outCat["decPM"], rtol=1e-4) 

207 npt.assert_allclose(self.parallaxes, outCat["parallax"], atol=1e-3) 

208 

209 npt.assert_allclose(predictedRADec["ra"], self.visitStars["ra"]) 

210 npt.assert_allclose(predictedRADec["dec"], self.visitStars["dec"]) 

211 

212 

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

214 unittest.main()