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 numpy as np 

26import pydantic 

27from astro_metadata_translator import ObservationInfo, VisitInfoTranslator 

28 

29from ._backgrounds import BackgroundMap, BackgroundMapSerializationModel 

30from ._concrete_bounds import BoundsSerializationModel 

31from ._geom import Bounds, Box 

32from ._image import Image, ImageSerializationModel 

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

34from ._masked_image import MaskedImage, MaskedImageSerializationModel 

35from ._observation_summary_stats import ObservationSummaryStats 

36from ._polygon import Polygon 

37from ._transforms import ( 

38 DetectorFrame, 

39 SkyProjection, 

40 SkyProjectionAstropyView, 

41 SkyProjectionSerializationModel, 

42) 

43from .aperture_corrections import ( 

44 ApertureCorrectionMap, 

45 ApertureCorrectionMapSerializationModel, 

46 aperture_corrections_from_legacy, 

47 aperture_corrections_to_legacy, 

48) 

49from .cameras import Detector, DetectorSerializationModel 

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

51from .fits import FitsOpaqueMetadata 

52from .psfs import ( 

53 GaussianPointSpreadFunction, 

54 GaussianPSFSerializationModel, 

55 LegacyPointSpreadFunction, 

56 PiffSerializationModel, 

57 PiffWrapper, 

58 PointSpreadFunction, 

59 PSFExSerializationModel, 

60 PSFExWrapper, 

61) 

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

63from .utils import is_none 

64 

65if TYPE_CHECKING: 

66 try: 

67 from lsst.afw.cameraGeom import Detector as LegacyDetector 

68 from lsst.afw.image import Exposure as LegacyExposure 

69 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

70 from lsst.afw.image import VisitInfo as LegacyVisitInfo 

71 except ImportError: 

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

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

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

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

76 

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

78 

79 

80class VisitImage(MaskedImage): 

81 """A calibrated single-visit image. 

82 

83 Parameters 

84 ---------- 

85 image 

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

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

88 mask 

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

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

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

92 variance 

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

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

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

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

97 (possibly by `None`). 

98 mask_schema 

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

100 not provided. 

101 sky_projection 

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

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

104 bounds 

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

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

107 components (``psf``, ``sky_projection``, ``aperture_corrections``, 

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

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

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

111 bounding box is used instead. 

112 obs_info 

113 General information about this visit in standardized form. 

114 summary_stats 

115 Summary statistics associated with this visit. Initialized to default 

116 values if not provided. 

117 photometric_scaling 

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

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

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

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

122 psf 

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

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

125 detector 

126 Geometry and electronic information for the detector attached to this 

127 image. 

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

129 Mapping from photometry algorithm name to the aperture correction for 

130 that algorithm. 

131 backgrounds 

132 Background models associated with this image. 

133 band 

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

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

136 metadata 

137 Arbitrary flexible metadata to associate with the image. 

138 """ 

139 

140 def __init__( 

141 self, 

142 image: Image, 

143 *, 

144 mask: Mask | None = None, 

145 variance: Image | None = None, 

146 mask_schema: MaskSchema | None = None, 

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

148 bounds: Bounds | None = None, 

149 obs_info: ObservationInfo | None = None, 

150 summary_stats: ObservationSummaryStats | None = None, 

151 photometric_scaling: Field | None = None, 

152 psf: PointSpreadFunction | ArchiveReadError, 

153 detector: Detector, 

154 aperture_corrections: ApertureCorrectionMap | None = None, 

155 backgrounds: BackgroundMap | None = None, 

156 band: str, 

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

158 ) -> None: 

159 super().__init__( 

160 image, 

161 mask=mask, 

162 variance=variance, 

163 mask_schema=mask_schema, 

164 sky_projection=sky_projection, 

165 metadata=metadata, 

166 ) 

167 if self.image.unit is None: 

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

169 if self.image.sky_projection is None: 

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

171 if obs_info is None: 

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

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

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

175 self._obs_info = obs_info 

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

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

178 if summary_stats is None: 

179 summary_stats = ObservationSummaryStats() 

180 self._summary_stats = summary_stats 

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

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

183 self._photometric_scaling = photometric_scaling 

184 self._psf = psf 

185 self._detector = detector 

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

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

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

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

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

191 self._band = band 

192 

193 @property 

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

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

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

197 

198 @property 

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

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

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

202 """ 

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

204 

205 @property 

206 def bounds(self) -> Bounds: 

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

208 return self._bounds 

209 

210 @property 

211 def obs_info(self) -> ObservationInfo: 

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

213 (`~astro_metadata_translator.ObservationInfo`). 

214 """ 

215 return self._obs_info 

216 

217 @property 

218 def physical_filter(self) -> str: 

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

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

221 return self._obs_info.physical_filter 

222 

223 @property 

224 def band(self) -> str: 

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

226 return self._band 

227 

228 @property 

229 def astropy_wcs(self) -> SkyProjectionAstropyView: 

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

231 

232 Notes 

233 ----- 

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

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

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

237 

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

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

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

241 """ 

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

243 

244 @property 

245 def summary_stats(self) -> ObservationSummaryStats: 

246 """Optional summary statistics for this observation 

247 (`ObservationSummaryStats`). 

248 """ 

249 return self._summary_stats 

250 

251 @property 

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

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

254 image (~`fields.BaseField`). 

255 """ 

256 return self._photometric_scaling 

257 

258 @photometric_scaling.setter 

259 def photometric_scaling(self, value: Field) -> None: 

260 if value.unit is None: 

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

262 self._photometric_scaling = value 

263 

264 @property 

265 def psf(self) -> PointSpreadFunction: 

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

267 (`.psfs.PointSpreadFunction`). 

268 """ 

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

270 raise self._psf 

271 return self._psf 

272 

273 @property 

274 def detector(self) -> Detector: 

275 """Geometry and electronic information about the detector 

276 (`.cameras.Detector`). 

277 """ 

278 return self._detector 

279 

280 @property 

281 def aperture_corrections(self) -> ApertureCorrectionMap: 

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

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

284 """ 

285 return self._aperture_corrections 

286 

287 @property 

288 def backgrounds(self) -> BackgroundMap: 

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

290 (`BackgroundMap`). 

291 """ 

292 return self._backgrounds 

293 

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

295 if bbox is ...: 

296 return self 

297 super().__getitem__(bbox) 

298 return self._transfer_metadata( 

299 VisitImage( 

300 self.image[bbox], 

301 mask=self.mask[bbox], 

302 variance=self.variance[bbox], 

303 sky_projection=self.sky_projection, 

304 psf=self.psf, 

305 obs_info=self.obs_info, 

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

307 summary_stats=self.summary_stats, 

308 detector=self._detector, 

309 photometric_scaling=self._photometric_scaling, 

310 aperture_corrections=self.aperture_corrections, 

311 backgrounds=self._backgrounds, 

312 band=self._band, 

313 ), 

314 bbox=bbox, 

315 ) 

316 

317 def __str__(self) -> str: 

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

319 

320 def __repr__(self) -> str: 

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

322 

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

324 """Deep-copy the visit image. 

325 

326 Parameters 

327 ---------- 

328 copy_detector 

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

330 """ 

331 return self._transfer_metadata( 

332 VisitImage( 

333 image=self._image.copy(), 

334 mask=self._mask.copy(), 

335 variance=self._variance.copy(), 

336 psf=self._psf, 

337 obs_info=self.obs_info, 

338 bounds=self._bounds, 

339 summary_stats=self.summary_stats.model_copy(), 

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

341 photometric_scaling=self._photometric_scaling, 

342 aperture_corrections=self.aperture_corrections.copy(), 

343 backgrounds=self._backgrounds.copy(), 

344 band=self.band, 

345 ), 

346 copy=True, 

347 ) 

348 

349 def convert_unit( 

350 self, 

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

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

353 copy_detector: bool = False, 

354 ) -> VisitImage: 

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

356 

357 Parameters 

358 ---------- 

359 unit 

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

361 

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

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

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

365 instrumental units but we know how to calibrate them) 

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

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

368 calibrated units and we want to revert back to instrumental 

369 units). 

370 copy 

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

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

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

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

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

376 already in the right units. 

377 copy_detector 

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

379 

380 Returns 

381 ------- 

382 `VisitImage` 

383 An image with the given units. 

384 """ 

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

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

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

388 if factor == 1.0: 

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

390 return self.copy() 

391 else: 

392 return self[...] 

393 elif copy is False: 

394 raise astropy.units.UnitConversionError( 

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

396 ) 

397 image = Image( 

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

399 ) 

400 variance = Image( 

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

402 bbox=self.bbox, 

403 unit=unit**2, 

404 ) 

405 elif self._photometric_scaling is None: 

406 raise astropy.units.UnitConversionError( 

407 "VisitImage.photometric_scaling is None, and there " 

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

409 ) 

410 else: 

411 if copy is False: 

412 raise astropy.units.UnitConversionError( 

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

414 ) 

415 scaling = self._photometric_scaling 

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

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

418 if constant_factor != 1.0: 

419 scaling = scaling * constant_factor 

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

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

422 if constant_factor != 1.0: 

423 scaling = scaling / constant_factor 

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

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

426 else: 

427 raise astropy.units.UnitConversionError( 

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

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

430 ) 

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

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

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

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

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

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

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

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

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

440 image.array *= self._image.array 

441 variance.array *= self._variance.array 

442 copy_components = copy is True 

443 return self._transfer_metadata( 

444 VisitImage( 

445 image=image, 

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

447 variance=variance, 

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

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

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

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

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

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

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

455 aperture_corrections=( 

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

457 ), 

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

459 band=self.band, 

460 ) 

461 ) 

462 

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

464 return self._serialize_impl(VisitImageSerializationModel, archive) 

465 

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

467 # VisitImageSerializationModel, not just a MaskedImageSerializationModel, 

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

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

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

471 ) -> M: 

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

473 match self._psf: 

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

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

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

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

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

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

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

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

482 case _: 

483 raise TypeError( 

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

485 ) 

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

487 result.obs_info = self.obs_info 

488 result.summary_stats = self.summary_stats 

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

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

491 result.band = self.band 

492 result.photometric_scaling = ( 

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

494 # method return types. 

495 archive.serialize_direct( 

496 "photometric_scaling", 

497 self._photometric_scaling.serialize, 

498 ) # type: ignore[assignment] 

499 if self._photometric_scaling is not None 

500 else None 

501 ) 

502 result.aperture_corrections = archive.serialize_direct( 

503 "aperture_corrections", 

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

505 ) 

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

507 return result 

508 

509 @staticmethod 

510 def _get_archive_tree_type[P: pydantic.BaseModel]( 

511 pointer_type: type[P], 

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

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

514 type that uses the given pointer type. 

515 """ 

516 return VisitImageSerializationModel[pointer_type] # type: ignore 

517 

518 @staticmethod 

519 def from_legacy( 

520 legacy: LegacyExposure, 

521 *, 

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

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

524 instrument: str | None = None, 

525 visit: int | None = None, 

526 ) -> VisitImage: 

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

528 

529 Parameters 

530 ---------- 

531 legacy 

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

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

534 unit 

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

536 key will be used, if available. 

537 plane_map 

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

539 description. If `None` (default) 

540 `get_legacy_visit_image_mask_planes` is used. 

541 instrument 

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

543 provided. 

544 visit 

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

546 """ 

547 if plane_map is None: 

548 plane_map = get_legacy_visit_image_mask_planes() 

549 md = legacy.getMetadata() 

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

551 instrument = _extract_or_check_header( 

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

553 ) 

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

555 legacy_wcs = legacy.getWcs() 

556 if legacy_wcs is None: 

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

558 legacy_detector = legacy.getDetector() 

559 if legacy_detector is None: 

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

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

562 

563 # Update the ObservationInfo from other components. 

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

565 

566 opaque_fits_metadata = FitsOpaqueMetadata() 

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

568 with warnings.catch_warnings(): 

569 # Silence warnings about long keys becoming HIERARCH. 

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

571 primary_header.update(md.toOrderedDict()) 

572 metadata = opaque_fits_metadata.extract_legacy_primary_header(primary_header) 

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

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

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

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

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

578 # Fix incorrect BUNIT='adu' in LSST 

579 # preliminary_visit_image. 

580 hdr_unit = astropy.units.electron 

581 if unit is None: 

582 unit = hdr_unit 

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

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

585 sky_projection = SkyProjection.from_legacy( 

586 legacy_wcs, 

587 DetectorFrame( 

588 instrument=instrument, 

589 visit=visit, 

590 detector=legacy_detector.getId(), 

591 bbox=detector_bbox, 

592 ), 

593 ) 

594 legacy_psf = legacy.getPsf() 

595 if legacy_psf is None: 

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

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

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

599 legacy_summary_stats = legacy.info.getSummaryStats() 

600 legacy_ap_corr_map = legacy.info.getApCorrMap() 

601 legacy_polygon = legacy.info.getValidPolygon() 

602 legacy_photo_calib = legacy.info.getPhotoCalib() 

603 detector = Detector.from_legacy( 

604 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True 

605 ) 

606 _reconcile_detector_serial(obs_info, detector) 

607 result = VisitImage( 

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

609 mask=masked_image.mask, 

610 variance=masked_image.variance, 

611 sky_projection=sky_projection, 

612 psf=psf, 

613 obs_info=obs_info, 

614 summary_stats=( 

615 ObservationSummaryStats.from_legacy(legacy_summary_stats) 

616 if legacy_summary_stats is not None 

617 else None 

618 ), 

619 detector=detector, 

620 aperture_corrections=( 

621 aperture_corrections_from_legacy(legacy_ap_corr_map) 

622 if legacy_ap_corr_map is not None 

623 else None 

624 ), 

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

626 photometric_scaling=( 

627 field_from_legacy_photo_calib( 

628 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

629 ) 

630 if legacy_photo_calib is not None 

631 else None 

632 ), 

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

634 metadata=metadata, 

635 ) 

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

637 result._opaque_metadata = opaque_fits_metadata 

638 return result 

639 

640 def to_legacy( 

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

642 ) -> LegacyExposure: 

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

644 

645 Parameters 

646 ---------- 

647 copy 

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

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

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

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

652 plane_map 

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

654 description. If `None` (default), 

655 `get_legacy_visit_image_mask_planes` is used. 

656 """ 

657 from lsst.afw.image import Exposure as LegacyExposure 

658 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

659 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo 

660 

661 if plane_map is None: 

662 plane_map = get_legacy_visit_image_mask_planes() 

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

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

665 result_info = result.info 

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

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

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

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

670 if self._photometric_scaling is not None: 

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

672 else: 

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

674 self._fill_legacy_metadata(result_info.getMetadata()) 

675 if isinstance(self._psf, LegacyPointSpreadFunction): 

676 result_info.setPsf(self._psf.legacy_psf) 

677 elif isinstance(self._psf, PiffWrapper): 

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

679 if isinstance(self.bounds, Polygon): 

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

681 if self.aperture_corrections: 

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

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

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

685 return result 

686 

687 @staticmethod 

688 def read_legacy( # type: ignore[override] 

689 filename: str, 

690 *, 

691 preserve_quantization: bool = False, 

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

693 instrument: str | None = None, 

694 visit: int | None = None, 

695 component: Literal[ 

696 "bbox", 

697 "image", 

698 "mask", 

699 "variance", 

700 "sky_projection", 

701 "psf", 

702 "detector", 

703 "photometric_scaling", 

704 "obs_info", 

705 "summary_stats", 

706 "aperture_corrections", 

707 ] 

708 | None = None, 

709 ) -> Any: 

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

711 

712 Parameters 

713 ---------- 

714 filename 

715 Full name of the file. 

716 preserve_quantization 

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

718 exactly preserve quantization-compressed pixel values. This causes 

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

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

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

722 not transferred to the copy. 

723 plane_map 

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

725 description. If `None` (default) 

726 `get_legacy_visit_image_mask_planes` is used. 

727 instrument 

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

729 provided. 

730 visit 

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

732 provided. 

733 component 

734 A component to read instead of the full image. 

735 """ 

736 from lsst.afw.image import ExposureFitsReader 

737 

738 reader = ExposureFitsReader(filename) 

739 if component == "bbox": 

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

741 legacy_detector = reader.readDetector() 

742 if legacy_detector is None: 

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

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

745 legacy_wcs = None 

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

747 legacy_wcs = reader.readWcs() 

748 if legacy_wcs is None: 

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

750 legacy_exposure_info = reader.readExposureInfo() 

751 summary_stats = None 

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

753 legacy_stats = legacy_exposure_info.getSummaryStats() 

754 if legacy_stats is not None: 

755 summary_stats = ObservationSummaryStats.from_legacy(legacy_stats) 

756 if component == "summary_stats": 

757 return summary_stats 

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

759 legacy_psf = reader.readPsf() 

760 if legacy_psf is None: 

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

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

763 if component == "psf": 

764 return psf 

765 aperture_corrections: ApertureCorrectionMap = {} 

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

767 legacy_ap_corr_map = reader.readApCorrMap() 

768 if legacy_ap_corr_map is not None: 

769 aperture_corrections = aperture_corrections_from_legacy(legacy_ap_corr_map) 

770 if component == "aperture_corrections": 

771 return aperture_corrections 

772 assert component in ( 

773 None, 

774 "image", 

775 "mask", 

776 "variance", 

777 "sky_projection", 

778 "obs_info", 

779 "detector", 

780 "photometric_scaling", 

781 ), component # for MyPy 

782 filter_label = reader.readFilter() 

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

784 primary_header = hdu_list[0].header 

785 obs_info = _obs_info_from_md(primary_header) 

786 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, filter_label) 

787 if component == "obs_info": 

788 return obs_info 

789 instrument = _extract_or_check_header( 

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

791 ) 

792 visit = _extract_or_check_header( 

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

794 ) 

795 opaque_metadata = FitsOpaqueMetadata() 

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

797 # this opaque_metadata down to MaskedImage._read_legacy_hdus 

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

799 metadata = opaque_metadata.extract_legacy_primary_header(primary_header) 

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

801 instrumental_unit = astropy.units.electron 

802 photometric_scaling: Field | None = None 

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

804 legacy_photo_calib = reader.readPhotoCalib() 

805 if legacy_photo_calib is not None: 

806 photometric_scaling = field_from_legacy_photo_calib( 

807 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

808 ) 

809 if component == "photometric_scaling": 

810 return photometric_scaling 

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

812 detector = Detector.from_legacy(