Coverage for python/lsst/images/_visit_image.py: 17%

1# This file is part of lsst-images. 

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# Use of this source code is governed by a 3-clause BSD-style 

10# license that can be found in the LICENSE file. 

11 

12from __future__ import annotations 

13 

14__all__ = ("VisitImage", "VisitImageSerializationModel") 

15 

16import functools 

17import logging 

18import warnings 

19from collections.abc import Callable, Mapping, MutableMapping 

20from types import EllipsisType 

21from typing import TYPE_CHECKING, Any, ClassVar, Literal, cast 

22 

23import astropy.io.fits 

24import astropy.units 

25import astropy.wcs 

26import numpy as np 

27import pydantic 

28from astro_metadata_translator import ObservationInfo, VisitInfoTranslator 

29 

30from ._backgrounds import BackgroundMap, BackgroundMapSerializationModel 

31from ._concrete_bounds import SerializableBounds 

32from ._geom import Bounds, Box 

33from ._image import Image, ImageSerializationModel 

34from ._mask import Mask, MaskPlane, MaskSchema, MaskSerializationModel, get_legacy_visit_image_mask_planes 

35from ._masked_image import MaskedImage, MaskedImageSerializationModel 

36from ._observation_summary_stats import ObservationSummaryStats 

37from ._polygon import Polygon 

38from ._transforms import DetectorFrame, Projection, ProjectionAstropyView, ProjectionSerializationModel 

39from .aperture_corrections import ( 

40 ApertureCorrectionMap, 

41 ApertureCorrectionMapSerializationModel, 

42 aperture_corrections_from_legacy, 

43 aperture_corrections_to_legacy, 

44) 

45from .cameras import Detector, DetectorSerializationModel 

46from .fields import BaseField, Field, FieldSerializationModel, field_from_legacy_photo_calib 

47from .fits import ExtensionKey, FitsOpaqueMetadata 

48from .psfs import ( 

49 GaussianPointSpreadFunction, 

50 GaussianPSFSerializationModel, 

51 LegacyPointSpreadFunction, 

52 PiffSerializationModel, 

53 PiffWrapper, 

54 PointSpreadFunction, 

55 PSFExSerializationModel, 

56 PSFExWrapper, 

57) 

58from .serialization import ArchiveReadError, InputArchive, InvalidParameterError, MetadataValue, OutputArchive 

59from .utils import is_none 

60 

61if TYPE_CHECKING: 

62 try: 

63 from lsst.afw.cameraGeom import Detector as LegacyDetector 

64 from lsst.afw.image import Exposure as LegacyExposure 

65 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

66 from lsst.afw.image import VisitInfo as LegacyVisitInfo 

67 except ImportError: 

68 type LegacyDetector = Any # type: ignore[no-redef] 

69 type LegacyExposure = Any # type: ignore[no-redef] 

70 type LegacyFilterLabel = Any # type: ignore[no-redef] 

71 type LegacyVisitInfo = Any # type: ignore[no-redef] 

72 

73_LOG = logging.getLogger("lsst.images") 

74 

75 

76class VisitImage(MaskedImage): 

77 """A calibrated single-visit image. 

78 

79 Parameters 

80 ---------- 

81 image 

82 The main image plane. If this has a `Projection`, it will be used 

83 for all planes unless a ``projection`` is passed separately. 

84 mask 

85 A bitmask image that annotates the main image plane. Must have the 

86 same bounding box as ``image`` if provided. Any attached projection 

87 is replaced (possibly by `None`). 

88 variance 

89 The per-pixel uncertainty of the main image as an image of variance 

90 values. Must have the same bounding box as ``image`` if provided, and 

91 its units must be the square of ``image.unit`` or `None`. 

92 Values default to ``1.0``. Any attached projection is replaced 

93 (possibly by `None`). 

94 mask_schema 

95 Schema for the mask plane. Must be provided if and only if ``mask`` is 

96 not provided. 

97 projection 

98 Projection that maps the pixel grid to the sky. Can only be `None` if 

99 a projection is already attached to ``image``. 

100 bounds 

101 The region where this image's pixels and other properties are valid. 

102 If not provided, the bounding box of the image is used. Other 

103 components (``psf``, ``projection``, ``aperture_corrections``, etc.) 

104 are assumed to have their own bounds which may or may not be the same 

105 as the image bounds. If ``bounds`` extends beyond the image bounding 

106 box, the intersection between ``bounds`` and the image bounding box 

107 is used instead. 

108 obs_info 

109 General information about this visit in standardized form. 

110 summary_stats 

111 Summary statistics associated with this visit. Initialized to default 

112 values if not provided. 

113 photometric_scaling 

114 Field that can be used to multiply a post-ISR image units to yield 

115 calibrated image units. This may be a scaling that was already 

116 applied (so dividing by it will recover the post-ISR units) or a 

117 scaling that has not been applied, depending on ``image.unit``. 

118 psf 

119 Point-spread function model for this image, or an exception explaining 

120 why it could not be read (to be raised if the PSF is requested later). 

121 detector 

122 Geometry and electronic information for the detector attached to this 

123 image. 

124 aperture_corrections : `dict` [`str`, `~fields.BaseField`] 

125 Mapping from photometry algorithm name to the aperture correction for 

126 that algorithm. 

127 backgrounds 

128 Background models associated with this image. 

129 band 

130 Name of the passband the image was observed with (this is a shorter, 

131 less specific version of ``obs_info.physical_filter``). 

132 metadata 

133 Arbitrary flexible metadata to associate with the image. 

134 """ 

135 

136 def __init__( 

137 self, 

138 image: Image, 

139 *, 

140 mask: Mask | None = None, 

141 variance: Image | None = None, 

142 mask_schema: MaskSchema | None = None, 

143 projection: Projection[DetectorFrame] | None = None, 

144 bounds: Bounds | None = None, 

145 obs_info: ObservationInfo | None = None, 

146 summary_stats: ObservationSummaryStats | None = None, 

147 photometric_scaling: Field | None = None, 

148 psf: PointSpreadFunction | ArchiveReadError, 

149 detector: Detector, 

150 aperture_corrections: ApertureCorrectionMap | None = None, 

151 backgrounds: BackgroundMap | None = None, 

152 band: str, 

153 metadata: dict[str, MetadataValue] | None = None, 

154 ): 

155 super().__init__( 

156 image, 

157 mask=mask, 

158 variance=variance, 

159 mask_schema=mask_schema, 

160 projection=projection, 

161 metadata=metadata, 

162 ) 

163 if self.image.unit is None: 

164 raise TypeError("The image component of a VisitImage must have units.") 

165 if self.image.projection is None: 

166 raise TypeError("The projection component of a VisitImage cannot be None.") 

167 if obs_info is None: 

168 raise TypeError("The observation info component of a VisitImage cannot be None.") 

169 if obs_info.physical_filter is None: 

170 raise ValueError("The obs_info.physical_filter attribute of a VisitImage cannot be None.") 

171 self._obs_info = obs_info 

172 if not isinstance(self.image.projection.pixel_frame, DetectorFrame): 

173 raise TypeError("The projection's pixel frame must be a DetectorFrame for VisitImage.") 

174 if summary_stats is None: 

175 summary_stats = ObservationSummaryStats() 

176 self._summary_stats = summary_stats 

177 if photometric_scaling is not None and photometric_scaling.unit is None: 

178 raise TypeError("If a photometric_scaling is provided, it must have units.") 

179 self._photometric_scaling = photometric_scaling 

180 self._psf = psf 

181 self._detector = detector 

182 self._aperture_corrections = aperture_corrections if aperture_corrections is not None else {} 

183 self._bounds = bounds if bounds is not None else self.bbox 

184 if not self.bbox.contains(self._bounds.bbox): 

185 self._bounds = self._bounds.intersection(self.bbox) 

186 self._backgrounds = backgrounds if backgrounds is not None else BackgroundMap() 

187 self._band = band 

188 

189 @property 

190 def unit(self) -> astropy.units.UnitBase: 

191 """The units of the image plane (`astropy.units.Unit`).""" 

192 return cast(astropy.units.UnitBase, super().unit) 

193 

194 @property 

195 def projection(self) -> Projection[DetectorFrame]: 

196 """The projection that maps the pixel grid to the sky 

197 (`Projection` [`DetectorFrame`]). 

198 """ 

199 return cast(Projection[DetectorFrame], super().projection) 

200 

201 @property 

202 def bounds(self) -> Bounds: 

203 """The region where pixels are valid (`Bounds`).""" 

204 return self._bounds 

205 

206 @property 

207 def obs_info(self) -> ObservationInfo: 

208 """General information about this observation in standard form. 

209 (`~astro_metadata_translator.ObservationInfo`). 

210 """ 

211 return self._obs_info 

212 

213 @property 

214 def physical_filter(self) -> str: 

215 """Full name of the physical bandpass filter (`str`).""" 

216 assert self._obs_info.physical_filter is not None, "Guaranteed at construction." 

217 return self._obs_info.physical_filter 

218 

219 @property 

220 def band(self) -> str: 

221 """Short name of the bandpass filter (`str`).""" 

222 return self._band 

223 

224 @property 

225 def astropy_wcs(self) -> ProjectionAstropyView: 

226 """An Astropy WCS for the pixel arrays (`ProjectionAstropyView`). 

227 

228 Notes 

229 ----- 

230 As expected for Astropy WCS objects, this defines pixel coordinates 

231 such that the first row and column in the arrays are ``(0, 0)``, not 

232 ``bbox.start``, as is the case for `projection`. 

233 

234 This object satisfies the `astropy.wcs.wcsapi.BaseHighLevelWCS` and 

235 `astropy.wcs.wcsapi.BaseLowLevelWCS` interfaces, but it is not an 

236 `astropy.wcs.WCS` (use `fits_wcs` for that). 

237 """ 

238 return cast(ProjectionAstropyView, super().astropy_wcs) 

239 

240 @property 

241 def summary_stats(self) -> ObservationSummaryStats: 

242 """Optional summary statistics for this observation 

243 (`ObservationSummaryStats`). 

244 """ 

245 return self._summary_stats 

246 

247 @property 

248 def photometric_scaling(self) -> Field | None: 

249 """Field that multiplies a post-ISR image to yield the calibrated 

250 image (~`fields.BaseField`). 

251 """ 

252 return self._photometric_scaling 

253 

254 @photometric_scaling.setter 

255 def photometric_scaling(self, value: Field): 

256 if value.unit is None: 

257 raise TypeError("The photometric_scaling for a VisitImage must have units.") 

258 self._photometric_scaling = value 

259 

260 @property 

261 def psf(self) -> PointSpreadFunction: 

262 """The point-spread function model for this image 

263 (`.psfs.PointSpreadFunction`). 

264 """ 

265 if isinstance(self._psf, ArchiveReadError): 

266 raise self._psf 

267 return self._psf 

268 

269 @property 

270 def detector(self) -> Detector: 

271 """Geometry and electronic information about the detector 

272 (`.cameras.Detector`). 

273 """ 

274 return self._detector 

275 

276 @property 

277 def aperture_corrections(self) -> ApertureCorrectionMap: 

278 """A mapping from photometry algorithm name to the aperture correction 

279 field for that algorithm (`dict` [`str`, `~.fields.BaseField`]). 

280 """ 

281 return self._aperture_corrections 

282 

283 @property 

284 def backgrounds(self) -> BackgroundMap: 

285 """A mapping of backgrounds associated with this image 

286 (`BackgroundMap`). 

287 """ 

288 return self._backgrounds 

289 

290 def __getitem__(self, bbox: Box | EllipsisType) -> VisitImage: 

291 if bbox is ...: 

292 return self 

293 super().__getitem__(bbox) 

294 return self._transfer_metadata( 

295 VisitImage( 

296 self.image[bbox], 

297 mask=self.mask[bbox], 

298 variance=self.variance[bbox], 

299 projection=self.projection, 

300 psf=self.psf, 

301 obs_info=self.obs_info, 

302 bounds=self._bounds, # don't need to intersect here, because __init__ will do that. 

303 summary_stats=self.summary_stats, 

304 detector=self._detector, 

305 photometric_scaling=self._photometric_scaling, 

306 aperture_corrections=self.aperture_corrections, 

307 backgrounds=self._backgrounds, 

308 band=self._band, 

309 ), 

310 bbox=bbox, 

311 ) 

312 

313 def __str__(self) -> str: 

314 return f"VisitImage({self.image!s}, {list(self.mask.schema.names)})" 

315 

316 def __repr__(self) -> str: 

317 return f"VisitImage({self.image!r}, mask_schema={self.mask.schema!r})" 

318 

319 def copy(self, *, copy_detector: bool = False) -> VisitImage: 

320 """Deep-copy the visit image. 

321 

322 Parameters 

323 ---------- 

324 copy_detector 

325 Whether to deep-copy the `detector` attribute. 

326 """ 

327 return self._transfer_metadata( 

328 VisitImage( 

329 image=self._image.copy(), 

330 mask=self._mask.copy(), 

331 variance=self._variance.copy(), 

332 psf=self._psf, 

333 obs_info=self.obs_info, 

334 bounds=self._bounds, 

335 summary_stats=self.summary_stats.model_copy(), 

336 detector=self._detector.copy() if copy_detector else self._detector, 

337 photometric_scaling=self._photometric_scaling, 

338 aperture_corrections=self.aperture_corrections.copy(), 

339 backgrounds=self._backgrounds.copy(), 

340 band=self.band, 

341 ), 

342 copy=True, 

343 ) 

344 

345 def convert_unit( 

346 self, 

347 unit: astropy.units.UnitBase = astropy.units.nJy, 

348 copy: Literal["as-needed"] | bool = True, 

349 copy_detector: bool = False, 

350 ) -> VisitImage: 

351 """Return an equivalent image with different pixel units. 

352 

353 Parameters 

354 ---------- 

355 unit 

356 The unit to transform to. This may be any of the following: 

357 

358 - any unit directly relatable to the current units via Astropy; 

359 - any unit relatable to the product of the current units with the 

360 `photometric_scaling` (i.e. if the current image is in 

361 instrumental units but we know how to calibrate them) 

362 - any unit relatable to the quotient of the current units with the 

363 `photometric_scaling` (i.e. if the current image is in 

364 calibrated units and we want to revert back to instrumental 

365 units). 

366 copy 

367 Whether to copy the images and other components. If `True`, all 

368 components that aren't controlled by some other argument will 

369 always be deep-copied. If `False`, the operation will fail if the 

370 image is not already in the right units. If ``as-needed``, only 

371 the image and variance will be copied, and only if they are not 

372 already in the right units. 

373 copy_detector 

374 Whether to deep-copy the `detector` attribute. 

375 

376 Returns 

377 ------- 

378 `VisitImage` 

379 An image with the given units. 

380 """ 

381 if copy not in (True, False, "as-needed"): 

382 raise TypeError(f"Invalid value for 'copy' parameter: {copy!r}.") 

383 if (factor := _get_unit_conversion_factor(self.unit, unit)) is not None: 

384 if factor == 1.0: 

385 if copy is True: # not "as-needed" 

386 return self.copy() 

387 else: 

388 return self[...] 

389 elif copy is False: 

390 raise astropy.units.UnitConversionError( 

391 f"Units must be converted ({self.unit} -> {unit}), but copy=False." 

392 ) 

393 image = Image(self._image.array * factor, bbox=self.bbox, projection=self.projection, unit=unit) 

394 variance = Image( 

395 self._variance.array * factor**2, 

396 bbox=self.bbox, 

397 unit=unit**2, 

398 ) 

399 elif self._photometric_scaling is None: 

400 raise astropy.units.UnitConversionError( 

401 "VisitImage.photometric_scaling is None, and there " 

402 f"is no constant conversion from {self.unit} to {unit}." 

403 ) 

404 else: 

405 if copy is False: 

406 raise astropy.units.UnitConversionError( 

407 f"Photometric scaling must be applied to go from ={self.unit} to {unit}, but copy=False." 

408 ) 

409 scaling = self._photometric_scaling 

410 assert scaling.unit is not None, "Checked at construction." 

411 if (constant_factor := _get_unit_conversion_factor(self.unit * scaling.unit, unit)) is not None: 

412 if constant_factor != 1.0: 

413 scaling = scaling * constant_factor 

414 scaling_array = scaling.render(self.bbox, dtype=self.image.array.dtype).array 

415 elif (constant_factor := _get_unit_conversion_factor(self.unit / scaling.unit, unit)) is not None: 

416 if constant_factor != 1.0: 

417 scaling = scaling / constant_factor 

418 scaling_array = scaling.render(self.bbox, dtype=self.image.array.dtype).array 

419 np.true_divide(1.0, scaling_array, out=scaling_array) 

420 else: 

421 raise astropy.units.UnitConversionError( 

422 f"photometric_scaling with units {scaling.unit} does not " 

423 f"provide a path from {self.unit} to {unit}." 

424 ) 

425 # We needed to allocate a new array to evaluate the scaling field, 

426 # and then we need to allocate another to hold its square for the 

427 # variance scaling. But then we can multiply those arrays in-place 

428 # to get the output image and variance to avoid yet more 

429 # allocations (note we can't instead multiply the visit image's 

430 # image and variance arrays in place because they might have other 

431 # references that are still associated with the old units). 

432 image = Image(scaling_array, bbox=self.bbox, unit=unit) 

433 variance = Image(np.square(scaling_array), bbox=self.bbox, unit=unit**2) 

434 image.array *= self._image.array 

435 variance.array *= self._variance.array 

436 copy_components = copy is True 

437 return self._transfer_metadata( 

438 VisitImage( 

439 image=image, 

440 mask=self._mask if not copy_components else self._mask.copy(), 

441 variance=variance, 

442 projection=self.projection, # never copied; immutable 

443 obs_info=self.obs_info if not copy_components else self.obs_info.model_copy(), 

444 psf=self._psf, # never copied; immutable 

445 bounds=self._bounds, # never copied; immutable 

446 summary_stats=self.summary_stats if not copy_components else self.summary_stats.model_copy(), 

447 detector=self._detector if not copy_detector else self._detector.copy(), 

448 photometric_scaling=self._photometric_scaling, # never copied; immutable 

449 aperture_corrections=( 

450 self.aperture_corrections if not copy_components else self.aperture_corrections.copy() 

451 ), 

452 backgrounds=self.backgrounds if not copy_components else self.backgrounds.copy(), 

453 band=self.band, 

454 ) 

455 ) 

456 

457 def serialize(self, archive: OutputArchive[Any]) -> VisitImageSerializationModel: 

458 return self._serialize_impl(VisitImageSerializationModel, archive) 

459 

460 # This is slightly bad Liskov substitution - we're demanding M be a 

461 # VisitImageSerializationModel, not just a MaskedImageSerializationModel, 

462 # but that's because we know only `serialize` will call it. 

463 def _serialize_impl[M: VisitImageSerializationModel]( # type: ignore[override] 

464 self, model_type: type[M], archive: OutputArchive[Any] 

465 ) -> M: 

466 result = super()._serialize_impl(model_type, archive) 

467 match self._psf: 

468 # MyPy is able to figure things out here with this match statement, 

469 # but not a single isinstance check on the three types. 

470 case PiffWrapper(): 

471 result.psf = archive.serialize_direct("psf", self._psf.serialize) 

472 case PSFExWrapper(): 

473 result.psf = archive.serialize_direct("psf", self._psf.serialize) 

474 case GaussianPointSpreadFunction(): 

475 result.psf = archive.serialize_direct("psf", self._psf.serialize) 

476 case _: 

477 raise TypeError( 

478 f"Cannot serialize VisitImage with unrecognized PSF type {type(self._psf).__name__}." 

479 ) 

480 assert result.projection is not None, "VisitImage always has a projection." 

481 result.obs_info = self.obs_info 

482 result.summary_stats = self.summary_stats 

483 result.bounds = self._bounds.serialize() if self._bounds != self.bbox else None 

484 result.detector = archive.serialize_direct("detector", self._detector.serialize) 

485 result.band = self.band 

486 result.photometric_scaling = ( 

487 # MyPy can't quite follow the type union through the serialize 

488 # method return types. 

489 archive.serialize_direct( 

490 "photometric_scaling", 

491 self._photometric_scaling.serialize, 

492 ) # type: ignore[assignment] 

493 if self._photometric_scaling is not None 

494 else None 

495 ) 

496 result.aperture_corrections = archive.serialize_direct( 

497 "aperture_corrections", 

498 functools.partial(ApertureCorrectionMapSerializationModel.serialize, self.aperture_corrections), 

499 ) 

500 result.backgrounds = archive.serialize_direct("backgrounds", self._backgrounds.serialize) 

501 return result 

502 

503 @staticmethod 

504 def _get_archive_tree_type[P: pydantic.BaseModel]( 

505 pointer_type: type[P], 

506 ) -> type[VisitImageSerializationModel[P]]: 

507 """Return the serialization model type for this object for an archive 

508 type that uses the given pointer type. 

509 """ 

510 return VisitImageSerializationModel[pointer_type] # type: ignore 

511 

512 @staticmethod 

513 def from_legacy( 

514 legacy: LegacyExposure, 

515 *, 

516 unit: astropy.units.UnitBase | None = None, 

517 plane_map: Mapping[str, MaskPlane] | None = None, 

518 instrument: str | None = None, 

519 visit: int | None = None, 

520 ) -> VisitImage: 

521 """Convert from an `lsst.afw.image.Exposure` instance. 

522 

523 Parameters 

524 ---------- 

525 legacy 

526 An `lsst.afw.image.Exposure` instance that will share image and 

527 variance (but not mask) pixel data with the returned object. 

528 unit 

529 Units of the image. If not provided, the ``BUNIT`` metadata 

530 key will be used, if available. 

531 plane_map 

532 A mapping from legacy mask plane name to the new plane name and 

533 description. If `None` (default) 

534 `get_legacy_visit_image_mask_planes` is used. 

535 instrument 

536 Name of the instrument. Extracted from the metadata if not 

537 provided. 

538 visit 

539 ID of the visit. Extracted from the metadata if not provided. 

540 """ 

541 if plane_map is None: 

542 plane_map = get_legacy_visit_image_mask_planes() 

543 md = legacy.getMetadata() 

544 obs_info = _obs_info_from_md(md, visit_info=legacy.info.getVisitInfo()) 

545 instrument = _extract_or_check_header( 

546 "LSST BUTLER DATAID INSTRUMENT", instrument, md, obs_info.instrument, str 

547 ) 

548 visit = _extract_or_check_header("LSST BUTLER DATAID VISIT", visit, md, obs_info.exposure_id, int) 

549 legacy_wcs = legacy.getWcs() 

550 if legacy_wcs is None: 

551 raise ValueError("Exposure does not have a SkyWcs.") 

552 legacy_detector = legacy.getDetector() 

553 if legacy_detector is None: 

554 raise ValueError("Exposure does not have a Detector.") 

555 detector_bbox = Box.from_legacy(legacy_detector.getBBox()) 

556 

557 # Update the ObservationInfo from other components. 

558 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, legacy.info.getFilter()) 

559 

560 opaque_fits_metadata = FitsOpaqueMetadata() 

561 primary_header = astropy.io.fits.Header() 

562 with warnings.catch_warnings(): 

563 # Silence warnings about long keys becoming HIERARCH. 

564 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning) 

565 primary_header.update(md.toOrderedDict()) 

566 metadata = opaque_fits_metadata.extract_legacy_primary_header(primary_header) 

567 instrumental_unit = opaque_fits_metadata.get_instrumental_unit() or astropy.units.electron 

568 hdr_unit: astropy.units.UnitBase | None = None 

569 if hdr_unit_str := md.get("BUNIT"): 

570 hdr_unit = astropy.units.Unit(hdr_unit_str, format="FITS") 

571 if hdr_unit == astropy.units.adu and instrumental_unit == astropy.units.electron: 

572 # Fix incorrect BUNIT='adu' in LSST 

573 # preliminary_visit_image. 

574 hdr_unit = astropy.units.electron 

575 if unit is None: 

576 unit = hdr_unit 

577 elif hdr_unit is not None and hdr_unit != unit: 

578 raise ValueError(f"BUNIT value {hdr_unit} disagrees with given unit {unit}.") 

579 projection = Projection.from_legacy( 

580 legacy_wcs, 

581 DetectorFrame( 

582 instrument=instrument, 

583 visit=visit, 

584 detector=legacy_detector.getId(), 

585 bbox=detector_bbox, 

586 ), 

587 ) 

588 legacy_psf = legacy.getPsf() 

589 if legacy_psf is None: 

590 raise ValueError("Exposure file does not have a Psf.") 

591 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox) 

592 masked_image = MaskedImage.from_legacy(legacy.getMaskedImage(), unit=unit, plane_map=plane_map) 

593 legacy_summary_stats = legacy.info.getSummaryStats() 

594 legacy_ap_corr_map = legacy.info.getApCorrMap() 

595 legacy_polygon = legacy.info.getValidPolygon() 

596 legacy_photo_calib = legacy.info.getPhotoCalib() 

597 detector = Detector.from_legacy( 

598 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True 

599 ) 

600 _reconcile_detector_serial(obs_info, detector) 

601 result = VisitImage( 

602 image=masked_image.image.view(unit=unit), 

603 mask=masked_image.mask, 

604 variance=masked_image.variance, 

605 projection=projection, 

606 psf=psf, 

607 obs_info=obs_info, 

608 summary_stats=( 

609 ObservationSummaryStats.from_legacy(legacy_summary_stats) 

610 if legacy_summary_stats is not None 

611 else None 

612 ), 

613 detector=detector, 

614 aperture_corrections=( 

615 aperture_corrections_from_legacy(legacy_ap_corr_map) 

616 if legacy_ap_corr_map is not None 

617 else None 

618 ), 

619 bounds=Polygon.from_legacy(legacy_polygon) if legacy_polygon is not None else None, 

620 photometric_scaling=( 

621 field_from_legacy_photo_calib( 

622 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

623 ) 

624 if legacy_photo_calib is not None 

625 else None 

626 ), 

627 band=legacy.info.getFilter().bandLabel, 

628 metadata=metadata, 

629 ) 

630 result.metadata["id"] = legacy.info.getId() 

631 result._opaque_metadata = opaque_fits_metadata 

632 return result 

633 

634 def to_legacy( 

635 self, *, copy: bool | None = None, plane_map: Mapping[str, MaskPlane] | None = None 

636 ) -> LegacyExposure: 

637 """Convert to an `lsst.afw.image.Exposure` instance. 

638 

639 Parameters 

640 ---------- 

641 copy 

642 If `True`, always copy the image and variance pixel data. 

643 If `False`, return a view, and raise `TypeError` if the pixel data 

644 is read-only (this is not supported by afw). If `None`, only copy 

645 if the pixel data is read-only. Mask pixel data is always copied. 

646 plane_map 

647 A mapping from legacy mask plane name to the new plane name and 

648 description. If `None` (default), 

649 `get_legacy_visit_image_mask_planes` is used. 

650 """ 

651 from lsst.afw.image import Exposure as LegacyExposure 

652 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

653 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo 

654 

655 if plane_map is None: 

656 plane_map = get_legacy_visit_image_mask_planes() 

657 legacy_masked_image = super().to_legacy(copy=copy, plane_map=plane_map) 

658 result = LegacyExposure(legacy_masked_image, dtype=self.image.array.dtype) 

659 result_info = result.info 

660 result_info.setId(self.metadata.get("id")) 

661 result_info.setWcs(self.projection.to_legacy()) 

662 result_info.setDetector(self.detector.to_legacy()) 

663 result_info.setFilter(LegacyFilterLabel.fromBandPhysical(self.band, self.obs_info.physical_filter)) 

664 if self._photometric_scaling is not None: 

665 result_info.setPhotoCalib(self._photometric_scaling.to_legacy_photo_calib(self.unit)) 

666 else: 

667 result_info.setPhotoCalib(BaseField.make_legacy_photo_calib(self.unit)) 

668 # We just dump all of the FITS headers and non-FITS metadata into the 

669 # legacy metadata component, to make sure we have everything. We dump 

670 # the latter into a pair of special cards to be able to full round-trip 

671 # them (including case preservation). 

672 result_md = result_info.getMetadata() 

673 try: 

674 result_md["BUNIT"] = self.unit.to_string(format="fits") 

675 except ValueError: 

676 # Write units that astropy doesn't think FITS will accept anyway; 

677 # FITS standard says "SHOULD" about using it's recommended units, 

678 # and coloring outside the lines is better than lying. 

679 result_md["BUNIT"] = self.unit.to_string() 

680 if isinstance(self._opaque_metadata, FitsOpaqueMetadata): 

681 result_md.update(self._opaque_metadata.headers[ExtensionKey()]) 

682 for n, (k, v) in enumerate(self.metadata.items()): 

683 result_md[f"LSST IMAGES KEY {n + 1}"] = k 

684 result_md[f"LSST IMAGES VALUE {n + 1}"] = v 

685 if isinstance(self._psf, LegacyPointSpreadFunction): 

686 result_info.setPsf(self._psf.legacy_psf) 

687 elif isinstance(self._psf, PiffWrapper): 

688 result_info.setPsf(self._psf.to_legacy()) 

689 if isinstance(self.bounds, Polygon): 

690 result_info.setValidPolygon(self.bounds.to_legacy()) 

691 if self.aperture_corrections: 

692 result_info.setApCorrMap(aperture_corrections_to_legacy(self.aperture_corrections)) 

693 result_info.setVisitInfo(MakeRawVisitInfoViaObsInfo.observationInfo2visitInfo(self.obs_info)) 

694 result_info.setSummaryStats(self.summary_stats.to_legacy()) 

695 return result 

696 

697 @staticmethod 

698 def read_legacy( # type: ignore[override] 

699 filename: str, 

700 *, 

701 preserve_quantization: bool = False, 

702 plane_map: Mapping[str, MaskPlane] | None = None, 

703 instrument: str | None = None, 

704 visit: int | None = None, 

705 component: Literal[ 

706 "bbox", 

707 "image", 

708 "mask", 

709 "variance", 

710 "projection", 

711 "psf", 

712 "detector", 

713 "photometric_scaling", 

714 "obs_info", 

715 "summary_stats", 

716 "aperture_corrections", 

717 ] 

718 | None = None, 

719 ) -> Any: 

720 """Read a FITS file written by `lsst.afw.image.Exposure.writeFits`. 

721 

722 Parameters 

723 ---------- 

724 filename 

725 Full name of the file. 

726 preserve_quantization 

727 If `True`, ensure that writing the masked image back out again will 

728 exactly preserve quantization-compressed pixel values. This causes 

729 the image and variance plane arrays to be marked as read-only and 

730 stores the original binary table data for those planes in memory. 

731 If the `MaskedImage` is copied, the precompressed pixel values are 

732 not transferred to the copy. 

733 plane_map 

734 A mapping from legacy mask plane name to the new plane name and 

735 description. If `None` (default) 

736 `get_legacy_visit_image_mask_planes` is used. 

737 instrument 

738 Name of the instrument. Read from the primary header if not 

739 provided. 

740 visit 

741 ID of the visit. Read from the primary header if not 

742 provided. 

743 component 

744 A component to read instead of the full image. 

745 """ 

746 from lsst.afw.image import ExposureFitsReader 

747 

748 reader = ExposureFitsReader(filename) 

749 if component == "bbox": 

750 return Box.from_legacy(reader.readBBox()) 

751 legacy_detector = reader.readDetector() 

752 if legacy_detector is None: 

753 raise ValueError(f"Exposure file {filename!r} does not have a Detector.") 

754 detector_bbox = Box.from_legacy(legacy_detector.getBBox()) 

755 legacy_wcs = None 

756 if component in (None, "image", "mask", "variance", "projection"): 

757 legacy_wcs = reader.readWcs() 

758 if legacy_wcs is None: 

759 raise ValueError(f"Exposure file {filename!r} does not have a SkyWcs.") 

760 legacy_exposure_info = reader.readExposureInfo() 

761 summary_stats = None 

762 if component in (None, "summary_stats"): 

763 legacy_stats = legacy_exposure_info.getSummaryStats() 

764 if legacy_stats is not None: 

765 summary_stats = ObservationSummaryStats.from_legacy(legacy_stats) 

766 if component == "summary_stats": 

767 return summary_stats 

768 if component in (None, "psf"): 

769 legacy_psf = reader.readPsf() 

770 if legacy_psf is None: 

771 raise ValueError(f"Exposure file {filename!r} does not have a Psf.") 

772 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds=detector_bbox) 

773 if component == "psf": 

774 return psf 

775 aperture_corrections: ApertureCorrectionMap = {} 

776 if component in (None, "aperture_corrections"): 

777 legacy_ap_corr_map = reader.readApCorrMap() 

778 if legacy_ap_corr_map is not None: 

779 aperture_corrections = aperture_corrections_from_legacy(legacy_ap_corr_map) 

780 if component == "aperture_corrections": 

781 return aperture_corrections 

782 assert component in ( 

783 None, 

784 "image", 

785 "mask", 

786 "variance", 

787 "projection", 

788 "obs_info", 

789 "detector", 

790 "photometric_scaling", 

791 ), component # for MyPy 

792 filter_label = reader.readFilter() 

793 with astropy.io.fits.open(filename) as hdu_list: 

794 primary_header = hdu_list[0].header 

795 obs_info = _obs_info_from_md(primary_header) 

796 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, filter_label) 

797 if component == "obs_info": 

798 return obs_info 

799 instrument = _extract_or_check_header( 

800 "LSST BUTLER DATAID INSTRUMENT", instrument, primary_header, obs_info.instrument, str 

801 ) 

802 visit = _extract_or_check_header( 

803 "LSST BUTLER DATAID VISIT", visit, primary_header, obs_info.exposure_id, int 

804 ) 

805 opaque_metadata = FitsOpaqueMetadata() 

806 # This extraction is destructive, so we need to be sure to pass 

807 # this opaque_metadata down to MaskedImage._read_legacy_hdus 

808 # so it doesn't try to extract it again. 

809 metadata = opaque_metadata.extract_legacy_primary_header(primary_header) 

810 if (instrumental_unit := opaque_metadata.get_instrumental_unit()) is None: