Coverage for python / lsst / analysis / tools / atools / astrometricRepeatability.py: 23%

1# This file is part of analysis_tools. 

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/>. 

21from __future__ import annotations 

22 

23__all__ = ( 

24 "AstrometricRepeatability", 

25 "AstrometricRelativeRepeatability", 

26 "StellarAstrometricResidualsRAFocalPlanePlot", 

27 "StellarAstrometricResidualsDecFocalPlanePlot", 

28 "StellarAstrometricResidualStdDevRAFocalPlanePlot", 

29 "StellarAstrometricResidualStdDevDecFocalPlanePlot", 

30 "StellarAstrometricResidualsRASkyPlot", 

31 "StellarAstrometricResidualsDecSkyPlot", 

32) 

33 

34from lsst.pex.config import ChoiceField, Field 

35 

36from ..actions.keyedData import CalcRelativeDistances 

37from ..actions.plot import FocalPlanePlot, HistPanel, HistPlot, SkyPlot 

38from ..actions.scalar import MedianAction 

39from ..actions.vector import ( 

40 BandSelector, 

41 ConvertFluxToMag, 

42 ConvertUnits, 

43 DownselectVector, 

44 LoadVector, 

45 PerGroupStatistic, 

46 RAcosDec, 

47 RangeSelector, 

48 ResidualWithPerGroupStatistic, 

49 SnSelector, 

50 ThresholdSelector, 

51) 

52from ..interfaces import AnalysisTool 

53 

54 

55class StellarAstrometricResidualsBase(AnalysisTool): 

56 """Plot mean astrometric residuals. 

57 

58 The individual source measurements are grouped together by object index 

59 and the per-group centroid is computed. The residuals between the 

60 individual sources and these centroids are then used to construct a plot 

61 showing the mean residual as a function of the focal-plane or sky position. 

62 """ 

63 

64 def setDefaults(self): 

65 super().setDefaults() 

66 

67 # Apply per-source selection criteria 

68 self.prep.selectors.bandSelector = BandSelector() 

69 

70 self.process.buildActions.mags = ConvertFluxToMag(vectorKey="psfFlux") 

71 self.process.buildActions.residual = ConvertUnits() 

72 self.process.buildActions.residual.inUnit = "degree" 

73 self.process.buildActions.residual.outUnit = "marcsec" 

74 

75 self.process.buildActions.residual.buildAction = ResidualWithPerGroupStatistic() 

76 

77 self.process.buildActions.x = LoadVector(vectorKey="x") 

78 self.process.buildActions.y = LoadVector(vectorKey="y") 

79 

80 self.process.buildActions.detector = LoadVector(vectorKey="detector") 

81 

82 self.process.filterActions.x = DownselectVector(vectorKey="x") 

83 self.process.filterActions.x.selector = ThresholdSelector( 

84 vectorKey="mags", 

85 op="le", 

86 threshold=24, 

87 ) 

88 self.process.filterActions.y = DownselectVector( 

89 vectorKey="y", selector=self.process.filterActions.x.selector 

90 ) 

91 self.process.filterActions.z = DownselectVector( 

92 vectorKey="residual", selector=self.process.filterActions.x.selector 

93 ) 

94 self.process.filterActions.detector = DownselectVector( 

95 vectorKey="detector", selector=self.process.filterActions.x.selector 

96 ) 

97 

98 self.process.buildActions.statMask = SnSelector() 

99 self.process.buildActions.statMask.threshold = 0 

100 self.process.buildActions.statMask.fluxType = "psfFlux" 

101 

102 

103class StellarAstrometricResidualsRASkyPlot(StellarAstrometricResidualsBase): 

104 """Plot mean astrometric residuals in RA as a function of the position in 

105 RA and Dec. 

106 """ 

107 

108 def setDefaults(self): 

109 super().setDefaults() 

110 

111 # Compute per-group quantities 

112 self.process.buildActions.residual.buildAction.buildAction = RAcosDec() 

113 self.process.buildActions.x = LoadVector(vectorKey="coord_ra") 

114 self.process.buildActions.y = LoadVector(vectorKey="coord_dec") 

115 

116 self.produce = SkyPlot() 

117 

118 self.produce.plotTypes = ["any"] 

119 self.produce.plotName = "ra_residuals" 

120 self.produce.xAxisLabel = "R.A. (deg)" 

121 self.produce.yAxisLabel = "Dec. (deg)" 

122 self.produce.zAxisLabel = "RAcos(Dec) - RAcos(Dec)$_{mean}$" 

123 

124 

125class StellarAstrometricResidualsDecSkyPlot(StellarAstrometricResidualsBase): 

126 """Plot mean astrometric residuals in RA as a function of the position in 

127 RA and Dec. 

128 """ 

129 

130 def setDefaults(self): 

131 super().setDefaults() 

132 

133 # Compute per-group quantities 

134 self.process.buildActions.residual.buildAction.buildAction.vectorKey = "coord_dec" 

135 self.process.buildActions.x = LoadVector(vectorKey="coord_ra") 

136 self.process.buildActions.y = LoadVector(vectorKey="coord_dec") 

137 

138 self.produce = SkyPlot() 

139 

140 self.produce.plotTypes = ["any"] 

141 self.produce.plotName = "ra_residuals" 

142 self.produce.xAxisLabel = "R.A. (deg)" 

143 self.produce.yAxisLabel = "Dec. (deg)" 

144 self.produce.zAxisLabel = "RAcos(Dec) - RAcos(Dec)$_{mean}$" 

145 

146 

147class StellarAstrometricResidualsRAFocalPlanePlot(StellarAstrometricResidualsBase): 

148 """Plot mean astrometric residuals in RA as a function of the focal plane 

149 position. 

150 """ 

151 

152 def setDefaults(self): 

153 super().setDefaults() 

154 

155 # Compute per-group quantities 

156 self.process.buildActions.residual.buildAction.buildAction = RAcosDec() 

157 

158 self.produce = FocalPlanePlot() 

159 self.produce.zAxisLabel = "RAcos(Dec) - RAcos(Dec)$_{mean}$ (mas)" 

160 

161 

162class StellarAstrometricResidualStdDevRAFocalPlanePlot(StellarAstrometricResidualsBase): 

163 """Plot mean astrometric residuals in RA as a function of the focal plane 

164 position. 

165 """ 

166 

167 def setDefaults(self): 

168 super().setDefaults() 

169 

170 # Compute per-group quantities 

171 self.process.buildActions.residual.buildAction.buildAction = RAcosDec() 

172 

173 self.produce = FocalPlanePlot() 

174 self.produce.statistic = "std" 

175 self.produce.zAxisLabel = "Std(RAcos(Dec) - RAcos(Dec)$_{mean}$) (mas)" 

176 

177 

178class StellarAstrometricResidualsDecFocalPlanePlot(StellarAstrometricResidualsBase): 

179 """Plot mean astrometric residuals in RA as a function of the focal plane 

180 position. 

181 """ 

182 

183 def setDefaults(self): 

184 super().setDefaults() 

185 

186 # Compute per-group quantities 

187 self.process.buildActions.residual.buildAction.buildAction.vectorKey = "coord_dec" 

188 

189 self.produce = FocalPlanePlot() 

190 self.produce.zAxisLabel = "Dec - Dec$_{mean}$ (mas)" 

191 

192 

193class StellarAstrometricResidualStdDevDecFocalPlanePlot(StellarAstrometricResidualsBase): 

194 """Plot mean astrometric residuals in RA as a function of the focal plane 

195 position. 

196 """ 

197 

198 def setDefaults(self): 

199 super().setDefaults() 

200 

201 # Compute per-group quantities 

202 self.process.buildActions.residual.buildAction.buildAction.vectorKey = "coord_dec" 

203 

204 self.produce = FocalPlanePlot() 

205 self.produce.statistic = "std" 

206 self.produce.zAxisLabel = "Std(Dec - Dec$_{mean}$) (mas)" 

207 

208 

209class AstrometricRelativeRepeatability(AnalysisTool): 

210 """Calculate the AMx, ADx, AFx metrics and make histograms showing the data 

211 used to compute the metrics. 

212 """ 

213 

214 fluxType = Field[str](doc="Flux type to calculate repeatability with", default="psfFlux") 

215 xValue = Field[int](doc="Metric suffix corresponding to annulus size (1, 2, or 3)", default=1) 

216 

217 def setDefaults(self): 

218 super().setDefaults() 

219 self.prep.selectors.bandSelector = BandSelector() 

220 # Following what was done in faro, only sources with S/N between 50 

221 # and 50000 are included. The other filtering that was done in faro 

222 # is now covered by only including sources from 

223 # isolated_star_presources. 

224 self.prep.selectors.snSelector = SnSelector() 

225 self.prep.selectors.snSelector.threshold = 50 

226 self.prep.selectors.snSelector.maxSN = 50000 

227 

228 # Select only sources with magnitude between 17 and 21.5 

229 self.process.filterActions.coord_ra = DownselectVector(vectorKey="coord_ra") 

230 self.process.filterActions.coord_ra.selector = RangeSelector( 

231 vectorKey="mags", minimum=17, maximum=21.5 

232 ) 

233 self.process.filterActions.coord_dec = DownselectVector( 

234 vectorKey="coord_dec", selector=self.process.filterActions.coord_ra.selector 

235 ) 

236 self.process.filterActions.obj_index = DownselectVector( 

237 vectorKey="obj_index", selector=self.process.filterActions.coord_ra.selector 

238 ) 

239 self.process.filterActions.visit = DownselectVector( 

240 vectorKey="visit", selector=self.process.filterActions.coord_ra.selector 

241 ) 

242 

243 self.process.calculateActions.rms = CalcRelativeDistances() 

244 

245 self.produce.metric.units = { 

246 "AMx": "mas", 

247 "AFx": "percent", 

248 "ADx": "mas", 

249 } 

250 

251 self.produce.plot = HistPlot() 

252 

253 self.produce.plot.panels["panel_sep"] = HistPanel() 

254 self.produce.plot.panels["panel_sep"].hists = dict(separationResiduals="Source separations") 

255 self.produce.plot.panels["panel_sep"].label = "Separation Distances (mas)" 

256 self.produce.plot.panels["panel_sep"].lowerRange = 0 

257 self.produce.plot.panels["panel_sep"].upperRange = 80 

258 self.produce.plot.panels["panel_sep"].rangeType = "fixed" 

259 

260 self.produce.plot.panels["panel_rms"] = HistPanel() 

261 self.produce.plot.panels["panel_rms"].hists = dict(rmsDistances="Object RMS") 

262 self.produce.plot.panels["panel_rms"].label = "Per-Object RMS (mas)" 

263 self.produce.plot.panels["panel_rms"].lowerRange = 0 

264 self.produce.plot.panels["panel_rms"].upperRange = 80 

265 self.produce.plot.panels["panel_rms"].rangeType = "fixed" 

266 # TODO: DM-39163 add reference lines for ADx, AMx, and AFx. 

267 

268 def finalize(self): 

269 super().finalize() 

270 self.prep.selectors.snSelector.fluxType = self.fluxType 

271 self.process.buildActions.mags = ConvertFluxToMag(vectorKey=self.fluxType) 

272 

273 self.produce.metric.newNames = { 

274 "AMx": f"{{band}}_AM{self.xValue}", 

275 "AFx": f"{{band}}_AF{self.xValue}", 

276 "ADx": f"{{band}}_AD{self.xValue}", 

277 } 

278 

279 

280class AstrometricRepeatability(AnalysisTool): 

281 """Calculate the median position RMS of point sources.""" 

282 

283 fluxType = Field[str](doc="Flux type to calculate repeatability with", default="psfFlux") 

284 level = Field[int]( 

285 doc="Set metric name for level 1 or 2 data product (1 or 2). Input connections must be set separately" 

286 " to correspond with this value.", 

287 default=1, 

288 ) 

289 coordinate = ChoiceField[str]( 

290 doc="RA or Dec", 

291 allowed={"RA": "Repeatability in RA direction", "Dec": "Repeatability in Dec direction"}, 

292 ) 

293 

294 def setDefaults(self): 

295 super().setDefaults() 

296 self.prep.selectors.bandSelector = BandSelector() 

297 

298 self.process.buildActions.perGroupStd = ConvertUnits() 

299 self.process.buildActions.perGroupStd.inUnit = "degree" 

300 self.process.buildActions.perGroupStd.outUnit = "marcsec" 

301 self.process.buildActions.perGroupStd.buildAction = PerGroupStatistic() 

302 self.process.buildActions.perGroupStd.buildAction.func = "std" 

303 self.process.buildActions.perGroupStd.buildAction.buildAction.vectorKey = "coord_dec" 

304 

305 self.process.buildActions.perGroupMag = PerGroupStatistic() 

306 self.process.buildActions.perGroupMag.func = "mean" 

307 self.process.buildActions.perGroupMag.buildAction = ConvertFluxToMag(vectorKey=self.fluxType) 

308 

309 # Select only sources with magnitude between 17 and 21.5 

310 self.process.filterActions.perGroupStdFiltered = DownselectVector(vectorKey="perGroupStd") 

311 self.process.filterActions.perGroupStdFiltered.selector = RangeSelector( 

312 vectorKey="perGroupMag", minimum=17, maximum=21.5 

313 ) 

314 

315 self.process.calculateActions.astromRepeatStdev = MedianAction(vectorKey="perGroupStdFiltered") 

316 

317 self.produce.metric.units = { 

318 "astromRepeatStdev": "mas", 

319 } 

320 

321 self.produce.plot = HistPlot() 

322 

323 self.produce.plot.panels["panel_rms"] = HistPanel() 

324 self.produce.plot.panels["panel_rms"].hists = dict(perGroupStdFiltered="Per-Object RMS") 

325 self.produce.plot.panels["panel_rms"].label = "Per-Object RMS (mas)" 

326 

327 def finalize(self): 

328 super().finalize() 

329 self.process.buildActions.perGroupMag.buildAction.vectorKey = self.fluxType 

330 

331 if self.coordinate == "RA": 

332 self.process.buildActions.perGroupStd.buildAction.buildAction = RAcosDec() 

333 

334 self.produce.metric.newNames = { 

335 "astromRepeatStdev": f"{{band}}_dmL{self.level}AstroErr_{self.coordinate}" 

336 }