Coverage for python/lsst/summit/utils/guiders/reading.py: 24%

344 statements  

« prev     ^ index     » next       coverage.py v7.11.0, created at 2025-11-20 02:57 -0800

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 

23import logging 

24from typing import TYPE_CHECKING, Any 

25 

26__all__ = [ 

27 "GuiderReader", 

28 "convertRawStampsToView", 

29 "GuiderData", 

30] 

31 

32from dataclasses import dataclass 

33from functools import cached_property 

34 

35import matplotlib.pyplot as plt 

36import numpy as np 

37from astropy.time import Time 

38 

39from lsst.afw import cameraGeom 

40from lsst.afw.image import ExposureF, ImageF, MaskedImageF, VisitInfo 

41from lsst.daf.butler import Butler, DatasetNotFoundError 

42from lsst.meas.algorithms.stamps import Stamp, Stamps 

43from lsst.obs.lsst import LsstCam 

44from lsst.utils.plotting.figures import make_figure 

45 

46from .plotting import GuiderPlotter, labelDetector 

47from .transformation import ( 

48 getCamRotAngle, 

49 makeCcdToDvcsTransform, 

50 makeInitGuiderWcs, 

51 makeRoiBbox, 

52 roiImageToDvcs, 

53) 

54 

55if TYPE_CHECKING: 

56 from lsst.afw.cameraGeom import Camera, Detector 

57 from lsst.geom import SkyWcs 

58 

59 

60def make_subplot(nrows=1, ncols=1, **fig_kwargs): 

61 """Return (fig, axs) using LSST's make_figure.""" 

62 fig = make_figure(**fig_kwargs) 

63 axs = fig.subplots(nrows=nrows, ncols=ncols, squeeze=True) 

64 return fig, axs 

65 

66 

67@dataclass 

68class GuiderData: 

69 """ 

70 LSST guider data container. 

71 

72 Holds raw guider Stamps, WCS objects, and exposure metadata. Provides 

73 cached properties for timestamps, header information, and convenience 

74 (dict-like) access. 

75 

76 Parameters 

77 ---------- 

78 seqNum : `int` 

79 Sequence number. 

80 dayObs : `int` 

81 Day of observation (YYYYMMDD). 

82 guiderNameMap : `dict[str, int]` 

83 Detector name to detector number map. 

84 rawStampsMap : `dict[str, Stamps]` 

85 Detector name to raw Stamps map. 

86 wcsMap : `dict[str, SkyWcs]` 

87 Detector name to SkyWcs map. 

88 camRotAngle : `float` 

89 Camera rotation angle in degrees. 

90 isMedianSubtracted : `bool`, optional 

91 If True, the raw stamps have been subtracted. 

92 view : `str`, optional 

93 Output view: 'dvcs', 'ccd', or 'roi'. 

94 

95 Cached properties 

96 ----------------- 

97 guiderNames : sorted detector names. 

98 detNameMax : detector with most stamps. 

99 metadata : metadata dict for detNameMax. 

100 header : exposure-level header summary. 

101 expid : DAYOBS * 100000 + SEQNUM. 

102 timestamps : masked astropy Time array (TAI scale). 

103 stampsMap : Stamps for all detectors in requested view. 

104 roiAmpNames : detector → amplifier name map. 

105 guiderFrequency : stamp cadence in Hz. 

106 

107 Public methods 

108 -------------- 

109 printHeaderInfo() : print header summary. 

110 getWcs(det) : return SkyWcs for detector. 

111 getStampArrayCoadd(det) : median stack of stamps. 

112 getGuiderAmpName(det) : ROI amplifier name. 

113 getGuiderDetNum(det) : detector number. 

114 plotMosaic(...) : plot all detectors as a mosaic. 

115 plotStamp(...) : plot one stamp from a single detector. 

116 makeGif(...) : create a GIF animation of guider frames. 

117 

118 Iteration and indexing 

119 ---------------------- 

120 for det in guiderData # iterates detector names 

121 guiderData.items() # yields (det, Stamps) 

122 guiderData['R44_SG0'] # full Stamps object 

123 guiderData['R44_SG0', 3] # stamp ndarray 

124 guiderData[3] # stamp from detNameMax 

125 """ 

126 

127 seqNum: int 

128 dayObs: int 

129 guiderNameMap: dict[str, int] 

130 rawStampsMap: dict[str, Stamps] 

131 wcsMap: dict[str, "SkyWcs"] 

132 camRotAngle: float 

133 isMedianSubtracted: bool = False 

134 view: str = "dvcs" # view type, either 'dvcs' or 'ccd' or 'roi' 

135 

136 def __repr__(self) -> str: 

137 """ 

138 Return a compact one-line representation. 

139 

140 Returns 

141 ------- 

142 reprStr : `str` 

143 String summary with exposure id, number of stamps, view, and guider 

144 names. 

145 """ 

146 return ( 

147 f"GuiderData(seqNum={self.seqNum}, " 

148 f"dayObs={self.dayObs}, " 

149 f"nStamps={len(self)}, view='{self.view}', " 

150 f"guiders={self.guiderNames})" 

151 ) 

152 

153 def __len__(self) -> int: 

154 """ 

155 Return number of stamps for the detector with the maximum count. 

156 

157 Returns 

158 ------- 

159 nStamps : `int` 

160 Number of stamps (time samples) in the exposure for the primary 

161 detector. 

162 """ 

163 return len(self.rawStampsMap[self.detNameMax]) 

164 

165 @cached_property 

166 def expid(self) -> int: 

167 """ 

168 Exposure identifier. 

169 

170 Returns 

171 ------- 

172 expid : `int` 

173 The exposure id. 

174 """ 

175 return self.dayObs * 100000 + self.seqNum 

176 

177 @cached_property 

178 def guiderNames(self) -> list[str]: 

179 """ 

180 Names of the guider detectors. 

181 

182 Returns 

183 ------- 

184 guiderNames : `list[str]` 

185 Sorted list of guider detector names. 

186 """ 

187 return sorted(list(self.guiderNameMap.keys())) 

188 

189 @cached_property 

190 def guiderIds(self) -> list[int]: 

191 """ 

192 IDs of the guider detectors. 

193 

194 Returns 

195 ------- 

196 guiderIds : `list[int]` 

197 Sorted list of guider detector numeric IDs. 

198 """ 

199 return sorted(list(self.guiderNameMap.values())) 

200 

201 @cached_property 

202 def detNameMax(self) -> str: 

203 """ 

204 Detector name with the maximum number of stamps. 

205 

206 Returns 

207 ------- 

208 detName : `str` 

209 Detector name possessing the most stamps. 

210 """ 

211 return max(self.rawStampsMap, key=lambda k: len(self.rawStampsMap[k])) 

212 

213 @cached_property 

214 def metadata(self) -> dict: 

215 """ 

216 Metadata for the detector with the most stamps. 

217 

218 Returns 

219 ------- 

220 metadata : `dict` 

221 Copy of the detector-level metadata dictionary. 

222 """ 

223 return self.rawStampsMap[self.detNameMax].metadata.toDict() 

224 

225 @cached_property 

226 def header(self) -> dict[str, str | float | None]: 

227 """ 

228 Dictionary of header metadata for this GuiderData. 

229 Fields: n_stamps, freq, expid, filter, cam_rot_angle, start_time, 

230 roi_cols, roi_rows, shuttime, az_start, el_start, az_end, el_end, 

231 seeing. 

232 

233 Returns 

234 ------- 

235 header : `dict[str, str | float | None]` 

236 The headers. 

237 """ 

238 md = self.metadata 

239 info: dict[str, str | float | None] = {} 

240 info["n_stamps"] = len(self) 

241 info["freq"] = self.guiderFrequency 

242 # Ensure DAYOBS and SEQNUM are ints before math, keep ≤79 chars 

243 info["expid"] = self.expid 

244 info["filter"] = md.get("FILTBAND", "Unknown") 

245 info["cam_rot_angle"] = self.camRotAngle 

246 info["start_time"] = md.get("GDSSTART", None) 

247 info["roi_cols"] = int(md.get("ROICOLS", 0)) 

248 info["roi_rows"] = int(md.get("ROIROWS", 0)) 

249 info["shuttime"] = metadata_to_float(md, "SHUTTIME", np.nan) 

250 info["guider_duration"] = self.guiderDurationSec 

251 info["az_start"] = metadata_to_float(md, "AZSTART", np.nan) 

252 info["el_start"] = metadata_to_float(md, "ELSTART", np.nan) 

253 info["az_end"] = metadata_to_float(md, "AZEND", np.nan) 

254 info["el_end"] = metadata_to_float(md, "ELEND", np.nan) 

255 return info 

256 

257 def printHeaderInfo(self) -> None: 

258 """Print a concise summary of key header fields.""" 

259 print(f"Data Id: {self.expid}, filter-band: {self.header['filter']}") 

260 print(f"ROI Shape (row, col): {self.header['roi_rows']}, " f"{self.header['roi_cols']}") 

261 print(f"With nStamps {len(self)}" f" at {self.guiderFrequency} Hz") 

262 print( 

263 f"Acq. Start Time: {self.header['start_time']} \n" 

264 f"with readout duration: {self.header['guider_duration']:.2f} sec" 

265 ) 

266 

267 @cached_property 

268 def weather(self) -> WeatherInfo: 

269 """ 

270 Weather information from metadata. 

271 

272 Returns 

273 ------- 

274 weather : `WeatherInfo` 

275 Parsed weather information. 

276 """ 

277 return WeatherInfo.fromMetadata(self.metadata) 

278 

279 @cached_property 

280 def timestampMap(self) -> dict[str, Time]: 

281 """ 

282 Aligned timestamp arrays for all detectors. 

283 

284 Returns 

285 ------- 

286 timestampMap : `dict[str, Time]` 

287 Mapping from detector name to masked Time array (TAI scale). 

288 """ 

289 nStamps = len(self.rawStampsMap[self.detNameMax]) 

290 timestampMap: dict[str, Time] = {} 

291 for detName, stamps in self.rawStampsMap.items(): 

292 timestampMap[detName] = standardizeGuiderTimestamps(stamps, nStamps) 

293 return timestampMap 

294 

295 @cached_property 

296 def guiderFrequency(self) -> float: 

297 """ 

298 Guider stamp cadence. 

299 

300 Returns 

301 ------- 

302 frequency : `float` 

303 Median stamp frequency in Hz. 

304 """ 

305 timestamps = self.timestampMap[self.detNameMax] 

306 # extract only the unmasked MJD values 

307 jd = timestamps[~timestamps.mask].jd 

308 period_sec = np.median(np.diff(jd)) * 86400.0 

309 return float(1.0 / period_sec) 

310 

311 @cached_property 

312 def guiderDurationSec(self) -> float: 

313 """ 

314 Total guider duration. 

315 

316 Returns 

317 ------- 

318 durationSec : `float` 

319 Total elapsed guider duration in seconds (including final period). 

320 """ 

321 timestamps = self.timestampMap[self.detNameMax] 

322 jd = timestamps[~timestamps.mask].jd 

323 duration_sec = (jd[-1] - jd[0]) * 86400.0 + 1 / self.guiderFrequency 

324 return float(duration_sec) 

325 

326 @cached_property 

327 def stampsMap(self) -> dict[str, Stamps]: 

328 """ 

329 Stamps converted to requested view. 

330 

331 Returns 

332 ------- 

333 stampsMap : `dict[str, Stamps]` 

334 Mapping of detector name to converted Stamps object. 

335 """ 

336 result: dict[str, Stamps] = {} 

337 if self.view == "roi": 

338 return self.rawStampsMap 

339 else: 

340 for detName, rawStamps in self.rawStampsMap.items(): 

341 result[detName] = convertRawStampsToView( 

342 rawStamps, 

343 detName, 

344 len(self), 

345 view=self.view, 

346 ) 

347 return result 

348 

349 @cached_property 

350 def roiAmpNames(self) -> dict[str, str]: 

351 """ 

352 Amplifier names used in the ROI. 

353 

354 Returns 

355 ------- 

356 ampNames : `dict[str, str]` 

357 Mapping from detector name to amplifier name active in ROI. 

358 """ 

359 ampNames: dict[str, str] = {} 

360 for detName in self.rawStampsMap.keys(): 

361 md = self.rawStampsMap[detName].metadata.toDict() 

362 segment = md["ROISEG"] 

363 ampName = "C" + segment[7:] 

364 ampNames[detName] = ampName 

365 return ampNames 

366 

367 def getStampArrayCoadd(self, detName: str) -> np.ndarray: 

368 """ 

369 Get the median-stacked stamp across time for the specified detector. 

370 

371 The stack is computed across the time axis after optional median-row 

372 bias removal (controlled by ``self.isMedianSubtracted``). 

373 

374 Parameters 

375 ---------- 

376 detName : `str` 

377 Guider detector name. 

378 

379 Returns 

380 ------- 

381 coadd : `np.ndarray` 

382 Median of all stamp arrays (time axis collapsed). 

383 """ 

384 if detName not in self.stampsMap: 

385 raise KeyError(f"{detName!r} not present in stampsMap") 

386 

387 stamps = self.stampsMap[detName] 

388 if len(stamps) == 0: 

389 raise ValueError(f"No stamps found for detector {detName!r}") 

390 

391 # Collect arrays, with optional bias subtraction 

392 arrList = [self[detName, idx] for idx in range(len(stamps))] 

393 stack = np.nanmedian(arrList, axis=0) 

394 return stack 

395 

396 def getGuiderAmpName(self, detName: str) -> str: 

397 """ 

398 Return amplifier name for a guider detector. 

399 

400 Parameters 

401 ---------- 

402 detName : `str` 

403 Guider detector name. 

404 

405 Returns 

406 ------- 

407 ampName : `str` 

408 Amplifier name string. 

409 """ 

410 if detName not in self.roiAmpNames: 

411 raise ValueError(f"Detector {detName} not found in roiAmpNames.") 

412 return self.roiAmpNames[detName] 

413 

414 def getGuiderDetNum(self, detName: str) -> int: 

415 """ 

416 Return detector number for a guider detector name. 

417 

418 Parameters 

419 ---------- 

420 detName : `str` 

421 Guider detector name. 

422 

423 Returns 

424 ------- 

425 detNum : `int` 

426 The numeric detector id. 

427 """ 

428 if detName not in self.guiderNameMap: 

429 raise ValueError(f"Detector {detName} not found in guiderNameMap.") 

430 return self.guiderNameMap[detName] 

431 

432 def getWcs(self, detName: str) -> SkyWcs: 

433 """ 

434 Return the wcs for a guider detector. 

435 

436 Parameters 

437 ---------- 

438 detName : `str` 

439 Guider detector name. 

440 

441 Returns 

442 ------- 

443 wcs : `SkyWcs` 

444 SkyWcs for the detector. 

445 """ 

446 try: 

447 return self.wcsMap[detName] 

448 except KeyError as e: 

449 raise KeyError(f"{detName!r} not found in wcsMap") from e 

450 

451 @cached_property 

452 def obsTime(self) -> Time: 

453 """ 

454 Observation start time in TAI. 

455 

456 Returns 

457 ------- 

458 obsTime : `Time` 

459 Start time parsed from metadata (TAI scale). 

460 """ 

461 gdstart = self.header["start_time"] 

462 return Time(gdstart, format="isot", scale="tai") 

463 

464 @cached_property 

465 def missingStampsMap(self) -> dict[str, int]: 

466 """ 

467 Number of missing stamps in the primary detector. 

468 

469 Returns 

470 ------- 

471 nMissing : dict[str, int] 

472 Count of missing stamps (NaN timestamps) in detNameMax. 

473 """ 

474 nMissing: dict[str, int] = {} 

475 for detName, timestamps in self.timestampMap.items(): 

476 nMissing[detName] = np.sum(timestamps.mask) 

477 return nMissing 

478 

479 @cached_property 

480 def nMissingStamps(self) -> int: 

481 """ 

482 Total number of missing stamps across all detectors. 

483 

484 Returns 

485 ------- 

486 nMissing : int 

487 Total count of missing stamps across all detectors. 

488 """ 

489 return sum(self.missingStampsMap.values()) 

490 

491 # Iterable / dict-like helpers 

492 def __iter__(self): 

493 """Iterate over detector names in guiderNames order.""" 

494 return iter(self.guiderNames) 

495 

496 def items(self): 

497 """Yield (detName, stamps) pairs like dict.items().""" 

498 for det in self.guiderNames: 

499 yield det, self.stampsMap[det] 

500 

501 def keys(self): 

502 """Iterate over detector names (dict-like .keys()).""" 

503 return iter(self.guiderNames) 

504 

505 def values(self): 

506 """Iterate over Stamps objects in guiderNames order.""" 

507 for det in self.guiderNames: 

508 yield self.stampsMap[det] 

509 

510 def __getitem__(self, key): 

511 """ 

512 Direct stamp access helper. 

513 

514 Example usage: 

515 gd["R44_SG0"] -> full `Stamps` object 

516 gd["R44_SG0", 3] -> stamp #3 ndarray 

517 gd[3] -> stamp #3 from `detNameMax` 

518 

519 Parameters 

520 ---------- 

521 key : `str | int | slice | tuple` 

522 Access pattern: 

523 - 'DET' -> Stamps 

524 - ('DET', i) -> single stamp ndarray 

525 - i (int) -> stamp i from primary detector 

526 - ('DET', slice) -> list of ndarrays 

527 

528 Returns 

529 ------- 

530 result : `Stamps | np.ndarray | list[np.ndarray]` 

531 Retrieved object according to key specification. 

532 """ 

533 # Single detector name -> full Stamps 

534 if isinstance(key, str): 

535 return self.stampsMap[key] 

536 

537 # idx only -> assume detNameMax 

538 if isinstance(key, (int, slice)): 

539 key = (self.detNameMax, key) 

540 

541 if isinstance(key, tuple): 

542 if len(key) != 2: 

543 raise TypeError("Key must be (detName, idx) or (idx,)") 

544 detName, idx = key[0], key[1] 

545 

546 if detName not in self.stampsMap: 

547 raise KeyError(f"{detName!r} not found in stampsMap") 

548 

549 stamps = self.stampsMap[detName] 

550 # slice returns list of ndarrays 

551 if isinstance(idx, slice): 

552 arrays = [self._processStampArray(stamp, detName) for stamp in stamps[idx]] 

553 return arrays 

554 # int -> single ndarray 

555 return self._processStampArray(stamps[idx], detName) 

556 

557 raise TypeError("Invalid key type for GuiderData indexing.") 

558 

559 def _processStampArray(self, stamp: Stamp, detName: str) -> np.ndarray: 

560 """ 

561 Convert a Stamp to an ndarray with optional median-row subtraction. 

562 

563 Parameters 

564 ---------- 

565 stamp : `Stamp` 

566 Input stamp object. 

567 detName : `str` 

568 Detector name (for row-axis logic). 

569 

570 Returns 

571 ------- 

572 array : `np.ndarray` 

573 2D image array (bias-corrected if configured). 

574 """ 

575 arr = stamp.stamp_im.image.array 

576 return arr 

577 

578 @cached_property 

579 def plotter(self): 

580 return GuiderPlotter(self) 

581 

582 def plotStamp( 

583 self, detName: str, stampNum: int, plo: float = 90, phi: float = 99.5, figsize=(10, 8) 

584 ) -> plt.Figure: 

585 """ 

586 Plot a single guider stamp. 

587 

588 Parameters 

589 ---------- 

590 detName : `str` 

591 Detector name. 

592 stampNum : `int` 

593 Stamp index. 

594 plo : `float`, optional 

595 Lower percentile stretch. 

596 phi : `float`, optional 

597 Upper percentile stretch. 

598 figsize : `tuple`, optional 

599 Figure size. 

600 

601 Returns 

602 ------- 

603 fig : `plt.Figure` 

604 The resulting figure. 

605 """ 

606 fig, axs = make_subplot(nrows=1, ncols=1, figsize=figsize) 

607 img = self[detName, stampNum] if stampNum > 0 else self.getStampArrayCoadd(detName) 

608 _ = axs.imshow( 

609 img, 

610 origin="lower", 

611 vmin=np.nanpercentile(img, plo), 

612 vmax=np.nanpercentile(img, phi), 

613 cmap="Greys", 

614 ) 

615 _ = labelDetector(axs, detName) 

616 axs.set_xlabel("X (pixels)", fontsize=11) 

617 axs.set_ylabel("Y (pixels)", fontsize=11) 

618 axs.set_title(f"{self.expid}") 

619 return fig 

620 

621 

622class GuiderReader: 

623 """ 

624 Utility to fetch LSST guider data via Butler. 

625 

626 Example: 

627 from lsst.summit.utils.guiders.reading import GuiderReader 

628 from lsst.daf.butler import Butler 

629 butler = Butler("embargo", collections="LSSTCam/raw/guider") 

630 

631 seqNum, dayObs = 461, 20250425 

632 reader = GuiderReader(butler, view="dvcs") 

633 guiderData = reader.get(dayObs=dayObs, seqNum=seqNum) 

634 """ 

635 

636 def __init__(self, butler: Butler, view: str = "dvcs"): 

637 self.butler = butler 

638 self.view = view 

639 self.log = logging.getLogger(__name__) 

640 # Define camera objects 

641 self.camera = LsstCam.getCamera() 

642 

643 # Build guiderNameMap 

644 self.guiderNameMap: dict[str, int] = {} 

645 for detector in self.camera: 

646 if detector.getType() == cameraGeom.DetectorType.GUIDER: 

647 detName = detector.getName() 

648 self.guiderNameMap[detName] = detector.getId() 

649 

650 self.guiderDetNames = list(self.guiderNameMap.keys()) 

651 self.nGuiders = len(self.guiderNameMap) 

652 

653 def get( 

654 self, dayObs: int, seqNum: int, doSubtractMedian: bool = True, scienceDetNum: int = 94 

655 ) -> GuiderData: 

656 """ 

657 Retrieve guider data for a given dayObs / seqNum. 

658 

659 Parameters 

660 ---------- 

661 dayObs : `int` 

662 Day of observation in YYYYMMDD format. 

663 seqNum : `int` 

664 Sequence number. 

665 doSubtractMedian : `bool`, optional 

666 If True, subtract median row bias from each stamp upon access. 

667 scienceDetNum : `int`, optional 

668 Science detector number for WCS reference. 

669 

670 Returns 

671 ------- 

672 guiderData : `GuiderData` 

673 Assembled guider data object. 

674 """ 

675 # check if the guider name is swaped (dayObs < 20250509) 

676 # modifies self.guiderNameMap in place if necessary 

677 self.applyGuiderNameSwapIfNeeded(dayObs) 

678 

679 # 1. Get guider_raw stamps all guiders 

680 rawStampsDict = self.getGuiderRawStamps(dayObs, seqNum) 

681 

682 # set the number of stamps 

683 # determine the maximum number of stamps among all guiders 

684 nStampsList = [len(stamps) for stamps in rawStampsDict.values()] 

685 self.nStamps = max(nStampsList) 

686 

687 if self.nStamps <= 1: 

688 raise RuntimeError( 

689 f"Only {self.nStamps} stamps found for dayObs {dayObs}, seqNum {seqNum}. " 

690 "At least 2 stamps are required to create GuiderData." 

691 ) 

692 

693 # Create a visitinfo dataId to get the WCS and camera rotation angle 

694 visitInfo = getVisitInfo(self.butler, dayObs, seqNum, scienceDetNum) 

695 

696 # 2. Get WCS map for all guiders 

697 wcsMapDict = makeInitGuiderWcs(self.camera, visitInfo) 

698 

699 # 3. Get camera rotation angle 

700 camRotAngle = getCamRotAngle(visitInfo) 

701 

702 # 4. subtract median row bias if requested 

703 if doSubtractMedian: 

704 rawStampsDict = self.subtractMedianRowBias(rawStampsDict) 

705 

706 # 4. Pack everything into a GuiderData object 

707 guiderData = GuiderData( 

708 seqNum=seqNum, 

709 dayObs=dayObs, 

710 view=self.view, 

711 rawStampsMap=rawStampsDict, 

712 guiderNameMap=self.guiderNameMap, 

713 wcsMap=wcsMapDict, 

714 camRotAngle=camRotAngle, 

715 isMedianSubtracted=doSubtractMedian, 

716 ) 

717 return guiderData 

718 

719 def getGuiderRawStamps(self, dayObs: int, seqNum: int) -> dict[str, Stamps]: 

720 """ 

721 Fetch raw guider Stamps for all detectors. 

722 

723 Parameters 

724 ---------- 

725 dayObs : `int` 

726 Observation day (YYYYMMDD). 

727 seqNum : `int` 

728 Sequence number. 

729 

730 Returns 

731 ------- 

732 rawStamps : `dict[str, Stamps]` 

733 Mapping from detector name to raw Stamps object. 

734 """ 

735 rawStamps: dict[str, Stamps] = {} 

736 for detName, detNum in self.guiderNameMap.items(): 

737 try: 

738 rawStamps[detName] = self.butler.get( 

739 "guider_raw", 

740 day_obs=dayObs, 

741 seq_num=seqNum, 

742 detector=detNum, 

743 instrument="LSSTCam", 

744 ) 

745 except DatasetNotFoundError as e: 

746 raise DatasetNotFoundError(f"No data for {detName} on {dayObs=} {seqNum=}") from e 

747 return rawStamps 

748 

749 def subtractMedianRowBias(self, rawStampsDict: dict[str, Stamps]) -> dict[str, Stamps]: 

750 """ 

751 Subtract median row bias from all raw stamps in-place. 

752 

753 This modifies the rawStampsMap of the GuiderReader instance. 

754 

755 Parameters 

756 ---------- 

757 rawStampsDict : `dict[str, Stamps]` 

758 Mapping from detector name to raw Stamps object. 

759 """ 

760 axisRowMap = self.getAxisRowMap() 

761 for detName, stamps in rawStampsDict.items(): 

762 # Subtract median along the row axis 

763 rowAxis = axisRowMap[detName] 

764 for i in range(len(stamps)): 

765 # Get the stamp image array 

766 data = stamps[i].stamp_im.image.array 

767 # Subtract median along the specified row axis 

768 medianRow = np.nanmedian(data, axis=rowAxis, keepdims=True) 

769 # Update the stamp with the new data 

770 stamps[i].stamp_im.image.array = data - medianRow 

771 # Update the rawStampsMap in-place 

772 rawStampsDict[detName] = stamps 

773 return rawStampsDict 

774 

775 def getAxisRowMap(self) -> dict[str, int]: 

776 """ 

777 Axis mapping (row axis) for guider detectors. 

778 

779 Returns 

780 ------- 

781 axisMap : `dict[str, int]` 

782 Dictionary with detector names as keys and axis mapping as values. 

783 The axis mapping is 1 for X and 0 for Y. 

784 """ 

785 camera = LsstCam.getCamera() 

786 axisMap: dict[str, int] = {} 

787 for detName, detNum in self.guiderNameMap.items(): 

788 # Get the detector object from the camera 

789 detector = camera[detNum] 

790 # Determine the row axis based on the orientation 

791 nq = detector.getOrientation().getNQuarter() 

792 axisMap[detName] = 0 if nq % 2 == 0 else 1 

793 return axisMap 

794 

795 def applyGuiderNameSwapIfNeeded(self, dayObs: int) -> None: 

796 """ 

797 Apply guider name swap (SG0/SG1) for early data if required. 

798 

799 Parameters 

800 ---------- 

801 dayObs : `int` 

802 Observation day (YYYYMMDD) used to decide whether to swap. 

803 """ 

804 if getattr(self, "_guiderNameMapSwapped", False): 

805 return # Already swapped; do nothing 

806 if dayObs < 20250509: 

807 newMap = {} 

808 for detName, detNum in self.guiderNameMap.items(): 

809 if detName.endswith("SG0"): 

810 swapped = detName.replace("SG0", "SG1") 

811 elif detName.endswith("SG1"): 

812 swapped = detName.replace("SG1", "SG0") 

813 else: 

814 swapped = detName 

815 swappedDetNum = self.camera[swapped].getId() 

816 newMap[swapped] = swappedDetNum 

817 self.guiderNameMap = newMap 

818 self._guiderNameMapSwapped = True 

819 

820 

821def getVisitInfo(butler: Butler, dayObs: int, seqNum: int, scienceDetNum: int) -> VisitInfo: 

822 """ 

823 Retrieve VisitInfo for a given dayObs / seqNum. 

824 

825 Parameters 

826 ---------- 

827 butler : `Butler` 

828 Active Butler instance. 

829 dayObs : `int` 

830 Observation day (YYYYMMDD). 

831 seqNum : `int` 

832 Sequence number. 

833 scienceDetNum : `int` 

834 Science detector number. 

835 

836 Returns 

837 ------- 

838 visitInfo : `VisitInfo` 

839 VisitInfo object retrieved from butler. 

840 """ 

841 dataId = { 

842 "instrument": "LSSTCam", 

843 "day_obs": dayObs, 

844 "seq_num": seqNum, 

845 "detector": scienceDetNum, 

846 } 

847 visitInfo = butler.get("raw.visitInfo", dataId) 

848 return visitInfo 

849 

850 

851# Helper functions for stamp conversion 

852def _makeRoiTransforms(metadata: dict, detector: Detector, camera: Camera) -> tuple[tuple, str]: 

853 """ 

854 Construct ROI transforms and derive amplifier name. 

855 

856 Parameters 

857 ---------- 

858 metadata : `dict` 

859 Exposure-level metadata. 

860 detector : `Detector` 

861 The detector. 

862 camera : `Camera` 

863 The camera object. 

864 

865 Returns 

866 ------- 

867 result : `tuple[tuple, str]` 

868 ((ccdViewBbox, fwd, back), ampName) 

869 """ 

870 ampName = "C" + metadata["ROISEG"][7:] 

871 ccdViewBbox = makeRoiBbox(metadata, camera) 

872 fwd, back = makeCcdToDvcsTransform(ccdViewBbox, detector.getOrientation().getNQuarter()) 

873 return (ccdViewBbox, fwd, back), ampName 

874 

875 

876def _convertMaskedImage( 

877 maskedImage: MaskedImageF, 

878 stampMetadata: dict[str, Any], 

879 metadata: dict[str, Any], 

880 transforms: tuple, 

881 detector: Detector, 

882 ampName: str, 

883 camera: Camera, 

884 view: str, 

885) -> Stamp: 

886 """ 

887 Convert one masked image ROI and build a Stamp. 

888 

889 Parameters 

890 ---------- 

891 maskedImage : `MaskedImageF` 

892 Input masked image for the ROI. 

893 stampMetadata : `dict` 

894 Stamp-level metadata (modified in-place). 

895 metadata : `dict` 

896 Exposure-level metadata. 

897 transforms : `tuple` 

898 Tuple containing (ccdViewBbox, forwardTransform, inverseTransform). 

899 detector : `Detector` 

900 Camera detector. 

901 ampName : `str` 

902 Amplifier name. 

903 camera : `Camera` 

904 Camera object. 

905 view : `str` 

906 Output view ('dvcs', 'ccd', or 'roi'). 

907 

908 Returns 

909 ------- 

910 stamp : `Stamp` 

911 Converted stamp object. 

912 """ 

913 ccdViewBbox, fwd, back = transforms 

914 rawArray = maskedImage.getImage().getArray() 

915 dvcsArray = roiImageToDvcs(rawArray, metadata, detector, ampName, camera, view=view) 

916 outImg = MaskedImageF(dvcsArray) 

917 archiveElement = [ccdViewBbox, fwd, back] 

918 return Stamp(outImg, archiveElement, metadata=stampMetadata) 

919 

920 

921def _blankStamp( 

922 stampIdx: int, 

923 metadata: dict[str, Any], 

924 transforms: tuple, 

925) -> Stamp: 

926 """ 

927 Create a blank (zero) stamp for a missing index. 

928 

929 Parameters 

930 ---------- 

931 stampIdx : `int` 

932 Missing stamp index. 

933 metadata : `dict` 

934 Original metadata container. 

935 transforms : `tuple` 

936 (ccdViewBbox, forwardTransform, inverseTransform). 

937 

938 Returns 

939 ------- 

940 stamp : `Stamp` 

941 Placeholder blank stamp with NaN timestamp. 

942 """ 

943 ccdViewBbox, fwd, back = transforms 

944 nRows, nCols = int(metadata["ROIROWS"]), int(metadata["ROICOLS"]) 

945 blankArray = np.zeros((nRows, nCols), dtype=np.float32) 

946 blankImg = ExposureF(MaskedImageF(ImageF(array=blankArray))) 

947 missingMetadata = metadata.copy() 

948 missingMetadata["DAQSTAMP"] = stampIdx 

949 missingMetadata["STMPTMJD"] = np.nan 

950 archiveElement = [ccdViewBbox, fwd, back] 

951 return Stamp(blankImg, archiveElement, metadata=missingMetadata) 

952 

953 

954def convertRawStampsToView( 

955 rawStamps: Stamps, 

956 detName: str, 

957 nStamps: int, 

958 view: str = "dvcs", 

959) -> Stamps: 

960 """ 

961 Convert guider stamps from raw ROI to a requested view ('dvcs', 'ccd', or 

962 'roi'). Handles missing stamps and preserves metadata and order. 

963 

964 Parameters 

965 ---------- 

966 rawStamps : `Stamps` 

967 Input raw ROI stamps. 

968 detName : `str` 

969 Detector name. 

970 nStamps : `int` 

971 Target total number of stamps (after filling gaps). 

972 view : `str`, optional 

973 Output view: 'dvcs', 'ccd', or 'roi'. 

974 

975 Returns 

976 ------- 

977 stampsOut : `Stamps` 

978 Converted stamps with gaps filled by blank stamps. 

979 """ 

980 camera = LsstCam.getCamera() 

981 detector = camera[detName] 

982 metadata = rawStamps.metadata 

983 metadataDict = metadata.toDict() 

984 

985 # Ensure CCDSLOT matches swapped name if applicable 

986 metadata["CCDSLOT"] = detName[4:7] 

987 

988 # Align timestamps and find valid/missing indices 

989 timestamps = standardizeGuiderTimestamps(rawStamps, nStamps) 

990 validIndices = np.where(~timestamps.mask)[0].tolist() 

991 missingIndices = np.where(timestamps.mask)[0].tolist() 

992 

993 # Pre‑compute transforms once 

994 transforms, ampName = _makeRoiTransforms(metadataDict, detector, camera) 

995 

996 stampsDict: dict[int, Stamp] = {} 

997 mIdx = 0 # index into masked images list 

998 

999 for idx in validIndices: 

1000 maskedImage = rawStamps.getMaskedImages()[mIdx] 

1001 stampMeta = rawStamps[mIdx].metadata 

1002 stampMeta["DAQSTAMP"] = stampMeta.get("DAQSTAMP", idx) 

1003 stampsDict[idx] = _convertMaskedImage( 

1004 maskedImage, 

1005 stampMeta, 

1006 metadataDict, 

1007 transforms, 

1008 detector, 

1009 ampName, 

1010 camera, 

1011 view, 

1012 ) 

1013 mIdx += 1 

1014 

1015 # Fill gaps with blanks 

1016 for idx in missingIndices: 

1017 stampsDict[idx] = _blankStamp(idx, metadataDict, transforms) 

1018 

1019 # Assemble in order 

1020 stampList = [stampsDict[i] for i in range(nStamps)] 

1021 return Stamps(stampList, metadata, use_archive=True) 

1022 

1023 

1024def standardizeGuiderTimestamps(rawStamps: Stamps, nStamps: int) -> Time: 

1025 """ 

1026 Return a masked `Time` array of length `nStamps` for one guider. 

1027 

1028 Missing stamps are filled with `NaN` and masked so every detector shares a 

1029 uniform, sortable timestamp vector. Result is in MJD, `scale="tai"`. 

1030 

1031 Parameters 

1032 ---------- 

1033 rawStamps : `Stamps` 

1034 Input stamps (may have missing indices). 

1035 nStamps : `int` 

1036 Desired total length (maximum among guiders). 

1037 

1038 Returns 

1039 ------- 

1040 fullTimestamps : `Time` 

1041 Masked Time array (scale='tai') with inferred cadence and gaps masked. 

1042 """ 

1043 timestampsList = [stamp.metadata.get("STMPTMJD", np.nan) for stamp in rawStamps] 

1044 mjdArray = np.ma.masked_invalid(timestampsList) 

1045 timestamps = Time(mjdArray, format="mjd", scale="tai") 

1046 

1047 # infer frequency from valid timestamps 

1048 jdDiff = np.diff(timestamps.jd) 

1049 freqDays = np.nanmedian(jdDiff) 

1050 startJd = timestamps[0].jd 

1051 

1052 if np.isnan(freqDays) or freqDays <= 0: 

1053 raise ValueError( 

1054 f"Invalid frequency {freqDays} derived from timestamps. " 

1055 "Ensure that the timestamps are valid and evenly spaced." 

1056 ) 

1057 

1058 timestampsIdeal = Time(startJd + np.arange(nStamps) * freqDays, format="jd", scale="tai") 

1059 

1060 tolerance = 0.5 * freqDays 

1061 actualJd = timestamps.jd 

1062 

1063 # Build aligned timestamp list with np.nan where missing 

1064 mjdList = [] 

1065 for idealJd in timestampsIdeal.jd: 

1066 # check if any actual timestamp is close to this ideal timestamp 

1067 if np.any(np.abs(actualJd - idealJd) < tolerance): 

1068 # match found — find the closest 

1069 idx = np.argmin(np.abs(actualJd - idealJd)) 

1070 mjdList.append(actualJd[idx]) 

1071 else: 

1072 # no match — this is a missing timestamp 

1073 mjdList.append(np.nan) 

1074 

1075 mjdArray = np.ma.masked_invalid(mjdList) 

1076 fullTimestamps = Time(mjdArray, format="mjd", scale="tai") 

1077 fullTimestamps.freq = freqDays 

1078 return fullTimestamps 

1079 

1080 

1081@dataclass 

1082class WeatherInfo: 

1083 """ 

1084 Container for weather information. 

1085 

1086 Parameters 

1087 ---------- 

1088 temperature : `float` 

1089 Temperature in degrees Celsius. 

1090 humidity : `float` 

1091 Relative humidity in percent. 

1092 pressure : `float` 

1093 Atmospheric pressure in hPa. 

1094 windSpeed : `float` 

1095 Wind speed in m/s. 

1096 windDir : `float` 

1097 Wind direction in degrees from North. 

1098 """ 

1099 

1100 temperature: float 

1101 humidity: float 

1102 pressure: float 

1103 windSpeed: float 

1104 windDir: float 

1105 seeing: float 

1106 

1107 @staticmethod 

1108 def fromMetadata(metadata: dict) -> WeatherInfo: 

1109 """ 

1110 Create a WeatherInfo object from metadata dictionary. 

1111 

1112 Parameters 

1113 ---------- 

1114 metadata : `dict` 

1115 Metadata dictionary containing weather keys. 

1116 

1117 Returns 

1118 ------- 

1119 weather : `WeatherInfo` 

1120 Parsed weather information. 

1121 """ 

1122 return WeatherInfo( 

1123 temperature=metadata_to_float(metadata, "AIRTEMP", np.nan), 

1124 humidity=metadata_to_float(metadata, "HUMIDITY", np.nan), 

1125 pressure=metadata_to_float(metadata, "PRESSURE", np.nan), 

1126 windSpeed=metadata_to_float(metadata, "WINDSPD", np.nan), 

1127 windDir=metadata_to_float(metadata, "WINDDIR", np.nan), 

1128 seeing=metadata_to_float(metadata, "SEEING", np.nan), 

1129 ) 

1130 

1131 

1132def metadata_to_float(metadata: dict, key: str, default: float = np.nan) -> float: 

1133 """ 

1134 Safely convert a metadata value to float. 

1135 

1136 Parameters 

1137 ---------- 

1138 metadata : `dict` 

1139 Metadata dictionary. 

1140 key : `str` 

1141 Key to look up in the metadata. 

1142 default : `float`, optional 

1143 Default value if key is missing or conversion fails. 

1144 

1145 Returns 

1146 ------- 

1147 value : `float` 

1148 Converted float value or default. 

1149 """ 

1150 value = metadata[key] 

1151 if value is None: 

1152 return default 

1153 else: 

1154 return float(value)