Coverage for python/lsst/summit/utils/guiders/reading.py: 26%
410 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-09-17 02:41 -0700
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 biasPercentile: float = 10.0,
728 scienceDetNum: int = 94,
729 ) -> GuiderData:
730 """
731 Retrieve guider data for a given dayObs / seqNum.
733 Parameters
734 ----------
735 dayObs : `int`
736 Day of observation in YYYYMMDD format.
737 seqNum : `int`
738 Sequence number.
739 doSubtractMedian : `bool`, optional
740 If True, subtract column bias from each stamp.
741 columnMaskK : `float`, optional
742 Threshold factor for column mask detection.
743 biasPercentile : `float`, optional
744 Percentile (0-100) for column bias estimation.
745 Lower values avoid star/trail contamination.
746 scienceDetNum : `int`, optional
747 Science detector number for WCS reference.
749 Returns
750 -------
751 guiderData : `GuiderData`
752 Assembled guider data object.
753 """
754 # Check if the guider name is swapped (dayObs < 20250509)
755 # modifies self.guiderNameMap in place if necessary
756 self.applyGuiderNameSwapIfNeeded(dayObs)
758 rawStampsDict = self.getGuiderRawStamps(dayObs, seqNum)
760 # Determine the maximum number of stamps among all guiders
761 nStampsList = [len(stamps) for stamps in rawStampsDict.values()]
762 nStamps = max(nStampsList)
764 if nStamps <= 1:
765 raise RuntimeError(
766 f"Only {nStamps} stamps found for dayObs {dayObs}, seqNum {seqNum}. "
767 "At least 2 stamps are required to create GuiderData."
768 )
770 visitInfo = getVisitInfo(self.butler, dayObs, seqNum, scienceDetNum)
771 wcsMapDict = makeInitGuiderWcs(self.camera, visitInfo)
772 camRotAngle = getCamRotAngle(visitInfo)
774 processedStampsDict = self.processStamps(
775 rawStampsDict,
776 doSubtractMedian,
777 columnMaskK,
778 biasPercentile,
779 )
780 guiderData = GuiderData(
781 seqNum=seqNum,
782 dayObs=dayObs,
783 view=self.view,
784 rawStampsMap=processedStampsDict,
785 guiderNameMap=self.guiderNameMap,
786 wcsMap=wcsMapDict,
787 camRotAngle=camRotAngle,
788 isMedianSubtracted=doSubtractMedian,
789 columnMaskK=columnMaskK,
790 )
791 return guiderData
793 def getGuiderRawStamps(self, dayObs: int, seqNum: int) -> dict[str, Stamps]:
794 """
795 Fetch raw guider Stamps for all detectors.
797 Parameters
798 ----------
799 dayObs : `int`
800 Observation day (YYYYMMDD).
801 seqNum : `int`
802 Sequence number.
804 Returns
805 -------
806 rawStamps : `dict[str, Stamps]`
807 Mapping from detector name to raw Stamps object.
808 """
809 rawStamps: dict[str, Stamps] = {}
810 for detName, detNum in self.guiderNameMap.items():
811 try:
812 rawStamps[detName] = self.butler.get(
813 "guider_raw",
814 day_obs=dayObs,
815 seq_num=seqNum,
816 detector=detNum,
817 instrument="LSSTCam",
818 )
819 except DatasetNotFoundError as e:
820 raise DatasetNotFoundError(f"No data for {detName} on {dayObs=} {seqNum=}") from e
821 return rawStamps
823 def processStamps(
824 self,
825 rawStampsDict: dict[str, Stamps],
826 doSubtractMedian: bool,
827 columnMaskK: float,
828 biasPercentile: float = 10.0,
829 ) -> dict[str, Stamps]:
830 """
831 Apply column bias subtraction and per-stamp column masking.
834 Parameters
835 ----------
836 rawStampsDict : `dict[str, Stamps]`
837 ROI view of stamps from butler `guider_raw`.
838 doSubtractMedian : `bool`, optional
839 If True, subtract column bias.
840 columnMaskK : `float`, optional
841 Threshold factor for column mask detection.
842 biasPercentile : `float`, optional
843 Percentile for column bias estimation.
845 Returns
846 -------
847 processedStamps : `dict[str, Stamps]`
848 New Stamps objects with processed (copied) data.
849 Per-stamp column masks are included in each stamp's MaskedImageF.
850 """
851 processedStamps: dict[str, Stamps] = {}
853 for detName, stamps in rawStampsDict.items():
854 if len(stamps) == 0:
855 processedStamps[detName] = stamps
856 continue
858 stampList: list[Stamp] = []
859 for i in range(len(stamps)):
860 # Work on a copy - never modify original
861 data = stamps[i].stamp_im.image.array.copy()
863 # Compute column bias using low percentile to avoid
864 # star contamination (median pulls dark dips at star cols).
865 percRows = np.nanpercentile(data, biasPercentile, axis=0)
866 rPerc = np.nanpercentile(data, biasPercentile)
867 medianValue = np.nanmedian(data)
868 # Replace NaN values (fully masked columns) with global median
869 percRows = np.where(np.isnan(percRows), medianValue, percRows)
871 # Compute column mask on bias-subtracted data
872 colMask = getColumnMask(
873 data - percRows[np.newaxis, :] + (rPerc - medianValue),
874 k=columnMaskK,
875 )
877 # Create mask array for MaskedImageF (BAD=1 for masked cols)
878 nRows, nCols = data.shape
879 maskArray = np.zeros((nRows, nCols), dtype=np.uint32)
880 maskArray[colMask] = 1 # Set BAD bit for masked columns
882 # Apply column bias subtraction
883 if doSubtractMedian:
884 data = data - percRows[np.newaxis, :] + (rPerc - medianValue)
886 # Fill masked columns with global median
887 if colMask.any():
888 medianValue = np.nanmedian(data[~colMask]) if (~colMask).any() else np.nan
889 data[colMask] = medianValue
891 # Create MaskedImageF with image and mask
892 image = ImageF(array=data.astype(np.float32))
893 mask = MaskX(nCols, nRows)
894 mask.array[:] = maskArray
895 outImg = MaskedImageF(image, mask)
897 stampList.append(
898 Stamp(
899 stamp_im=outImg,
900 archive_element=stamps[i].archive_element,
901 metadata=stamps[i].metadata,
902 )
903 )
905 processedStamps[detName] = Stamps(stampList, stamps.metadata, use_mask=True, use_archive=True)
907 return processedStamps
909 def getAxisRowMap(self) -> dict[str, int]:
910 """
911 Axis mapping (row axis) for guider detectors.
913 Returns
914 -------
915 axisMap : `dict[str, int]`
916 Dictionary with detector names as keys and axis mapping as values.
917 The axis mapping is 1 for X and 0 for Y.
918 """
919 camera = LsstCam.getCamera()
920 axisMap: dict[str, int] = {}
921 for detName, detNum in self.guiderNameMap.items():
922 # Get the detector object from the camera
923 detector = camera[detNum]
924 # Determine the row axis based on the orientation
925 nq = detector.getOrientation().getNQuarter()
926 axisMap[detName] = 0 if nq % 2 == 0 else 1
927 return axisMap
929 def applyGuiderNameSwapIfNeeded(self, dayObs: int) -> None:
930 """
931 Apply guider name swap (SG0/SG1) for early data if required.
933 Parameters
934 ----------
935 dayObs : `int`
936 Observation day (YYYYMMDD) used to decide whether to swap.
937 """
938 if getattr(self, "_guiderNameMapSwapped", False):
939 return # Already swapped; do nothing
940 if dayObs < 20250509:
941 newMap = {}
942 for detName, detNum in self.guiderNameMap.items():
943 if detName.endswith("SG0"):
944 swapped = detName.replace("SG0", "SG1")
945 elif detName.endswith("SG1"):
946 swapped = detName.replace("SG1", "SG0")
947 else:
948 swapped = detName
949 swappedDetNum = self.camera[swapped].getId()
950 newMap[swapped] = swappedDetNum
951 self.guiderNameMap = newMap
952 self._guiderNameMapSwapped = True
955def getVisitInfo(butler: Butler, dayObs: int, seqNum: int, scienceDetNum: int) -> VisitInfo:
956 """
957 Retrieve VisitInfo for a given dayObs / seqNum.
959 Parameters
960 ----------
961 butler : `Butler`
962 Active Butler instance.
963 dayObs : `int`
964 Observation day (YYYYMMDD).
965 seqNum : `int`
966 Sequence number.
967 scienceDetNum : `int`
968 Science detector number.
970 Returns
971 -------
972 visitInfo : `VisitInfo`
973 VisitInfo object retrieved from butler.
974 """
975 dataId = {
976 "instrument": "LSSTCam",
977 "day_obs": dayObs,
978 "seq_num": seqNum,
979 "detector": scienceDetNum,
980 }
981 visitInfo = butler.get("raw.visitInfo", dataId)
982 return visitInfo
985# Helper functions for stamp conversion
986def _makeRoiTransforms(metadata: dict, detector: Detector, camera: Camera) -> tuple[tuple, str]:
987 """
988 Construct ROI transforms and derive amplifier name.
990 Parameters
991 ----------
992 metadata : `dict`
993 Exposure-level metadata.
994 detector : `Detector`
995 The detector.
996 camera : `Camera`
997 The camera object.
999 Returns
1000 -------
1001 result : `tuple[tuple, str]`
1002 ((ccdViewBbox, fwd, back), ampName)
1003 """
1004 ampName = "C" + metadata["ROISEG"][7:]
1005 ccdViewBbox = makeRoiBbox(metadata, camera)
1006 fwd, back = makeCcdToDvcsTransform(ccdViewBbox, detector.getOrientation().getNQuarter())
1007 return (ccdViewBbox, fwd, back), ampName
1010def _convertMaskedImage(
1011 maskedImage: MaskedImageF,
1012 stampMetadata: dict[str, Any],
1013 metadata: dict[str, Any],
1014 transforms: tuple,
1015 detector: Detector,
1016 ampName: str,
1017 camera: Camera,
1018 view: str,
1019) -> Stamp:
1020 """
1021 Convert one masked image ROI and build a Stamp.
1023 Parameters
1024 ----------
1025 maskedImage : `MaskedImageF`
1026 Input masked image for the ROI.
1027 stampMetadata : `dict`
1028 Stamp-level metadata (modified in-place).
1029 metadata : `dict`
1030 Exposure-level metadata.
1031 transforms : `tuple`
1032 Tuple containing (ccdViewBbox, forwardTransform, inverseTransform).
1033 detector : `Detector`
1034 Camera detector.
1035 ampName : `str`
1036 Amplifier name.
1037 camera : `Camera`
1038 Camera object.
1039 view : `str`
1040 Output view ('dvcs', 'ccd', or 'roi').
1042 Returns
1043 -------
1044 stamp : `Stamp`
1045 Converted stamp object.
1046 """
1047 ccdViewBbox, fwd, back = transforms
1048 rawArray = maskedImage.getImage().getArray()
1049 dvcsArray = roiImageToDvcs(rawArray, metadata, detector, ampName, camera, view=view)
1050 outImg = MaskedImageF(dvcsArray)
1051 archiveElement = [ccdViewBbox, fwd, back]
1052 return Stamp(outImg, archiveElement, metadata=stampMetadata)
1055def _blankStamp(
1056 stampIdx: int,
1057 metadata: dict[str, Any],
1058 transforms: tuple,
1059) -> Stamp:
1060 """
1061 Create a blank (zero) stamp for a missing index.
1063 Parameters
1064 ----------
1065 stampIdx : `int`
1066 Missing stamp index.
1067 metadata : `dict`
1068 Original metadata container.
1069 transforms : `tuple`
1070 (ccdViewBbox, forwardTransform, inverseTransform).
1072 Returns
1073 -------
1074 stamp : `Stamp`
1075 Placeholder blank stamp with NaN timestamp.
1076 """
1077 ccdViewBbox, fwd, back = transforms
1078 nRows, nCols = int(metadata["ROIROWS"]), int(metadata["ROICOLS"])
1079 blankArray = np.zeros((nRows, nCols), dtype=np.float32)
1080 blankImg = ExposureF(MaskedImageF(ImageF(array=blankArray)))
1081 missingMetadata = metadata.copy()
1082 missingMetadata["DAQSTAMP"] = stampIdx
1083 missingMetadata["STMPTMJD"] = np.nan
1084 archiveElement = [ccdViewBbox, fwd, back]
1085 return Stamp(blankImg, archiveElement, metadata=missingMetadata)
1088def convertRawStampsToView(
1089 rawStamps: Stamps,
1090 detName: str,
1091 nStamps: int,
1092 view: str = "dvcs",
1093) -> Stamps:
1094 """
1095 Convert guider stamps from raw ROI to a requested view ('dvcs', 'ccd', or
1096 'roi'). Handles missing stamps and preserves metadata and order.
1098 Parameters
1099 ----------
1100 rawStamps : `Stamps`
1101 Input raw ROI stamps.
1102 detName : `str`
1103 Detector name.
1104 nStamps : `int`
1105 Target total number of stamps (after filling gaps).
1106 view : `str`, optional
1107 Output view: 'dvcs', 'ccd', or 'roi'.
1109 Returns
1110 -------
1111 stampsOut : `Stamps`
1112 Converted stamps with gaps filled by blank stamps.
1113 """
1114 camera = LsstCam.getCamera()
1115 detector = camera[detName]
1116 metadata = rawStamps.metadata
1117 metadataDict = metadata.toDict()
1119 # Ensure CCDSLOT matches swapped name if applicable
1120 metadata["CCDSLOT"] = detName[4:7]
1122 # Align timestamps and find valid/missing indices
1123 timestamps = standardizeGuiderTimestamps(rawStamps, nStamps)
1124 validIndices = np.where(~timestamps.mask)[0].tolist()
1125 missingIndices = np.where(timestamps.mask)[0].tolist()
1127 # Pre‑compute transforms once
1128 transforms, ampName = _makeRoiTransforms(metadataDict, detector, camera)
1130 stampsDict: dict[int, Stamp] = {}
1131 mIdx = 0 # index into masked images list
1133 # Column masking is now done in GuiderReader.processStamps() before this
1134 # function is called, so we no longer need to compute or apply it here.
1136 for idx in validIndices:
1137 maskedImage = rawStamps.getMaskedImages()[mIdx]
1138 stampMeta = rawStamps[mIdx].metadata
1139 stampMeta["DAQSTAMP"] = stampMeta.get("DAQSTAMP", idx)
1140 stampsDict[idx] = _convertMaskedImage(
1141 maskedImage, stampMeta, metadataDict, transforms, detector, ampName, camera, view
1142 )
1143 mIdx += 1
1145 # Fill gaps with blanks
1146 for idx in missingIndices:
1147 stampsDict[idx] = _blankStamp(idx, metadataDict, transforms)
1149 # Assemble in order
1150 stampList = [stampsDict[i] for i in range(nStamps)]
1151 return Stamps(stampList, metadata, use_archive=True)
1154def standardizeGuiderTimestamps(rawStamps: Stamps, nStamps: int) -> Time:
1155 """
1156 Return a masked `Time` array of length `nStamps` for one guider.
1158 Missing stamps are filled with `NaN` and masked so every detector shares a
1159 uniform, sortable timestamp vector. Result is in MJD, `scale="tai"`.
1161 Parameters
1162 ----------
1163 rawStamps : `Stamps`
1164 Input stamps (may have missing indices).
1165 nStamps : `int`
1166 Desired total length (maximum among guiders).
1168 Returns
1169 -------
1170 fullTimestamps : `Time`
1171 Masked Time array (scale='tai') with inferred cadence and gaps masked.
1172 """
1173 timestampsList = [stamp.metadata.get("STMPTMJD", np.nan) for stamp in rawStamps]
1174 mjdArray = np.ma.masked_invalid(timestampsList)
1175 timestamps = Time(mjdArray, format="mjd", scale="tai")
1177 # infer frequency from valid timestamps
1178 jdDiff = np.diff(timestamps.jd)
1179 freqDays = np.nanmedian(jdDiff)
1180 startJd = timestamps[0].jd
1182 if np.isnan(freqDays) or freqDays <= 0:
1183 raise ValueError(
1184 f"Invalid frequency {freqDays} derived from timestamps. "
1185 "Ensure that the timestamps are valid and evenly spaced."
1186 )
1188 timestampsIdeal = Time(startJd + np.arange(nStamps) * freqDays, format="jd", scale="tai")
1190 tolerance = 0.5 * freqDays
1191 actualJd = timestamps.jd
1193 # Build aligned timestamp list with np.nan where missing
1194 mjdList = []
1195 for idealJd in timestampsIdeal.jd:
1196 # check if any actual timestamp is close to this ideal timestamp
1197 if np.any(np.abs(actualJd - idealJd) < tolerance):
1198 # match found — find the closest
1199 idx = np.argmin(np.abs(actualJd - idealJd))
1200 mjdList.append(actualJd[idx])
1201 else:
1202 # no match — this is a missing timestamp
1203 mjdList.append(np.nan)
1205 mjdArray = np.ma.masked_invalid(mjdList)
1206 fullTimestamps = Time(mjdArray, format="mjd", scale="tai")
1207 fullTimestamps.freq = freqDays
1208 return fullTimestamps
1211@dataclass
1212class WeatherInfo:
1213 """
1214 Container for weather information.
1216 Parameters
1217 ----------
1218 temperature : `float`
1219 Temperature in degrees Celsius.
1220 humidity : `float`
1221 Relative humidity in percent.
1222 pressure : `float`
1223 Atmospheric pressure in hPa.
1224 windSpeed : `float`
1225 Wind speed in m/s.
1226 windDir : `float`
1227 Wind direction in degrees from North.
1228 """
1230 temperature: float
1231 humidity: float
1232 pressure: float
1233 windSpeed: float
1234 windDir: float
1235 seeing: float
1237 @staticmethod
1238 def fromMetadata(metadata: dict) -> WeatherInfo:
1239 """
1240 Create a WeatherInfo object from metadata dictionary.
1242 Parameters
1243 ----------
1244 metadata : `dict`
1245 Metadata dictionary containing weather keys.
1247 Returns
1248 -------
1249 weather : `WeatherInfo`
1250 Parsed weather information.
1251 """
1252 return WeatherInfo(
1253 temperature=metadata_to_float(metadata, "AIRTEMP", np.nan),
1254 humidity=metadata_to_float(metadata, "HUMIDITY", np.nan),
1255 pressure=metadata_to_float(metadata, "PRESSURE", np.nan),
1256 windSpeed=metadata_to_float(metadata, "WINDSPD", np.nan),
1257 windDir=metadata_to_float(metadata, "WINDDIR", np.nan),
1258 seeing=metadata_to_float(metadata, "SEEING", np.nan),
1259 )
1262def metadata_to_float(metadata: dict, key: str, default: float = np.nan) -> float:
1263 """
1264 Safely convert a metadata value to float.
1266 Parameters
1267 ----------
1268 metadata : `dict`
1269 Metadata dictionary.
1270 key : `str`
1271 Key to look up in the metadata.
1272 default : `float`, optional
1273 Default value if key is missing or conversion fails.
1275 Returns
1276 -------
1277 value : `float`
1278 Converted float value or default.
1279 """
1280 value = metadata[key]
1281 if value is None:
1282 return default
1283 else:
1284 return float(value)
1287def mad(x: np.ndarray) -> float:
1288 """
1289 Median absolute deviation metric for std.
1291 Parameters
1292 ----------
1293 x : `ndarray`
1294 Input array.
1296 Returns
1297 -------
1298 value : `float`
1299 Robust standard deviation
1300 """
1301 med = np.nanmedian(x)
1302 return 1.4826 * np.nanmedian(np.abs(x - med))
1305def getColumnMask(img: NDArray[np.floating], k: float = 6.0) -> NDArray[np.bool_]:
1306 """Return a boolean mask for bad columns based on robust median statistics.
1308 Parameters
1309 ----------
1310 img : `ndarray`
1311 2D image array.
1312 k : `float`, optional
1313 Threshold factor for MAD.
1315 Returns
1316 -------
1317 mask : `ndarray`
1318 Boolean mask of same shape as img, ``True`` for bad columns.
1319 """
1320 columnMedian = np.nanmedian(img, axis=0)
1321 center = np.nanmedian(columnMedian)
1323 # MAD of column medians
1324 diffs = np.abs(columnMedian - center)
1325 s = np.nanmedian(diffs)
1327 if s == 0 or not np.isfinite(s):
1328 return np.zeros(img.shape, dtype=bool)
1330 badCols = diffs > k * s
1331 return np.broadcast_to(badCols, img.shape).copy()