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

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 .describe import DescribeOptions, FieldRole, Report, ReportField 

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 ) -> None: 

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) -> None: 

261 if value.unit is None: 261 ↛ 262line 261 didn't jump to line 262 because the condition on line 261 was never true

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

297 return self._transfer_metadata( 

298 VisitImage( 

299 self.image[bbox], 

300 mask=self.mask[bbox], 

301 variance=self.variance[bbox], 

302 sky_projection=self.sky_projection, 

303 psf=self.psf, 

304 obs_info=self.obs_info, 

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

306 summary_stats=self.summary_stats, 

307 detector=self._detector, 

308 photometric_scaling=self._photometric_scaling, 

309 aperture_corrections=self.aperture_corrections, 

310 backgrounds=self._backgrounds, 

311 band=self._band, 

312 ), 

313 bbox=bbox, 

314 ) 

315 

316 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report: 

317 """Return a `Report` describing this visit image. 

318 

319 Parameters 

320 ---------- 

321 options : `DescribeOptions`, optional 

322 Rendering options; forwarded to all children. Child construction 

323 can be expensive or raise for an unreadable component, so 

324 `DescribeOptions.brief` skips it. 

325 """ 

326 # The image and mask schema are rendered as children below, and the 

327 # image renders as ``Image(bbox, dtype)``, which would restate the 

328 # shared bbox. Both are REPR_ONLY so only repr sees them. 

329 fields = [ 

330 ReportField( 

331 label="image", 

332 value=self.image, 

333 repr_value=repr(self.image), 

334 positional=True, 

335 role=FieldRole.REPR_ONLY, 

336 ), 

337 ReportField( 

338 label="mask_schema", 

339 value=self.mask.schema, 

340 repr_value=repr(self.mask.schema), 

341 role=FieldRole.REPR_ONLY, 

342 ), 

343 ReportField(label="band", value=self.band, role=FieldRole.DERIVED), 

344 ReportField(label="physical_filter", value=self.physical_filter, role=FieldRole.DERIVED), 

345 ReportField(label="bbox", value=self.bbox, repr_value=repr(self.bbox), role=FieldRole.DERIVED), 

346 ] 

347 summary = f"VisitImage({self.image!s}, {list(self.mask.schema.names)})" 

348 if options.brief: 

349 return Report(type_name="VisitImage", summary=summary, fields=fields) 

350 child = options.for_child() 

351 plane = options.for_child("sky_projection", "bbox") 

352 children: dict[str, Report] = { 

353 "image": self.image._describe(plane), 

354 "mask": self.mask._describe(plane), 

355 "variance": self.variance._describe(plane), 

356 "sky_projection": self.sky_projection._describe(child, bbox=self.bbox), 

357 "psf": self.psf._describe(child), 

358 "detector": self.detector._describe(child), 

359 "summary_stats": self.summary_stats._describe(child), 

360 "backgrounds": self.backgrounds._describe(child), 

361 } 

362 if self.photometric_scaling is not None: 362 ↛ 363line 362 didn't jump to line 363 because the condition on line 362 was never true

363 children["photometric_scaling"] = self.photometric_scaling._describe(child) 

364 return Report( 

365 type_name="VisitImage", 

366 summary=summary, 

367 fields=fields, 

368 children=children, 

369 ) 

370 

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

372 """Deep-copy the visit image. 

373 

374 Parameters 

375 ---------- 

376 copy_detector 

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

378 """ 

379 return self._transfer_metadata( 

380 VisitImage( 

381 image=self._image.copy(), 

382 mask=self._mask.copy(), 

383 variance=self._variance.copy(), 

384 psf=self._psf, 

385 obs_info=self.obs_info, 

386 bounds=self._bounds, 

387 summary_stats=self.summary_stats.model_copy(), 

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

389 photometric_scaling=self._photometric_scaling, 

390 aperture_corrections=self.aperture_corrections.copy(), 

391 backgrounds=self._backgrounds.copy(), 

392 band=self.band, 

393 ), 

394 copy=True, 

395 ) 

396 

397 def convert_unit( 

398 self, 

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

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

401 copy_detector: bool = False, 

402 ) -> VisitImage: 

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

404 

405 Parameters 

406 ---------- 

407 unit 

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

409 

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

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

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

413 instrumental units but we know how to calibrate them) 

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

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

416 calibrated units and we want to revert back to instrumental 

417 units). 

418 copy 

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

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

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

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

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

424 already in the right units. 

425 copy_detector 

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

427 

428 Returns 

429 ------- 

430 `VisitImage` 

431 An image with the given units. 

432 """ 

433 if copy not in (True, False, "as-needed"): 433 ↛ 434line 433 didn't jump to line 434 because the condition on line 433 was never true

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

435 if (factor := _get_unit_conversion_factor(self.unit, unit)) is not None: 435 ↛ 436line 435 didn't jump to line 436 because the condition on line 435 was never true

436 if factor == 1.0: 

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

438 return self.copy() 

439 else: 

440 return self[...] 

441 elif copy is False: 

442 raise astropy.units.UnitConversionError( 

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

444 ) 

445 image = Image( 

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

447 ) 

448 variance = Image( 

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

450 bbox=self.bbox, 

451 unit=unit**2, 

452 ) 

453 elif self._photometric_scaling is None: 453 ↛ 454line 453 didn't jump to line 454 because the condition on line 453 was never true

454 raise astropy.units.UnitConversionError( 

455 "VisitImage.photometric_scaling is None, and there " 

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

457 ) 

458 else: 

459 if copy is False: 459 ↛ 460line 459 didn't jump to line 460 because the condition on line 459 was never true

460 raise astropy.units.UnitConversionError( 

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

462 ) 

463 scaling = self._photometric_scaling 

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

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

466 if constant_factor != 1.0: 466 ↛ 467line 466 didn't jump to line 467 because the condition on line 466 was never true

467 scaling = scaling * constant_factor 

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

469 elif (constant_factor := _get_unit_conversion_factor(self.unit / scaling.unit, unit)) is not None: 469 ↛ 475line 469 didn't jump to line 475 because the condition on line 469 was always true

470 if constant_factor != 1.0: 470 ↛ 471line 470 didn't jump to line 471 because the condition on line 470 was never true

471 scaling = scaling / constant_factor 

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

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

474 else: 

475 raise astropy.units.UnitConversionError( 

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

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

478 ) 

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

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

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

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

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

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

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

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

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

488 image.array *= self._image.array 

489 variance.array *= self._variance.array 

490 copy_components = copy is True 

491 return self._transfer_metadata( 

492 VisitImage( 

493 image=image, 

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

495 variance=variance, 

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

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

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

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

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

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

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

503 aperture_corrections=( 

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

505 ), 

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

507 band=self.band, 

508 ) 

509 ) 

510 

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

512 return self._serialize_impl(VisitImageSerializationModel, archive) 

513 

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

515 # VisitImageSerializationModel, not just a MaskedImageSerializationModel, 

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

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

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

519 ) -> M: 

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

521 match self._psf: 

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

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

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

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

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

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

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

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

530 case _: 

531 raise TypeError( 

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

533 ) 

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

535 result.obs_info = self.obs_info 

536 result.summary_stats = self.summary_stats 

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

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

539 result.band = self.band 

540 result.photometric_scaling = ( 

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

542 # method return types. 

543 archive.serialize_direct( 

544 "photometric_scaling", 

545 self._photometric_scaling.serialize, 

546 ) # type: ignore[assignment] 

547 if self._photometric_scaling is not None 

548 else None 

549 ) 

550 result.aperture_corrections = archive.serialize_direct( 

551 "aperture_corrections", 

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

553 ) 

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

555 return result 

556 

557 @staticmethod 

558 def _get_archive_tree_type[P: pydantic.BaseModel]( 

559 pointer_type: type[P], 

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

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

562 type that uses the given pointer type. 

563 """ 

564 return VisitImageSerializationModel[pointer_type] # type: ignore 

565 

566 @staticmethod 

567 def from_legacy( # type: ignore[override] 

568 legacy: LegacyExposure, 

569 *, 

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

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

572 instrument: str | None = None, 

573 visit: int | None = None, 

574 ) -> VisitImage: 

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

576 

577 Parameters 

578 ---------- 

579 legacy 

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

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

582 unit 

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

584 key will be used, if available. 

585 plane_map 

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

587 description. If `None` (default) 

588 `get_legacy_visit_image_mask_planes` is used. 

589 instrument 

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

591 provided. 

592 visit 

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

594 """ 

595 if plane_map is None: 

596 plane_map = get_legacy_visit_image_mask_planes() 

597 md = legacy.getMetadata() 

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

599 instrument = _extract_or_check_header( 

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

601 ) 

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

603 legacy_wcs = legacy.getWcs() 

604 if legacy_wcs is None: 

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

606 legacy_detector = legacy.getDetector() 

607 if legacy_detector is None: 

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

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

610 

611 # Update the ObservationInfo from other components. 

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

613 

614 opaque_fits_metadata = FitsOpaqueMetadata() 

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

616 with warnings.catch_warnings(): 

617 # Silence warnings about long keys becoming HIERARCH. 

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

619 primary_header.update(md.toOrderedDict()) 

620 metadata = opaque_fits_metadata.extract_legacy_primary_header(primary_header) 

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

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

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

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

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

626 # Fix incorrect BUNIT='adu' in LSST 

627 # preliminary_visit_image. 

628 hdr_unit = astropy.units.electron 

629 if unit is None: 

630 unit = hdr_unit 

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

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

633 sky_projection = SkyProjection.from_legacy( 

634 legacy_wcs, 

635 DetectorFrame( 

636 instrument=instrument, 

637 visit=visit, 

638 detector=legacy_detector.getId(), 

639 bbox=detector_bbox, 

640 ), 

641 ) 

642 legacy_psf = legacy.getPsf() 

643 if legacy_psf is None: 

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

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

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

647 legacy_summary_stats = legacy.info.getSummaryStats() 

648 legacy_ap_corr_map = legacy.info.getApCorrMap() 

649 legacy_polygon = legacy.info.getValidPolygon() 

650 legacy_photo_calib = legacy.info.getPhotoCalib() 

651 detector = Detector.from_legacy( 

652 legacy_detector, instrument=instrument, visit=visit, is_raw_assembled=True 

653 ) 

654 _reconcile_detector_serial(obs_info, detector) 

655 result = VisitImage( 

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

657 mask=masked_image.mask, 

658 variance=masked_image.variance, 

659 sky_projection=sky_projection, 

660 psf=psf, 

661 obs_info=obs_info, 

662 summary_stats=( 

663 ObservationSummaryStats.from_legacy(legacy_summary_stats) 

664 if legacy_summary_stats is not None 

665 else None 

666 ), 

667 detector=detector, 

668 aperture_corrections=( 

669 aperture_corrections_from_legacy(legacy_ap_corr_map) 

670 if legacy_ap_corr_map is not None 

671 else None 

672 ), 

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

674 photometric_scaling=( 

675 field_from_legacy_photo_calib( 

676 legacy_photo_calib, bounds=detector_bbox, instrumental_unit=instrumental_unit 

677 ) 

678 if legacy_photo_calib is not None 

679 else None 

680 ), 

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

682 metadata=metadata, 

683 ) 

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

685 result._opaque_metadata = opaque_fits_metadata 

686 return result 

687 

688 def to_legacy( 

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

690 ) -> LegacyExposure: 

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

692 

693 Parameters 

694 ---------- 

695 copy 

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

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

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

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

700 plane_map 

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

702 description. If `None` (default), 

703 `get_legacy_visit_image_mask_planes` is used. 

704 """ 

705 from lsst.afw.image import Exposure as LegacyExposure 

706 from lsst.afw.image import FilterLabel as LegacyFilterLabel 

707 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo 

708 

709 if plane_map is None: 

710 plane_map = get_legacy_visit_image_mask_planes() 

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

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

713 result_info = result.info 

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

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

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

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

718 if self._photometric_scaling is not None: 

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

720 else: 

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

722 self._fill_legacy_metadata(result_info.getMetadata()) 

723 if isinstance(self._psf, LegacyPointSpreadFunction): 

724 result_info.setPsf(self._psf.legacy_psf) 

725 elif isinstance(self._psf, PiffWrapper): 

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

727 if isinstance(self.bounds, Polygon): 

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

729 if self.aperture_corrections: 

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

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

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

733 return result 

734 

735 @staticmethod 

736 def read_legacy( # type: ignore[override] 

737 filename: str, 

738 *, 

739 preserve_quantization: bool = False, 

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

741 instrument: str | None = None, 

742 visit: int | None = None, 

743 component: Literal[ 

744 "bbox", 

745 "image", 

746 "mask", 

747 "variance", 

748 "sky_projection", 

749 "psf", 

750 "detector", 

751 "photometric_scaling", 

752 "obs_info", 

753 "summary_stats", 

754 "aperture_corrections", 

755 ] 

756 | None = None, 

757 ) -> Any: 

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

759 

760 Parameters 

761 ---------- 

762 filename 

763 Full name of the file. 

764 preserve_quantization 

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

766 exactly preserve quantization-compressed pixel values. This causes 

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

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

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

770 not transferred to the copy. 

771 plane_map 

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

773 description. If `None` (default) 

774 `get_legacy_visit_image_mask_planes` is used. 

775 instrument 

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

777 provided. 

778 visit 

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

780 provided. 

781 component 

782 A component to read instead of the full image. 

783 """ 

784 from lsst.afw.image import ExposureFitsReader 

785 

786 reader = ExposureFitsReader(filename) 

787 if component == "bbox": 

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

789 legacy_detector = reader.readDetector() 

790 if legacy_detector is None: 

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

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

793 legacy_wcs = None 

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

795 legacy_wcs = reader.readWcs() 

796 if legacy_wcs is None: 

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

798 legacy_exposure_info = reader.readExposureInfo() 

799 summary_stats = None 

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

801 legacy_stats = legacy_exposure_info.getSummaryStats() 

802 if legacy_stats is not None: 

803 summary_stats = ObservationSummaryStats.from_legacy(legacy_stats) 

804 if component == "summary_stats": 

805 return summary_stats