Coverage for python / lsst / summit / utils / guiders / reading.py: 23%
373 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:17 +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 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 @cached_property
493 def alt(self) -> float:
494 """Return altitude (el_start) from the guider data header."""
495 raw = self.header.get("el_start")
496 if raw is None:
497 raise KeyError("Header key 'el_start' is missing or None")
498 return float(raw)
500 @cached_property
501 def az(self) -> float:
502 """Return azimuth (az_start) from the guider data header."""
503 raw = self.header.get("az_start")
504 if raw is None:
505 raise KeyError("Header key 'az_start' is missing or None")
506 return float(raw)
508 # Iterable / dict-like helpers
509 def __iter__(self):
510 """Iterate over detector names in guiderNames order."""
511 return iter(self.guiderNames)
513 def items(self):
514 """Yield (detName, stamps) pairs like dict.items()."""
515 for det in self.guiderNames:
516 yield det, self.stampsMap[det]
518 def keys(self):
519 """Iterate over detector names (dict-like .keys())."""
520 return iter(self.guiderNames)
522 def values(self):
523 """Iterate over Stamps objects in guiderNames order."""
524 for det in self.guiderNames:
525 yield self.stampsMap[det]
527 def __getitem__(self, key):
528 """
529 Direct stamp access helper.
531 Example usage:
532 gd["R44_SG0"] -> full `Stamps` object
533 gd["R44_SG0", 3] -> stamp #3 ndarray
534 gd[3] -> stamp #3 from `detNameMax`
536 Parameters
537 ----------
538 key : `str | int | slice | tuple`
539 Access pattern:
540 - 'DET' -> Stamps
541 - ('DET', i) -> single stamp ndarray
542 - i (int) -> stamp i from primary detector
543 - ('DET', slice) -> list of ndarrays
545 Returns
546 -------
547 result : `Stamps | np.ndarray | list[np.ndarray]`
548 Retrieved object according to key specification.
549 """
550 # Single detector name -> full Stamps
551 if isinstance(key, str):
552 return self.stampsMap[key]
554 # idx only -> assume detNameMax
555 if isinstance(key, (int, slice)):
556 key = (self.detNameMax, key)
558 if isinstance(key, tuple):
559 if len(key) != 2:
560 raise TypeError("Key must be (detName, idx) or (idx,)")
561 detName, idx = key[0], key[1]
563 if detName not in self.stampsMap:
564 raise KeyError(f"{detName!r} not found in stampsMap")
566 stamps = self.stampsMap[detName]
567 # slice returns list of ndarrays
568 if isinstance(idx, slice):
569 arrays = [self._processStampArray(stamp, detName) for stamp in stamps[idx]]
570 return arrays
571 # int -> single ndarray
572 return self._processStampArray(stamps[idx], detName)
574 raise TypeError("Invalid key type for GuiderData indexing.")
576 def _processStampArray(self, stamp: Stamp, detName: str) -> np.ndarray:
577 """
578 Convert a Stamp to an ndarray with optional median-row subtraction.
580 Parameters
581 ----------
582 stamp : `Stamp`
583 Input stamp object.
584 detName : `str`
585 Detector name (for row-axis logic).
587 Returns
588 -------
589 array : `np.ndarray`
590 2D image array (bias-corrected if configured).
591 """
592 arr = stamp.stamp_im.image.array
593 return arr
595 @cached_property
596 def plotter(self):
597 return GuiderPlotter(self)
599 def plotStamp(
600 self, detName: str, stampNum: int, plo: float = 90, phi: float = 99.5, figsize=(10, 8)
601 ) -> plt.Figure:
602 """
603 Plot a single guider stamp.
605 Parameters
606 ----------
607 detName : `str`
608 Detector name.
609 stampNum : `int`
610 Stamp index.
611 plo : `float`, optional
612 Lower percentile stretch.
613 phi : `float`, optional
614 Upper percentile stretch.
615 figsize : `tuple`, optional
616 Figure size.
618 Returns
619 -------
620 fig : `plt.Figure`
621 The resulting figure.
622 """
623 fig, axs = make_subplot(nrows=1, ncols=1, figsize=figsize)
624 img = self[detName, stampNum] if stampNum > 0 else self.getStampArrayCoadd(detName)
625 _ = axs.imshow(
626 img,
627 origin="lower",
628 vmin=np.nanpercentile(img, plo),
629 vmax=np.nanpercentile(img, phi),
630 cmap="Greys",
631 )
632 _ = labelDetector(axs, detName)
633 axs.set_xlabel("X (pixels)", fontsize=11)
634 axs.set_ylabel("Y (pixels)", fontsize=11)
635 axs.set_title(f"{self.expid}")
636 return fig
639class GuiderReader:
640 """
641 Utility to fetch LSST guider data via Butler.
643 Example:
644 from lsst.summit.utils.guiders.reading import GuiderReader
645 from lsst.daf.butler import Butler
646 butler = Butler("embargo", collections="LSSTCam/raw/guider")
648 seqNum, dayObs = 461, 20250425
649 reader = GuiderReader(butler, view="dvcs")
650 guiderData = reader.get(dayObs=dayObs, seqNum=seqNum)
651 """
653 def __init__(self, butler: Butler, view: str = "dvcs"):
654 self.butler = butler
655 self.view = view
656 self.log = logging.getLogger(__name__)
657 # Define camera objects
658 self.camera = LsstCam.getCamera()
660 # Build guiderNameMap
661 self.guiderNameMap: dict[str, int] = {}
662 for detector in self.camera:
663 if detector.getType() == cameraGeom.DetectorType.GUIDER:
664 detName = detector.getName()
665 self.guiderNameMap[detName] = detector.getId()
667 self.guiderDetNames = list(self.guiderNameMap.keys())
668 self.nGuiders = len(self.guiderNameMap)
670 def get(
671 self, dayObs: int, seqNum: int, doSubtractMedian: bool = True, scienceDetNum: int = 94
672 ) -> GuiderData:
673 """
674 Retrieve guider data for a given dayObs / seqNum.
676 Parameters
677 ----------
678 dayObs : `int`
679 Day of observation in YYYYMMDD format.
680 seqNum : `int`
681 Sequence number.
682 doSubtractMedian : `bool`, optional
683 If True, subtract median row bias from each stamp upon access.
684 scienceDetNum : `int`, optional
685 Science detector number for WCS reference.
687 Returns
688 -------
689 guiderData : `GuiderData`
690 Assembled guider data object.
691 """
692 # check if the guider name is swaped (dayObs < 20250509)
693 # modifies self.guiderNameMap in place if necessary
694 self.applyGuiderNameSwapIfNeeded(dayObs)
696 # 1. Get guider_raw stamps all guiders
697 rawStampsDict = self.getGuiderRawStamps(dayObs, seqNum)
699 # set the number of stamps
700 # determine the maximum number of stamps among all guiders
701 nStampsList = [len(stamps) for stamps in rawStampsDict.values()]
702 self.nStamps = max(nStampsList)
704 if self.nStamps <= 1:
705 raise RuntimeError(
706 f"Only {self.nStamps} stamps found for dayObs {dayObs}, seqNum {seqNum}. "
707 "At least 2 stamps are required to create GuiderData."
708 )
710 # Create a visitinfo dataId to get the WCS and camera rotation angle
711 visitInfo = getVisitInfo(self.butler, dayObs, seqNum, scienceDetNum)
713 # 2. Get WCS map for all guiders
714 wcsMapDict = makeInitGuiderWcs(self.camera, visitInfo)
716 # 3. Get camera rotation angle
717 camRotAngle = getCamRotAngle(visitInfo)
719 # 4. subtract median row bias if requested
720 if doSubtractMedian:
721 rawStampsDict = self.subtractMedianRowBias(rawStampsDict)
723 # 4. Pack everything into a GuiderData object
724 guiderData = GuiderData(
725 seqNum=seqNum,
726 dayObs=dayObs,
727 view=self.view,
728 rawStampsMap=rawStampsDict,
729 guiderNameMap=self.guiderNameMap,
730 wcsMap=wcsMapDict,
731 camRotAngle=camRotAngle,
732 isMedianSubtracted=doSubtractMedian,
733 )
734 return guiderData
736 def getGuiderRawStamps(self, dayObs: int, seqNum: int) -> dict[str, Stamps]:
737 """
738 Fetch raw guider Stamps for all detectors.
740 Parameters
741 ----------
742 dayObs : `int`
743 Observation day (YYYYMMDD).
744 seqNum : `int`
745 Sequence number.
747 Returns
748 -------
749 rawStamps : `dict[str, Stamps]`
750 Mapping from detector name to raw Stamps object.
751 """
752 rawStamps: dict[str, Stamps] = {}
753 for detName, detNum in self.guiderNameMap.items():
754 try:
755 rawStamps[detName] = self.butler.get(
756 "guider_raw",
757 day_obs=dayObs,
758 seq_num=seqNum,
759 detector=detNum,
760 instrument="LSSTCam",
761 )
762 except DatasetNotFoundError as e:
763 raise DatasetNotFoundError(f"No data for {detName} on {dayObs=} {seqNum=}") from e
764 return rawStamps
766 def subtractMedianRowBias(self, rawStampsDict: dict[str, Stamps]) -> dict[str, Stamps]:
767 """
768 Subtract median row bias from all raw stamps in-place.
770 This modifies the rawStampsMap of the GuiderReader instance.
772 Parameters
773 ----------
774 rawStampsDict : `dict[str, Stamps]`
775 ROI view of stamps that come from butler `guider_raw`.
776 """
777 for detName, stamps in rawStampsDict.items():
778 # Subtract median along the row axis
779 rowAxis = 0
780 for i in range(len(stamps)):
781 # Get the stamp image array
782 data = stamps[i].stamp_im.image.array
783 # Subtract median along the specified row axis
784 medianRow = np.nanmedian(data, axis=rowAxis, keepdims=True)
785 # Update the stamp with the new data
786 stamps[i].stamp_im.image.array = data - medianRow
787 # Update the rawStampsMap in-place
788 rawStampsDict[detName] = stamps
789 return rawStampsDict
791 def getAxisRowMap(self) -> dict[str, int]:
792 """
793 Axis mapping (row axis) for guider detectors.
795 Returns
796 -------
797 axisMap : `dict[str, int]`
798 Dictionary with detector names as keys and axis mapping as values.
799 The axis mapping is 1 for X and 0 for Y.
800 """
801 camera = LsstCam.getCamera()
802 axisMap: dict[str, int] = {}
803 for detName, detNum in self.guiderNameMap.items():
804 # Get the detector object from the camera
805 detector = camera[detNum]
806 # Determine the row axis based on the orientation
807 nq = detector.getOrientation().getNQuarter()
808 axisMap[detName] = 0 if nq % 2 == 0 else 1
809 return axisMap
811 def applyGuiderNameSwapIfNeeded(self, dayObs: int) -> None:
812 """
813 Apply guider name swap (SG0/SG1) for early data if required.
815 Parameters
816 ----------
817 dayObs : `int`
818 Observation day (YYYYMMDD) used to decide whether to swap.
819 """
820 if getattr(self, "_guiderNameMapSwapped", False):
821 return # Already swapped; do nothing
822 if dayObs < 20250509:
823 newMap = {}
824 for detName, detNum in self.guiderNameMap.items():
825 if detName.endswith("SG0"):
826 swapped = detName.replace("SG0", "SG1")
827 elif detName.endswith("SG1"):
828 swapped = detName.replace("SG1", "SG0")
829 else:
830 swapped = detName
831 swappedDetNum = self.camera[swapped].getId()
832 newMap[swapped] = swappedDetNum
833 self.guiderNameMap = newMap
834 self._guiderNameMapSwapped = True
837def getVisitInfo(butler: Butler, dayObs: int, seqNum: int, scienceDetNum: int) -> VisitInfo:
838 """
839 Retrieve VisitInfo for a given dayObs / seqNum.
841 Parameters
842 ----------
843 butler : `Butler`
844 Active Butler instance.
845 dayObs : `int`
846 Observation day (YYYYMMDD).
847 seqNum : `int`
848 Sequence number.
849 scienceDetNum : `int`
850 Science detector number.
852 Returns
853 -------
854 visitInfo : `VisitInfo`
855 VisitInfo object retrieved from butler.
856 """
857 dataId = {
858 "instrument": "LSSTCam",
859 "day_obs": dayObs,
860 "seq_num": seqNum,
861 "detector": scienceDetNum,
862 }
863 visitInfo = butler.get("raw.visitInfo", dataId)
864 return visitInfo
867# Helper functions for stamp conversion
868def _makeRoiTransforms(metadata: dict, detector: Detector, camera: Camera) -> tuple[tuple, str]:
869 """
870 Construct ROI transforms and derive amplifier name.
872 Parameters
873 ----------
874 metadata : `dict`
875 Exposure-level metadata.
876 detector : `Detector`
877 The detector.
878 camera : `Camera`
879 The camera object.
881 Returns
882 -------
883 result : `tuple[tuple, str]`
884 ((ccdViewBbox, fwd, back), ampName)
885 """
886 ampName = "C" + metadata["ROISEG"][7:]
887 ccdViewBbox = makeRoiBbox(metadata, camera)
888 fwd, back = makeCcdToDvcsTransform(ccdViewBbox, detector.getOrientation().getNQuarter())
889 return (ccdViewBbox, fwd, back), ampName
892def _convertMaskedImage(
893 maskedImage: MaskedImageF,
894 stampMetadata: dict[str, Any],
895 metadata: dict[str, Any],
896 transforms: tuple,
897 detector: Detector,
898 ampName: str,
899 camera: Camera,
900 view: str,
901 mask: None | np.ndarray = None,
902) -> Stamp:
903 """
904 Convert one masked image ROI and build a Stamp.
906 Parameters
907 ----------
908 maskedImage : `MaskedImageF`
909 Input masked image for the ROI.
910 stampMetadata : `dict`
911 Stamp-level metadata (modified in-place).
912 metadata : `dict`
913 Exposure-level metadata.
914 transforms : `tuple`
915 Tuple containing (ccdViewBbox, forwardTransform, inverseTransform).
916 detector : `Detector`
917 Camera detector.
918 ampName : `str`
919 Amplifier name.
920 camera : `Camera`
921 Camera object.
922 view : `str`
923 Output view ('dvcs', 'ccd', or 'roi').
924 mask: `ndarray or None`
925 Masked bad columns.
927 Returns
928 -------
929 stamp : `Stamp`
930 Converted stamp object.
931 """
932 ccdViewBbox, fwd, back = transforms
933 rawArray = maskedImage.getImage().getArray()
934 if mask is not None:
935 median = np.nanmedian(rawArray[~mask].flatten())
936 rawArray[mask] = median
937 dvcsArray = roiImageToDvcs(rawArray, metadata, detector, ampName, camera, view=view)
938 outImg = MaskedImageF(dvcsArray)
939 archiveElement = [ccdViewBbox, fwd, back]
940 return Stamp(outImg, archiveElement, metadata=stampMetadata)
943def _blankStamp(
944 stampIdx: int,
945 metadata: dict[str, Any],
946 transforms: tuple,
947) -> Stamp:
948 """
949 Create a blank (zero) stamp for a missing index.
951 Parameters
952 ----------
953 stampIdx : `int`
954 Missing stamp index.
955 metadata : `dict`
956 Original metadata container.
957 transforms : `tuple`
958 (ccdViewBbox, forwardTransform, inverseTransform).
960 Returns
961 -------
962 stamp : `Stamp`
963 Placeholder blank stamp with NaN timestamp.
964 """
965 ccdViewBbox, fwd, back = transforms
966 nRows, nCols = int(metadata["ROIROWS"]), int(metadata["ROICOLS"])
967 blankArray = np.zeros((nRows, nCols), dtype=np.float32)
968 blankImg = ExposureF(MaskedImageF(ImageF(array=blankArray)))
969 missingMetadata = metadata.copy()
970 missingMetadata["DAQSTAMP"] = stampIdx
971 missingMetadata["STMPTMJD"] = np.nan
972 archiveElement = [ccdViewBbox, fwd, back]
973 return Stamp(blankImg, archiveElement, metadata=missingMetadata)
976def convertRawStampsToView(
977 rawStamps: Stamps,
978 detName: str,
979 nStamps: int,
980 view: str = "dvcs",
981) -> Stamps:
982 """
983 Convert guider stamps from raw ROI to a requested view ('dvcs', 'ccd', or
984 'roi'). Handles missing stamps and preserves metadata and order.
986 Parameters
987 ----------
988 rawStamps : `Stamps`
989 Input raw ROI stamps.
990 detName : `str`
991 Detector name.
992 nStamps : `int`
993 Target total number of stamps (after filling gaps).
994 view : `str`, optional
995 Output view: 'dvcs', 'ccd', or 'roi'.
997 Returns
998 -------
999 stampsOut : `Stamps`
1000 Converted stamps with gaps filled by blank stamps.
1001 """
1002 camera = LsstCam.getCamera()
1003 detector = camera[detName]
1004 metadata = rawStamps.metadata
1005 metadataDict = metadata.toDict()
1007 # Ensure CCDSLOT matches swapped name if applicable
1008 metadata["CCDSLOT"] = detName[4:7]
1010 # Align timestamps and find valid/missing indices
1011 timestamps = standardizeGuiderTimestamps(rawStamps, nStamps)
1012 validIndices = np.where(~timestamps.mask)[0].tolist()
1013 missingIndices = np.where(timestamps.mask)[0].tolist()
1015 # Pre‑compute transforms once
1016 transforms, ampName = _makeRoiTransforms(metadataDict, detector, camera)
1018 stampsDict: dict[int, Stamp] = {}
1019 mIdx = 0 # index into masked images list
1021 # get bad columns
1022 img0 = rawStamps.getMaskedImages()[mIdx].getImage().getArray()
1023 mask = getColumnMask(img0, k=6)
1025 for idx in validIndices:
1026 maskedImage = rawStamps.getMaskedImages()[mIdx]
1027 stampMeta = rawStamps[mIdx].metadata
1028 stampMeta["DAQSTAMP"] = stampMeta.get("DAQSTAMP", idx)
1029 stampsDict[idx] = _convertMaskedImage(
1030 maskedImage, stampMeta, metadataDict, transforms, detector, ampName, camera, view, mask
1031 )
1032 mIdx += 1
1034 # Fill gaps with blanks
1035 for idx in missingIndices:
1036 stampsDict[idx] = _blankStamp(idx, metadataDict, transforms)
1038 # Assemble in order
1039 stampList = [stampsDict[i] for i in range(nStamps)]
1040 return Stamps(stampList, metadata, use_archive=True)
1043def standardizeGuiderTimestamps(rawStamps: Stamps, nStamps: int) -> Time:
1044 """
1045 Return a masked `Time` array of length `nStamps` for one guider.
1047 Missing stamps are filled with `NaN` and masked so every detector shares a
1048 uniform, sortable timestamp vector. Result is in MJD, `scale="tai"`.
1050 Parameters
1051 ----------
1052 rawStamps : `Stamps`
1053 Input stamps (may have missing indices).
1054 nStamps : `int`
1055 Desired total length (maximum among guiders).
1057 Returns
1058 -------
1059 fullTimestamps : `Time`
1060 Masked Time array (scale='tai') with inferred cadence and gaps masked.
1061 """
1062 timestampsList = [stamp.metadata.get("STMPTMJD", np.nan) for stamp in rawStamps]
1063 mjdArray = np.ma.masked_invalid(timestampsList)
1064 timestamps = Time(mjdArray, format="mjd", scale="tai")
1066 # infer frequency from valid timestamps
1067 jdDiff = np.diff(timestamps.jd)
1068 freqDays = np.nanmedian(jdDiff)
1069 startJd = timestamps[0].jd
1071 if np.isnan(freqDays) or freqDays <= 0:
1072 raise ValueError(
1073 f"Invalid frequency {freqDays} derived from timestamps. "
1074 "Ensure that the timestamps are valid and evenly spaced."
1075 )
1077 timestampsIdeal = Time(startJd + np.arange(nStamps) * freqDays, format="jd", scale="tai")
1079 tolerance = 0.5 * freqDays
1080 actualJd = timestamps.jd
1082 # Build aligned timestamp list with np.nan where missing
1083 mjdList = []
1084 for idealJd in timestampsIdeal.jd:
1085 # check if any actual timestamp is close to this ideal timestamp
1086 if np.any(np.abs(actualJd - idealJd) < tolerance):
1087 # match found — find the closest
1088 idx = np.argmin(np.abs(actualJd - idealJd))
1089 mjdList.append(actualJd[idx])
1090 else:
1091 # no match — this is a missing timestamp
1092 mjdList.append(np.nan)
1094 mjdArray = np.ma.masked_invalid(mjdList)
1095 fullTimestamps = Time(mjdArray, format="mjd", scale="tai")
1096 fullTimestamps.freq = freqDays
1097 return fullTimestamps
1100@dataclass
1101class WeatherInfo:
1102 """
1103 Container for weather information.
1105 Parameters
1106 ----------
1107 temperature : `float`
1108 Temperature in degrees Celsius.
1109 humidity : `float`
1110 Relative humidity in percent.
1111 pressure : `float`
1112 Atmospheric pressure in hPa.
1113 windSpeed : `float`
1114 Wind speed in m/s.
1115 windDir : `float`
1116 Wind direction in degrees from North.
1117 """
1119 temperature: float
1120 humidity: float
1121 pressure: float
1122 windSpeed: float
1123 windDir: float
1124 seeing: float
1126 @staticmethod
1127 def fromMetadata(metadata: dict) -> WeatherInfo:
1128 """
1129 Create a WeatherInfo object from metadata dictionary.
1131 Parameters
1132 ----------
1133 metadata : `dict`
1134 Metadata dictionary containing weather keys.
1136 Returns
1137 -------
1138 weather : `WeatherInfo`
1139 Parsed weather information.
1140 """
1141 return WeatherInfo(
1142 temperature=metadata_to_float(metadata, "AIRTEMP", np.nan),
1143 humidity=metadata_to_float(metadata, "HUMIDITY", np.nan),
1144 pressure=metadata_to_float(metadata, "PRESSURE", np.nan),
1145 windSpeed=metadata_to_float(metadata, "WINDSPD", np.nan),
1146 windDir=metadata_to_float(metadata, "WINDDIR", np.nan),
1147 seeing=metadata_to_float(metadata, "SEEING", np.nan),
1148 )
1151def metadata_to_float(metadata: dict, key: str, default: float = np.nan) -> float:
1152 """
1153 Safely convert a metadata value to float.
1155 Parameters
1156 ----------
1157 metadata : `dict`
1158 Metadata dictionary.
1159 key : `str`
1160 Key to look up in the metadata.
1161 default : `float`, optional
1162 Default value if key is missing or conversion fails.
1164 Returns
1165 -------
1166 value : `float`
1167 Converted float value or default.
1168 """
1169 value = metadata[key]
1170 if value is None:
1171 return default
1172 else:
1173 return float(value)
1176def mad(x: np.ndarray) -> float:
1177 """
1178 Median absolute deviation metric for std.
1180 Parameters
1181 ----------
1182 x : `ndarray`
1183 Input array.
1185 Returns
1186 -------
1187 value : `float`
1188 Robust standard deviation
1189 """
1190 med = np.nanmedian(x)
1191 return 1.4826 * np.nanmedian(np.abs(x - med))
1194def getColumnMask(img: NDArray[np.floating], k: int = 6) -> NDArray[np.bool_]:
1195 """Return a boolean mask for bad columns based on robust median statistics.
1197 Parameters
1198 ----------
1199 img : `ndarray`
1200 2D image array.
1201 k : `int`, optional
1202 Threshold factor for MAD.
1204 Returns
1205 -------
1206 mask : `ndarray`
1207 Boolean mask of same shape as img, ``True`` for bad columns.
1208 """
1209 columnMedian = np.nanmedian(img, axis=0)
1210 center = np.nanmedian(columnMedian)
1212 # MAD of column medians
1213 diffs = np.abs(columnMedian - center)
1214 s = np.nanmedian(diffs)
1216 if s == 0 or not np.isfinite(s):
1217 return np.zeros(img.shape, dtype=bool)
1219 badCols = diffs > k * s
1220 return np.broadcast_to(badCols, img.shape).copy()