Coverage for python/lsst/summit/utils/guiders/reading.py: 24%
344 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-14 02:25 -0800
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-14 02:25 -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
23import logging
24from typing import TYPE_CHECKING, Any
26__all__ = [
27 "GuiderReader",
28 "convertRawStampsToView",
29 "GuiderData",
30]
32from dataclasses import dataclass
33from functools import cached_property
35import matplotlib.pyplot as plt
36import numpy as np
37from astropy.time import Time
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
46from .plotting import GuiderPlotter, labelDetector
47from .transformation import (
48 getCamRotAngle,
49 makeCcdToDvcsTransform,
50 makeInitGuiderWcs,
51 makeRoiBbox,
52 roiImageToDvcs,
53)
55if TYPE_CHECKING:
56 from lsst.afw.cameraGeom import Camera, Detector
57 from lsst.geom import SkyWcs
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
67@dataclass
68class GuiderData:
69 """
70 LSST guider data container.
72 Holds raw guider Stamps, WCS objects, and exposure metadata. Provides
73 cached properties for timestamps, header information, and convenience
74 (dict-like) access.
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'.
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.
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.
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 """
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'
136 def __repr__(self) -> str:
137 """
138 Return a compact one-line representation.
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 )
153 def __len__(self) -> int:
154 """
155 Return number of stamps for the detector with the maximum count.
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])
165 @cached_property
166 def expid(self) -> int:
167 """
168 Exposure identifier.
170 Returns
171 -------
172 expid : `int`
173 The exposure id.
174 """
175 return self.dayObs * 100000 + self.seqNum
177 @cached_property
178 def guiderNames(self) -> list[str]:
179 """
180 Names of the guider detectors.
182 Returns
183 -------
184 guiderNames : `list[str]`
185 Sorted list of guider detector names.
186 """
187 return sorted(list(self.guiderNameMap.keys()))
189 @cached_property
190 def guiderIds(self) -> list[int]:
191 """
192 IDs of the guider detectors.
194 Returns
195 -------
196 guiderIds : `list[int]`
197 Sorted list of guider detector numeric IDs.
198 """
199 return sorted(list(self.guiderNameMap.values()))
201 @cached_property
202 def detNameMax(self) -> str:
203 """
204 Detector name with the maximum number of stamps.
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]))
213 @cached_property
214 def metadata(self) -> dict:
215 """
216 Metadata for the detector with the most stamps.
218 Returns
219 -------
220 metadata : `dict`
221 Copy of the detector-level metadata dictionary.
222 """
223 return self.rawStampsMap[self.detNameMax].metadata.toDict()
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.
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
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 )
267 @cached_property
268 def weather(self) -> WeatherInfo:
269 """
270 Weather information from metadata.
272 Returns
273 -------
274 weather : `WeatherInfo`
275 Parsed weather information.
276 """
277 return WeatherInfo.fromMetadata(self.metadata)
279 @cached_property
280 def timestampMap(self) -> dict[str, Time]:
281 """
282 Aligned timestamp arrays for all detectors.
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
295 @cached_property
296 def guiderFrequency(self) -> float:
297 """
298 Guider stamp cadence.
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)
311 @cached_property
312 def guiderDurationSec(self) -> float:
313 """
314 Total guider duration.
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)
326 @cached_property
327 def stampsMap(self) -> dict[str, Stamps]:
328 """
329 Stamps converted to requested view.
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
349 @cached_property
350 def roiAmpNames(self) -> dict[str, str]:
351 """
352 Amplifier names used in the ROI.
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
367 def getStampArrayCoadd(self, detName: str) -> np.ndarray:
368 """
369 Get the median-stacked stamp across time for the specified detector.
371 The stack is computed across the time axis after optional median-row
372 bias removal (controlled by ``self.isMedianSubtracted``).
374 Parameters
375 ----------
376 detName : `str`
377 Guider detector name.
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")
387 stamps = self.stampsMap[detName]
388 if len(stamps) == 0:
389 raise ValueError(f"No stamps found for detector {detName!r}")
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
396 def getGuiderAmpName(self, detName: str) -> str:
397 """
398 Return amplifier name for a guider detector.
400 Parameters
401 ----------
402 detName : `str`
403 Guider detector name.
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]
414 def getGuiderDetNum(self, detName: str) -> int:
415 """
416 Return detector number for a guider detector name.
418 Parameters
419 ----------
420 detName : `str`
421 Guider detector name.
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]
432 def getWcs(self, detName: str) -> SkyWcs:
433 """
434 Return the wcs for a guider detector.
436 Parameters
437 ----------
438 detName : `str`
439 Guider detector name.
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
451 @cached_property
452 def obsTime(self) -> Time:
453 """
454 Observation start time in TAI.
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")
464 @cached_property
465 def missingStampsMap(self) -> dict[str, int]:
466 """
467 Number of missing stamps in the primary detector.
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
479 @cached_property
480 def nMissingStamps(self) -> int:
481 """
482 Total number of missing stamps across all detectors.
484 Returns
485 -------
486 nMissing : int
487 Total count of missing stamps across all detectors.
488 """
489 return sum(self.missingStampsMap.values())
491 # Iterable / dict-like helpers
492 def __iter__(self):
493 """Iterate over detector names in guiderNames order."""
494 return iter(self.guiderNames)
496 def items(self):
497 """Yield (detName, stamps) pairs like dict.items()."""
498 for det in self.guiderNames:
499 yield det, self.stampsMap[det]
501 def keys(self):
502 """Iterate over detector names (dict-like .keys())."""
503 return iter(self.guiderNames)
505 def values(self):
506 """Iterate over Stamps objects in guiderNames order."""
507 for det in self.guiderNames:
508 yield self.stampsMap[det]
510 def __getitem__(self, key):
511 """
512 Direct stamp access helper.
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`
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
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]
537 # idx only -> assume detNameMax
538 if isinstance(key, (int, slice)):
539 key = (self.detNameMax, key)
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]
546 if detName not in self.stampsMap:
547 raise KeyError(f"{detName!r} not found in stampsMap")
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)
557 raise TypeError("Invalid key type for GuiderData indexing.")
559 def _processStampArray(self, stamp: Stamp, detName: str) -> np.ndarray:
560 """
561 Convert a Stamp to an ndarray with optional median-row subtraction.
563 Parameters
564 ----------
565 stamp : `Stamp`
566 Input stamp object.
567 detName : `str`
568 Detector name (for row-axis logic).
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
578 @cached_property
579 def plotter(self):
580 return GuiderPlotter(self)
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.
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.
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
622class GuiderReader:
623 """
624 Utility to fetch LSST guider data via Butler.
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")
631 seqNum, dayObs = 461, 20250425
632 reader = GuiderReader(butler, view="dvcs")
633 guiderData = reader.get(dayObs=dayObs, seqNum=seqNum)
634 """
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()
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()
650 self.guiderDetNames = list(self.guiderNameMap.keys())
651 self.nGuiders = len(self.guiderNameMap)
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.
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.
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)
679 # 1. Get guider_raw stamps all guiders
680 rawStampsDict = self.getGuiderRawStamps(dayObs, seqNum)
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)
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 )
693 # Create a visitinfo dataId to get the WCS and camera rotation angle
694 visitInfo = getVisitInfo(self.butler, dayObs, seqNum, scienceDetNum)
696 # 2. Get WCS map for all guiders
697 wcsMapDict = makeInitGuiderWcs(self.camera, visitInfo)
699 # 3. Get camera rotation angle
700 camRotAngle = getCamRotAngle(visitInfo)
702 # 4. subtract median row bias if requested
703 if doSubtractMedian:
704 rawStampsDict = self.subtractMedianRowBias(rawStampsDict)
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
719 def getGuiderRawStamps(self, dayObs: int, seqNum: int) -> dict[str, Stamps]:
720 """
721 Fetch raw guider Stamps for all detectors.
723 Parameters
724 ----------
725 dayObs : `int`
726 Observation day (YYYYMMDD).
727 seqNum : `int`
728 Sequence number.
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
749 def subtractMedianRowBias(self, rawStampsDict: dict[str, Stamps]) -> dict[str, Stamps]:
750 """
751 Subtract median row bias from all raw stamps in-place.
753 This modifies the rawStampsMap of the GuiderReader instance.
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
775 def getAxisRowMap(self) -> dict[str, int]:
776 """
777 Axis mapping (row axis) for guider detectors.
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
795 def applyGuiderNameSwapIfNeeded(self, dayObs: int) -> None:
796 """
797 Apply guider name swap (SG0/SG1) for early data if required.
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
821def getVisitInfo(butler: Butler, dayObs: int, seqNum: int, scienceDetNum: int) -> VisitInfo:
822 """
823 Retrieve VisitInfo for a given dayObs / seqNum.
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.
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
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.
856 Parameters
857 ----------
858 metadata : `dict`
859 Exposure-level metadata.
860 detector : `Detector`
861 The detector.
862 camera : `Camera`
863 The camera object.
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
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.
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').
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)
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.
929 Parameters
930 ----------
931 stampIdx : `int`
932 Missing stamp index.
933 metadata : `dict`
934 Original metadata container.
935 transforms : `tuple`
936 (ccdViewBbox, forwardTransform, inverseTransform).
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)
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.
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'.
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()
985 # Ensure CCDSLOT matches swapped name if applicable
986 metadata["CCDSLOT"] = detName[4:7]
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()
993 # Pre‑compute transforms once
994 transforms, ampName = _makeRoiTransforms(metadataDict, detector, camera)
996 stampsDict: dict[int, Stamp] = {}
997 mIdx = 0 # index into masked images list
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
1015 # Fill gaps with blanks
1016 for idx in missingIndices:
1017 stampsDict[idx] = _blankStamp(idx, metadataDict, transforms)
1019 # Assemble in order
1020 stampList = [stampsDict[i] for i in range(nStamps)]
1021 return Stamps(stampList, metadata, use_archive=True)
1024def standardizeGuiderTimestamps(rawStamps: Stamps, nStamps: int) -> Time:
1025 """
1026 Return a masked `Time` array of length `nStamps` for one guider.
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"`.
1031 Parameters
1032 ----------
1033 rawStamps : `Stamps`
1034 Input stamps (may have missing indices).
1035 nStamps : `int`
1036 Desired total length (maximum among guiders).
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")
1047 # infer frequency from valid timestamps
1048 jdDiff = np.diff(timestamps.jd)
1049 freqDays = np.nanmedian(jdDiff)
1050 startJd = timestamps[0].jd
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 )
1058 timestampsIdeal = Time(startJd + np.arange(nStamps) * freqDays, format="jd", scale="tai")
1060 tolerance = 0.5 * freqDays
1061 actualJd = timestamps.jd
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)
1075 mjdArray = np.ma.masked_invalid(mjdList)
1076 fullTimestamps = Time(mjdArray, format="mjd", scale="tai")
1077 fullTimestamps.freq = freqDays
1078 return fullTimestamps
1081@dataclass
1082class WeatherInfo:
1083 """
1084 Container for weather information.
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 """
1100 temperature: float
1101 humidity: float
1102 pressure: float
1103 windSpeed: float
1104 windDir: float
1105 seeing: float
1107 @staticmethod
1108 def fromMetadata(metadata: dict) -> WeatherInfo:
1109 """
1110 Create a WeatherInfo object from metadata dictionary.
1112 Parameters
1113 ----------
1114 metadata : `dict`
1115 Metadata dictionary containing weather keys.
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 )
1132def metadata_to_float(metadata: dict, key: str, default: float = np.nan) -> float:
1133 """
1134 Safely convert a metadata value to float.
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.
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)