Coverage for python/lsst/summit/utils/guiders/reading.py: 26%
409 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-14 11:15 +0000
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 collections.abc import Iterator
25from typing import TYPE_CHECKING, Any, overload
27__all__ = [
28 "GuiderReader",
29 "convertRawStampsToView",
30 "GuiderData",
31]
33from dataclasses import dataclass
34from functools import cached_property
36import matplotlib.pyplot as plt
37import numpy as np
38from astropy.time import Time
39from numpy.typing import NDArray
40from pydantic import BaseModel, ConfigDict, computed_field
42from lsst.afw import cameraGeom
43from lsst.afw.image import ExposureF, ImageF, MaskedImageF, MaskX, VisitInfo
44from lsst.daf.butler import Butler, DatasetNotFoundError
45from lsst.meas.algorithms.stamps import Stamp, Stamps
46from lsst.obs.lsst import LsstCam
47from lsst.utils.plotting.figures import make_figure
49from .plotting import GuiderPlotter, labelDetector
50from .transformation import (
51 getCamRotAngle,
52 makeCcdToDvcsTransform,
53 makeInitGuiderWcs,
54 makeRoiBbox,
55 roiImageToDvcs,
56)
58if TYPE_CHECKING:
59 from lsst.afw.cameraGeom import Camera, Detector
60 from lsst.geom import SkyWcs
63def make_subplot(nrows: int = 1, ncols: int = 1, **fig_kwargs: Any) -> tuple[plt.Figure, Any]:
64 """Return (fig, axs) using LSST's make_figure."""
65 fig = make_figure(**fig_kwargs)
66 axs = fig.subplots(nrows=nrows, ncols=ncols, squeeze=True)
67 return fig, axs
70class GuiderData(BaseModel):
71 """
72 LSST guider data container.
74 Holds raw guider Stamps, WCS objects, and exposure metadata. Provides
75 cached properties for timestamps, header information, and convenience
76 (dict-like) access.
78 Parameters
79 ----------
80 seqNum : `int`
81 Sequence number.
82 dayObs : `int`
83 Day of observation (YYYYMMDD).
84 guiderNameMap : `dict[str, int]`
85 Detector name to detector number map.
86 rawStampsMap : `dict[str, Stamps]`
87 Detector name to raw Stamps map.
88 wcsMap : `dict[str, SkyWcs]`
89 Detector name to SkyWcs map.
90 camRotAngle : `float`
91 Camera rotation angle in degrees.
92 isMedianSubtracted : `bool`
93 If True, the raw stamps have been median-subtracted.
94 columnMaskK : `float`
95 Threshold factor used for column mask detection.
96 view : `str`
97 Output view: 'dvcs', 'ccd', or 'roi'.
99 Cached properties
100 -----------------
101 guiderNames : sorted detector names.
102 detNameMax : detector with most stamps.
103 metadata : metadata dict for detNameMax.
104 header : exposure-level header summary.
105 expid : DAYOBS * 100000 + SEQNUM.
106 timestamps : masked astropy Time array (TAI scale).
107 stampsMap : Stamps for all detectors in requested view.
108 roiAmpNames : detector → amplifier name map.
109 guiderFrequency : stamp cadence in Hz.
111 Public methods
112 --------------
113 printHeaderInfo() : print header summary.
114 getWcs(det) : return SkyWcs for detector.
115 getStampArrayCoadd(det) : median stack of stamps.
116 getGuiderAmpName(det) : ROI amplifier name.
117 getGuiderDetNum(det) : detector number.
118 plotMosaic(...) : plot all detectors as a mosaic.
119 plotStamp(...) : plot one stamp from a single detector.
120 makeGif(...) : create a GIF animation of guider frames.
122 Iteration and indexing
123 ----------------------
124 for det in guiderData # iterates detector names
125 guiderData.items() # yields (det, Stamps)
126 guiderData['R44_SG0'] # full Stamps object
127 guiderData['R44_SG0', 3] # stamp ndarray
128 guiderData[3] # stamp from detNameMax
129 """
131 model_config = ConfigDict(frozen=True, arbitrary_types_allowed=True)
133 seqNum: int
134 dayObs: int
135 guiderNameMap: dict[str, int]
136 rawStampsMap: dict[str, Stamps]
137 wcsMap: dict[str, Any]
138 camRotAngle: float
139 isMedianSubtracted: bool = False
140 columnMaskK: float = 50.0
141 view: str = "dvcs"
143 def __repr__(self) -> str:
144 """
145 Return a compact one-line representation.
147 Returns
148 -------
149 reprStr : `str`
150 String summary with exposure id, number of stamps, view, and guider
151 names.
152 """
153 return (
154 f"GuiderData(seqNum={self.seqNum}, "
155 f"dayObs={self.dayObs}, "
156 f"nStamps={len(self)}, view='{self.view}', "
157 f"guiders={self.guiderNames})"
158 )
160 def __len__(self) -> int:
161 """
162 Return number of stamps for the detector with the maximum count.
164 Returns
165 -------
166 nStamps : `int`
167 Number of stamps (time samples) in the exposure for the primary
168 detector.
169 """
170 return len(self.rawStampsMap[self.detNameMax])
172 @computed_field # type: ignore[prop-decorator]
173 @cached_property
174 def expid(self) -> int:
175 """
176 Exposure identifier.
178 Returns
179 -------
180 expid : `int`
181 The exposure id.
182 """
183 return self.dayObs * 100000 + self.seqNum
185 @computed_field # type: ignore[prop-decorator]
186 @cached_property
187 def guiderNames(self) -> list[str]:
188 """
189 Names of the guider detectors.
191 Returns
192 -------
193 guiderNames : `list[str]`
194 Sorted list of guider detector names.
195 """
196 return sorted(list(self.guiderNameMap.keys()))
198 @computed_field # type: ignore[prop-decorator]
199 @cached_property
200 def guiderIds(self) -> list[int]:
201 """
202 IDs of the guider detectors.
204 Returns
205 -------
206 guiderIds : `list[int]`
207 Sorted list of guider detector numeric IDs.
208 """
209 return sorted(list(self.guiderNameMap.values()))
211 @computed_field # type: ignore[prop-decorator]
212 @cached_property
213 def detNameMax(self) -> str:
214 """
215 Detector name with the maximum number of stamps.
217 Returns
218 -------
219 detName : `str`
220 Detector name possessing the most stamps.
221 """
222 return max(self.rawStampsMap, key=lambda k: len(self.rawStampsMap[k]))
224 @computed_field # type: ignore[prop-decorator]
225 @cached_property
226 def metadata(self) -> dict:
227 """
228 Metadata for the detector with the most stamps.
230 Returns
231 -------
232 metadata : `dict`
233 Copy of the detector-level metadata dictionary.
234 """
235 return self.rawStampsMap[self.detNameMax].metadata.toDict()
237 @computed_field # type: ignore[prop-decorator]
238 @cached_property
239 def header(self) -> dict[str, str | float | None]:
240 """
241 Dictionary of header metadata for this GuiderData.
242 Fields: n_stamps, freq, expid, filter, cam_rot_angle, start_time,
243 roi_cols, roi_rows, shuttime, az_start, el_start, az_end, el_end,
244 seeing.
246 Returns
247 -------
248 header : `dict[str, str | float | None]`
249 The headers.
250 """
251 md = self.metadata
252 info: dict[str, str | float | None] = {}
253 info["n_stamps"] = len(self)
254 info["freq"] = self.guiderFrequency
255 # Ensure DAYOBS and SEQNUM are ints before math, keep ≤79 chars
256 info["expid"] = self.expid
257 info["filter"] = md.get("FILTBAND", "Unknown")
258 info["cam_rot_angle"] = self.camRotAngle
259 info["start_time"] = md.get("GDSSTART", None)
260 info["roi_cols"] = int(md.get("ROICOLS", 0))
261 info["roi_rows"] = int(md.get("ROIROWS", 0))
262 info["shuttime"] = metadata_to_float(md, "SHUTTIME", np.nan)
263 info["guider_duration"] = self.guiderDurationSec
264 info["az_start"] = metadata_to_float(md, "AZSTART", np.nan)
265 info["el_start"] = metadata_to_float(md, "ELSTART", np.nan)
266 info["az_end"] = metadata_to_float(md, "AZEND", np.nan)
267 info["el_end"] = metadata_to_float(md, "ELEND", np.nan)
268 return info
270 def printHeaderInfo(self) -> None:
271 """Print a concise summary of key header fields."""
272 print(f"Data Id: {self.expid}, filter-band: {self.header['filter']}")
273 print(f"ROI Shape (row, col): {self.header['roi_rows']}, {self.header['roi_cols']}")
274 print(f"With nStamps {len(self)} at {self.guiderFrequency} Hz")
275 print(
276 f"Acq. Start Time: {self.header['start_time']} \n"
277 f"with readout duration: {self.header['guider_duration']:.2f} sec"
278 )
280 @computed_field # type: ignore[prop-decorator]
281 @cached_property
282 def weather(self) -> WeatherInfo:
283 """
284 Weather information from metadata.
286 Returns
287 -------
288 weather : `WeatherInfo`
289 Parsed weather information.
290 """
291 return WeatherInfo.fromMetadata(self.metadata)
293 @computed_field # type: ignore[prop-decorator]
294 @cached_property
295 def timestampMap(self) -> dict[str, Time]:
296 """
297 Aligned timestamp arrays for all detectors.
299 Returns
300 -------
301 timestampMap : `dict[str, Time]`
302 Mapping from detector name to masked Time array (TAI scale).
303 """
304 nStamps = len(self.rawStampsMap[self.detNameMax])
305 timestampMap: dict[str, Time] = {}
306 for detName, stamps in self.rawStampsMap.items():
307 timestampMap[detName] = standardizeGuiderTimestamps(stamps, nStamps)
308 return timestampMap
310 @computed_field # type: ignore[prop-decorator]
311 @cached_property
312 def guiderFrequency(self) -> float:
313 """
314 Guider stamp cadence.
316 Returns
317 -------
318 frequency : `float`
319 Median stamp frequency in Hz.
320 """
321 timestamps = self.timestampMap[self.detNameMax]
322 # extract only the unmasked MJD values
323 jd = timestamps[~timestamps.mask].jd
324 period_sec = np.median(np.diff(jd)) * 86400.0
325 return float(1.0 / period_sec)
327 @computed_field # type: ignore[prop-decorator]
328 @cached_property
329 def guiderDurationSec(self) -> float:
330 """
331 Total guider duration.
333 Returns
334 -------
335 durationSec : `float`
336 Total elapsed guider duration in seconds (including final period).
337 """
338 timestamps = self.timestampMap[self.detNameMax]
339 jd = timestamps[~timestamps.mask].jd
340 duration_sec = (jd[-1] - jd[0]) * 86400.0 + 1 / self.guiderFrequency
341 return float(duration_sec)
343 @computed_field # type: ignore[prop-decorator]
344 @cached_property
345 def stampsMap(self) -> dict[str, Stamps]:
346 """
347 Stamps converted to requested view.
349 Returns
350 -------
351 stampsMap : `dict[str, Stamps]`
352 Mapping of detector name to converted Stamps object.
353 """
354 result: dict[str, Stamps] = {}
355 if self.view == "roi":
356 return self.rawStampsMap
357 else:
358 for detName, rawStamps in self.rawStampsMap.items():
359 result[detName] = convertRawStampsToView(
360 rawStamps,
361 detName,
362 len(self),
363 view=self.view,
364 )
365 return result
367 @computed_field # type: ignore[prop-decorator]
368 @cached_property
369 def roiAmpNames(self) -> dict[str, str]:
370 """
371 Amplifier names used in the ROI.
373 Returns
374 -------
375 ampNames : `dict[str, str]`
376 Mapping from detector name to amplifier name active in ROI.
377 """
378 ampNames: dict[str, str] = {}
379 for detName in self.rawStampsMap.keys():
380 md = self.rawStampsMap[detName].metadata.toDict()
381 segment = md["ROISEG"]
382 ampName = "C" + segment[7:]
383 ampNames[detName] = ampName
384 return ampNames
386 def getStampArrayCoadd(self, detName: str) -> np.ndarray:
387 """
388 Get the median-stacked stamp across time for the specified detector.
390 The stack is computed across the time axis after optional median-row
391 bias removal (controlled by ``self.isMedianSubtracted``).
393 Parameters
394 ----------
395 detName : `str`
396 Guider detector name.
398 Returns
399 -------
400 coadd : `np.ndarray`
401 Median of all stamp arrays (time axis collapsed).
402 """
403 if detName not in self.stampsMap:
404 raise KeyError(f"{detName!r} not present in stampsMap")
406 stamps = self.stampsMap[detName]
407 if len(stamps) == 0:
408 raise ValueError(f"No stamps found for detector {detName!r}")
410 # Collect arrays, with optional bias subtraction.
411 arrList = [self[detName, idx] for idx in range(len(stamps))]
412 # Add a uniform [0, 1) dither per stamp before the median to break
413 # integer quantization. Without this, the coadd pixel distribution
414 # can be so narrow that STDEVCLIP returns 0. (See DM-54263.)
415 stack = np.array(arrList, dtype=np.float32)
416 rng = np.random.default_rng(seed=0)
417 stack += rng.uniform(0, 1, size=stack.shape).astype(np.float32)
418 return np.nanmedian(stack, axis=0)
420 def getGuiderAmpName(self, detName: str) -> str:
421 """
422 Return amplifier name for a guider detector.
424 Parameters
425 ----------
426 detName : `str`
427 Guider detector name.
429 Returns
430 -------
431 ampName : `str`
432 Amplifier name string.
433 """
434 if detName not in self.roiAmpNames:
435 raise ValueError(f"Detector {detName} not found in roiAmpNames.")
436 return self.roiAmpNames[detName]
438 def getGuiderDetNum(self, detName: str) -> int:
439 """
440 Return detector number for a guider detector name.
442 Parameters
443 ----------
444 detName : `str`
445 Guider detector name.
447 Returns
448 -------
449 detNum : `int`
450 The numeric detector id.
451 """
452 if detName not in self.guiderNameMap:
453 raise ValueError(f"Detector {detName} not found in guiderNameMap.")
454 return self.guiderNameMap[detName]
456 def getWcs(self, detName: str) -> SkyWcs:
457 """
458 Return the wcs for a guider detector.
460 Parameters
461 ----------
462 detName : `str`
463 Guider detector name.
465 Returns
466 -------
467 wcs : `SkyWcs`
468 SkyWcs for the detector.
469 """
470 try:
471 return self.wcsMap[detName]
472 except KeyError as e:
473 raise KeyError(f"{detName!r} not found in wcsMap") from e
475 @computed_field # type: ignore[prop-decorator]
476 @cached_property
477 def obsTime(self) -> Time:
478 """
479 Observation start time in TAI.
481 Returns
482 -------
483 obsTime : `Time`
484 Start time parsed from metadata (TAI scale).
485 """
486 gdstart = self.header["start_time"]
487 return Time(gdstart, format="isot", scale="tai")
489 @computed_field # type: ignore[prop-decorator]
490 @cached_property
491 def missingStampsMap(self) -> dict[str, int]:
492 """
493 Number of missing stamps in the primary detector.
495 Returns
496 -------
497 nMissing : dict[str, int]
498 Count of missing stamps (NaN timestamps) in detNameMax.
499 """
500 nMissing: dict[str, int] = {}
501 for detName, timestamps in self.timestampMap.items():
502 nMissing[detName] = np.sum(timestamps.mask)
503 return nMissing
505 @computed_field # type: ignore[prop-decorator]
506 @cached_property
507 def nMissingStamps(self) -> int:
508 """
509 Total number of missing stamps across all detectors.
511 Returns
512 -------
513 nMissing : int
514 Total count of missing stamps across all detectors.
515 """
516 return sum(self.missingStampsMap.values())
518 @computed_field # type: ignore[prop-decorator]
519 @cached_property
520 def alt(self) -> float:
521 """Return altitude (el_start) from the guider data header."""
522 raw = self.header.get("el_start")
523 if raw is None:
524 raise KeyError("Header key 'el_start' is missing or None")
525 return float(raw)
527 @computed_field # type: ignore[prop-decorator]
528 @cached_property
529 def az(self) -> float:
530 """Return azimuth (az_start) from the guider data header."""
531 raw = self.header.get("az_start")
532 if raw is None:
533 raise KeyError("Header key 'az_start' is missing or None")
534 return float(raw)
536 # Iterable / dict-like helpers
537 def __iter__(self) -> Iterator[str]: # type: ignore[override]
538 """Iterate over detector names in guiderNames order."""
539 return iter(self.guiderNames)
541 def items(self) -> Iterator[tuple[str, Stamps]]:
542 """Yield (detName, stamps) pairs like dict.items()."""
543 for det in self.guiderNames:
544 yield det, self.stampsMap[det]
546 def keys(self) -> Iterator[str]:
547 """Iterate over detector names (dict-like .keys())."""
548 return iter(self.guiderNames)
550 def values(self) -> Iterator[Stamps]:
551 """Iterate over Stamps objects in guiderNames order."""
552 for det in self.guiderNames:
553 yield self.stampsMap[det]
555 @overload
556 def __getitem__(self, key: str) -> Stamps: ...
558 @overload
559 def __getitem__(self, key: int) -> np.ndarray: ...
561 @overload
562 def __getitem__(self, key: slice) -> list[np.ndarray]: ...
564 @overload
565 def __getitem__(self, key: tuple[str, int]) -> np.ndarray: ...
567 @overload
568 def __getitem__(self, key: tuple[str, slice]) -> list[np.ndarray]: ...
570 def __getitem__(
571 self, key: str | int | slice | tuple[str, int | slice]
572 ) -> Stamps | np.ndarray | list[np.ndarray]:
573 """
574 Direct stamp access helper.
576 Example usage:
577 gd["R44_SG0"] -> full `Stamps` object
578 gd["R44_SG0", 3] -> stamp #3 ndarray
579 gd[3] -> stamp #3 from `detNameMax`
581 Parameters
582 ----------
583 key : `str | int | slice | tuple`
584 Access pattern:
585 - 'DET' -> Stamps
586 - ('DET', i) -> single stamp ndarray
587 - i (int) -> stamp i from primary detector
588 - ('DET', slice) -> list of ndarrays
590 Returns
591 -------
592 result : `Stamps | np.ndarray | list[np.ndarray]`
593 Retrieved object according to key specification.
594 """
595 # Single detector name -> full Stamps
596 if isinstance(key, str):
597 return self.stampsMap[key]
599 # idx only -> assume detNameMax
600 if isinstance(key, (int, slice)):
601 key = (self.detNameMax, key)
603 if isinstance(key, tuple):
604 if len(key) != 2:
605 raise TypeError("Key must be (detName, idx) or (idx,)")
606 detName, idx = key[0], key[1]
608 if detName not in self.stampsMap:
609 raise KeyError(f"{detName!r} not found in stampsMap")
611 stamps = self.stampsMap[detName]
612 # slice returns list of ndarrays
613 if isinstance(idx, slice):
614 arrays = [self._processStampArray(stamp, detName) for stamp in stamps[idx]]
615 return arrays
616 # int -> single ndarray
617 return self._processStampArray(stamps[idx], detName)
619 raise TypeError("Invalid key type for GuiderData indexing.")
621 def _processStampArray(self, stamp: Stamp, detName: str) -> np.ndarray:
622 """
623 Convert a Stamp to an ndarray with optional median-row subtraction.
625 Parameters
626 ----------
627 stamp : `Stamp`
628 Input stamp object.
629 detName : `str`
630 Detector name (for row-axis logic).
632 Returns
633 -------
634 array : `np.ndarray`
635 2D image array (bias-corrected if configured).
636 """
637 arr = stamp.stamp_im.image.array
638 return arr
640 @computed_field # type: ignore[prop-decorator]
641 @cached_property
642 def plotter(self) -> GuiderPlotter:
643 return GuiderPlotter(self)
645 def plotStamp(
646 self,
647 detName: str,
648 stampNum: int,
649 plo: float = 90,
650 phi: float = 99.5,
651 figsize: tuple[float, float] = (10, 8),
652 ) -> plt.Figure:
653 """
654 Plot a single guider stamp.
656 Parameters
657 ----------
658 detName : `str`
659 Detector name.
660 stampNum : `int`
661 Stamp index.
662 plo : `float`, optional
663 Lower percentile stretch.
664 phi : `float`, optional
665 Upper percentile stretch.
666 figsize : `tuple`, optional
667 Figure size.
669 Returns
670 -------
671 fig : `plt.Figure`
672 The resulting figure.
673 """
674 fig, axs = make_subplot(nrows=1, ncols=1, figsize=figsize)
675 img = self[detName, stampNum] if stampNum > 0 else self.getStampArrayCoadd(detName)
676 _ = axs.imshow(
677 img,
678 origin="lower",
679 vmin=np.nanpercentile(img, plo),
680 vmax=np.nanpercentile(img, phi),
681 cmap="Greys",
682 )
683 _ = labelDetector(axs, detName)
684 axs.set_xlabel("X (pixels)", fontsize=11)
685 axs.set_ylabel("Y (pixels)", fontsize=11)
686 axs.set_title(f"{self.expid}")
687 return fig
690class GuiderReader:
691 """
692 Utility to fetch LSST guider data via Butler.
694 Example:
695 from lsst.summit.utils.guiders.reading import GuiderReader
696 from lsst.daf.butler import Butler
697 butler = Butler("embargo", collections="LSSTCam/raw/guider")
699 seqNum, dayObs = 461, 20250425
700 reader = GuiderReader(butler, view="dvcs")
701 guiderData = reader.get(dayObs=dayObs, seqNum=seqNum)
702 """
704 def __init__(self, butler: Butler, view: str = "dvcs"):
705 self.butler = butler
706 self.view = view
707 self.log = logging.getLogger(__name__)
708 # Define camera objects
709 self.camera = LsstCam.getCamera()
711 # Build guiderNameMap
712 self.guiderNameMap: dict[str, int] = {}
713 for detector in self.camera:
714 if detector.getType() == cameraGeom.DetectorType.GUIDER:
715 detName = detector.getName()
716 self.guiderNameMap[detName] = detector.getId()
718 self.guiderDetNames = list(self.guiderNameMap.keys())
719 self.nGuiders = len(self.guiderNameMap)
721 def get(
722 self,
723 dayObs: int,
724 seqNum: int,
725 doSubtractMedian: bool = True,
726 columnMaskK: float = 50.0,
727 scienceDetNum: int = 94,
728 ) -> GuiderData:
729 """
730 Retrieve guider data for a given dayObs / seqNum.
732 Parameters
733 ----------
734 dayObs : `int`
735 Day of observation in YYYYMMDD format.
736 seqNum : `int`
737 Sequence number.
738 doSubtractMedian : `bool`, optional
739 If True, subtract median row bias from each stamp.
740 columnMaskK : `float`, optional
741 Threshold factor for column mask detection.
742 scienceDetNum : `int`, optional
743 Science detector number for WCS reference.
745 Returns
746 -------
747 guiderData : `GuiderData`
748 Assembled guider data object.
749 """
750 # Check if the guider name is swapped (dayObs < 20250509)
751 # modifies self.guiderNameMap in place if necessary
752 self.applyGuiderNameSwapIfNeeded(dayObs)
754 rawStampsDict = self.getGuiderRawStamps(dayObs, seqNum)
756 # Determine the maximum number of stamps among all guiders
757 nStampsList = [len(stamps) for stamps in rawStampsDict.values()]
758 nStamps = max(nStampsList)
760 if nStamps <= 1:
761 raise RuntimeError(
762 f"Only {nStamps} stamps found for dayObs {dayObs}, seqNum {seqNum}. "
763 "At least 2 stamps are required to create GuiderData."
764 )
766 visitInfo = getVisitInfo(self.butler, dayObs, seqNum, scienceDetNum)
767 wcsMapDict = makeInitGuiderWcs(self.camera, visitInfo)
768 camRotAngle = getCamRotAngle(visitInfo)
770 processedStampsDict = self.processStamps(
771 rawStampsDict,
772 doSubtractMedian,
773 columnMaskK,
774 )
775 guiderData = GuiderData(
776 seqNum=seqNum,
777 dayObs=dayObs,
778 view=self.view,
779 rawStampsMap=processedStampsDict,
780 guiderNameMap=self.guiderNameMap,
781 wcsMap=wcsMapDict,
782 camRotAngle=camRotAngle,
783 isMedianSubtracted=doSubtractMedian,
784 columnMaskK=columnMaskK,
785 )
786 return guiderData
788 def getGuiderRawStamps(self, dayObs: int, seqNum: int) -> dict[str, Stamps]:
789 """
790 Fetch raw guider Stamps for all detectors.
792 Parameters
793 ----------
794 dayObs : `int`
795 Observation day (YYYYMMDD).
796 seqNum : `int`
797 Sequence number.
799 Returns
800 -------
801 rawStamps : `dict[str, Stamps]`
802 Mapping from detector name to raw Stamps object.
803 """
804 rawStamps: dict[str, Stamps] = {}
805 for detName, detNum in self.guiderNameMap.items():
806 try:
807 rawStamps[detName] = self.butler.get(
808 "guider_raw",
809 day_obs=dayObs,
810 seq_num=seqNum,
811 detector=detNum,
812 instrument="LSSTCam",
813 )
814 except DatasetNotFoundError as e:
815 raise DatasetNotFoundError(f"No data for {detName} on {dayObs=} {seqNum=}") from e
816 return rawStamps
818 def processStamps(
819 self,
820 rawStampsDict: dict[str, Stamps],
821 doSubtractMedian: bool,
822 columnMaskK: float,
823 ) -> dict[str, Stamps]:
824 """
825 Apply median row bias subtraction and per-stamp column masking.
828 Parameters
829 ----------
830 rawStampsDict : `dict[str, Stamps]`
831 ROI view of stamps from butler `guider_raw`.
832 doSubtractMedian : `bool`, optional
833 If True, subtract median row bias.
834 columnMaskK : `float`, optional
835 Threshold factor for column mask detection.
837 Returns
838 -------
839 processedStamps : `dict[str, Stamps]`
840 New Stamps objects with processed (copied) data.
841 Per-stamp column masks are included in each stamp's MaskedImageF.
842 """
843 processedStamps: dict[str, Stamps] = {}
845 for detName, stamps in rawStampsDict.items():
846 if len(stamps) == 0:
847 processedStamps[detName] = stamps
848 continue
850 stampList: list[Stamp] = []
851 for i in range(len(stamps)):
852 # Work on a copy - never modify original
853 data = stamps[i].stamp_im.image.array.copy()
855 # Compute median row bias
856 medianRows = np.nanmedian(data, axis=0)
857 medianValue = np.nanmedian(data)
858 # Replace NaN medians (fully masked columns) with global median
859 medianRows = np.where(np.isnan(medianRows), medianValue, medianRows)
861 # Compute column mask on bias-subtracted data
862 colMask = getColumnMask(data - medianRows[np.newaxis, :], k=columnMaskK)
864 # Create mask array for MaskedImageF (BAD=1 for masked cols)
865 nRows, nCols = data.shape
866 maskArray = np.zeros((nRows, nCols), dtype=np.uint32)
867 maskArray[colMask] = 1 # Set BAD bit for masked columns
869 # Apply median row bias subtraction
870 if doSubtractMedian:
871 data = data - medianRows[np.newaxis, :]
873 # Fill masked columns with global median
874 if colMask.any():
875 medianValue = np.nanmedian(data[~colMask]) if (~colMask).any() else np.nan
876 data[colMask] = medianValue
878 # Create MaskedImageF with image and mask
879 image = ImageF(array=data.astype(np.float32))
880 mask = MaskX(nCols, nRows)
881 mask.array[:] = maskArray
882 outImg = MaskedImageF(image, mask)
884 stampList.append(
885 Stamp(
886 stamp_im=outImg,
887 archive_element=stamps[i].archive_element,
888 metadata=stamps[i].metadata,
889 )
890 )
892 processedStamps[detName] = Stamps(stampList, stamps.metadata, use_mask=True, use_archive=True)
894 return processedStamps
896 def getAxisRowMap(self) -> dict[str, int]:
897 """
898 Axis mapping (row axis) for guider detectors.
900 Returns
901 -------
902 axisMap : `dict[str, int]`
903 Dictionary with detector names as keys and axis mapping as values.
904 The axis mapping is 1 for X and 0 for Y.
905 """
906 camera = LsstCam.getCamera()
907 axisMap: dict[str, int] = {}
908 for detName, detNum in self.guiderNameMap.items():
909 # Get the detector object from the camera
910 detector = camera[detNum]
911 # Determine the row axis based on the orientation
912 nq = detector.getOrientation().getNQuarter()
913 axisMap[detName] = 0 if nq % 2 == 0 else 1
914 return axisMap
916 def applyGuiderNameSwapIfNeeded(self, dayObs: int) -> None:
917 """
918 Apply guider name swap (SG0/SG1) for early data if required.
920 Parameters
921 ----------
922 dayObs : `int`
923 Observation day (YYYYMMDD) used to decide whether to swap.
924 """
925 if getattr(self, "_guiderNameMapSwapped", False):
926 return # Already swapped; do nothing
927 if dayObs < 20250509:
928 newMap = {}
929 for detName, detNum in self.guiderNameMap.items():
930 if detName.endswith("SG0"):
931 swapped = detName.replace("SG0", "SG1")
932 elif detName.endswith("SG1"):
933 swapped = detName.replace("SG1", "SG0")
934 else:
935 swapped = detName
936 swappedDetNum = self.camera[swapped].getId()
937 newMap[swapped] = swappedDetNum
938 self.guiderNameMap = newMap
939 self._guiderNameMapSwapped = True
942def getVisitInfo(butler: Butler, dayObs: int, seqNum: int, scienceDetNum: int) -> VisitInfo:
943 """
944 Retrieve VisitInfo for a given dayObs / seqNum.
946 Parameters
947 ----------
948 butler : `Butler`
949 Active Butler instance.
950 dayObs : `int`
951 Observation day (YYYYMMDD).
952 seqNum : `int`
953 Sequence number.
954 scienceDetNum : `int`
955 Science detector number.
957 Returns
958 -------
959 visitInfo : `VisitInfo`
960 VisitInfo object retrieved from butler.
961 """
962 dataId = {
963 "instrument": "LSSTCam",
964 "day_obs": dayObs,
965 "seq_num": seqNum,
966 "detector": scienceDetNum,
967 }
968 visitInfo = butler.get("raw.visitInfo", dataId)
969 return visitInfo
972# Helper functions for stamp conversion
973def _makeRoiTransforms(metadata: dict, detector: Detector, camera: Camera) -> tuple[tuple, str]:
974 """
975 Construct ROI transforms and derive amplifier name.
977 Parameters
978 ----------
979 metadata : `dict`
980 Exposure-level metadata.
981 detector : `Detector`
982 The detector.
983 camera : `Camera`
984 The camera object.
986 Returns
987 -------
988 result : `tuple[tuple, str]`
989 ((ccdViewBbox, fwd, back), ampName)
990 """
991 ampName = "C" + metadata["ROISEG"][7:]
992 ccdViewBbox = makeRoiBbox(metadata, camera)
993 fwd, back = makeCcdToDvcsTransform(ccdViewBbox, detector.getOrientation().getNQuarter())
994 return (ccdViewBbox, fwd, back), ampName
997def _convertMaskedImage(
998 maskedImage: MaskedImageF,
999 stampMetadata: dict[str, Any],
1000 metadata: dict[str, Any],
1001 transforms: tuple,
1002 detector: Detector,
1003 ampName: str,
1004 camera: Camera,
1005 view: str,
1006) -> Stamp:
1007 """
1008 Convert one masked image ROI and build a Stamp.
1010 Parameters
1011 ----------
1012 maskedImage : `MaskedImageF`
1013 Input masked image for the ROI.
1014 stampMetadata : `dict`
1015 Stamp-level metadata (modified in-place).
1016 metadata : `dict`
1017 Exposure-level metadata.
1018 transforms : `tuple`
1019 Tuple containing (ccdViewBbox, forwardTransform, inverseTransform).
1020 detector : `Detector`
1021 Camera detector.
1022 ampName : `str`
1023 Amplifier name.
1024 camera : `Camera`
1025 Camera object.
1026 view : `str`
1027 Output view ('dvcs', 'ccd', or 'roi').
1029 Returns
1030 -------
1031 stamp : `Stamp`
1032 Converted stamp object.
1033 """
1034 ccdViewBbox, fwd, back = transforms
1035 rawArray = maskedImage.getImage().getArray()
1036 dvcsArray = roiImageToDvcs(rawArray, metadata, detector, ampName, camera, view=view)
1037 outImg = MaskedImageF(dvcsArray)
1038 archiveElement = [ccdViewBbox, fwd, back]
1039 return Stamp(outImg, archiveElement, metadata=stampMetadata)
1042def _blankStamp(
1043 stampIdx: int,
1044 metadata: dict[str, Any],
1045 transforms: tuple,
1046) -> Stamp:
1047 """
1048 Create a blank (zero) stamp for a missing index.
1050 Parameters
1051 ----------
1052 stampIdx : `int`
1053 Missing stamp index.
1054 metadata : `dict`
1055 Original metadata container.
1056 transforms : `tuple`
1057 (ccdViewBbox, forwardTransform, inverseTransform).
1059 Returns
1060 -------
1061 stamp : `Stamp`
1062 Placeholder blank stamp with NaN timestamp.
1063 """
1064 ccdViewBbox, fwd, back = transforms
1065 nRows, nCols = int(metadata["ROIROWS"]), int(metadata["ROICOLS"])
1066 blankArray = np.zeros((nRows, nCols), dtype=np.float32)
1067 blankImg = ExposureF(MaskedImageF(ImageF(array=blankArray)))
1068 missingMetadata = metadata.copy()
1069 missingMetadata["DAQSTAMP"] = stampIdx
1070 missingMetadata["STMPTMJD"] = np.nan
1071 archiveElement = [ccdViewBbox, fwd, back]
1072 return Stamp(blankImg, archiveElement, metadata=missingMetadata)
1075def convertRawStampsToView(
1076 rawStamps: Stamps,
1077 detName: str,
1078 nStamps: int,
1079 view: str = "dvcs",
1080) -> Stamps:
1081 """
1082 Convert guider stamps from raw ROI to a requested view ('dvcs', 'ccd', or
1083 'roi'). Handles missing stamps and preserves metadata and order.
1085 Parameters
1086 ----------
1087 rawStamps : `Stamps`
1088 Input raw ROI stamps.
1089 detName : `str`
1090 Detector name.
1091 nStamps : `int`
1092 Target total number of stamps (after filling gaps).
1093 view : `str`, optional
1094 Output view: 'dvcs', 'ccd', or 'roi'.
1096 Returns
1097 -------
1098 stampsOut : `Stamps`
1099 Converted stamps with gaps filled by blank stamps.
1100 """
1101 camera = LsstCam.getCamera()
1102 detector = camera[detName]
1103 metadata = rawStamps.metadata
1104 metadataDict = metadata.toDict()
1106 # Ensure CCDSLOT matches swapped name if applicable
1107 metadata["CCDSLOT"] = detName[4:7]
1109 # Align timestamps and find valid/missing indices
1110 timestamps = standardizeGuiderTimestamps(rawStamps, nStamps)
1111 validIndices = np.where(~timestamps.mask)[0].tolist()
1112 missingIndices = np.where(timestamps.mask)[0].tolist()
1114 # Pre‑compute transforms once
1115 transforms, ampName = _makeRoiTransforms(metadataDict, detector, camera)
1117 stampsDict: dict[int, Stamp] = {}
1118 mIdx = 0 # index into masked images list
1120 # Column masking is now done in GuiderReader.processStamps() before this
1121 # function is called, so we no longer need to compute or apply it here.
1123 for idx in validIndices:
1124 maskedImage = rawStamps.getMaskedImages()[mIdx]
1125 stampMeta = rawStamps[mIdx].metadata
1126 stampMeta["DAQSTAMP"] = stampMeta.get("DAQSTAMP", idx)
1127 stampsDict[idx] = _convertMaskedImage(
1128 maskedImage, stampMeta, metadataDict, transforms, detector, ampName, camera, view
1129 )
1130 mIdx += 1
1132 # Fill gaps with blanks
1133 for idx in missingIndices:
1134 stampsDict[idx] = _blankStamp(idx, metadataDict, transforms)
1136 # Assemble in order
1137 stampList = [stampsDict[i] for i in range(nStamps)]
1138 return Stamps(stampList, metadata, use_archive=True)
1141def standardizeGuiderTimestamps(rawStamps: Stamps, nStamps: int) -> Time:
1142 """
1143 Return a masked `Time` array of length `nStamps` for one guider.
1145 Missing stamps are filled with `NaN` and masked so every detector shares a
1146 uniform, sortable timestamp vector. Result is in MJD, `scale="tai"`.
1148 Parameters
1149 ----------
1150 rawStamps : `Stamps`
1151 Input stamps (may have missing indices).
1152 nStamps : `int`
1153 Desired total length (maximum among guiders).
1155 Returns
1156 -------
1157 fullTimestamps : `Time`
1158 Masked Time array (scale='tai') with inferred cadence and gaps masked.
1159 """
1160 timestampsList = [stamp.metadata.get("STMPTMJD", np.nan) for stamp in rawStamps]
1161 mjdArray = np.ma.masked_invalid(timestampsList)
1162 timestamps = Time(mjdArray, format="mjd", scale="tai")
1164 # infer frequency from valid timestamps
1165 jdDiff = np.diff(timestamps.jd)
1166 freqDays = np.nanmedian(jdDiff)
1167 startJd = timestamps[0].jd
1169 if np.isnan(freqDays) or freqDays <= 0:
1170 raise ValueError(
1171 f"Invalid frequency {freqDays} derived from timestamps. "
1172 "Ensure that the timestamps are valid and evenly spaced."
1173 )
1175 timestampsIdeal = Time(startJd + np.arange(nStamps) * freqDays, format="jd", scale="tai")
1177 tolerance = 0.5 * freqDays
1178 actualJd = timestamps.jd
1180 # Build aligned timestamp list with np.nan where missing
1181 mjdList = []
1182 for idealJd in timestampsIdeal.jd:
1183 # check if any actual timestamp is close to this ideal timestamp
1184 if np.any(np.abs(actualJd - idealJd) < tolerance):
1185 # match found — find the closest
1186 idx = np.argmin(np.abs(actualJd - idealJd))
1187 mjdList.append(actualJd[idx])
1188 else:
1189 # no match — this is a missing timestamp
1190 mjdList.append(np.nan)
1192 mjdArray = np.ma.masked_invalid(mjdList)
1193 fullTimestamps = Time(mjdArray, format="mjd", scale="tai")
1194 fullTimestamps.freq = freqDays
1195 return fullTimestamps
1198@dataclass
1199class WeatherInfo:
1200 """
1201 Container for weather information.
1203 Parameters
1204 ----------
1205 temperature : `float`
1206 Temperature in degrees Celsius.
1207 humidity : `float`
1208 Relative humidity in percent.
1209 pressure : `float`
1210 Atmospheric pressure in hPa.
1211 windSpeed : `float`
1212 Wind speed in m/s.
1213 windDir : `float`
1214 Wind direction in degrees from North.
1215 """
1217 temperature: float
1218 humidity: float
1219 pressure: float
1220 windSpeed: float
1221 windDir: float
1222 seeing: float
1224 @staticmethod
1225 def fromMetadata(metadata: dict) -> WeatherInfo:
1226 """
1227 Create a WeatherInfo object from metadata dictionary.
1229 Parameters
1230 ----------
1231 metadata : `dict`
1232 Metadata dictionary containing weather keys.
1234 Returns
1235 -------
1236 weather : `WeatherInfo`
1237 Parsed weather information.
1238 """
1239 return WeatherInfo(
1240 temperature=metadata_to_float(metadata, "AIRTEMP", np.nan),
1241 humidity=metadata_to_float(metadata, "HUMIDITY", np.nan),
1242 pressure=metadata_to_float(metadata, "PRESSURE", np.nan),
1243 windSpeed=metadata_to_float(metadata, "WINDSPD", np.nan),
1244 windDir=metadata_to_float(metadata, "WINDDIR", np.nan),
1245 seeing=metadata_to_float(metadata, "SEEING", np.nan),
1246 )
1249def metadata_to_float(metadata: dict, key: str, default: float = np.nan) -> float:
1250 """
1251 Safely convert a metadata value to float.
1253 Parameters
1254 ----------
1255 metadata : `dict`
1256 Metadata dictionary.
1257 key : `str`
1258 Key to look up in the metadata.
1259 default : `float`, optional
1260 Default value if key is missing or conversion fails.
1262 Returns
1263 -------
1264 value : `float`
1265 Converted float value or default.
1266 """
1267 value = metadata[key]
1268 if value is None:
1269 return default
1270 else:
1271 return float(value)
1274def mad(x: np.ndarray) -> float:
1275 """
1276 Median absolute deviation metric for std.
1278 Parameters
1279 ----------
1280 x : `ndarray`
1281 Input array.
1283 Returns
1284 -------
1285 value : `float`
1286 Robust standard deviation
1287 """
1288 med = np.nanmedian(x)
1289 return 1.4826 * np.nanmedian(np.abs(x - med))
1292def getColumnMask(img: NDArray[np.floating], k: float = 6.0) -> NDArray[np.bool_]:
1293 """Return a boolean mask for bad columns based on robust median statistics.
1295 Parameters
1296 ----------
1297 img : `ndarray`
1298 2D image array.
1299 k : `float`, optional
1300 Threshold factor for MAD.
1302 Returns
1303 -------
1304 mask : `ndarray`
1305 Boolean mask of same shape as img, ``True`` for bad columns.
1306 """
1307 columnMedian = np.nanmedian(img, axis=0)
1308 center = np.nanmedian(columnMedian)
1310 # MAD of column medians
1311 diffs = np.abs(columnMedian - center)
1312 s = np.nanmedian(diffs)
1314 if s == 0 or not np.isfinite(s):
1315 return np.zeros(img.shape, dtype=bool)
1317 badCols = diffs > k * s
1318 return np.broadcast_to(badCols, img.shape).copy()