Coverage for python / lsst / meas / algorithms / brightStarStamps.py: 26%
165 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-14 23:55 +0000
1# This file is part of meas_algorithms.
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/>.
22"""Collection of small images (postage stamps) centered on bright stars."""
24from __future__ import annotations
26__all__ = ["BrightStarStamp", "BrightStarStamps"]
28from collections.abc import Sequence
29from dataclasses import dataclass
31import numpy as np
33from lsst.afw.detection import Psf
34from lsst.afw.fits import Fits, readMetadata
35from lsst.afw.geom import SkyWcs
36from lsst.afw.image import ImageFitsReader, MaskedImageF, MaskFitsReader
37from lsst.afw.table.io import InputArchive, OutputArchive
38from lsst.daf.base import PropertyList
39from lsst.geom import Angle, Point2D, degrees
40from lsst.meas.algorithms.stamps import AbstractStamp
41from lsst.utils.introspection import get_full_type_name
44@dataclass
45class BrightStarStamp(AbstractStamp):
46 """A single postage stamp centered on a bright star.
48 Attributes
49 ----------
50 stamp_im : `~lsst.afw.image.MaskedImageF`
51 The pixel data for this stamp.
52 psf : `~lsst.afw.detection.Psf`, optional
53 The point-spread function for this star.
54 wcs : `~lsst.afw.geom.SkyWcs`, optional
55 World coordinate system associated with the stamp.
56 visit : `int`, optional
57 Visit number of the observation.
58 detector : `int`, optional
59 Detector ID within the visit.
60 ref_id : `int`, optional
61 Reference catalog ID of the star.
62 ref_mag : `float`, optional
63 Reference catalog magnitude of the star.
64 position : `~lsst.geom.Point2D`, optional
65 Center position of the star on the detector in pixel coordinates.
66 focal_plane_radius : `float`, optional
67 Radial distance from the focal plane center in tangent-plane pixels.
68 focal_plane_angle : `~lsst.geom.Angle`, optional
69 Azimuthal angle on the focal plane (counterclockwise from +X).
70 scale : `float`, optional
71 Flux scaling factor applied to the PSF model.
72 scale_err : `float`, optional
73 Error in the flux scale.
74 pedestal : `float`, optional
75 Background pedestal level.
76 pedestal_err : `float`, optional
77 Error on the pedestal.
78 pedestal_scale_cov : `float`, optional
79 Covariance between pedestal and scale.
80 gradient_x : `float`, optional
81 Background gradient in the X direction.
82 gradient_y : `float`, optional
83 Background gradient in the Y direction.
84 global_reduced_chi_squared : `float`, optional
85 Reduced chi-squared for the global model fit.
86 global_degrees_of_freedom : `int`, optional
87 Degrees of freedom for the global model fit.
88 psf_reduced_chi_squared : `float`, optional
89 Reduced chi-squared for the PSF fit.
90 psf_degrees_of_freedom : `int`, optional
91 Degrees of freedom for the PSF fit.
92 psf_masked_flux_fraction : `float`, optional
93 Fraction of flux masked in the PSF.
95 Notes
96 -----
97 This class is designed to be used with `BrightStarStamps`, which manages
98 collections of these stamps and handles reading/writing them to FITS files.
99 The `factory` class method provides a standard interface to construct
100 instances from image data and metadata, while the `_getMetadata` method
101 extracts metadata for storage in FITS headers.
102 """
104 stamp_im: MaskedImageF
105 psf: Psf | None
106 wcs: SkyWcs | None
107 visit: int | None
108 detector: int | None
109 ref_id: int | None
110 ref_mag: float | None
111 position: Point2D | None
112 focal_plane_radius: float | None
113 focal_plane_angle: Angle | None
114 scale: float | None
115 scale_err: float | None
116 pedestal: float | None
117 pedestal_err: float | None
118 pedestal_scale_cov: float | None
119 gradient_x: float | None
120 gradient_y: float | None
121 global_reduced_chi_squared: float | None
122 global_degrees_of_freedom: int | None
123 psf_reduced_chi_squared: float | None
124 psf_degrees_of_freedom: int | None
125 psf_masked_flux_fraction: float | None
127 # Mapping of metadata keys to attribute names
128 _metadata_attribute_map = {
129 "VISIT": "visit",
130 "DETECTOR": "detector",
131 "REF_ID": "ref_id",
132 "REF_MAG": "ref_mag",
133 "POSITION_X": "position.x",
134 "POSITION_Y": "position.y",
135 "FOCAL_PLANE_RADIUS": "focal_plane_radius",
136 "FOCAL_PLANE_ANGLE_DEGREES": "focal_plane_angle",
137 "SCALE": "scale",
138 "SCALE_ERR": "scale_err",
139 "PEDESTAL": "pedestal",
140 "PEDESTAL_ERR": "pedestal_err",
141 "PEDESTAL_SCALE_COV": "pedestal_scale_cov",
142 "GRADIENT_X": "gradient_x",
143 "GRADIENT_Y": "gradient_y",
144 "GLOBAL_REDUCED_CHI_SQUARED": "global_reduced_chi_squared",
145 "GLOBAL_DEGREES_OF_FREEDOM": "global_degrees_of_freedom",
146 "PSF_REDUCED_CHI_SQUARED": "psf_reduced_chi_squared",
147 "PSF_DEGREES_OF_FREEDOM": "psf_degrees_of_freedom",
148 "PSF_MASKED_FLUX_FRACTION": "psf_masked_flux_fraction",
149 }
151 def _getMetadata(self) -> PropertyList:
152 """Extract metadata from the stamp's attributes.
154 This method constructs a `PropertyList` containing metadata
155 extracted from the stamp's attributes. It is used when writing the
156 stamp to a FITS file to store relevant metadata in the FITS headers.
158 Returns
159 -------
160 metadata : `PropertyList`
161 A `PropertyList` containing the metadata, or `None` if no
162 metadata attributes are defined.
163 """
164 metadata = PropertyList()
165 for metadata_key, attribute_name in self._metadata_attribute_map.items():
166 if "." in attribute_name:
167 top_attr, sub_attr = attribute_name.split(".")
168 value = getattr(getattr(self, top_attr), sub_attr)
169 elif metadata_key == "FOCAL_PLANE_ANGLE_DEGREES":
170 value = getattr(self, attribute_name).asDegrees()
171 else:
172 value = getattr(self, attribute_name)
173 metadata[metadata_key] = value
174 return metadata
176 @property
177 def metadata(self) -> PropertyList:
178 """Return the stamp's metadata as a PropertyList."""
179 return self._getMetadata()
181 @classmethod
182 def factory(
183 cls,
184 stamp_im: MaskedImageF,
185 psf: Psf | None,
186 wcs: SkyWcs | None,
187 metadata: PropertyList,
188 ) -> BrightStarStamp:
189 """Construct a `BrightStarStamp` from image data and metadata.
191 This method provides a standard interface to create a `BrightStarStamp`
192 from its image data, PSF, WCS, and associated metadata.
193 It is used by the `BrightStarStamps.readFits` method to construct
194 individual bright star stamps from FITS files.
196 Parameters
197 ----------
198 stamp_im : `~lsst.afw.image.MaskedImageF`
199 Masked image for the stamp.
200 psf : `~lsst.afw.detection.Psf`, optional
201 Point-spread function for the stamp.
202 wcs : `~lsst.afw.geom.SkyWcs`, optional
203 World coordinate system for the stamp.
204 metadata : `PropertyList`
205 Metadata associated with the stamp, containing keys for all
206 required attributes.
208 Returns
209 -------
210 brightStarStamp : `BrightStarStamp`
211 The constructed `BrightStarStamp` instance.
212 """
213 kwargs = {}
215 for metadata_key, attribute_name in cls._metadata_attribute_map.items():
216 if "." in attribute_name: # for nested attributes like position.x
217 top_attr, sub_attr = attribute_name.split(".")
218 if top_attr not in kwargs: # avoid overwriting position
219 if top_attr == "position": # make an initial Point2D
220 kwargs[top_attr] = Point2D(0, 0)
221 setattr(kwargs[top_attr], sub_attr, metadata[metadata_key])
222 elif attribute_name == "focal_plane_angle":
223 kwargs[attribute_name] = Angle(metadata[metadata_key], degrees)
224 else:
225 kwargs[attribute_name] = metadata[metadata_key]
227 return cls(stamp_im=stamp_im, psf=psf, wcs=wcs, **kwargs)
230class BrightStarStamps(Sequence[BrightStarStamp]):
231 """A collection of bright star stamps.
233 Parameters
234 ----------
235 brightStarStamps : `Iterable` [`BrightStarStamp`]
236 Collection of `BrightStarStamp` instances.
237 metadata : `~lsst.daf.base.PropertyList`, optional
238 Global metadata associated with the collection.
239 """
241 def __init__(
242 self,
243 brightStarStamps: Sequence[BrightStarStamp],
244 metadata: PropertyList | None = None,
245 ):
246 self._stamps = list(brightStarStamps)
247 self._metadata = PropertyList() if metadata is None else metadata.deepCopy()
248 self.by_ref_id = {stamp.ref_id: stamp for stamp in self}
250 def __len__(self):
251 return len(self._stamps)
253 def __getitem__(self, index):
254 if isinstance(index, slice):
255 return BrightStarStamps(self._stamps[index], metadata=self._metadata)
256 return self._stamps[index]
258 def __iter__(self):
259 return iter(self._stamps)
261 @property
262 def metadata(self):
263 """Return the collection's global metadata as a PropertyList."""
264 return self._metadata
266 @classmethod
267 def readFits(cls, filename: str) -> BrightStarStamps:
268 """Make a `BrightStarStamps` object from a FITS file.
270 Parameters
271 ----------
272 filename : `str`
273 Name of the FITS file to read.
275 Returns
276 -------
277 brightStarStamps : `BrightStarStamps`
278 The constructed `BrightStarStamps` instance.
279 """
280 return cls.readFitsWithOptions(filename, None)
282 @classmethod
283 def readFitsWithOptions(cls, filename: str, options: PropertyList | None) -> BrightStarStamps:
284 """Make a `BrightStarStamps` object from a FITS file, with options.
286 Parameters
287 ----------
288 filename : `str`
289 Name of the FITS file to read.
290 options : `~lsst.daf.base.PropertyList`, optional
291 Options for reading the FITS file. Not currently used.
293 Returns
294 -------
295 brightStarStamps : `BrightStarStamps`
296 The constructed `BrightStarStamps` instance.
297 """
298 with Fits(filename, "r") as fits_file:
299 stamp_planes = {}
300 stamp_psf_ids = {}
301 stamp_wcs_ids = {}
302 stamp_metadata = {}
303 archive = None
305 for hdu_num in range(1, fits_file.countHdus()): # Skip primary HDU
306 metadata = readMetadata(filename, hdu=hdu_num)
307 extname = metadata["EXTNAME"]
308 stamp_id: int | None = metadata.get("EXTVER", None)
310 # Skip non-image BINTABLEs (except ARCHIVE_INDEX)
311 if metadata["XTENSION"] == "BINTABLE" and not metadata.get("ZIMAGE", False):
312 if extname != "ARCHIVE_INDEX":
313 continue
315 # Handle the archive index separately
316 if extname == "ARCHIVE_INDEX":
317 fits_file.setHdu(hdu_num)
318 archive = InputArchive.readFits(fits_file)
319 continue
320 elif metadata.get("EXTTYPE") == "ARCHIVE_DATA":
321 continue
323 # Select reader and dtype
324 if extname == "IMAGE":
325 reader = ImageFitsReader(filename, hdu=hdu_num)
326 dtype = np.dtype(MaskedImageF.dtype)
327 stamp_psf_ids[stamp_id] = metadata.pop("PSF", None)
328 stamp_wcs_ids[stamp_id] = metadata.pop("WCS", None)
329 stamp_metadata[stamp_id] = metadata
330 elif extname == "MASK":
331 reader = MaskFitsReader(filename, hdu=hdu_num)
332 dtype = None
333 elif extname == "VARIANCE":
334 reader = ImageFitsReader(filename, hdu=hdu_num)
335 dtype = np.dtype("float32")
336 else:
337 raise ValueError(f"Unknown extension type: {extname}")
339 if stamp_id is not None:
340 stamp_planes.setdefault(stamp_id, {})[extname.lower()] = reader.read(dtype=dtype)
342 primary_metadata = readMetadata(filename, hdu=0)
343 num_stamps = primary_metadata["N_STAMPS"]
345 if len(stamp_planes) != num_stamps:
346 raise ValueError(
347 f"Number of stamps read ({len(stamp_planes)}) does not agree with the "
348 f"number of stamps recorded in the primary HDU metadata ({num_stamps})."
349 )
350 if archive is None:
351 raise ValueError("No archive index was found in the FITS file; cannot read PSF or WCS.")
353 brightStarStamps = []
354 for stamp_id in range(1, num_stamps + 1): # Need to increment by one as EXTVER starts at 1
355 stamp = MaskedImageF(**stamp_planes[stamp_id])
356 psf = archive.get(stamp_psf_ids[stamp_id])
357 wcs = archive.get(stamp_wcs_ids[stamp_id])
358 brightStarStamps.append(BrightStarStamp.factory(stamp, psf, wcs, stamp_metadata[stamp_id]))
360 return cls(brightStarStamps, primary_metadata)
362 def writeFits(self, filename: str):
363 """Write this `BrightStarStamps` object to a FITS file.
365 Parameters
366 ----------
367 filename : `str`
368 Name of the FITS file to write.
369 """
370 metadata = self._metadata.deepCopy()
372 # Store metadata in the primary HDU
373 metadata["N_STAMPS"] = len(self._stamps)
374 metadata["VERSION"] = 2 # Record version number in case of future code changes
375 metadata["STAMPCLS"] = get_full_type_name(self)
377 # Create and write to the FITS file within a context manager
378 with Fits(filename, "w") as fits_file:
379 fits_file.createEmpty()
381 # Store Persistables in an OutputArchive
382 output_archive = OutputArchive()
383 stamp_psf_ids = []
384 stamp_wcs_ids = []
385 for stamp in self._stamps:
386 stamp_psf_ids.append(output_archive.put(stamp.psf))
387 stamp_wcs_ids.append(output_archive.put(stamp.wcs))
389 # Write to the FITS file
390 fits_file.writeMetadata(metadata)
391 del metadata
392 output_archive.writeFits(fits_file)
394 # Add all pixel data to extension HDUs; note: EXTVER should be 1-based
395 for stamp_id, (stamp, stamp_psf_id, stamp_wcs_id) in enumerate(
396 zip(self._stamps, stamp_psf_ids, stamp_wcs_ids),
397 start=1,
398 ):
399 metadata = PropertyList()
400 metadata.update({"EXTVER": stamp_id, "EXTNAME": "IMAGE"})
401 if stamp_metadata := stamp._getMetadata():
402 metadata.update(stamp_metadata)
403 metadata["PSF"] = stamp_psf_id
404 metadata["WCS"] = stamp_wcs_id
405 stamp.stamp_im.getImage().writeFits(filename, metadata=metadata, mode="a")
407 metadata = PropertyList()
408 metadata.update({"EXTVER": stamp_id, "EXTNAME": "MASK"})
409 stamp.stamp_im.getMask().writeFits(filename, metadata=metadata, mode="a")
411 metadata = PropertyList()
412 metadata.update({"EXTVER": stamp_id, "EXTNAME": "VARIANCE"})
413 stamp.stamp_im.getVariance().writeFits(filename, metadata=metadata, mode="a")