Coverage for python/lsst/summit/utils/guiders/reading.py: 23%
362 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-25 01:54 -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
38from numpy.typing import NDArray
40from lsst.afw import cameraGeom
41from lsst.afw.image import ExposureF, ImageF, MaskedImageF, VisitInfo
42from lsst.daf.butler import Butler, DatasetNotFoundError
43from lsst.meas.algorithms.stamps import Stamp, Stamps
44from lsst.obs.lsst import LsstCam
45from lsst.utils.plotting.figures import make_figure
47from .plotting import GuiderPlotter, labelDetector
48from .transformation import (
49 getCamRotAngle,
50 makeCcdToDvcsTransform,
51 makeInitGuiderWcs,
52 makeRoiBbox,
53 roiImageToDvcs,
54)
56if TYPE_CHECKING:
57 from lsst.afw.cameraGeom import Camera, Detector
58 from lsst.geom import SkyWcs
61def make_subplot(nrows=1, ncols=1, **fig_kwargs):
62 """Return (fig, axs) using LSST's make_figure."""
63 fig = make_figure(**fig_kwargs)
64 axs = fig.subplots(nrows=nrows, ncols=ncols, squeeze=True)
65 return fig, axs
68@dataclass
69class GuiderData:
70 """
71 LSST guider data container.
73 Holds raw guider Stamps, WCS objects, and exposure metadata. Provides
74 cached properties for timestamps, header information, and convenience
75 (dict-like) access.
77 Parameters
78 ----------
79 seqNum : `int`
80 Sequence number.
81 dayObs : `int`
82 Day of observation (YYYYMMDD).
83 guiderNameMap : `dict[str, int]`
84 Detector name to detector number map.
85 rawStampsMap : `dict[str, Stamps]`
86 Detector name to raw Stamps map.
87 wcsMap : `dict[str, SkyWcs]`
88 Detector name to SkyWcs map.
89 camRotAngle : `float`
90 Camera rotation angle in degrees.
91 isMedianSubtracted : `bool`, optional
92 If True, the raw stamps have been subtracted.
93 view : `str`, optional
94 Output view: 'dvcs', 'ccd', or 'roi'.
96 Cached properties
97 -----------------
98 guiderNames : sorted detector names.
99 detNameMax : detector with most stamps.
100 metadata : metadata dict for detNameMax.
101 header : exposure-level header summary.
102 expid : DAYOBS * 100000 + SEQNUM.
103 timestamps : masked astropy Time array (TAI scale).
104 stampsMap : Stamps for all detectors in requested view.
105 roiAmpNames : detector → amplifier name map.
106 guiderFrequency : stamp cadence in Hz.
108 Public methods
109 --------------
110 printHeaderInfo() : print header summary.
111 getWcs(det) : return SkyWcs for detector.
112 getStampArrayCoadd(det) : median stack of stamps.
113 getGuiderAmpName(det) : ROI amplifier name.
114 getGuiderDetNum(det) : detector number.
115 plotMosaic(...) : plot all detectors as a mosaic.
116 plotStamp(...) : plot one stamp from a single detector.
117 makeGif(...) : create a GIF animation of guider frames.
119 Iteration and indexing
120 ----------------------
121 for det in guiderData # iterates detector names
122 guiderData.items() # yields (det, Stamps)
123 guiderData['R44_SG0'] # full Stamps object
124 guiderData['R44_SG0', 3] # stamp ndarray
125 guiderData[3] # stamp from detNameMax
126 """
128 seqNum: int
129 dayObs: int
130 guiderNameMap: dict[str, int]
131 rawStampsMap: dict[str, Stamps]
132 wcsMap: dict[str, "SkyWcs"]
133 camRotAngle: float
134 isMedianSubtracted: bool = False
135 view: str = "dvcs" # view type, either 'dvcs' or 'ccd' or 'roi'
137 def __repr__(self) -> str:
138 """
139 Return a compact one-line representation.
141 Returns
142 -------
143 reprStr : `str`
144 String summary with exposure id, number of stamps, view, and guider
145 names.
146 """
147 return (
148 f"GuiderData(seqNum={self.seqNum}, "
149 f"dayObs={self.dayObs}, "
150 f"nStamps={len(self)}, view='{self.view}', "
151 f"guiders={self.guiderNames})"
152 )
154 def __len__(self) -> int:
155 """
156 Return number of stamps for the detector with the maximum count.
158 Returns
159 -------
160 nStamps : `int`
161 Number of stamps (time samples) in the exposure for the primary
162 detector.
163 """
164 return len(self.rawStampsMap[self.detNameMax])
166 @cached_property
167 def expid(self) -> int:
168 """
169 Exposure identifier.
171 Returns
172 -------
173 expid : `int`
174 The exposure id.
175 """
176 return self.dayObs * 100000 + self.seqNum
178 @cached_property
179 def guiderNames(self) -> list[str]:
180 """
181 Names of the guider detectors.
183 Returns
184 -------
185 guiderNames : `list[str]`
186 Sorted list of guider detector names.
187 """
188 return sorted(list(self.guiderNameMap.keys()))
190 @cached_property
191 def guiderIds(self) -> list[int]:
192 """
193 IDs of the guider detectors.
195 Returns
196 -------
197 guiderIds : `list[int]`
198 Sorted list of guider detector numeric IDs.
199 """
200 return sorted(list(self.guiderNameMap.values()))
202 @cached_property
203 def detNameMax(self) -> str:
204 """
205 Detector name with the maximum number of stamps.
207 Returns
208 -------
209 detName : `str`
210 Detector name possessing the most stamps.
211 """
212 return max(self.rawStampsMap, key=lambda k: len(self.rawStampsMap[k]))
214 @cached_property
215 def metadata(self) -> dict:
216 """
217 Metadata for the detector with the most stamps.
219 Returns
220 -------
221 metadata : `dict`
222 Copy of the detector-level metadata dictionary.
223 """
224 return self.rawStampsMap[self.detNameMax].metadata.toDict()
226 @cached_property
227 def header(self) -> dict[str, str | float | None]:
228 """
229 Dictionary of header metadata for this GuiderData.
230 Fields: n_stamps, freq, expid, filter, cam_rot_angle, start_time,
231 roi_cols, roi_rows, shuttime, az_start, el_start, az_end, el_end,
232 seeing.
234 Returns
235 -------
236 header : `dict[str, str | float | None]`
237 The headers.
238 """
239 md = self.metadata
240 info: dict[str, str | float | None] = {}
241 info["n_stamps"] = len(self)
242 info["freq"] = self.guiderFrequency
243 # Ensure DAYOBS and SEQNUM are ints before math, keep ≤79 chars
244 info["expid"] = self.expid
245 info["filter"] = md.get("FILTBAND", "Unknown")
246 info["cam_rot_angle"] = self.camRotAngle
247 info["start_time"] = md.get("GDSSTART", None)
248 info["roi_cols"] = int(md.get("ROICOLS", 0))
249 info["roi_rows"] = int(md.get("ROIROWS", 0))
250 info["shuttime"] = metadata_to_float(md, "SHUTTIME", np.nan)
251 info["guider_duration"] = self.guiderDurationSec
252 info["az_start"] = metadata_to_float(md, "AZSTART", np.nan)
253 info["el_start"] = metadata_to_float(md, "ELSTART", np.nan)
254 info["az_end"] = metadata_to_float(md, "AZEND", np.nan)
255 info["el_end"] = metadata_to_float(md, "ELEND", np.nan)
256 return info
258 def printHeaderInfo(self) -> None:
259 """Print a concise summary of key header fields."""
260 print(f"Data Id: {self.expid}, filter-band: {self.header['filter']}")
261 print(f"ROI Shape (row, col): {self.header['roi_rows']}, " f"{self.header['roi_cols']}")
262 print(f"With nStamps {len(self)}" f" at {self.guiderFrequency} Hz")
263 print(
264 f"Acq. Start Time: {self.header['start_time']} \n"
265 f"with readout duration: {self.header['guider_duration']:.2f} sec"
266 )
268 @cached_property
269 def weather(self) -> WeatherInfo:
270 """
271 Weather information from metadata.
273 Returns
274 -------
275 weather : `WeatherInfo`
276 Parsed weather information.
277 """
278 return WeatherInfo.fromMetadata(self.metadata)
280 @cached_property
281 def timestampMap(self) -> dict[str, Time]:
282 """
283 Aligned timestamp arrays for all detectors.
285 Returns
286 -------
287 timestampMap : `dict[str, Time]`
288 Mapping from detector name to masked Time array (TAI scale).
289 """
290 nStamps = len(self.rawStampsMap[self.detNameMax])
291 timestampMap: dict[str, Time] = {}
292 for detName, stamps in self.rawStampsMap.items():
293 timestampMap[detName] = standardizeGuiderTimestamps(stamps, nStamps)
294 return timestampMap
296 @cached_property
297 def guiderFrequency(self) -> float:
298 """
299 Guider stamp cadence.
301 Returns
302 -------
303 frequency : `float`
304 Median stamp frequency in Hz.
305 """
306 timestamps = self.timestampMap[self.detNameMax]
307 # extract only the unmasked MJD values
308 jd = timestamps[~timestamps.mask].jd
309 period_sec = np.median(np.diff(jd)) * 86400.0
310 return float(1.0 / period_sec)
312 @cached_property
313 def guiderDurationSec(self) -> float:
314 """
315 Total guider duration.
317 Returns
318 -------
319 durationSec : `float`
320 Total elapsed guider duration in seconds (including final period).
321 """
322 timestamps = self.timestampMap[self.detNameMax]
323 jd = timestamps[~timestamps.mask].jd
324 duration_sec = (jd[-1] - jd[0]) * 86400.0 + 1 / self.guiderFrequency
325 return float(duration_sec)
327 @cached_property
328 def stampsMap(self) -> dict[str, Stamps]:
329 """
330 Stamps converted to requested view.
332 Returns
333 -------
334 stampsMap : `dict[str, Stamps]`
335 Mapping of detector name to converted Stamps object.
336 """
337 result: dict[str, Stamps] = {}
338 if self.view == "roi":
339 return self.rawStampsMap
340 else:
341 for detName, rawStamps in self.rawStampsMap.items():
342 result[detName] = convertRawStampsToView(
343 rawStamps,
344 detName,
345 len(self),
346 view=self.view,
347 )
348 return result
350 @cached_property
351 def roiAmpNames(self) -> dict[str, str]:
352 """
353 Amplifier names used in the ROI.
355 Returns
356 -------
357 ampNames : `dict[str, str]`
358 Mapping from detector name to amplifier name active in ROI.
359 """
360 ampNames: dict[str, str] = {}
361 for detName in self.rawStampsMap.keys():
362 md = self.rawStampsMap[detName].metadata.toDict()
363 segment = md["ROISEG"]
364 ampName = "C" + segment[7:]
365 ampNames[detName] = ampName
366 return ampNames
368 def getStampArrayCoadd(self, detName: str) -> np.ndarray:
369 """
370 Get the median-stacked stamp across time for the specified detector.
372 The stack is computed across the time axis after optional median-row
373 bias removal (controlled by ``self.isMedianSubtracted``).
375 Parameters
376 ----------
377 detName : `str`
378 Guider detector name.
380 Returns
381 -------
382 coadd : `np.ndarray`
383 Median of all stamp arrays (time axis collapsed).
384 """
385 if detName not in self.stampsMap:
386 raise KeyError(f"{detName!r} not present in stampsMap")
388 stamps = self.stampsMap[detName]
389 if len(stamps) == 0:
390 raise ValueError(f"No stamps found for detector {detName!r}")
392 # Collect arrays, with optional bias subtraction
393 arrList = [self[detName, idx] for idx in range(len(stamps))]
394 stack = np.nanmedian(arrList, axis=0)
395 return stack
397 def getGuiderAmpName(self, detName: str) -> str:
398 """
399 Return amplifier name for a guider detector.
401 Parameters
402 ----------
403 detName : `str`
404 Guider detector name.
406 Returns
407 -------
408 ampName : `str`
409 Amplifier name string.
410 """
411 if detName not in self.roiAmpNames:
412 raise ValueError(f"Detector {detName} not found in roiAmpNames.")
413 return self.roiAmpNames[detName]
415 def getGuiderDetNum(self, detName: str) -> int:
416 """
417 Return detector number for a guider detector name.
419 Parameters
420 ----------
421 detName : `str`
422 Guider detector name.
424 Returns
425 -------
426 detNum : `int`
427 The numeric detector id.
428 """
429 if detName not in self.guiderNameMap:
430 raise ValueError(f"Detector {detName} not found in guiderNameMap.")
431 return self.guiderNameMap[detName]
433 def getWcs(self, detName: str) -> SkyWcs:
434 """
435 Return the wcs for a guider detector.
437 Parameters
438 ----------
439 detName : `str`
440 Guider detector name.
442 Returns
443 -------
444 wcs : `SkyWcs`
445 SkyWcs for the detector.
446 """
447 try:
448 return self.wcsMap[detName]
449 except KeyError as e:
450 raise KeyError(f"{detName!r} not found in wcsMap") from e
452 @cached_property
453 def obsTime(self) -> Time:
454 """
455 Observation start time in TAI.
457 Returns
458 -------
459 obsTime : `Time`
460 Start time parsed from metadata (TAI scale).
461 """
462 gdstart = self.header["start_time"]
463 return Time(gdstart, format="isot", scale="tai")
465 @cached_property
466 def missingStampsMap(self) -> dict[str, int]:
467 """
468 Number of missing stamps in the primary detector.
470 Returns
471 -------
472 nMissing : dict[str, int]
473 Count of missing stamps (NaN timestamps) in detNameMax.
474 """
475 nMissing: dict[str, int] = {}
476 for detName, timestamps in self.timestampMap.items():
477 nMissing[detName] = np.sum(timestamps.mask)
478 return nMissing
480 @cached_property
481 def nMissingStamps(self) -> int:
482 """
483 Total number of missing stamps across all detectors.
485 Returns
486 -------
487 nMissing : int
488 Total count of missing stamps across all detectors.
489 """
490 return sum(self.missingStampsMap.values())
492 # Iterable / dict-like helpers
493 def __iter__(self):
494 """Iterate over detector names in guiderNames order."""
495 return iter(self.guiderNames)
497 def items(self):
498 """Yield (detName, stamps) pairs like dict.items()."""
499 for det in self.guiderNames:
500 yield det, self.stampsMap[det]
502 def keys(self):
503 """Iterate over detector names (dict-like .keys())."""
504 return iter(self.guiderNames)
506 def values(self):
507 """Iterate over Stamps objects in guiderNames order."""
508 for det in self.guiderNames:
509 yield self.stampsMap[det]
511 def __getitem__(self, key):
512 """
513 Direct stamp access helper.
515 Example usage:
516 gd["R44_SG0"] -> full `Stamps` object
517 gd["R44_SG0", 3] -> stamp #3 ndarray
518 gd[3] -> stamp #3 from `detNameMax`
520 Parameters
521 ----------
522 key : `str | int | slice | tuple`
523 Access pattern:
524 - 'DET' -> Stamps
525 - ('DET', i) -> single stamp ndarray
526 - i (int) -> stamp i from primary detector
527 - ('DET', slice) -> list of ndarrays
529 Returns
530 -------
531 result : `Stamps | np.ndarray | list[np.ndarray]`
532 Retrieved object according to key specification.
533 """
534 # Single detector name -> full Stamps
535 if isinstance(key, str):
536 return self.stampsMap[key]
538 # idx only -> assume detNameMax
539 if isinstance(key, (int, slice)):
540 key = (self.detNameMax, key)
542 if isinstance(key, tuple):
543 if len(key) != 2:
544 raise TypeError("Key must be (detName, idx) or (idx,)")
545 detName, idx = key[0], key[1]
547 if detName not in self.stampsMap:
548 raise KeyError(f"{detName!r} not found in stampsMap")
550 stamps = self.stampsMap[detName]
551 # slice returns list of ndarrays
552 if isinstance(idx, slice):
553 arrays = [self._processStampArray(stamp, detName) for stamp in stamps[idx]]
554 return arrays
555 # int -> single ndarray
556 return self._processStampArray(stamps[idx], detName)
558 raise TypeError("Invalid key type for GuiderData indexing.")
560 def _processStampArray(self, stamp: Stamp, detName: str) -> np.ndarray:
561 """
562 Convert a Stamp to an ndarray with optional median-row subtraction.
564 Parameters
565 ----------
566 stamp : `Stamp`
567 Input stamp object.
568 detName : `str`
569 Detector name (for row-axis logic).
571 Returns
572 -------
573 array : `np.ndarray`
574 2D image array (bias-corrected if configured).
575 """
576 arr = stamp.stamp_im.image.array
577 return arr
579 @cached_property
580 def plotter(self):
581 return GuiderPlotter(self)
583 def plotStamp(
584 self, detName: str, stampNum: int, plo: float = 90, phi: float = 99.5, figsize=(10, 8)
585 ) -> plt.Figure:
586 """
587 Plot a single guider stamp.
589 Parameters
590 ----------
591 detName : `str`
592 Detector name.
593 stampNum : `int`
594 Stamp index.
595 plo : `float`, optional
596 Lower percentile stretch.
597 phi : `float`, optional
598 Upper percentile stretch.
599 figsize : `tuple`, optional
600 Figure size.
602 Returns
603 -------
604 fig : `plt.Figure`
605 The resulting figure.
606 """
607 fig, axs = make_subplot(nrows=1, ncols=1, figsize=figsize)
608 img = self[detName, stampNum] if stampNum > 0 else self.getStampArrayCoadd(detName)
609 _ = axs.imshow(
610 img,
611 origin="lower",
612 vmin=np.nanpercentile(img, plo),
613 vmax=np.nanpercentile(img, phi),
614 cmap="Greys",
615 )
616 _ = labelDetector(axs, detName)
617 axs.set_xlabel("X (pixels)", fontsize=11)
618 axs.set_ylabel("Y (pixels)", fontsize=11)
619 axs.set_title(f"{self.expid}")
620 return fig
623class GuiderReader:
624 """
625 Utility to fetch LSST guider data via Butler.
627 Example:
628 from lsst.summit.utils.guiders.reading import GuiderReader
629 from lsst.daf.butler import Butler
630 butler = Butler("embargo", collections="LSSTCam/raw/guider")
632 seqNum, dayObs = 461, 20250425
633 reader = GuiderReader(butler, view="dvcs")
634 guiderData = reader.get(dayObs=dayObs, seqNum=seqNum)
635 """
637 def __init__(self, butler: Butler, view: str = "dvcs"):
638 self.butler = butler
639 self.view = view
640 self.log = logging.getLogger(__name__)
641 # Define camera objects
642 self.camera = LsstCam.getCamera()
644 # Build guiderNameMap
645 self.guiderNameMap: dict[str, int] = {}
646 for detector in self.camera:
647 if detector.getType() == cameraGeom.DetectorType.GUIDER:
648 detName = detector.getName()
649 self.guiderNameMap[detName] = detector.getId()
651 self.guiderDetNames = list(self.guiderNameMap.keys())
652 self.nGuiders = len(self.guiderNameMap)
654 def get(
655 self, dayObs: int, seqNum: int, doSubtractMedian: bool = True, scienceDetNum: int = 94
656 ) -> GuiderData:
657 """
658 Retrieve guider data for a given dayObs / seqNum.
660 Parameters
661 ----------
662 dayObs : `int`
663 Day of observation in YYYYMMDD format.
664 seqNum : `int`
665 Sequence number.
666 doSubtractMedian : `bool`, optional
667 If True, subtract median row bias from each stamp upon access.
668 scienceDetNum : `int`, optional
669 Science detector number for WCS reference.
671 Returns
672 -------
673 guiderData : `GuiderData`
674 Assembled guider data object.
675 """
676 # check if the guider name is swaped (dayObs < 20250509)
677 # modifies self.guiderNameMap in place if necessary
678 self.applyGuiderNameSwapIfNeeded(dayObs)
680 # 1. Get guider_raw stamps all guiders
681 rawStampsDict = self.getGuiderRawStamps(dayObs, seqNum)
683 # set the number of stamps
684 # determine the maximum number of stamps among all guiders
685 nStampsList = [len(stamps) for stamps in rawStampsDict.values()]
686 self.nStamps = max(nStampsList)
688 if self.nStamps <= 1:
689 raise RuntimeError(
690 f"Only {self.nStamps} stamps found for dayObs {dayObs}, seqNum {seqNum}. "
691 "At least 2 stamps are required to create GuiderData."
692 )
694 # Create a visitinfo dataId to get the WCS and camera rotation angle
695 visitInfo = getVisitInfo(self.butler, dayObs, seqNum, scienceDetNum)
697 # 2. Get WCS map for all guiders
698 wcsMapDict = makeInitGuiderWcs(self.camera, visitInfo)
700 # 3. Get camera rotation angle
701 camRotAngle = getCamRotAngle(visitInfo)
703 # 4. subtract median row bias if requested
704 if doSubtractMedian:
705 rawStampsDict = self.subtractMedianRowBias(rawStampsDict)
707 # 4. Pack everything into a GuiderData object
708 guiderData = GuiderData(
709 seqNum=seqNum,
710 dayObs=dayObs,
711 view=self.view,
712 rawStampsMap=rawStampsDict,
713 guiderNameMap=self.guiderNameMap,
714 wcsMap=wcsMapDict,
715 camRotAngle=camRotAngle,
716 isMedianSubtracted=doSubtractMedian,
717 )
718 return guiderData
720 def getGuiderRawStamps(self, dayObs: int, seqNum: int) -> dict[str, Stamps]:
721 """
722 Fetch raw guider Stamps for all detectors.
724 Parameters
725 ----------
726 dayObs : `int`
727 Observation day (YYYYMMDD).
728 seqNum : `int`
729 Sequence number.
731 Returns
732 -------
733 rawStamps : `dict[str, Stamps]`
734 Mapping from detector name to raw Stamps object.
735 """
736 rawStamps: dict[str, Stamps] = {}
737 for detName, detNum in self.guiderNameMap.items():
738 try:
739 rawStamps[detName] = self.butler.get(
740 "guider_raw",
741 day_obs=dayObs,
742 seq_num=seqNum,
743 detector=detNum,
744 instrument="LSSTCam",
745 )
746 except DatasetNotFoundError as e:
747 raise DatasetNotFoundError(f"No data for {detName} on {dayObs=} {seqNum=}") from e
748 return rawStamps
750 def subtractMedianRowBias(self, rawStampsDict: dict[str, Stamps]) -> dict[str, Stamps]:
751 """
752 Subtract median row bias from all raw stamps in-place.
754 This modifies the rawStampsMap of the GuiderReader instance.
756 Parameters
757 ----------
758 rawStampsDict : `dict[str, Stamps]`
759 Mapping from detector name to raw Stamps object.
760 """
761 axisRowMap = self.getAxisRowMap()
762 for detName, stamps in rawStampsDict.items():
763 # Subtract median along the row axis
764 rowAxis = axisRowMap[detName]
765 for i in range(len(stamps)):
766 # Get the stamp image array
767 data = stamps[i].stamp_im.image.array
768 # Subtract median along the specified row axis
769 medianRow = np.nanmedian(data, axis=rowAxis, keepdims=True)
770 # Update the stamp with the new data
771 stamps[i].stamp_im.image.array = data - medianRow
772 # Update the rawStampsMap in-place
773 rawStampsDict[detName] = stamps
774 return rawStampsDict
776 def getAxisRowMap(self) -> dict[str, int]:
777 """
778 Axis mapping (row axis) for guider detectors.
780 Returns
781 -------
782 axisMap : `dict[str, int]`
783 Dictionary with detector names as keys and axis mapping as values.
784 The axis mapping is 1 for X and 0 for Y.
785 """
786 camera = LsstCam.getCamera()
787 axisMap: dict[str, int] = {}
788 for detName, detNum in self.guiderNameMap.items():
789 # Get the detector object from the camera
790 detector = camera[detNum]
791 # Determine the row axis based on the orientation
792 nq = detector.getOrientation().getNQuarter()
793 axisMap[detName] = 0 if nq % 2 == 0 else 1
794 return axisMap
796 def applyGuiderNameSwapIfNeeded(self, dayObs: int) -> None:
797 """
798 Apply guider name swap (SG0/SG1) for early data if required.
800 Parameters
801 ----------
802 dayObs : `int`
803 Observation day (YYYYMMDD) used to decide whether to swap.
804 """
805 if getattr(self, "_guiderNameMapSwapped", False):
806 return # Already swapped; do nothing
807 if dayObs < 20250509:
808 newMap = {}
809 for detName, detNum in self.guiderNameMap.items():
810 if detName.endswith("SG0"):
811 swapped = detName.replace("SG0", "SG1")
812 elif detName.endswith("SG1"):
813 swapped = detName.replace("SG1", "SG0")
814 else:
815 swapped = detName
816 swappedDetNum = self.camera[swapped].getId()
817 newMap[swapped] = swappedDetNum
818 self.guiderNameMap = newMap
819 self._guiderNameMapSwapped = True
822def getVisitInfo(butler: Butler, dayObs: int, seqNum: int, scienceDetNum: int) -> VisitInfo:
823 """
824 Retrieve VisitInfo for a given dayObs / seqNum.
826 Parameters
827 ----------
828 butler : `Butler`
829 Active Butler instance.
830 dayObs : `int`
831 Observation day (YYYYMMDD).
832 seqNum : `int`
833 Sequence number.
834 scienceDetNum : `int`
835 Science detector number.
837 Returns
838 -------
839 visitInfo : `VisitInfo`
840 VisitInfo object retrieved from butler.
841 """
842 dataId = {
843 "instrument": "LSSTCam",
844 "day_obs": dayObs,
845 "seq_num": seqNum,
846 "detector": scienceDetNum,
847 }
848 visitInfo = butler.get("raw.visitInfo", dataId)
849 return visitInfo
852# Helper functions for stamp conversion
853def _makeRoiTransforms(metadata: dict, detector: Detector, camera: Camera) -> tuple[tuple, str]:
854 """
855 Construct ROI transforms and derive amplifier name.
857 Parameters
858 ----------
859 metadata : `dict`
860 Exposure-level metadata.
861 detector : `Detector`
862 The detector.
863 camera : `Camera`
864 The camera object.
866 Returns
867 -------
868 result : `tuple[tuple, str]`
869 ((ccdViewBbox, fwd, back), ampName)
870 """
871 ampName = "C" + metadata["ROISEG"][7:]
872 ccdViewBbox = makeRoiBbox(metadata, camera)
873 fwd, back = makeCcdToDvcsTransform(ccdViewBbox, detector.getOrientation().getNQuarter())
874 return (ccdViewBbox, fwd, back), ampName
877def _convertMaskedImage(
878 maskedImage: MaskedImageF,
879 stampMetadata: dict[str, Any],
880 metadata: dict[str, Any],
881 transforms: tuple,
882 detector: Detector,
883 ampName: str,
884 camera: Camera,
885 view: str,
886 mask: None | np.ndarray = None,
887) -> Stamp:
888 """
889 Convert one masked image ROI and build a Stamp.
891 Parameters
892 ----------
893 maskedImage : `MaskedImageF`
894 Input masked image for the ROI.
895 stampMetadata : `dict`
896 Stamp-level metadata (modified in-place).
897 metadata : `dict`
898 Exposure-level metadata.
899 transforms : `tuple`
900 Tuple containing (ccdViewBbox, forwardTransform, inverseTransform).
901 detector : `Detector`
902 Camera detector.
903 ampName : `str`
904 Amplifier name.
905 camera : `Camera`
906 Camera object.
907 view : `str`
908 Output view ('dvcs', 'ccd', or 'roi').
909 mask: `ndarray or None`
910 Masked bad columns.
912 Returns
913 -------
914 stamp : `Stamp`
915 Converted stamp object.
916 """
917 ccdViewBbox, fwd, back = transforms
918 rawArray = maskedImage.getImage().getArray()
919 if mask is not None:
920 median = np.nanmedian(rawArray[~mask].flatten())
921 rawArray[mask] = median
922 dvcsArray = roiImageToDvcs(rawArray, metadata, detector, ampName, camera, view=view)
923 outImg = MaskedImageF(dvcsArray)
924 archiveElement = [ccdViewBbox, fwd, back]
925 return Stamp(outImg, archiveElement, metadata=stampMetadata)
928def _blankStamp(
929 stampIdx: int,
930 metadata: dict[str, Any],
931 transforms: tuple,
932) -> Stamp:
933 """
934 Create a blank (zero) stamp for a missing index.
936 Parameters
937 ----------
938 stampIdx : `int`
939 Missing stamp index.
940 metadata : `dict`
941 Original metadata container.
942 transforms : `tuple`
943 (ccdViewBbox, forwardTransform, inverseTransform).
945 Returns
946 -------
947 stamp : `Stamp`
948 Placeholder blank stamp with NaN timestamp.
949 """
950 ccdViewBbox, fwd, back = transforms
951 nRows, nCols = int(metadata["ROIROWS"]), int(metadata["ROICOLS"])
952 blankArray = np.zeros((nRows, nCols), dtype=np.float32)
953 blankImg = ExposureF(MaskedImageF(ImageF(array=blankArray)))
954 missingMetadata = metadata.copy()
955 missingMetadata["DAQSTAMP"] = stampIdx
956 missingMetadata["STMPTMJD"] = np.nan
957 archiveElement = [ccdViewBbox, fwd, back]
958 return Stamp(blankImg, archiveElement, metadata=missingMetadata)
961def convertRawStampsToView(
962 rawStamps: Stamps,
963 detName: str,
964 nStamps: int,
965 view: str = "dvcs",
966) -> Stamps:
967 """
968 Convert guider stamps from raw ROI to a requested view ('dvcs', 'ccd', or
969 'roi'). Handles missing stamps and preserves metadata and order.
971 Parameters
972 ----------
973 rawStamps : `Stamps`
974 Input raw ROI stamps.
975 detName : `str`
976 Detector name.
977 nStamps : `int`
978 Target total number of stamps (after filling gaps).
979 view : `str`, optional
980 Output view: 'dvcs', 'ccd', or 'roi'.
982 Returns
983 -------
984 stampsOut : `Stamps`
985 Converted stamps with gaps filled by blank stamps.
986 """
987 camera = LsstCam.getCamera()
988 detector = camera[detName]
989 metadata = rawStamps.metadata
990 metadataDict = metadata.toDict()
992 # Ensure CCDSLOT matches swapped name if applicable
993 metadata["CCDSLOT"] = detName[4:7]
995 # Align timestamps and find valid/missing indices
996 timestamps = standardizeGuiderTimestamps(rawStamps, nStamps)
997 validIndices = np.where(~timestamps.mask)[0].tolist()
998 missingIndices = np.where(timestamps.mask)[0].tolist()
1000 # Pre‑compute transforms once
1001 transforms, ampName = _makeRoiTransforms(metadataDict, detector, camera)
1003 stampsDict: dict[int, Stamp] = {}
1004 mIdx = 0 # index into masked images list
1006 # get bad columns
1007 img0 = rawStamps.getMaskedImages()[mIdx].getImage().getArray()
1008 mask = getColumnMask(img0, k=6)
1010 for idx in validIndices:
1011 maskedImage = rawStamps.getMaskedImages()[mIdx]
1012 stampMeta = rawStamps[mIdx].metadata
1013 stampMeta["DAQSTAMP"] = stampMeta.get("DAQSTAMP", idx)
1014 stampsDict[idx] = _convertMaskedImage(
1015 maskedImage, stampMeta, metadataDict, transforms, detector, ampName, camera, view, mask
1016 )
1017 mIdx += 1
1019 # Fill gaps with blanks
1020 for idx in missingIndices:
1021 stampsDict[idx] = _blankStamp(idx, metadataDict, transforms)
1023 # Assemble in order
1024 stampList = [stampsDict[i] for i in range(nStamps)]
1025 return Stamps(stampList, metadata, use_archive=True)
1028def standardizeGuiderTimestamps(rawStamps: Stamps, nStamps: int) -> Time:
1029 """
1030 Return a masked `Time` array of length `nStamps` for one guider.
1032 Missing stamps are filled with `NaN` and masked so every detector shares a
1033 uniform, sortable timestamp vector. Result is in MJD, `scale="tai"`.
1035 Parameters
1036 ----------
1037 rawStamps : `Stamps`
1038 Input stamps (may have missing indices).
1039 nStamps : `int`
1040 Desired total length (maximum among guiders).
1042 Returns
1043 -------
1044 fullTimestamps : `Time`
1045 Masked Time array (scale='tai') with inferred cadence and gaps masked.
1046 """
1047 timestampsList = [stamp.metadata.get("STMPTMJD", np.nan) for stamp in rawStamps]
1048 mjdArray = np.ma.masked_invalid(timestampsList)
1049 timestamps = Time(mjdArray, format="mjd", scale="tai")
1051 # infer frequency from valid timestamps
1052 jdDiff = np.diff(timestamps.jd)
1053 freqDays = np.nanmedian(jdDiff)
1054 startJd = timestamps[0].jd
1056 if np.isnan(freqDays) or freqDays <= 0:
1057 raise ValueError(
1058 f"Invalid frequency {freqDays} derived from timestamps. "
1059 "Ensure that the timestamps are valid and evenly spaced."
1060 )
1062 timestampsIdeal = Time(startJd + np.arange(nStamps) * freqDays, format="jd", scale="tai")
1064 tolerance = 0.5 * freqDays
1065 actualJd = timestamps.jd
1067 # Build aligned timestamp list with np.nan where missing
1068 mjdList = []
1069 for idealJd in timestampsIdeal.jd:
1070 # check if any actual timestamp is close to this ideal timestamp
1071 if np.any(np.abs(actualJd - idealJd) < tolerance):
1072 # match found — find the closest
1073 idx = np.argmin(np.abs(actualJd - idealJd))
1074 mjdList.append(actualJd[idx])
1075 else:
1076 # no match — this is a missing timestamp
1077 mjdList.append(np.nan)
1079 mjdArray = np.ma.masked_invalid(mjdList)
1080 fullTimestamps = Time(mjdArray, format="mjd", scale="tai")
1081 fullTimestamps.freq = freqDays
1082 return fullTimestamps
1085@dataclass
1086class WeatherInfo:
1087 """
1088 Container for weather information.
1090 Parameters
1091 ----------
1092 temperature : `float`
1093 Temperature in degrees Celsius.
1094 humidity : `float`
1095 Relative humidity in percent.
1096 pressure : `float`
1097 Atmospheric pressure in hPa.
1098 windSpeed : `float`
1099 Wind speed in m/s.
1100 windDir : `float`
1101 Wind direction in degrees from North.
1102 """
1104 temperature: float
1105 humidity: float
1106 pressure: float
1107 windSpeed: float
1108 windDir: float
1109 seeing: float
1111 @staticmethod
1112 def fromMetadata(metadata: dict) -> WeatherInfo:
1113 """
1114 Create a WeatherInfo object from metadata dictionary.
1116 Parameters
1117 ----------
1118 metadata : `dict`
1119 Metadata dictionary containing weather keys.
1121 Returns
1122 -------
1123 weather : `WeatherInfo`
1124 Parsed weather information.
1125 """
1126 return WeatherInfo(
1127 temperature=metadata_to_float(metadata, "AIRTEMP", np.nan),
1128 humidity=metadata_to_float(metadata, "HUMIDITY", np.nan),
1129 pressure=metadata_to_float(metadata, "PRESSURE", np.nan),
1130 windSpeed=metadata_to_float(metadata, "WINDSPD", np.nan),
1131 windDir=metadata_to_float(metadata, "WINDDIR", np.nan),
1132 seeing=metadata_to_float(metadata, "SEEING", np.nan),
1133 )
1136def metadata_to_float(metadata: dict, key: str, default: float = np.nan) -> float:
1137 """
1138 Safely convert a metadata value to float.
1140 Parameters
1141 ----------
1142 metadata : `dict`
1143 Metadata dictionary.
1144 key : `str`
1145 Key to look up in the metadata.
1146 default : `float`, optional
1147 Default value if key is missing or conversion fails.
1149 Returns
1150 -------
1151 value : `float`
1152 Converted float value or default.
1153 """
1154 value = metadata[key]
1155 if value is None:
1156 return default
1157 else:
1158 return float(value)
1161def mad(x: np.ndarray) -> float:
1162 """
1163 Median absolute deviation metric for std.
1165 Parameters
1166 ----------
1167 x : `ndarray`
1168 Input array.
1170 Returns
1171 -------
1172 value : `float`
1173 Robust standard deviation
1174 """
1175 med = np.nanmedian(x)
1176 return 1.4826 * np.nanmedian(np.abs(x - med))
1179def getColumnMask(img: NDArray[np.floating], k: int = 6) -> NDArray[np.bool_]:
1180 """Return a boolean mask for bad columns based on robust median statistics.
1182 Parameters
1183 ----------
1184 img : `ndarray`
1185 2D image array.
1186 k : `int`, optional
1187 Threshold factor for MAD.
1189 Returns
1190 -------
1191 mask : `ndarray`
1192 Boolean mask of same shape as img, ``True`` for bad columns.
1193 """
1194 columnMedian = np.nanmedian(img, axis=0)
1195 center = np.nanmedian(columnMedian)
1197 # MAD of column medians
1198 diffs = np.abs(columnMedian - center)
1199 s = np.nanmedian(diffs)
1201 if s == 0 or not np.isfinite(s):
1202 return np.zeros(img.shape, dtype=bool)
1204 badCols = diffs > k * s
1205 return np.broadcast_to(badCols, img.shape).copy()