Coverage for python / lsst / meas / algorithms / brightStarStamps.py: 27%

168 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-04-17 09:07 +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/>. 

21 

22"""Collection of small images (postage stamps) centered on bright stars.""" 

23 

24from __future__ import annotations 

25 

26__all__ = ["BrightStarStamp", "BrightStarStamps"] 

27 

28from collections.abc import Sequence 

29from dataclasses import dataclass 

30 

31import numpy as np 

32 

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 

42 

43 

44@dataclass 

45class BrightStarStamp(AbstractStamp): 

46 """A single postage stamp centered on a bright star. 

47 

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 curvature_x: `float`, optional 

85 Background curvature in the X direction. 

86 curvature_y: `float`, optional 

87 Background curvature in the Y direction. 

88 curvature_xy: `float`, optional 

89 The xy component of the second order fit. 

90 global_reduced_chi_squared : `float`, optional 

91 Reduced chi-squared for the global model fit. 

92 global_degrees_of_freedom : `int`, optional 

93 Degrees of freedom for the global model fit. 

94 psf_reduced_chi_squared : `float`, optional 

95 Reduced chi-squared for the PSF fit. 

96 psf_degrees_of_freedom : `int`, optional 

97 Degrees of freedom for the PSF fit. 

98 psf_masked_flux_fraction : `float`, optional 

99 Fraction of flux masked in the PSF. 

100 

101 Notes 

102 ----- 

103 This class is designed to be used with `BrightStarStamps`, which manages 

104 collections of these stamps and handles reading/writing them to FITS files. 

105 The `factory` class method provides a standard interface to construct 

106 instances from image data and metadata, while the `_getMetadata` method 

107 extracts metadata for storage in FITS headers. 

108 """ 

109 

110 stamp_im: MaskedImageF 

111 psf: Psf | None 

112 wcs: SkyWcs | None 

113 visit: int | None 

114 detector: int | None 

115 ref_id: int | None 

116 ref_mag: float | None 

117 position: Point2D | None 

118 focal_plane_radius: float | None 

119 focal_plane_angle: Angle | None 

120 scale: float | None 

121 scale_err: float | None 

122 pedestal: float | None 

123 pedestal_err: float | None 

124 pedestal_scale_cov: float | None 

125 gradient_x: float | None 

126 gradient_y: float | None 

127 curvature_x: float | None 

128 curvature_y: float | None 

129 curvature_xy: float | None 

130 global_reduced_chi_squared: float | None 

131 global_degrees_of_freedom: int | None 

132 psf_reduced_chi_squared: float | None 

133 psf_degrees_of_freedom: int | None 

134 psf_masked_flux_fraction: float | None 

135 

136 # Mapping of metadata keys to attribute names 

137 _metadata_attribute_map = { 

138 "VISIT": "visit", 

139 "DETECTOR": "detector", 

140 "REF_ID": "ref_id", 

141 "REF_MAG": "ref_mag", 

142 "POSITION_X": "position.x", 

143 "POSITION_Y": "position.y", 

144 "FOCAL_PLANE_RADIUS": "focal_plane_radius", 

145 "FOCAL_PLANE_ANGLE_DEGREES": "focal_plane_angle", 

146 "SCALE": "scale", 

147 "SCALE_ERR": "scale_err", 

148 "PEDESTAL": "pedestal", 

149 "PEDESTAL_ERR": "pedestal_err", 

150 "PEDESTAL_SCALE_COV": "pedestal_scale_cov", 

151 "GRADIENT_X": "gradient_x", 

152 "GRADIENT_Y": "gradient_y", 

153 "CURVATURE_X": "curvature_x", 

154 "CURVATURE_Y": "curvature_y", 

155 "CURVATURE_XY": "curvature_xy", 

156 "GLOBAL_REDUCED_CHI_SQUARED": "global_reduced_chi_squared", 

157 "GLOBAL_DEGREES_OF_FREEDOM": "global_degrees_of_freedom", 

158 "PSF_REDUCED_CHI_SQUARED": "psf_reduced_chi_squared", 

159 "PSF_DEGREES_OF_FREEDOM": "psf_degrees_of_freedom", 

160 "PSF_MASKED_FLUX_FRACTION": "psf_masked_flux_fraction", 

161 } 

162 

163 def _getMetadata(self) -> PropertyList: 

164 """Extract metadata from the stamp's attributes. 

165 

166 This method constructs a `PropertyList` containing metadata 

167 extracted from the stamp's attributes. It is used when writing the 

168 stamp to a FITS file to store relevant metadata in the FITS headers. 

169 

170 Returns 

171 ------- 

172 metadata : `PropertyList` 

173 A `PropertyList` containing the metadata, or `None` if no 

174 metadata attributes are defined. 

175 """ 

176 metadata = PropertyList() 

177 for metadata_key, attribute_name in self._metadata_attribute_map.items(): 

178 if "." in attribute_name: 

179 top_attr, sub_attr = attribute_name.split(".") 

180 value = getattr(getattr(self, top_attr), sub_attr) 

181 elif metadata_key == "FOCAL_PLANE_ANGLE_DEGREES": 

182 value = getattr(self, attribute_name).asDegrees() 

183 else: 

184 value = getattr(self, attribute_name) 

185 metadata[metadata_key] = value 

186 return metadata 

187 

188 @property 

189 def metadata(self) -> PropertyList: 

190 """Return the stamp's metadata as a PropertyList.""" 

191 return self._getMetadata() 

192 

193 @classmethod 

194 def factory( 

195 cls, 

196 stamp_im: MaskedImageF, 

197 psf: Psf | None, 

198 wcs: SkyWcs | None, 

199 metadata: PropertyList, 

200 ) -> BrightStarStamp: 

201 """Construct a `BrightStarStamp` from image data and metadata. 

202 

203 This method provides a standard interface to create a `BrightStarStamp` 

204 from its image data, PSF, WCS, and associated metadata. 

205 It is used by the `BrightStarStamps.readFits` method to construct 

206 individual bright star stamps from FITS files. 

207 

208 Parameters 

209 ---------- 

210 stamp_im : `~lsst.afw.image.MaskedImageF` 

211 Masked image for the stamp. 

212 psf : `~lsst.afw.detection.Psf`, optional 

213 Point-spread function for the stamp. 

214 wcs : `~lsst.afw.geom.SkyWcs`, optional 

215 World coordinate system for the stamp. 

216 metadata : `PropertyList` 

217 Metadata associated with the stamp, containing keys for all 

218 required attributes. 

219 

220 Returns 

221 ------- 

222 brightStarStamp : `BrightStarStamp` 

223 The constructed `BrightStarStamp` instance. 

224 """ 

225 kwargs = {} 

226 

227 for metadata_key, attribute_name in cls._metadata_attribute_map.items(): 

228 if "." in attribute_name: # for nested attributes like position.x 

229 top_attr, sub_attr = attribute_name.split(".") 

230 if top_attr not in kwargs: # avoid overwriting position 

231 if top_attr == "position": # make an initial Point2D 

232 kwargs[top_attr] = Point2D(0, 0) 

233 setattr(kwargs[top_attr], sub_attr, metadata[metadata_key]) 

234 elif attribute_name == "focal_plane_angle": 

235 kwargs[attribute_name] = Angle(metadata[metadata_key], degrees) 

236 else: 

237 kwargs[attribute_name] = metadata[metadata_key] 

238 

239 return cls(stamp_im=stamp_im, psf=psf, wcs=wcs, **kwargs) 

240 

241 

242class BrightStarStamps(Sequence[BrightStarStamp]): 

243 """A collection of bright star stamps. 

244 

245 Parameters 

246 ---------- 

247 brightStarStamps : `Iterable` [`BrightStarStamp`] 

248 Collection of `BrightStarStamp` instances. 

249 metadata : `~lsst.daf.base.PropertyList`, optional 

250 Global metadata associated with the collection. 

251 """ 

252 

253 def __init__( 

254 self, 

255 brightStarStamps: Sequence[BrightStarStamp], 

256 metadata: PropertyList | None = None, 

257 ): 

258 self._stamps = list(brightStarStamps) 

259 self._metadata = PropertyList() if metadata is None else metadata.deepCopy() 

260 self.by_ref_id = {stamp.ref_id: stamp for stamp in self} 

261 

262 def __len__(self): 

263 return len(self._stamps) 

264 

265 def __getitem__(self, index): 

266 if isinstance(index, slice): 

267 return BrightStarStamps(self._stamps[index], metadata=self._metadata) 

268 return self._stamps[index] 

269 

270 def __iter__(self): 

271 return iter(self._stamps) 

272 

273 @property 

274 def metadata(self): 

275 """Return the collection's global metadata as a PropertyList.""" 

276 return self._metadata 

277 

278 @classmethod 

279 def readFits(cls, filename: str) -> BrightStarStamps: 

280 """Make a `BrightStarStamps` object from a FITS file. 

281 

282 Parameters 

283 ---------- 

284 filename : `str` 

285 Name of the FITS file to read. 

286 

287 Returns 

288 ------- 

289 brightStarStamps : `BrightStarStamps` 

290 The constructed `BrightStarStamps` instance. 

291 """ 

292 return cls.readFitsWithOptions(filename, None) 

293 

294 @classmethod 

295 def readFitsWithOptions(cls, filename: str, options: PropertyList | None) -> BrightStarStamps: 

296 """Make a `BrightStarStamps` object from a FITS file, with options. 

297 

298 Parameters 

299 ---------- 

300 filename : `str` 

301 Name of the FITS file to read. 

302 options : `~lsst.daf.base.PropertyList`, optional 

303 Options for reading the FITS file. Not currently used. 

304 

305 Returns 

306 ------- 

307 brightStarStamps : `BrightStarStamps` 

308 The constructed `BrightStarStamps` instance. 

309 """ 

310 with Fits(filename, "r") as fits_file: 

311 stamp_planes = {} 

312 stamp_psf_ids = {} 

313 stamp_wcs_ids = {} 

314 stamp_metadata = {} 

315 archive = None 

316 

317 for hdu_num in range(1, fits_file.countHdus()): # Skip primary HDU 

318 metadata = readMetadata(filename, hdu=hdu_num) 

319 extname = metadata["EXTNAME"] 

320 stamp_id: int | None = metadata.get("EXTVER", None) 

321 

322 # Skip non-image BINTABLEs (except ARCHIVE_INDEX) 

323 if metadata["XTENSION"] == "BINTABLE" and not metadata.get("ZIMAGE", False): 

324 if extname != "ARCHIVE_INDEX": 

325 continue 

326 

327 # Handle the archive index separately 

328 if extname == "ARCHIVE_INDEX": 

329 fits_file.setHdu(hdu_num) 

330 archive = InputArchive.readFits(fits_file) 

331 continue 

332 elif metadata.get("EXTTYPE") == "ARCHIVE_DATA": 

333 continue 

334 

335 # Select reader and dtype 

336 if extname == "IMAGE": 

337 reader = ImageFitsReader(filename, hdu=hdu_num) 

338 dtype = np.dtype(MaskedImageF.dtype) 

339 stamp_psf_ids[stamp_id] = metadata.pop("PSF", None) 

340 stamp_wcs_ids[stamp_id] = metadata.pop("WCS", None) 

341 stamp_metadata[stamp_id] = metadata 

342 elif extname == "MASK": 

343 reader = MaskFitsReader(filename, hdu=hdu_num) 

344 dtype = None 

345 elif extname == "VARIANCE": 

346 reader = ImageFitsReader(filename, hdu=hdu_num) 

347 dtype = np.dtype("float32") 

348 else: 

349 raise ValueError(f"Unknown extension type: {extname}") 

350 

351 if stamp_id is not None: 

352 stamp_planes.setdefault(stamp_id, {})[extname.lower()] = reader.read(dtype=dtype) 

353 

354 primary_metadata = readMetadata(filename, hdu=0) 

355 num_stamps = primary_metadata["N_STAMPS"] 

356 

357 if len(stamp_planes) != num_stamps: 

358 raise ValueError( 

359 f"Number of stamps read ({len(stamp_planes)}) does not agree with the " 

360 f"number of stamps recorded in the primary HDU metadata ({num_stamps})." 

361 ) 

362 if archive is None: 

363 raise ValueError("No archive index was found in the FITS file; cannot read PSF or WCS.") 

364 

365 brightStarStamps = [] 

366 for stamp_id in range(1, num_stamps + 1): # Need to increment by one as EXTVER starts at 1 

367 stamp = MaskedImageF(**stamp_planes[stamp_id]) 

368 psf = archive.get(stamp_psf_ids[stamp_id]) 

369 wcs = archive.get(stamp_wcs_ids[stamp_id]) 

370 brightStarStamps.append(BrightStarStamp.factory(stamp, psf, wcs, stamp_metadata[stamp_id])) 

371 

372 return cls(brightStarStamps, primary_metadata) 

373 

374 def writeFits(self, filename: str): 

375 """Write this `BrightStarStamps` object to a FITS file. 

376 

377 Parameters 

378 ---------- 

379 filename : `str` 

380 Name of the FITS file to write. 

381 """ 

382 metadata = self._metadata.deepCopy() 

383 

384 # Store metadata in the primary HDU 

385 metadata["N_STAMPS"] = len(self._stamps) 

386 metadata["VERSION"] = 2 # Record version number in case of future code changes 

387 metadata["STAMPCLS"] = get_full_type_name(self) 

388 

389 # Create and write to the FITS file within a context manager 

390 with Fits(filename, "w") as fits_file: 

391 fits_file.createEmpty() 

392 

393 # Store Persistables in an OutputArchive 

394 output_archive = OutputArchive() 

395 stamp_psf_ids = [] 

396 stamp_wcs_ids = [] 

397 for stamp in self._stamps: 

398 stamp_psf_ids.append(output_archive.put(stamp.psf)) 

399 stamp_wcs_ids.append(output_archive.put(stamp.wcs)) 

400 

401 # Write to the FITS file 

402 fits_file.writeMetadata(metadata) 

403 del metadata 

404 output_archive.writeFits(fits_file) 

405 

406 # Add all pixel data to extension HDUs; note: EXTVER should be 1-based 

407 for stamp_id, (stamp, stamp_psf_id, stamp_wcs_id) in enumerate( 

408 zip(self._stamps, stamp_psf_ids, stamp_wcs_ids), 

409 start=1, 

410 ): 

411 metadata = PropertyList() 

412 metadata.update({"EXTVER": stamp_id, "EXTNAME": "IMAGE"}) 

413 if stamp_metadata := stamp._getMetadata(): 

414 metadata.update(stamp_metadata) 

415 metadata["PSF"] = stamp_psf_id 

416 metadata["WCS"] = stamp_wcs_id 

417 stamp.stamp_im.getImage().writeFits(filename, metadata=metadata, mode="a") 

418 

419 metadata = PropertyList() 

420 metadata.update({"EXTVER": stamp_id, "EXTNAME": "MASK"}) 

421 stamp.stamp_im.getMask().writeFits(filename, metadata=metadata, mode="a") 

422 

423 metadata = PropertyList() 

424 metadata.update({"EXTVER": stamp_id, "EXTNAME": "VARIANCE"}) 

425 stamp.stamp_im.getVariance().writeFits(filename, metadata=metadata, mode="a")