Coverage for python / lsst / summit / extras / slewTimingAuxTel.py: 0%

117 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-05-05 08:54 +0000

1# This file is part of summit_extras. 

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 itertools 

23 

24import astropy 

25import matplotlib 

26import matplotlib.pyplot as plt 

27import pandas as pd 

28from astropy.time import TimeDelta 

29from lsst_efd_client import EfdClient 

30from lsst_efd_client import merge_packed_time_series as mpts 

31from matplotlib.lines import Line2D 

32from matplotlib.patches import Patch 

33 

34import lsst.daf.butler as dafButler 

35import lsst.summit.utils.butlerUtils as butlerUtils 

36from lsst.summit.utils.efdUtils import getCommands, getEfdData 

37 

38__all__ = ["plotExposureTiming"] 

39 

40READOUT_TIME = TimeDelta(2.3, format="sec") 

41 

42COMMANDS_TO_QUERY = [ 

43 # at the time of writing this was all the commands that existed for ATPtg 

44 # and ATMCS. We explicitly exclude the 20Hz ATMCS.command_trackTarget 

45 # command, and include all others. Perhaps this should be done dynamically 

46 # by using the findTopics function and removing command_trackTarget from 

47 # the list? 

48 "lsst.sal.ATPtg.command_azElTarget", 

49 "lsst.sal.ATPtg.command_disable", 

50 "lsst.sal.ATPtg.command_enable", 

51 "lsst.sal.ATPtg.command_exitControl", 

52 "lsst.sal.ATPtg.command_offsetAbsorb", 

53 "lsst.sal.ATPtg.command_offsetAzEl", 

54 "lsst.sal.ATPtg.command_offsetClear", 

55 "lsst.sal.ATPtg.command_offsetPA", 

56 "lsst.sal.ATPtg.command_offsetRADec", 

57 "lsst.sal.ATPtg.command_pointAddData", 

58 "lsst.sal.ATPtg.command_pointLoadModel", 

59 "lsst.sal.ATPtg.command_pointNewFile", 

60 "lsst.sal.ATPtg.command_poriginAbsorb", 

61 "lsst.sal.ATPtg.command_poriginClear", 

62 "lsst.sal.ATPtg.command_poriginOffset", 

63 "lsst.sal.ATPtg.command_poriginXY", 

64 "lsst.sal.ATPtg.command_raDecTarget", 

65 "lsst.sal.ATPtg.command_rotOffset", 

66 "lsst.sal.ATPtg.command_standby", 

67 "lsst.sal.ATPtg.command_start", 

68 "lsst.sal.ATPtg.command_startTracking", 

69 "lsst.sal.ATPtg.command_stopTracking", 

70 "lsst.sal.ATMCS.command_disable", 

71 "lsst.sal.ATMCS.command_enable", 

72 "lsst.sal.ATMCS.command_exitControl", 

73 "lsst.sal.ATMCS.command_setInstrumentPort", 

74 "lsst.sal.ATMCS.command_standby", 

75 "lsst.sal.ATMCS.command_start", 

76 "lsst.sal.ATMCS.command_startTracking", 

77 "lsst.sal.ATMCS.command_stopTracking", 

78 # 'lsst.sal.ATMCS.command_trackTarget', # exclude the 20Hz data 

79] 

80 

81 

82def getMountPositionData( 

83 client: EfdClient, 

84 begin: astropy.time.Time, 

85 end: astropy.time.Time, 

86 prePadding: float = 0, 

87 postPadding: float = 0, 

88) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]: 

89 """Retrieve the mount position data from the EFD. 

90 

91 Parameters 

92 ---------- 

93 client : `EfdClient` 

94 The EFD client used to retrieve the data. 

95 begin : `astropy.time.Time` 

96 The start time of the data retrieval window. 

97 end : `astropy.time.Time` 

98 The end time of the data retrieval window. 

99 prePadding : `float`, optional 

100 The amount of time to pad before the begin time, in seconds. 

101 postPadding : `float`, optional 

102 The amount of time to pad after the end time, in seconds. 

103 

104 Returns 

105 ------- 

106 az : `pandas.DataFrame` 

107 Azimuth encoder time series. 

108 el : `pandas.DataFrame` 

109 Elevation encoder time series. 

110 rot : `pandas.DataFrame` 

111 Nasmyth (rotation) encoder time series. 

112 """ 

113 mountPosition = getEfdData( 

114 client, 

115 "lsst.sal.ATMCS.mount_AzEl_Encoders", 

116 begin=begin, 

117 end=end, 

118 prePadding=prePadding, 

119 postPadding=postPadding, 

120 ) 

121 nasmythPosition = getEfdData( 

122 client, 

123 "lsst.sal.ATMCS.mount_Nasmyth_Encoders", 

124 begin=begin, 

125 end=end, 

126 prePadding=prePadding, 

127 postPadding=postPadding, 

128 ) 

129 

130 az = mpts(mountPosition, "azimuthCalculatedAngle", stride=1) 

131 el = mpts(mountPosition, "elevationCalculatedAngle", stride=1) 

132 rot = mpts(nasmythPosition, "nasmyth2CalculatedAngle", stride=1) 

133 return az, el, rot 

134 

135 

136def getAxesInPosition( 

137 client: EfdClient, 

138 begin: astropy.time.Time, 

139 end: astropy.time.Time, 

140 prePadding: float = 0, 

141 postPadding: float = 0, 

142) -> pd.DataFrame: 

143 """Return the ATMCS ``allAxesInPosition`` transitions for a window. 

144 

145 Parameters 

146 ---------- 

147 client : `lsst_efd_client.EfdClient` 

148 The EFD client used to retrieve the data. 

149 begin : `astropy.time.Time` 

150 The start time of the query window. 

151 end : `astropy.time.Time` 

152 The end time of the query window. 

153 prePadding : `float`, optional 

154 Seconds of padding before ``begin``. 

155 postPadding : `float`, optional 

156 Seconds of padding after ``end``. 

157 

158 Returns 

159 ------- 

160 transitions : `pandas.DataFrame` 

161 DataFrame of ``lsst.sal.ATMCS.logevent_allAxesInPosition`` 

162 events, indexed by time with a boolean ``inPosition`` column. 

163 """ 

164 return getEfdData( 

165 client, 

166 "lsst.sal.ATMCS.logevent_allAxesInPosition", 

167 begin=begin, 

168 end=end, 

169 prePadding=prePadding, 

170 postPadding=postPadding, 

171 ) 

172 

173 

174def plotExposureTiming( 

175 client: EfdClient, 

176 expRecords: list[dafButler.DimensionRecord], 

177 plotHexapod: bool = False, 

178 prePadding: float = 1, 

179 postPadding: float = 3, 

180) -> matplotlib.figure.Figure: 

181 """Plot the AuxTel mount command timings for a set of exposures. 

182 

183 Plots the mount position data for the entire time range of the 

184 exposures, regardless of whether the exposures are contiguous. 

185 Exposure integration and readout windows are shaded, and any 

186 commands and in-position transitions within the time range are 

187 overplotted as vertical lines. 

188 

189 Parameters 

190 ---------- 

191 client : `lsst_efd_client.EfdClient` 

192 The client object used to retrieve EFD data. 

193 expRecords : `list` [`lsst.daf.butler.DimensionRecord`] 

194 Exposure records to plot. The time axis spans from the start 

195 of the first exposure to the end of the last exposure, 

196 regardless of whether intermediate exposures are included. 

197 plotHexapod : `bool`, optional 

198 If `True`, overplot ``ATAOS.logevent_hexapodCorrectionStarted`` 

199 and ``ATAOS.logevent_hexapodCorrectionCompleted`` transitions. 

200 prePadding : `float`, optional 

201 Seconds of padding before the start of the first exposure. 

202 postPadding : `float`, optional 

203 Seconds of padding after the end of the last exposure. 

204 

205 Returns 

206 ------- 

207 fig : `matplotlib.figure.Figure` 

208 The figure containing the plot. 

209 """ 

210 inPositionAlpha = 0.5 

211 commandAlpha = 0.5 

212 integrationColor = "grey" 

213 readoutColor = "blue" 

214 

215 legendHandles: list = [] 

216 

217 expRecords = sorted(expRecords, key=lambda x: (x.day_obs, x.seq_num)) # ensure we're sorted 

218 

219 startSeqNum = expRecords[0].seq_num 

220 endSeqNum = expRecords[-1].seq_num 

221 title = f"Mount command timings for seqNums {startSeqNum} - {endSeqNum}" 

222 

223 begin = expRecords[0].timespan.begin 

224 end = expRecords[-1].timespan.end 

225 

226 az, el, rot = getMountPositionData(client, begin, end, prePadding=prePadding, postPadding=postPadding) 

227 

228 # Create a figure with a grid specification and have axes share x 

229 # and have no room between each 

230 fig = plt.figure(figsize=(12, 6)) 

231 gs = fig.add_gridspec(3, 1, hspace=0) 

232 azimuth_ax = fig.add_subplot(gs[0, 0]) 

233 elevation_ax = fig.add_subplot(gs[1, 0], sharex=azimuth_ax) 

234 rotation_ax = fig.add_subplot(gs[2, 0], sharex=azimuth_ax) 

235 axes = {"az": azimuth_ax, "el": elevation_ax, "rot": rotation_ax} 

236 

237 # plot the telemetry 

238 axes["az"].plot(az["azimuthCalculatedAngle"]) 

239 axes["el"].plot(el["elevationCalculatedAngle"]) 

240 axes["rot"].plot(rot["nasmyth2CalculatedAngle"]) 

241 

242 # shade the expRecords' regions including the readout time 

243 for i, record in enumerate(expRecords): 

244 # these need to be in UTC because matplotlib magic turns all the axis 

245 # timings into UTC when plotting from a dataframe. 

246 startExposing = record.timespan.begin.utc.datetime 

247 endExposing = record.timespan.end.utc.datetime 

248 

249 readoutEnd = (record.timespan.end + READOUT_TIME).utc.datetime 

250 seqNum = record.seq_num 

251 for axName, ax in axes.items(): 

252 ax.axvspan(startExposing, endExposing, color=integrationColor, alpha=0.3) 

253 ax.axvspan(endExposing, readoutEnd, color=readoutColor, alpha=0.1) 

254 if axName == "el": # only add seqNum annotation to bottom axis 

255 label = f"seqNum = {seqNum}" 

256 midpoint = startExposing + (endExposing - startExposing) / 2 

257 ax.annotate( 

258 label, 

259 xy=(midpoint, 0.5), 

260 xycoords=("data", "axes fraction"), 

261 ha="center", 

262 va="center", 

263 fontsize=10, 

264 color="black", 

265 ) 

266 

267 # place vertical lines at the times when axes transition in/out of position 

268 inPostionTransitions = getAxesInPosition(client, begin, end, prePadding, postPadding) 

269 for time, data in inPostionTransitions.iterrows(): 

270 inPosition = data["inPosition"] 

271 if inPosition: 

272 axes["az"].axvline(time, color="green", linestyle="--", alpha=inPositionAlpha) 

273 axes["el"].axvline(time, color="green", linestyle="--", alpha=inPositionAlpha) 

274 axes["rot"].axvline(time, color="green", linestyle="--", alpha=inPositionAlpha) 

275 else: 

276 axes["az"].axvline(time, color="red", linestyle="-", alpha=inPositionAlpha) 

277 axes["el"].axvline(time, color="red", linestyle="-", alpha=inPositionAlpha) 

278 axes["rot"].axvline(time, color="red", linestyle="-", alpha=inPositionAlpha) 

279 handle = Line2D( 

280 [0], [0], color="green", linestyle="--", label="allAxesInPosition=True", alpha=inPositionAlpha 

281 ) 

282 legendHandles.append(handle) 

283 handle = Line2D( 

284 [0], [0], color="red", linestyle="-", label="allAxesInPosition=False", alpha=inPositionAlpha 

285 ) 

286 legendHandles.append(handle) 

287 

288 # place vertical lines at the times when commands were issued 

289 commandTimes = getCommands( 

290 client, COMMANDS_TO_QUERY, begin, end, prePadding, postPadding, timeFormat="python" 

291 ) 

292 if plotHexapod: 

293 hexMoveStarts = getEfdData( 

294 client, 

295 "lsst.sal.ATAOS.logevent_hexapodCorrectionStarted", 

296 begin=begin, 

297 end=end, 

298 prePadding=prePadding, 

299 postPadding=postPadding, 

300 ) 

301 hexMoveEnds = getEfdData( 

302 client, 

303 "lsst.sal.ATAOS.logevent_hexapodCorrectionCompleted", 

304 begin=begin, 

305 end=end, 

306 prePadding=prePadding, 

307 postPadding=postPadding, 

308 ) 

309 newCommands = {} 

310 for time, data in hexMoveStarts.iterrows(): 

311 newCommands[time] = "lsst.sal.ATAOS.logevent_hexapodCorrectionStarted" 

312 for time, data in hexMoveEnds.iterrows(): 

313 newCommands[time] = "lsst.sal.ATAOS.logevent_hexapodCorrectionCompleted" 

314 commandTimes.update(newCommands) 

315 

316 uniqueCommands = list(set(commandTimes.values())) 

317 colorCycle = itertools.cycle(["b", "g", "r", "c", "m", "y", "k"]) 

318 commandColors = {command: next(colorCycle) for command in uniqueCommands} 

319 for time, command in commandTimes.items(): 

320 color = commandColors[command] 

321 axes["az"].axvline(time, linestyle="-.", alpha=commandAlpha, color=color) # type: ignore[arg-type] 

322 axes["el"].axvline(time, linestyle="-.", alpha=commandAlpha, color=color) # type: ignore[arg-type] 

323 axes["rot"].axvline(time, linestyle="-.", alpha=commandAlpha, color=color) # type: ignore[arg-type] 

324 

325 # manually build the legend to avoid duplicating the labels due to multiple 

326 # commands of the same name 

327 handles = [ 

328 Line2D([0], [0], color=color, linestyle="-.", label=label, alpha=commandAlpha) 

329 for label, color in commandColors.items() 

330 ] 

331 legendHandles.extend(handles) 

332 

333 axes["az"].set_ylabel("Azimuth (deg)") 

334 axes["el"].set_ylabel("Elevation (deg)") 

335 axes["rot"].set_ylabel("Rotation (deg)") 

336 axes["rot"].set_xlabel("Time (UTC)") # this is UTC because of the magic matplotlib does on time indices 

337 fig.suptitle(title) 

338 

339 shaded_handle = Patch(facecolor=integrationColor, alpha=0.3, label="Shutter open period") 

340 legendHandles.append(shaded_handle) 

341 shaded_handle = Patch(facecolor=readoutColor, alpha=0.1, label="Readout period") 

342 legendHandles.append(shaded_handle) 

343 # put the legend under the plot itself 

344 axes["rot"].legend(handles=legendHandles, loc="upper center", bbox_to_anchor=(0.5, -0.3), ncol=2) 

345 

346 fig.tight_layout() 

347 plt.show() 

348 return fig 

349 

350 

351if __name__ == "__main__": 

352 # example usage 

353 import lsst.summit.utils.butlerUtils as butlerUtils # noqa: F811 

354 from lsst.summit.utils.efdUtils import makeEfdClient 

355 

356 client = makeEfdClient() 

357 butler = butlerUtils.makeDefaultLatissButler(embargo=True) 

358 

359 where = "exposure.day_obs=20240215" 

360 records = list(butler.registry.queryDimensionRecords("exposure", where=where)) 

361 records = sorted(records, key=lambda x: (x.day_obs, x.seq_num)) 

362 print(f"Found {len(records)} records from {len(set(r.day_obs for r in records))} days") 

363 

364 expRecords = [records[61], records[62]] 

365 az = plotExposureTiming(client, expRecords)