Coverage for python / lsst / summit / utils / guiders / tracking.py: 12%

1# This file is part of summit_utils. 

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__ = ["GuiderStarTracker"] 

24 

25import logging 

26from typing import TYPE_CHECKING 

27 

28import numpy as np 

29import pandas as pd 

30 

31from .transformation import convertRoiToCcd, convertToAltaz, convertToFocalPlane 

32 

33if TYPE_CHECKING: 

34 from .reading import GuiderData 

35 

36from .detection import GuiderStarTrackerConfig, buildReferenceCatalog, makeBlankCatalog, trackStarAcrossStamp 

37 

38 

39def _selectBrighestStar(refCatalog: pd.DataFrame, guiderName: str) -> tuple[float, float]: 

40 """ 

41 Select the brightest star from the reference catalog for a given guider. 

42 

43 Parameters 

44 ---------- 

45 refCatalog : `pd.DataFrame` 

46 Reference catalog with star positions and SNR. 

47 guiderName : `str` 

48 Name of the guider. 

49 

50 Returns 

51 ------- 

52 (xroi, yroi) : `tuple[float, float]` 

53 Coordinates of the brightest star in the ROI. 

54 """ 

55 ref = refCatalog[refCatalog["detector"] == guiderName].copy() 

56 ref.sort_values(by=["snr"], ascending=False, inplace=True) 

57 return (ref["xroi"].values[0], ref["yroi"].values[0]) 

58 

59 

60class GuiderStarTracker: 

61 """ 

62 Class to track stars in the Guider data. 

63 

64 Parameters 

65 ---------- 

66 guiderData : `GuiderData` 

67 GuiderData instance containing guider data and metadata. 

68 config : `GuiderStarTrackerConfig`, optional 

69 Config object with setup and quality control parameters. If None, a new 

70 default GuiderStarTrackerConfig is created per GuiderStarTracker 

71 instance. 

72 """ 

73 

74 def __init__( 

75 self, 

76 guiderData: GuiderData, 

77 config: GuiderStarTrackerConfig | None = None, 

78 ) -> None: 

79 """ 

80 Initialize the GuiderStarTracker with guider data and configuration. 

81 

82 Parameters 

83 ---------- 

84 guiderData : `GuiderData` 

85 GuiderData instance containing guider data and metadata. 

86 config : `GuiderStarTrackerConfig`, optional 

87 Config object with setup and quality control parameters. 

88 """ 

89 self.log = logging.getLogger(__name__) 

90 self.guiderData = guiderData 

91 self.nStamps = len(guiderData) 

92 self.expid = guiderData.expid 

93 self.shape = guiderData[guiderData.detNameMax, 0].shape 

94 

95 # detection and QC parameters from config 

96 if config is None: 

97 config = GuiderStarTrackerConfig() 

98 self.config = config 

99 

100 # initialize outputs 

101 self.blankStars = makeBlankCatalog() 

102 self.columns = list(self.blankStars.columns) 

103 

104 def trackGuiderStars(self, refCatalog: None | pd.DataFrame = None) -> pd.DataFrame: 

105 """ 

106 Track stars across guider exposures using a reference catalog. 

107 

108 Parameters 

109 ---------- 

110 refCatalog : `pd.DataFrame`, optional 

111 Reference catalog with known star positions per detector. 

112 

113 Returns 

114 ------- 

115 stars : `pd.DataFrame` 

116 DataFrame with tracked stars and their properties, 

117 including positions, fluxes, and residual offsets. 

118 """ 

119 if refCatalog is None: 

120 self.log.info("Using self-generated refcat") 

121 _refCatalog = buildReferenceCatalog(self.guiderData, self.log, self.config) 

122 refCatalog = applyQualityCuts(_refCatalog, self.shape, self.config) 

123 

124 if refCatalog.empty: 

125 self.log.warning(f"Reference catalog is empty for {self.expid}. No stars to track.") 

126 return self.blankStars 

127 

128 trackedStarTables = [] 

129 for guiderName in self.guiderData.guiderNames: 

130 refDet = refCatalog.query(f"detector == '{guiderName}'") 

131 if refDet.empty: 

132 self.log.warning(f"No stars found in reference catalog for {guiderName} in {self.expid}.") 

133 continue 

134 

135 table = self._trackStarForOneGuider(refDet, guiderName) 

136 if not table.empty: 

137 trackedStarTables.append(table) 

138 

139 if not trackedStarTables: 

140 self.log.warning(f"No stars tracked for {self.expid} for any guider.") 

141 return self.blankStars 

142 

143 # build the final catalog 

144 trackedStarCatalog = pd.concat(trackedStarTables, ignore_index=True) 

145 

146 # Set unique IDs for all tracked stars 

147 trackedStarCatalog = setUniqueId(self.guiderData, trackedStarCatalog) 

148 

149 # Make the final DataFrame with selected columns 

150 return trackedStarCatalog[self.columns] 

151 

152 def _trackStarForOneGuider(self, refCatalog: pd.DataFrame, guiderName: str) -> pd.DataFrame: 

153 """ 

154 Track stars for a single guider using the reference catalog. 

155 

156 Parameters 

157 ---------- 

158 refCatalog : `pd.DataFrame` 

159 Reference catalog with known star positions per detector. 

160 guiderName : `str` 

161 Name of the guider to process. 

162 

163 Returns 

164 ------- 

165 stars : `pd.DataFrame` 

166 DataFrame with tracked stars for the specified guider. 

167 """ 

168 gd = self.guiderData 

169 shape = self.shape 

170 cfg = self.config 

171 

172 # Select the brightest star from the reference catalog 

173 refCenter = _selectBrighestStar(refCatalog, guiderName) 

174 

175 # Measure the stars across all stamps for this guider 

176 starStamps = trackStarAcrossStamp(refCenter, gd, guiderName, cfg) 

177 

178 # Apply quality cuts to the tracked stars 

179 stars = applyQualityCuts(starStamps, shape, cfg) 

180 

181 # Filter by minimum number of detections 

182 minStampDetections = int(len(gd) * cfg.minValidStampFraction) 

183 mask = stars["stamp"].groupby(stars["detector"]).transform("count") >= minStampDetections 

184 stars = stars[mask].copy() 

185 

186 if stars.empty: 

187 self.log.warning( 

188 f"No tracked stars passed the quality cuts" f" for {guiderName} in {self.expid}." 

189 ) 

190 return self.blankStars 

191 

192 # convert to CCD, focal plane and Alt/Az coordinates 

193 stars = convertToCcdFocalPlaneAltAz(stars, gd, guiderName) 

194 

195 # Compute the rotator angle: theta = np.arctan2(yfp, xfp) 

196 stars = computeRotatorAngle(stars) 

197 

198 # Convert e1, e2 to Alt/Az coordinates 

199 stars = convertEllipticity(stars, gd.camRotAngle) 

200 

201 # Add timestamp and elapsed time 

202 stars = addTimeStamp(stars, gd, guiderName) 

203 

204 # Compute Offsets 

205 stars = computeOffsets(stars) 

206 

207 return stars 

208 

209 

210def addTimeStamp(stars: pd.DataFrame, guiderData: GuiderData, guiderName: str) -> pd.DataFrame: 

211 """ 

212 Add timestamp and elapsed time to the star DataFrame. 

213 

214 Parameters 

215 ---------- 

216 stars : `pd.DataFrame` 

217 DataFrame with star measurements. 

218 guiderData : `GuiderData` 

219 GuiderData instance containing guider data and metadata. 

220 guiderName : `str` 

221 Name of the guider. 

222 

223 Returns 

224 ------- 

225 stars : `pd.DataFrame` 

226 DataFrame with added 'timestamp' and 'elapsed_time' columns. 

227 """ 

228 gd = guiderData 

229 # the stamp are aligned with the index of the timestamps 

230 stampIndex = stars["stamp"].to_numpy() 

231 indices = np.array([ix for ix in stampIndex], dtype=int) 

232 

233 # get the timestamp for each stamp 

234 timeStamp = gd.timestampMap[guiderName][indices] 

235 

236 # inital time is the time of the first stamp of the guider with most stamps 

237 t0 = gd.timestampMap[gd.detNameMax][0] 

238 

239 stars["timestamp"] = timeStamp 

240 stars["elapsed_time"] = np.where(timeStamp.mask, np.nan, (timeStamp - t0).sec) 

241 return stars 

242 

243 

244def convertToCcdFocalPlaneAltAz(stars: pd.DataFrame, guiderData: GuiderData, guiderName: str) -> pd.DataFrame: 

245 """ 

246 Convert star positions to CCD, focal plane, and Alt/Az coordinates. 

247 

248 Parameters 

249 ---------- 

250 stars : `pd.DataFrame` 

251 DataFrame with star measurements in ROI coordinates. 

252 guiderData : `GuiderData` 

253 GuiderData instance containing guider data and metadata. 

254 guiderName : `str` 

255 Name of the guider. 

256 

257 Returns 

258 ------- 

259 stars : `pd.DataFrame` 

260 DataFrame with added columns for CCD, focal plane, 

261 and Alt/Az coordinates. 

262 """ 

263 gd = guiderData 

264 obsTime = gd.obsTime 

265 detNum = gd.getGuiderDetNum(guiderName) 

266 wcs = gd.getWcs(guiderName) 

267 pixelScale = wcs.getPixelScale().asArcseconds() 

268 

269 # Convert to CCD coordinates 

270 stars["xccd"], stars["yccd"] = convertRoiToCcd( 

271 stars["xroi"].to_numpy(), stars["yroi"].to_numpy(), gd, guiderName 

272 ) 

273 

274 # Convert to focal plane coordinates 

275 stars["xfp"], stars["yfp"] = convertToFocalPlane( 

276 stars["xccd"].to_numpy(), stars["yccd"].to_numpy(), detNum 

277 ) 

278 

279 # Convert to Alt/Az coordinates 

280 stars["alt"], stars["az"] = convertToAltaz( 

281 stars["xccd"].to_numpy(), stars["yccd"].to_numpy(), wcs, obsTime 

282 ) 

283 

284 # Convert fwhm to arcsec 

285 stars["fwhm"] = stars["fwhm"] * pixelScale 

286 

287 return stars 

288 

289 

290def applyQualityCuts( 

291 stars: pd.DataFrame, shape: tuple[float, float], config: GuiderStarTrackerConfig 

292) -> pd.DataFrame: 

293 """ 

294 Apply cuts according to min SNR, maximum ellipticity and edge margin. 

295 

296 Parameters 

297 ---------- 

298 stars : `pd.DataFrame` 

299 DataFrame with star measurements. 

300 shape : `tuple[float, float]` 

301 Shape of the ROI (cols, rows). 

302 config : `GuiderStarTrackerConfig` 

303 Configuration object with quality cut parameters. 

304 

305 Returns 

306 ------- 

307 stars : `pd.DataFrame` 

308 DataFrame with stars that passed the quality cuts. 

309 """ 

310 minSnr = config.minSnr 

311 maxEllipticity = config.maxEllipticity 

312 edgeMargin = config.edgeMargin 

313 roiCols, roiRows = shape 

314 

315 # Filter by minimum SNR 

316 mask1 = (stars["snr"] >= minSnr) & (stars["flux"] > 0) & (stars["flux_err"] > 0) 

317 

318 # Filter by minimum number of detections and maximum ellipticity 

319 eabs = np.hypot(stars["e1"], stars["e2"]) 

320 mask3 = (stars["e1"].abs() <= maxEllipticity) & (stars["e1"].abs() <= maxEllipticity) 

321 mask3 &= eabs <= maxEllipticity 

322 

323 # Add edge margin mask 

324 mask4 = ( 

325 (stars["xroi"] >= edgeMargin) 

326 & (stars["xroi"] <= roiRows - edgeMargin) 

327 & (stars["yroi"] >= edgeMargin) 

328 & (stars["yroi"] <= roiCols - edgeMargin) 

329 ) 

330 # Combine all masks 

331 mask = mask1 & mask3 & mask4 

332 return stars[mask].copy() 

333 

334 

335def setUniqueId(guiderData: GuiderData, stars: pd.DataFrame) -> pd.DataFrame: 

336 """ 

337 Assign unique IDs to tracked stars. 

338 

339 Parameters 

340 ---------- 

341 guiderData : `GuiderData` 

342 GuiderData instance containing guider data and metadata. 

343 stars : `pd.DataFrame` 

344 DataFrame with star measurements. 

345 

346 Returns 

347 ------- 

348 stars : `pd.DataFrame` 

349 DataFrame with added 'detid' and 'trackid' columns. 

350 

351 """ 

352 # Create a numeric “global” starid: 

353 detMap = guiderData.guiderNameMap 

354 stars["detid"] = stars["detector"].map(detMap) 

355 stars["trackid"] = stars["detid"] * 1000 + stars["stamp"] 

356 return stars 

357 

358 

359def computeOffsets(stars: pd.DataFrame) -> pd.DataFrame: 

360 """ 

361 Compute the offsets for each star in the catalog. 

362 

363 Parameters 

364 ---------- 

365 stars : `pd.DataFrame` 

366 DataFrame with star measurements. 

367 

368 Returns 

369 ------- 

370 stars : `pd.DataFrame` 

371 DataFrame with added offset columns: 

372 - dx, dy : offsets in CCD coordinates (pixels) 

373 - dxfp, dyfp : offsets in focal plane coordinates (mm) 

374 - dalt, daz : offsets in Alt/Az coordinates (arcsec) 

375 - dtheta : offset in rotator angle (arcsec) 

376 - magoffset : magnitude offset (mmag) 

377 

378 """ 

379 # make reference positions 

380 stars["xroi_ref"] = stars.groupby("detector")["xroi"].transform("median") 

381 stars["yroi_ref"] = stars.groupby("detector")["yroi"].transform("median") 

382 stars["xccd_ref"] = stars.groupby("detector")["xccd"].transform("median") 

383 stars["yccd_ref"] = stars.groupby("detector")["yccd"].transform("median") 

384 stars["xfp_ref"] = stars.groupby("detector")["xfp"].transform("median") 

385 stars["yfp_ref"] = stars.groupby("detector")["yfp"].transform("median") 

386 stars["alt_ref"] = stars.groupby("detector")["alt"].transform("median") 

387 stars["az_ref"] = stars.groupby("detector")["az"].transform("median") 

388 stars["theta_ref"] = stars.groupby("detector")["theta"].transform("median") 

389 

390 # Compute all your offsets 

391 stars["dx"] = stars["xccd"] - stars["xccd_ref"] 

392 stars["dy"] = stars["yccd"] - stars["yccd_ref"] 

393 stars["dxfp"] = stars["xfp"] - stars["xfp_ref"] 

394 stars["dyfp"] = stars["yfp"] - stars["yfp_ref"] 

395 stars["dalt"] = (stars["alt"] - stars["alt_ref"]) * 3600 

396 stars["daz"] = (stars["az"] - stars["az_ref"]) * 3600 

397 stars["dtheta"] = (stars["theta"] - stars["theta_ref"]) * 3600 

398 

399 # Correct for cos(alt) in daz 

400 stars["daz"] = np.cos(stars["alt_ref"] * np.pi / 180) * stars["daz"] 

401 

402 # compute mag offset 

403 stars["flux_ref"] = stars.groupby("detector")["flux"].transform("median") 

404 stars["flux_ref"] = pd.to_numeric(stars["flux_ref"], errors="coerce") 

405 stars["flux"] = pd.to_numeric(stars["flux"], errors="coerce") 

406 stars["magoffset"] = -2.5 * np.log10((stars["flux"] + 1e-12) / (stars["flux_ref"] + 1e-12)) 

407 # convert to mmag and replace inf with nan 

408 stars["magoffset"] = 1000 * stars["magoffset"].replace([np.inf, -np.inf], np.nan) 

409 

410 # not detections have nan offsets 

411 offsetCols = ["dx", "dy", "dxfp", "dyfp", "dalt", "daz", "dtheta", "magoffset"] 

412 stars.loc[stars["flux"].isna(), offsetCols] = np.nan 

413 

414 return stars 

415 

416 

417def computeRotatorAngle(stars: pd.DataFrame) -> pd.DataFrame: 

418 """ 

419 Compute the rotator angle (theta) and its propagated 1-sigma uncertainty 

420 for each row in `stars`. 

421 

422 Parameters 

423 ---------- 

424 stars : `pd.DataFrame` 

425 DataFrame with star measurements. Must contain 'xfp', 'yfp', 'xerr', 

426 and 'yerr' columns. 

427 

428 Returns 

429 ------- 

430 stars : `pd.DataFrame` 

431 DataFrame with added 'theta' and 'theta_err' columns. 

432 """ 

433 mm = 0.001 # convert microns to mm 

434 xfp = stars["xfp"].to_numpy() 

435 yfp = stars["yfp"].to_numpy() 

436 sigmaX = stars["xerr"].to_numpy() * 10 * mm 

437 sigmaY = stars["yerr"].to_numpy() * 10 * mm 

438 

439 # Angle in degrees 

440 theta = np.degrees(np.arctan2(yfp, xfp)) 

441 

442 # Denominator r^2 = x^2 + y^2 

443 denom = xfp**2 + yfp**2 

444 

445 # Suppress warnings for rows where denom == 0; set sigmaTheta to NaN there. 

446 with np.errstate(divide="ignore", invalid="ignore"): 

447 sigmaThetaRad = np.sqrt((yfp**2 * sigmaX**2 + xfp**2 * sigmaY**2) / denom**2) 

448 sigmaTheta = np.degrees(sigmaThetaRad) 

449 sigmaTheta = np.where(denom == 0, np.nan, sigmaTheta) 

450 

451 stars["theta"] = theta 

452 stars["theta_err"] = sigmaTheta * 3600 # convert to arcsec 

453 return stars 

454 

455 

456def convertEllipticity(stars: pd.DataFrame, camRotAngleDeg: float) -> pd.DataFrame: 

457 """ 

458 Rotate ellipticity components (e1, e2) to Alt/Az. 

459 

460 Parameters 

461 ---------- 

462 stars : `pd.DataFrame` 

463 DataFrame with star measurements. 

464 camRotAngleDeg : `float` 

465 Camera rotator angle in degrees. 

466 

467 Returns 

468 ------- 

469 stars : `pd.DataFrame` 

470 DataFrame with added columns 'e1_altaz' and 'e2_altaz'. 

471 """ 

472 e1_rot, e2_rot = rotateEllipticity(stars["e1"], stars["e2"], camRotAngleDeg) 

473 stars["e1_altaz"] = e1_rot 

474 stars["e2_altaz"] = e2_rot 

475 return stars 

476 

477 

478def rotateEllipticity(e1: float, e2: float, theta_deg: float) -> tuple[float, float]: 

479 """ 

480 Rotate ellipticity components (e1, e2) by theta_deg degrees. 

481 

482 Parameters 

483 ---------- 

484 e1 : `float` 

485 First ellipticity component. 

486 e2 : `float` 

487 Second ellipticity component. 

488 theta_deg : `float` 

489 Rotation angle in degrees. 

490 

491 Returns 

492 ------- 

493 e1_rot, e2_rot : `tuple[float, float]` 

494 Rotated ellipticity components. 

495 """ 

496 theta = np.deg2rad(theta_deg) 

497 cos2t = np.cos(2 * theta) 

498 sin2t = np.sin(2 * theta) 

499 e1_rot = e1 * cos2t + e2 * sin2t 

500 e2_rot = -e1 * sin2t + e2 * cos2t 

501 return e1_rot, e2_rot