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

1# This file is part of summit_utils.

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.

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

23__all__ = ["GuiderStarTracker"]

25import logging

26from typing import TYPE_CHECKING

28import numpy as np

29import pandas as pd

31from .transformation import convertRoiToCcd, convertToAltaz, convertToFocalPlane

33if TYPE_CHECKING:

34 from .reading import GuiderData

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

39class GuiderStarTracker:

40 """

41 Class to track stars in the Guider data.

43 Parameters

44 ----------

45 guiderData : `GuiderData`

46 GuiderData instance containing guider data and metadata.

47 config : `GuiderStarTrackerConfig`, optional

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

49 default GuiderStarTrackerConfig is created per GuiderStarTracker

50 instance.

51 """

53 def __init__(

54 self,

55 guiderData: GuiderData,

56 config: GuiderStarTrackerConfig | None = None,

57 ) -> None:

58 """

59 Initialize the GuiderStarTracker with guider data and configuration.

61 Parameters

62 ----------

63 guiderData : `GuiderData`

64 GuiderData instance containing guider data and metadata.

65 config : `GuiderStarTrackerConfig`, optional

66 Config object with setup and quality control parameters.

67 """

68 self.log = logging.getLogger(__name__)

69 self.guiderData = guiderData

70 self.nStamps = len(guiderData)

71 self.expid = guiderData.expid

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

74 # detection and QC parameters from config

75 if config is None:

76 config = GuiderStarTrackerConfig()

77 self.config = config

79 # initialize outputs

80 self.blankStars = makeBlankCatalog()

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

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

84 """

85 Track stars across guider exposures using a reference catalog.

87 Parameters

88 ----------

89 refCatalog : `pd.DataFrame`, optional

90 Reference catalog with known star positions per detector.

92 Returns

93 -------

94 stars : `pd.DataFrame`

95 DataFrame with tracked stars and their properties,

96 including positions, fluxes, and residual offsets.

97 """

98 _refCatalog = None

99 if refCatalog is None:

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

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

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

103

104 if refCatalog.empty:

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

106 return self.blankStars

107

108 trackedStarTables = []

109 for guiderName in self.guiderData.guiderNames:

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

111 if refDet.empty:

112 # Diagnose: check what was in the unfiltered refcat

113 preCut = _refCatalog.query(f"detector == '{guiderName}'") if _refCatalog is not None else None

114 if preCut is not None and not preCut.empty:

115 best = preCut.sort_values("snr", ascending=False).iloc[0]

116 reason = _diagnoseQualityCutRejections(

117 preCut.sort_values("snr", ascending=False).iloc[:1], self.shape, self.config

118 )[0]

119 self.log.warning(

120 f"No stars in refcat for {guiderName} in {self.expid}. "

121 f"Best candidate: flux={best['flux']:.0f} snr={best['snr']:.1f} "

122 f"e={np.hypot(best['e1'], best['e2']):.2f} fwhm={best['fwhm']:.1f} "

123 f"pos=({best['xroi']:.0f},{best['yroi']:.0f}). Rejected: {reason}"

124 )

125 else:

126 self.log.warning(

127 f"No stars detected for {guiderName} in {self.expid} " f"(empty before quality cuts)."

128 )

129 continue

130

131 table = self._trackStarForOneGuider(refDet, guiderName)

132 if not table.empty:

133 trackedStarTables.append(table)

134

135 if not trackedStarTables:

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

137 return self.blankStars

138

139 # build the final catalog

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

141

142 # Set unique IDs for all tracked stars

143 trackedStarCatalog = setUniqueId(self.guiderData, trackedStarCatalog)

144

145 # Make the final DataFrame with selected columns

146 return trackedStarCatalog[self.columns]

147

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

149 """

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

151

152 Iteratively tries to track stars starting with the highest SNR,

153 and falls back to the next best candidates if quality cuts fail.

154

155 Parameters

156 ----------

157 refCatalog : `pd.DataFrame`

158 Reference catalog with known star positions per detector.

159 guiderName : `str`

160 Name of the guider to process.

161

162 Returns

163 -------

164 stars : `pd.DataFrame`

165 DataFrame with tracked stars for the specified guider.

166 """

167 gd = self.guiderData

168 shape = self.shape

169 cfg = self.config

170

171 # Sort by SNR (brightest first) for iterative fallback

172 refCatalogSorted = refCatalog.sort_values(by="snr", ascending=False).reset_index(drop=True)

173

174 stars = None

175 for idx, row in refCatalogSorted.iterrows():

176 refCenter = (row["xroi"], row["yroi"])

177 starSnr = row["snr"]

178

179 self.log.debug(f"Attempting to track star {idx + 1} (SNR={starSnr:.2f}) for {guiderName}")

180

181 # Measure the stars across all stamps for this guider

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

183

184 if starStamps.empty:

185 self.log.debug(f"Star {idx + 1} produced no detections")

186 continue

187

188 # Apply quality cuts to the tracked stars

189 stars = applyQualityCuts(starStamps, shape, cfg)

190

191 # Filter by minimum number of detections

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

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

194 stars = stars[mask].copy()

195

196 if not stars.empty:

197 self.log.debug(

198 f"Successfully tracked star {idx + 1} (SNR={starSnr:.2f}) "

199 f"for {guiderName} with {len(stars)} detections"

200 )

201 break

202 else:

203 self.log.debug(f"Star {idx + 1} (SNR={starSnr:.2f}) failed quality cuts for {guiderName}")

204 stars = None

205

206 if stars is None or stars.empty:

207 self.log.warning(

208 f"No stars passed quality cuts for {guiderName} in {self.expid}. "

209 f"Tried {len(refCatalogSorted)} candidates."

210 )

211 return self.blankStars

212

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

214 stars = convertToCcdFocalPlaneAltAz(stars, gd, guiderName)

215

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

217 stars = computeRotatorAngle(stars)

218

219 # Convert e1, e2 to Alt/Az coordinates

220 stars = convertEllipticity(stars, gd.camRotAngle)

221

222 # Add timestamp and elapsed time

223 stars = addTimeStamp(stars, gd, guiderName)

224

225 # Compute Offsets

226 stars = computeOffsets(stars)

227

228 return stars

229

230

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

232 """

233 Add timestamp and elapsed time to the star DataFrame.

234

235 Parameters

236 ----------

237 stars : `pd.DataFrame`

238 DataFrame with star measurements.

239 guiderData : `GuiderData`

240 GuiderData instance containing guider data and metadata.

241 guiderName : `str`

242 Name of the guider.

243

244 Returns

245 -------

246 stars : `pd.DataFrame`

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

248 """

249 gd = guiderData

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

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

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

253

254 # get the timestamp for each stamp

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

256

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

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

259

260 stars["timestamp"] = timeStamp

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

262 return stars

263

264

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

266 """

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

268

269 Parameters

270 ----------

271 stars : `pd.DataFrame`

272 DataFrame with star measurements in ROI coordinates.

273 guiderData : `GuiderData`

274 GuiderData instance containing guider data and metadata.

275 guiderName : `str`

276 Name of the guider.

277

278 Returns

279 -------

280 stars : `pd.DataFrame`

281 DataFrame with added columns for CCD, focal plane,

282 and Alt/Az coordinates.

283 """

284 gd = guiderData

285 obsTime = gd.obsTime

286 detNum = gd.getGuiderDetNum(guiderName)

287 wcs = gd.getWcs(guiderName)

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

289

290 # Convert to CCD coordinates

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

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

293 )

294

295 # Convert to focal plane coordinates

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

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

298 )

299

300 # Convert to Alt/Az coordinates

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

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

303 )

304

305 # Convert fwhm to arcsec

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

307

308 return stars

309

310

311def applyQualityCuts(

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

313) -> pd.DataFrame:

314 """

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

316

317 Parameters

318 ----------

319 stars : `pd.DataFrame`

320 DataFrame with star measurements.

321 shape : `tuple[float, float]`

322 Shape of the ROI (cols, rows).

323 config : `GuiderStarTrackerConfig`

324 Configuration object with quality cut parameters.

325

326 Returns

327 -------

328 stars : `pd.DataFrame`

329 DataFrame with stars that passed the quality cuts.

330 """

331 minSnr = config.minSnr

332 maxEllipticity = config.maxEllipticity

333 edgeMargin = config.edgeMargin

334 roiCols, roiRows = shape

335

336 # Filter by minimum SNR

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

338

339 # Filter by minimum number of detections and maximum ellipticity

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

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

342 mask3 &= eabs <= maxEllipticity

343

344 # Add edge margin mask

345 mask4 = (

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

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

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

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

350 )

351 # Combine all masks

352 mask = mask1 & mask3 & mask4

353 return stars[mask].copy()

354

355

356def _diagnoseQualityCutRejections(

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

358) -> list[str]:

359 """Return per-row rejection reasons using the same masks as

360 applyQualityCuts.

361

362 Parameters

363 ----------

364 stars : `pd.DataFrame`

365 Pre-cut DataFrame (rows that were all rejected).

366 shape : `tuple[float, float]`

367 ROI shape (rows, cols).

368 config : `GuiderStarTrackerConfig`

369 Configuration with quality cut thresholds.

370

371 Returns

372 -------

373 reasons : `list[str]`

374 One string per row describing which cuts failed.

375 """

376 minSnr = config.minSnr

377 maxEllipticity = config.maxEllipticity

378 edgeMargin = config.edgeMargin

379 roiCols, roiRows = shape

380

381 # Same masks as applyQualityCuts

382 snrOk = (stars["snr"] >= minSnr) & (stars["flux"] > 0) & (stars["flux_err"] > 0)

383

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

385 ellipOk = (stars["e1"].abs() <= maxEllipticity) & (stars["e1"].abs() <= maxEllipticity)

386 ellipOk &= eabs <= maxEllipticity

387

388 edgeOk = (

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

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

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

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

393 )

394

395 reasons = []

396 for i, (_, row) in enumerate(stars.iterrows()):

397 r = []

398 if not snrOk.iloc[i]:

399 snr = row["snr"]

400 r.append(f"snr={snr:.1f}" if np.isfinite(snr) else "snr=NaN")

401 if not ellipOk.iloc[i]:

402 e = np.hypot(row["e1"], row["e2"])

403 r.append(f"e={e:.2f}" if np.isfinite(e) else "e=NaN")

404 if not edgeOk.iloc[i]:

405 r.append(f"edge(x={row['xroi']:.0f},y={row['yroi']:.0f})")

406 reasons.append("|".join(r) if r else "?")

407 return reasons

408

409

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

411 """

412 Assign unique IDs to tracked stars.

413

414 Parameters

415 ----------

416 guiderData : `GuiderData`

417 GuiderData instance containing guider data and metadata.

418 stars : `pd.DataFrame`

419 DataFrame with star measurements.

420

421 Returns

422 -------

423 stars : `pd.DataFrame`

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

425

426 """

427 # Create a numeric “global” starid:

428 detMap = guiderData.guiderNameMap

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

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

431 return stars

432

433

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

435 """

436 Compute the offsets for each star in the catalog.

437

438 Parameters

439 ----------

440 stars : `pd.DataFrame`

441 DataFrame with star measurements.

442

443 Returns

444 -------

445 stars : `pd.DataFrame`

446 DataFrame with added offset columns:

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

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

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

450 - dtheta : offset in rotator angle (arcsec)

451 - magoffset : magnitude offset (mmag)

452

453 """

454 # make reference positions

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

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

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

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

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

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

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

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

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

464

465 # Compute all your offsets

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

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

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

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

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

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

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

473

474 # Correct for cos(alt) in daz

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

476

477 # compute mag offset

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

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

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

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

482 # convert to mmag and replace inf with nan

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

484

485 # not detections have nan offsets

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

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

488

489 return stars

490

491

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

493 """

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

495 for each row in `stars`.

496

497 Parameters

498 ----------

499 stars : `pd.DataFrame`

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

501 and 'yerr' columns.

502

503 Returns

504 -------

505 stars : `pd.DataFrame`

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

507 """

508 mm = 0.001 # convert microns to mm

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

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

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

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

513

514 # Angle in degrees

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

516

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

518 denom = xfp**2 + yfp**2

519

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

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

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

523 sigmaTheta = np.degrees(sigmaThetaRad)

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

525

526 stars["theta"] = theta

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

528 return stars

529

530

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

532 """

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

534

535 Parameters

536 ----------

537 stars : `pd.DataFrame`

538 DataFrame with star measurements.

539 camRotAngleDeg : `float`

540 Camera rotator angle in degrees.

541

542 Returns

543 -------

544 stars : `pd.DataFrame`

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

546 """

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

548 stars["e1_altaz"] = e1_rot

549 stars["e2_altaz"] = e2_rot

550 return stars

551

552

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

554 """

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

556

557 Parameters

558 ----------

559 e1 : `float`

560 First ellipticity component.

561 e2 : `float`

562 Second ellipticity component.

563 theta_deg : `float`

564 Rotation angle in degrees.

565

566 Returns

567 -------

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

569 Rotated ellipticity components.

570 """

571 theta = np.deg2rad(theta_deg)

572 cos2t = np.cos(2 * theta)

573 sin2t = np.sin(2 * theta)

574 e1_rot = e1 * cos2t + e2 * sin2t

575 e2_rot = -e1 * sin2t + e2 * cos2t

576 return e1_rot, e2_rot