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

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 ( 

39 DetectorFrame, 

40 SkyProjection, 

41 SkyProjectionAstropyView, 

42 SkyProjectionSerializationModel, 

43) 

44from .aperture_corrections import ( 

45 ApertureCorrectionMap, 

46 ApertureCorrectionMapSerializationModel, 

47 aperture_corrections_from_legacy, 

48 aperture_corrections_to_legacy, 

49) 

50from .cameras import Detector, DetectorSerializationModel 

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

52from .fits import FitsOpaqueMetadata 

53from .psfs import ( 

54 GaussianPointSpreadFunction, 

55 GaussianPSFSerializationModel, 

56 LegacyPointSpreadFunction, 

57 PiffSerializationModel, 

58 PiffWrapper, 

59 PointSpreadFunction, 

60 PSFExSerializationModel, 

61 PSFExWrapper, 

62) 

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

64from .utils import is_none 

65 

66if TYPE_CHECKING: 

67 try: 

68 from lsst.afw.cameraGeom import Detector as LegacyDetector 

69 from lsst.afw.image import Exposure as LegacyExposure 

70 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

71 from lsst.afw.image import VisitInfo as LegacyVisitInfo 

72 except ImportError: 

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

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

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

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

77 

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

79 

80 

81class VisitImage(MaskedImage): 

82 """A calibrated single-visit image. 

83 

84 Parameters 

85 ---------- 

86 image 

87 The main image plane. If this has a `SkyProjection`, it will be used 

88 for all planes unless a ``sky_projection`` is passed separately. 

89 mask 

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

91 same bounding box as ``image`` if provided. Any attached 

92 ``sky_projection`` is replaced (possibly by `None`). 

93 variance 

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

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

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

97 Values default to ``1.0``. Any attached ``sky_projection`` is replaced 

98 (possibly by `None`). 

99 mask_schema 

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

101 not provided. 

102 sky_projection 

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

104 a ``sky_projection`` is already attached to ``image``. 

105 bounds 

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

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

108 components (``psf``, ``sky_projection``, ``aperture_corrections``, 

109 etc.) are assumed to have their own bounds which may or may not be the 

110 same as the image bounds. If ``bounds`` extends beyond the image 

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

112 bounding box is used instead. 

113 obs_info 

114 General information about this visit in standardized form. 

115 summary_stats 

116 Summary statistics associated with this visit. Initialized to default 

117 values if not provided. 

118 photometric_scaling 

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

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

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

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

123 psf 

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

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

126 detector 

127 Geometry and electronic information for the detector attached to this 

128 image. 

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

130 Mapping from photometry algorithm name to the aperture correction for 

131 that algorithm. 

132 backgrounds 

133 Background models associated with this image. 

134 band 

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

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

137 metadata 

138 Arbitrary flexible metadata to associate with the image. 

139 """ 

140 

141 def __init__( 

142 self, 

143 image: Image, 

144 *, 

145 mask: Mask | None = None, 

146 variance: Image | None = None, 

147 mask_schema: MaskSchema | None = None, 

148 sky_projection: SkyProjection[DetectorFrame] | None = None, 

149 bounds: Bounds | None = None, 

150 obs_info: ObservationInfo | None = None, 

151 summary_stats: ObservationSummaryStats | None = None, 

152 photometric_scaling: Field | None = None, 

153 psf: PointSpreadFunction | ArchiveReadError, 

154 detector: Detector, 

155 aperture_corrections: ApertureCorrectionMap | None = None, 

156 backgrounds: BackgroundMap | None = None, 

157 band: str, 

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

159 ): 

160 super().__init__( 

161 image, 

162 mask=mask, 

163 variance=variance, 

164 mask_schema=mask_schema, 

165 sky_projection=sky_projection, 

166 metadata=metadata, 

167 ) 

168 if self.image.unit is None: 

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

170 if self.image.sky_projection is None: 

171 raise TypeError("The sky_projection component of a VisitImage cannot be None.") 

172 if obs_info is None: 

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

174 if obs_info.physical_filter is None: 174 ↛ 175line 174 didn't jump to line 175 because the condition on line 174 was never true

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

176 self._obs_info = obs_info 

177 if not isinstance(self.image.sky_projection.pixel_frame, DetectorFrame): 

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

179 if summary_stats is None: 

180 summary_stats = ObservationSummaryStats() 

181 self._summary_stats = summary_stats 

182 if photometric_scaling is not None and photometric_scaling.unit is None: 182 ↛ 183line 182 didn't jump to line 183 because the condition on line 182 was never true

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

184 self._photometric_scaling = photometric_scaling 

185 self._psf = psf 

186 self._detector = detector 

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

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

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

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

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

192 self._band = band 

193 

194 @property 

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

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

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

198 

199 @property 

200 def sky_projection(self) -> SkyProjection[DetectorFrame]: 

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

202 (`SkyProjection` [`DetectorFrame`]). 

203 """ 

204 return cast(SkyProjection[DetectorFrame], super().sky_projection) 

205 

206 @property 

207 def bounds(self) -> Bounds: 

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

209 return self._bounds 

210 

211 @property 

212 def obs_info(self) -> ObservationInfo: 

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

214 (`~astro_metadata_translator.ObservationInfo`). 

215 """ 

216 return self._obs_info 

217 

218 @property 

219 def physical_filter(self) -> str: 

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

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

222 return self._obs_info.physical_filter 

223 

224 @property 

225 def band(self) -> str: 

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

227 return self._band 

228 

229 @property 

230 def astropy_wcs(self) -> SkyProjectionAstropyView: 

231 """An Astropy WCS for the pixel arrays (`SkyProjectionAstropyView`). 

232 

233 Notes 

234 ----- 

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

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

237 ``bbox.start``, as is the case for `sky_projection`. 

238 

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

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

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

242 """ 

243 return cast(SkyProjectionAstropyView, super().astropy_wcs) 

244 

245 @property 

246 def summary_stats(self) -> ObservationSummaryStats: 

247 """Optional summary statistics for this observation 

248 (`ObservationSummaryStats`). 

249 """ 

250 return self._summary_stats 

251 

252 @property 

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

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

255 image (~`fields.BaseField`). 

256 """ 

257 return self._photometric_scaling 

258 

259 @photometric_scaling.setter 

260 def photometric_scaling(self, value: Field): 

261 if value.unit is None: 

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

263 self._photometric_scaling = value 

264 

265 @property 

266 def psf(self) -> PointSpreadFunction: 

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

268 (`.psfs.PointSpreadFunction`). 

269 """ 

270 if isinstance(self._psf, ArchiveReadError): 270 ↛ 271line 270 didn't jump to line 271 because the condition on line 270 was never true

271 raise self._psf 

272 return self._psf 

273 

274 @property 

275 def detector(self) -> Detector: 

276 """Geometry and electronic information about the detector 

277 (`.cameras.Detector`). 

278 """ 

279 return self._detector 

280 

281 @property 

282 def aperture_corrections(self) -> ApertureCorrectionMap: 

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

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

285 """ 

286 return self._aperture_corrections 

287 

288 @property 

289 def backgrounds(self) -> BackgroundMap: 

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

291 (`BackgroundMap`). 

292 """ 

293 return self._backgrounds 

294 

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

296 if bbox is ...: 

297 return self 

298 super().__getitem__(bbox) 

299 return self._transfer_metadata( 

300 VisitImage( 

301 self.image[bbox], 

302 mask=self.mask[bbox], 

303 variance=self.variance[bbox], 

304 sky_projection=self.sky_projection, 

305 psf=self.psf, 

306 obs_info=self.obs_info, 

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

308 summary_stats=self.summary_stats, 

309 detector=self._detector, 

310 photometric_scaling=self._photometric_scaling, 

311 aperture_corrections=self.aperture_corrections, 

312 backgrounds=self._backgrounds, 

313 band=self._band, 

314 ), 

315 bbox=bbox, 

316 ) 

317 

318 def __str__(self) -> str: 

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

320 

321 def __repr__(self) -> str: 

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

323 

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

325 """Deep-copy the visit image. 

326 

327 Parameters 

328 ---------- 

329 copy_detector 

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

331 """ 

332 return self._transfer_metadata( 

333 VisitImage( 

334 image=self._image.copy(), 

335 mask=self._mask.copy(), 

336 variance=self._variance.copy(), 

337 psf=self._psf, 

338 obs_info=self.obs_info, 

339 bounds=self._bounds, 

340 summary_stats=self.summary_stats.model_copy(), 

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

342 photometric_scaling=self._photometric_scaling, 

343 aperture_corrections=self.aperture_corrections.copy(), 

344 backgrounds=self._backgrounds.copy(), 

345 band=self.band, 

346 ), 

347 copy=True, 

348 ) 

349 

350 def convert_unit( 

351 self, 

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

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

354 copy_detector: bool = False, 

355 ) -> VisitImage: 

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

357 

358 Parameters 

359 ---------- 

360 unit 

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

362 

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

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

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

366 instrumental units but we know how to calibrate them) 

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

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

369 calibrated units and we want to revert back to instrumental 

370 units). 

371 copy 

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

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

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

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

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

377 already in the right units. 

378 copy_detector 

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

380 

381 Returns 

382 ------- 

383 `VisitImage` 

384 An image with the given units. 

385 """ 

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

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

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

389 if factor == 1.0: 

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

391 return self.copy() 

392 else: 

393 return self[...] 

394 elif copy is False: 

395 raise astropy.units.UnitConversionError( 

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

397 ) 

398 image = Image( 

399 self._image.array * factor, bbox=self.bbox, sky_projection=self.sky_projection, unit=unit 

400 ) 

401 variance = Image( 

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

403 bbox=self.bbox, 

404 unit=unit**2, 

405 ) 

406 elif self._photometric_scaling is None: 

407 raise astropy.units.UnitConversionError( 

408 "VisitImage.photometric_scaling is None, and there " 

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

410 ) 

411 else: 

412 if copy is False: 

413 raise astropy.units.UnitConversionError( 

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

415 ) 

416 scaling = self._photometric_scaling 

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

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

419 if constant_factor != 1.0: 

420 scaling = scaling * constant_factor 

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

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

423 if constant_factor != 1.0: 

424 scaling = scaling / constant_factor 

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

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

427 else: 

428 raise astropy.units.UnitConversionError( 

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

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

431 ) 

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

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

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

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

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

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

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

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

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

441 image.array *= self._image.array 

442 variance.array *= self._variance.array 

443 copy_components = copy is True 

444 return self._transfer_metadata( 

445 VisitImage( 

446 image=image, 

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

448 variance=variance, 

449 sky_projection=self.sky_projection, # never copied; immutable 

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

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

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

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

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

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

456 aperture_corrections=( 

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

458 ), 

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

460 band=self.band, 

461 ) 

462 ) 

463 

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

465 return self._serialize_impl(VisitImageSerializationModel, archive) 

466 

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

468 # VisitImageSerializationModel, not just a MaskedImageSerializationModel, 

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

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

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

472 ) -> M: 

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

474 match self._psf: 

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

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

477 case PiffWrapper(): 477 ↛ 478line 477 didn't jump to line 478 because the pattern on line 477 never matched

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

479 case PSFExWrapper(): 479 ↛ 480line 479 didn't jump to line 480 because the pattern on line 479 never matched

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

481 case GaussianPointSpreadFunction(): 481 ↛ 483line 481 didn't jump to line 483 because the pattern on line 481 always matched

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

483 case _: 

484 raise TypeError( 

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

486 ) 

487 assert result.sky_projection is not None, "VisitImage always has a sky_projection." 

488 result.obs_info = self.obs_info 

489 result.summary_stats = self.summary_stats 

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

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

492 result.band = self.band 

493 result.photometric_scaling = ( 

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

495 # method return types. 

496 archive.serialize_direct( 

497 "photometric_scaling", 

498 self._photometric_scaling.serialize, 

499 ) # type: ignore[assignment] 

500 if self._photometric_scaling is not None 

501 else None 

502 ) 

503 result.aperture_corrections = archive.serialize_direct( 

504 "aperture_corrections", 

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

506 ) 

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

508 return result 

509 

510 @staticmethod 

511 def _get_archive_tree_type[P: pydantic.BaseModel]( 

512 pointer_type: type[P], 

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

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

515 type that uses the given pointer type. 

516 """ 

517 return VisitImageSerializationModel[pointer_type] # type: ignore 

518 

519 @staticmethod 

520 def from_legacy( 

521 legacy: LegacyExposure, 

522 *, 

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

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

525 instrument: str | None = None, 

526 visit: int | None = None, 

527 ) -> VisitImage: 

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

529 

530 Parameters 

531 ---------- 

532 legacy 

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

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

535 unit 

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

537 key will be used, if available. 

538 plane_map 

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

540 description. If `None` (default) 

541 `get_legacy_visit_image_mask_planes` is used. 

542 instrument 

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

544 provided. 

545 visit 

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

547 """ 

548 if plane_map is None: 

549 plane_map = get_legacy_visit_image_mask_planes() 

550 md = legacy.getMetadata() 

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

552 instrument = _extract_or_check_header( 

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

554 ) 

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

556 legacy_wcs = legacy.getWcs() 

557 if legacy_wcs is None: 

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

559 legacy_detector = legacy.getDetector() 

560 if legacy_detector is None: 

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

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

563 

564 # Update the ObservationInfo from other components. 

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

566 

567 opaque_fits_metadata = FitsOpaqueMetadata() 

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

569 with warnings.catch_warnings(): 

570 # Silence warnings about long keys becoming HIERARCH. 

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

572 primary_header.update(md.toOrderedDict()) 

573 metadata = opaque_fits_metadata.extract_legacy_primary_header(primary_header) 

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

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

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

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

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

579 # Fix incorrect BUNIT='adu' in LSST 

580 # preliminary_visit_image. 

581 hdr_unit = astropy.units.electron 

582 if unit is None: 

583 unit = hdr_unit 

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

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

586 sky_projection = SkyProjection.from_legacy( 

587 legacy_wcs, 

588 DetectorFrame( 

589 instrument=instrument, 

590 visit=visit, 

591 detector=legacy_detector.getId(), 

592 bbox=detector_bbox, 

593 ), 

594 ) 

595 legacy_psf = legacy.getPsf() 

596 if legacy_psf is None: 

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

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

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

600 legacy_summary_stats = legacy.info.getSummaryStats() 

601 legacy_ap_corr_map = legacy.info.getApCorrMap() 

602 legacy_polygon = legacy.info.getValidPolygon() 

603 legacy_photo_calib = legacy.info.getPhotoCalib() 

604 detector = Detector.from_legacy( 

605 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True 

606 ) 

607 _reconcile_detector_serial(obs_info, detector) 

608 result = VisitImage( 

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

610 mask=masked_image.mask, 

611 variance=masked_image.variance, 

612 sky_projection=sky_projection, 

613 psf=psf, 

614 obs_info=obs_info, 

615 summary_stats=( 

616 ObservationSummaryStats.from_legacy(legacy_summary_stats) 

617 if legacy_summary_stats is not None 

618 else None 

619 ), 

620 detector=detector, 

621 aperture_corrections=( 

622 aperture_corrections_from_legacy(legacy_ap_corr_map) 

623 if legacy_ap_corr_map is not None 

624 else None 

625 ), 

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

627 photometric_scaling=( 

628 field_from_legacy_photo_calib( 

629 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

630 ) 

631 if legacy_photo_calib is not None 

632 else None 

633 ), 

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

635 metadata=metadata, 

636 ) 

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

638 result._opaque_metadata = opaque_fits_metadata 

639 return result 

640 

641 def to_legacy( 

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

643 ) -> LegacyExposure: 

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

645 

646 Parameters 

647 ---------- 

648 copy 

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

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

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

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

653 plane_map 

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

655 description. If `None` (default), 

656 `get_legacy_visit_image_mask_planes` is used. 

657 """ 

658 from lsst.afw.image import Exposure as LegacyExposure 

659 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

660 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo 

661 

662 if plane_map is None: 

663 plane_map = get_legacy_visit_image_mask_planes() 

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

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

666 result_info = result.info 

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

668 result_info.setWcs(self.sky_projection.to_legacy()) 

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

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

671 if self._photometric_scaling is not None: 

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

673 else: 

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

675 self._fill_legacy_metadata(result_info.getMetadata()) 

676 if isinstance(self._psf, LegacyPointSpreadFunction): 

677 result_info.setPsf(self._psf.legacy_psf) 

678 elif isinstance(self._psf, PiffWrapper): 

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

680 if isinstance(self.bounds, Polygon): 

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

682 if self.aperture_corrections: 

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

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

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

686 return result 

687 

688 @staticmethod 

689 def read_legacy( # type: ignore[override] 

690 filename: str, 

691 *, 

692 preserve_quantization: bool = False, 

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

694 instrument: str | None = None, 

695 visit: int | None = None, 

696 component: Literal[ 

697 "bbox", 

698 "image", 

699 "mask", 

700 "variance", 

701 "sky_projection", 

702 "psf", 

703 "detector", 

704 "photometric_scaling", 

705 "obs_info", 

706 "summary_stats", 

707 "aperture_corrections", 

708 ] 

709 | None = None, 

710 ) -> Any: 

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

712 

713 Parameters 

714 ---------- 

715 filename 

716 Full name of the file. 

717 preserve_quantization 

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

719 exactly preserve quantization-compressed pixel values. This causes 

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

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

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

723 not transferred to the copy. 

724 plane_map 

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

726 description. If `None` (default) 

727 `get_legacy_visit_image_mask_planes` is used. 

728 instrument 

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

730 provided. 

731 visit 

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

733 provided. 

734 component 

735 A component to read instead of the full image. 

736 """ 

737 from lsst.afw.image import ExposureFitsReader 

738 

739 reader = ExposureFitsReader(filename) 

740 if component == "bbox": 

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

742 legacy_detector = reader.readDetector() 

743 if legacy_detector is None: 

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

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

746 legacy_wcs = None 

747 if component in (None, "image", "mask", "variance", "sky_projection"): 

748 legacy_wcs = reader.readWcs() 

749 if legacy_wcs is None: 

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

751 legacy_exposure_info = reader.readExposureInfo() 

752 summary_stats = None 

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

754 legacy_stats = legacy_exposure_info.getSummaryStats() 

755 if legacy_stats is not None: 

756 summary_stats = ObservationSummaryStats.from_legacy(legacy_stats) 

757 if component == "summary_stats": 

758 return summary_stats 

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

760 legacy_psf = reader.readPsf() 

761 if legacy_psf is None: 

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

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

764 if component == "psf": 

765 return psf 

766 aperture_corrections: ApertureCorrectionMap = {} 

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

768 legacy_ap_corr_map = reader.readApCorrMap() 

769 if legacy_ap_corr_map is not None: 

770 aperture_corrections = aperture_corrections_from_legacy(legacy_ap_corr_map) 

771 if component == "aperture_corrections": 

772 return aperture_corrections 

773 assert component in ( 

774 None, 

775 "image", 

776 "mask", 

777 "variance", 

778 "sky_projection", 

779 "obs_info", 

780 "detector", 

781 "photometric_scaling", 

782 ), component # for MyPy 

783 filter_label = reader.readFilter() 

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

785 primary_header = hdu_list[0].header 

786 obs_info = _obs_info_from_md(primary_header) 

787 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, filter_label) 

788 if component == "obs_info": 

789 return obs_info 

790 instrument = _extract_or_check_header( 

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

792 ) 

793 visit = _extract_or_check_header( 

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

795 ) 

796 opaque_metadata = FitsOpaqueMetadata() 

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

798 # this opaque_metadata down to MaskedImage._read_legacy_hdus 

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

800 metadata = opaque_metadata.extract_legacy_primary_header(primary_header) 

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

802 instrumental_unit = astropy.units.electron 

803 photometric_scaling: Field | None = None 

804 if component in (None, "photometric_scaling"): 

805 legacy_photo_calib = reader.readPhotoCalib() 

806 if legacy_photo_calib is not None: 

807 photometric_scaling = field_from_legacy_photo_calib( 

808 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

809 ) 

810 if component == "photometric_scaling": 

811 return photometric_scaling 

812 if component in ("detector", None): 

813 detector = Detector.from_legacy(