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 bbox, _ = self._handle_getitem_args(bbox) 

296 return self._transfer_metadata( 

297 VisitImage( 

298 self.image[bbox], 

299 mask=self.mask[bbox], 

300 variance=self.variance[bbox], 

301 sky_projection=self.sky_projection, 

302 psf=self.psf, 

303 obs_info=self.obs_info, 

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

305 summary_stats=self.summary_stats, 

306 detector=self._detector, 

307 photometric_scaling=self._photometric_scaling, 

308 aperture_corrections=self.aperture_corrections, 

309 backgrounds=self._backgrounds, 

310 band=self._band, 

311 ), 

312 bbox=bbox, 

313 ) 

314 

315 def __str__(self) -> str: 

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

317 

318 def __repr__(self) -> str: 

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

320 

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

322 """Deep-copy the visit image. 

323 

324 Parameters 

325 ---------- 

326 copy_detector 

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

328 """ 

329 return self._transfer_metadata( 

330 VisitImage( 

331 image=self._image.copy(), 

332 mask=self._mask.copy(), 

333 variance=self._variance.copy(), 

334 psf=self._psf, 

335 obs_info=self.obs_info, 

336 bounds=self._bounds, 

337 summary_stats=self.summary_stats.model_copy(), 

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

339 photometric_scaling=self._photometric_scaling, 

340 aperture_corrections=self.aperture_corrections.copy(), 

341 backgrounds=self._backgrounds.copy(), 

342 band=self.band, 

343 ), 

344 copy=True, 

345 ) 

346 

347 def convert_unit( 

348 self, 

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

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

351 copy_detector: bool = False, 

352 ) -> VisitImage: 

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

354 

355 Parameters 

356 ---------- 

357 unit 

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

359 

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

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

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

363 instrumental units but we know how to calibrate them) 

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

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

366 calibrated units and we want to revert back to instrumental 

367 units). 

368 copy 

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

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

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

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

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

374 already in the right units. 

375 copy_detector 

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

377 

378 Returns 

379 ------- 

380 `VisitImage` 

381 An image with the given units. 

382 """ 

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

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

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

386 if factor == 1.0: 

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

388 return self.copy() 

389 else: 

390 return self[...] 

391 elif copy is False: 

392 raise astropy.units.UnitConversionError( 

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

394 ) 

395 image = Image( 

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

397 ) 

398 variance = Image( 

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

400 bbox=self.bbox, 

401 unit=unit**2, 

402 ) 

403 elif self._photometric_scaling is None: 

404 raise astropy.units.UnitConversionError( 

405 "VisitImage.photometric_scaling is None, and there " 

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

407 ) 

408 else: 

409 if copy is False: 

410 raise astropy.units.UnitConversionError( 

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

412 ) 

413 scaling = self._photometric_scaling 

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

415 if (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 elif (constant_factor := _get_unit_conversion_factor(self.unit / scaling.unit, unit)) is not None: 

420 if constant_factor != 1.0: 

421 scaling = scaling / constant_factor 

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

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

424 else: 

425 raise astropy.units.UnitConversionError( 

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

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

428 ) 

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

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

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

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

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

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

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

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

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

438 image.array *= self._image.array 

439 variance.array *= self._variance.array 

440 copy_components = copy is True 

441 return self._transfer_metadata( 

442 VisitImage( 

443 image=image, 

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

445 variance=variance, 

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

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

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

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

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

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

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

453 aperture_corrections=( 

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

455 ), 

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

457 band=self.band, 

458 ) 

459 ) 

460 

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

462 return self._serialize_impl(VisitImageSerializationModel, archive) 

463 

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

465 # VisitImageSerializationModel, not just a MaskedImageSerializationModel, 

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

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

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

469 ) -> M: 

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

471 match self._psf: 

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

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

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

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

476 case PSFExWrapper(): 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 GaussianPointSpreadFunction(): 478 ↛ 480line 478 didn't jump to line 480 because the pattern on line 478 always matched

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

480 case _: 

481 raise TypeError( 

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

483 ) 

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

485 result.obs_info = self.obs_info 

486 result.summary_stats = self.summary_stats 

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

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

489 result.band = self.band 

490 result.photometric_scaling = ( 

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

492 # method return types. 

493 archive.serialize_direct( 

494 "photometric_scaling", 

495 self._photometric_scaling.serialize, 

496 ) # type: ignore[assignment] 

497 if self._photometric_scaling is not None 

498 else None 

499 ) 

500 result.aperture_corrections = archive.serialize_direct( 

501 "aperture_corrections", 

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

503 ) 

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

505 return result 

506 

507 @staticmethod 

508 def _get_archive_tree_type[P: pydantic.BaseModel]( 

509 pointer_type: type[P], 

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

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

512 type that uses the given pointer type. 

513 """ 

514 return VisitImageSerializationModel[pointer_type] # type: ignore 

515 

516 @staticmethod 

517 def from_legacy( 

518 legacy: LegacyExposure, 

519 *, 

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

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

522 instrument: str | None = None, 

523 visit: int | None = None, 

524 ) -> VisitImage: 

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

526 

527 Parameters 

528 ---------- 

529 legacy 

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

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

532 unit 

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

534 key will be used, if available. 

535 plane_map 

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

537 description. If `None` (default) 

538 `get_legacy_visit_image_mask_planes` is used. 

539 instrument 

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

541 provided. 

542 visit 

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

544 """ 

545 if plane_map is None: 

546 plane_map = get_legacy_visit_image_mask_planes() 

547 md = legacy.getMetadata() 

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

549 instrument = _extract_or_check_header( 

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

551 ) 

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

553 legacy_wcs = legacy.getWcs() 

554 if legacy_wcs is None: 

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

556 legacy_detector = legacy.getDetector() 

557 if legacy_detector is None: 

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

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

560 

561 # Update the ObservationInfo from other components. 

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

563 

564 opaque_fits_metadata = FitsOpaqueMetadata() 

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

566 with warnings.catch_warnings(): 

567 # Silence warnings about long keys becoming HIERARCH. 

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

569 primary_header.update(md.toOrderedDict()) 

570 metadata = opaque_fits_metadata.extract_legacy_primary_header(primary_header) 

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

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

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

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

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

576 # Fix incorrect BUNIT='adu' in LSST 

577 # preliminary_visit_image. 

578 hdr_unit = astropy.units.electron 

579 if unit is None: 

580 unit = hdr_unit 

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

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

583 sky_projection = SkyProjection.from_legacy( 

584 legacy_wcs, 

585 DetectorFrame( 

586 instrument=instrument, 

587 visit=visit, 

588 detector=legacy_detector.getId(), 

589 bbox=detector_bbox, 

590 ), 

591 ) 

592 legacy_psf = legacy.getPsf() 

593 if legacy_psf is None: 

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

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

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

597 legacy_summary_stats = legacy.info.getSummaryStats() 

598 legacy_ap_corr_map = legacy.info.getApCorrMap() 

599 legacy_polygon = legacy.info.getValidPolygon() 

600 legacy_photo_calib = legacy.info.getPhotoCalib() 

601 detector = Detector.from_legacy( 

602 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True 

603 ) 

604 _reconcile_detector_serial(obs_info, detector) 

605 result = VisitImage( 

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

607 mask=masked_image.mask, 

608 variance=masked_image.variance, 

609 sky_projection=sky_projection, 

610 psf=psf, 

611 obs_info=obs_info, 

612 summary_stats=( 

613 ObservationSummaryStats.from_legacy(legacy_summary_stats) 

614 if legacy_summary_stats is not None 

615 else None 

616 ), 

617 detector=detector, 

618 aperture_corrections=( 

619 aperture_corrections_from_legacy(legacy_ap_corr_map) 

620 if legacy_ap_corr_map is not None 

621 else None 

622 ), 

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

624 photometric_scaling=( 

625 field_from_legacy_photo_calib( 

626 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

627 ) 

628 if legacy_photo_calib is not None 

629 else None 

630 ), 

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

632 metadata=metadata, 

633 ) 

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

635 result._opaque_metadata = opaque_fits_metadata 

636 return result 

637 

638 def to_legacy( 

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

640 ) -> LegacyExposure: 

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

642 

643 Parameters 

644 ---------- 

645 copy 

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

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

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

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

650 plane_map 

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

652 description. If `None` (default), 

653 `get_legacy_visit_image_mask_planes` is used. 

654 """ 

655 from lsst.afw.image import Exposure as LegacyExposure 

656 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

657 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo 

658 

659 if plane_map is None: 

660 plane_map = get_legacy_visit_image_mask_planes() 

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

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

663 result_info = result.info 

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

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

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

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

668 if self._photometric_scaling is not None: 

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

670 else: 

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

672 self._fill_legacy_metadata(result_info.getMetadata()) 

673 if isinstance(self._psf, LegacyPointSpreadFunction): 

674 result_info.setPsf(self._psf.legacy_psf) 

675 elif isinstance(self._psf, PiffWrapper): 

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

677 if isinstance(self.bounds, Polygon): 

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

679 if self.aperture_corrections: 

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

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

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

683 return result 

684 

685 @staticmethod 

686 def read_legacy( # type: ignore[override] 

687 filename: str, 

688 *, 

689 preserve_quantization: bool = False, 

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

691 instrument: str | None = None, 

692 visit: int | None = None, 

693 component: Literal[ 

694 "bbox", 

695 "image", 

696 "mask", 

697 "variance", 

698 "sky_projection", 

699 "psf", 

700 "detector", 

701 "photometric_scaling", 

702 "obs_info", 

703 "summary_stats", 

704 "aperture_corrections", 

705 ] 

706 | None = None, 

707 ) -> Any: 

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

709 

710 Parameters 

711 ---------- 

712 filename 

713 Full name of the file. 

714 preserve_quantization 

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

716 exactly preserve quantization-compressed pixel values. This causes 

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

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

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

720 not transferred to the copy. 

721 plane_map 

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

723 description. If `None` (default) 

724 `get_legacy_visit_image_mask_planes` is used. 

725 instrument 

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

727 provided. 

728 visit 

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

730 provided. 

731 component 

732 A component to read instead of the full image. 

733 """ 

734 from lsst.afw.image import ExposureFitsReader 

735 

736 reader = ExposureFitsReader(filename) 

737 if component == "bbox": 

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

739 legacy_detector = reader.readDetector() 

740 if legacy_detector is None: 

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

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

743 legacy_wcs = None 

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

745 legacy_wcs = reader.readWcs() 

746 if legacy_wcs is None: 

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

748 legacy_exposure_info = reader.readExposureInfo() 

749 summary_stats = None 

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

751 legacy_stats = legacy_exposure_info.getSummaryStats() 

752 if legacy_stats is not None: 

753 summary_stats = ObservationSummaryStats.from_legacy(legacy_stats) 

754 if component == "summary_stats": 

755 return summary_stats 

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

757 legacy_psf = reader.readPsf() 

758 if legacy_psf is None: 

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

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

761 if component == "psf": 

762 return psf 

763 aperture_corrections: ApertureCorrectionMap = {} 

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

765 legacy_ap_corr_map = reader.readApCorrMap() 

766 if legacy_ap_corr_map is not None: 

767 aperture_corrections = aperture_corrections_from_legacy(legacy_ap_corr_map) 

768 if component == "aperture_corrections": 

769 return aperture_corrections 

770 assert component in ( 

771 None, 

772 "image", 

773 "mask", 

774 "variance", 

775 "sky_projection", 

776 "obs_info", 

777 "detector", 

778 "photometric_scaling", 

779 ), component # for MyPy 

780 filter_label = reader.readFilter() 

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

782 primary_header = hdu_list[0].header 

783 obs_info = _obs_info_from_md(primary_header) 

784 obs_info = _update_obs_info_from_legacy(obs_info, legacy_detector, filter_label) 

785 if component == "obs_info": 

786 return obs_info 

787 instrument = _extract_or_check_header( 

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

789 ) 

790 visit = _extract_or_check_header( 

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

792 ) 

793 opaque_metadata = FitsOpaqueMetadata() 

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

795 # this opaque_metadata down to MaskedImage._read_legacy_hdus 

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

797 metadata = opaque_metadata.extract_legacy_primary_header(primary_header) 

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

799 instrumental_unit = astropy.units.electron 

800 photometric_scaling: Field | None = None 

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

802 legacy_photo_calib = reader.readPhotoCalib() 

803 if legacy_photo_calib is not None: 

804 photometric_scaling = field_from_legacy_photo_calib( 

805 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

806 ) 

807 if component == "photometric_scaling": 

808 return photometric_scaling 

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

810 detector = Detector.from_legacy( 

811 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True