Coverage for python/lsst/images/_visit

1# This file is part of lsst-images.

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.

9# Use of this source code is governed by a 3-clause BSD-style

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

12from __future__ import annotations

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

16import functools

17import warnings

18from collections.abc import Callable, Mapping, MutableMapping

19from types import EllipsisType

20from typing import Any, Literal, cast, overload

22import astropy.io.fits

23import astropy.units

24import astropy.wcs

25import pydantic

26from astro_metadata_translator import ObservationInfo, VisitInfoTranslator

28from ._geom import Box

29from ._image import Image, ImageSerializationModel

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

31from ._masked_image import MaskedImage, MaskedImageSerializationModel

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

33from .fits import FitsOpaqueMetadata

34from .psfs import (

35 GaussianPointSpreadFunction,

36 GaussianPSFSerializationModel,

37 PiffSerializationModel,

38 PiffWrapper,

39 PointSpreadFunction,

40 PSFExSerializationModel,

41 PSFExWrapper,

42)

43from .serialization import ArchiveReadError, InputArchive, MetadataValue, OutputArchive

46def _obs_info_from_md(md: MutableMapping[str, Any], visit_info: Any = None) -> ObservationInfo:

47 # Try to get an ObservationInfo from the primary header as if

48 # it's a raw header. Else fallback.

49 try:

50 obs_info = ObservationInfo.from_header(md, quiet=True)

51 except ValueError:

52 # Not known translator. Must fall back to visit info. If we have

53 # an actual VisitInfo, serialize it since we know that it will be

54 # complete.

55 if visit_info is not None:

56 from lsst.afw.image import setVisitInfoMetadata

57 from lsst.daf.base import PropertyList

59 pl = PropertyList()

60 setVisitInfoMetadata(pl, visit_info)

61 # Merge so that we still have access to butler provenance.

62 md.update(pl)

64 # Try the given header looking for VisitInfo hints.

65 # We get lots of warnings if nothing can be found. Currently

66 # no way to disable those without capturing them.

67 obs_info = ObservationInfo.from_header(md, translator_class=VisitInfoTranslator, quiet=True)

68 return obs_info

71class VisitImage(MaskedImage):

72 """A calibrated single-visit image.

74 Parameters

75 ----------

76 image

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

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

79 mask

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

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

82 is replaced (possibly by `None`).

83 variance

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

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

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

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

88 (possibly by `None`).

89 mask_schema

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

91 not provided.

92 opaque_metadata

93 Opaque metadata obtained from reading this object from storage. It may

94 be provided when writing to storage to propagate that metadata and/or

95 preserve file-format-specific options (e.g. compression parameters).

96 projection

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

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

99 psf

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

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

102 obs_info

103 General information about this visit in standardized form.

104 metadata

105 Arbitrary flexible metadata to associate with the image.

106 """

107

108 def __init__(

109 self,

110 image: Image,

111 *,

112 mask: Mask | None = None,

113 variance: Image | None = None,

114 mask_schema: MaskSchema | None = None,

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

116 obs_info: ObservationInfo | None = None,

117 psf: PointSpreadFunction | ArchiveReadError,

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

119 ):

120 super().__init__(

121 image,

122 mask=mask,

123 variance=variance,

124 mask_schema=mask_schema,

125 projection=projection,

126 obs_info=obs_info,

127 metadata=metadata,

128 )

129 if self.image.unit is None:

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

131 if self.image.projection is None:

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

133 if self.image.obs_info is None:

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

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

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

137 self._psf = psf

138

139 @property

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

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

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

143

144 @property

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

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

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

148 """

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

150

151 @property

152 def obs_info(self) -> ObservationInfo:

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

154 (`~astro_metadata_translator.ObservationInfo`).

155 """

156 obs_info = self.image.obs_info

157 assert obs_info is not None

158 return obs_info

159

160 @property

161 def astropy_wcs(self) -> ProjectionAstropyView:

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

163

164 Notes

165 -----

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

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

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

169

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

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

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

173 """

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

175

176 @property

177 def psf(self) -> PointSpreadFunction:

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

179 (`.psfs.PointSpreadFunction`).

180 """

181 if isinstance(self._psf, ArchiveReadError):

182 raise self._psf

183 return self._psf

184

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

186 if bbox is ...:

187 return self

188 super().__getitem__(bbox)

189 return self._transfer_metadata(

190 VisitImage(

191 self.image[bbox],

192 mask=self.mask[bbox],

193 variance=self.variance[bbox],

194 projection=self.projection,

195 psf=self.psf,

196 obs_info=self.obs_info,

197 ),

198 bbox=bbox,

199 )

200

201 def __str__(self) -> str:

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

203

204 def __repr__(self) -> str:

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

206

207 def copy(self) -> VisitImage:

208 """Deep-copy the visit image."""

209 return self._transfer_metadata(

210 VisitImage(

211 image=self._image.copy(),

212 mask=self._mask.copy(),

213 variance=self._variance.copy(),

214 psf=self._psf,

215 obs_info=self.obs_info,

216 ),

217 copy=True,

218 )

219

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

221 serialized_image = archive.serialize_direct(

222 "image", functools.partial(self.image.serialize, save_projection=False)

223 )

224 serialized_mask = archive.serialize_direct(

225 "mask", functools.partial(self.mask.serialize, save_projection=False)

226 )

227 serialized_variance = archive.serialize_direct(

228 "variance", functools.partial(self.variance.serialize, save_projection=False)

229 )

230 serialized_projection = (

231 archive.serialize_direct("projection", self.projection.serialize)

232 if self.projection is not None

233 else None

234 )

235 serialized_psf: PiffSerializationModel | PSFExSerializationModel | GaussianPSFSerializationModel

236 match self._psf:

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

238 # but not a single isinstance check on both types.

239 case PiffWrapper():

240 serialized_psf = archive.serialize_direct("psf", self._psf.serialize)

241 case PSFExWrapper():

242 serialized_psf = archive.serialize_direct("psf", self._psf.serialize)

243 case GaussianPointSpreadFunction():

244 serialized_psf = archive.serialize_direct("psf", self._psf.serialize)

245 case _:

246 raise TypeError(

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

248 )

249 return VisitImageSerializationModel(

250 image=serialized_image,

251 mask=serialized_mask,

252 variance=serialized_variance,

253 projection=serialized_projection,

254 psf=serialized_psf,

255 obs_info=self.obs_info,

256 metadata=self.metadata,

257 )

258

259 # Type-checkers want the model argument to only require

260 # MaskedImageSerializationModel[Any], and they'd be absolutely right if

261 # this were a regular instance method. But whether Liskov substitution

262 # applies to classmethods and staticmethods is sort of context-dependent,

263 # and here we do not want it to.

264 @staticmethod

265 def deserialize(

266 model: VisitImageSerializationModel[Any], # type: ignore[override]

267 archive: InputArchive[Any],

268 *,

269 bbox: Box | None = None,

270 ) -> VisitImage:

271 masked_image = MaskedImage.deserialize(model, archive, bbox=bbox)

272 psf = model.deserialize_psf(archive)

273 projection = Projection.deserialize(model.projection, archive)

274 return VisitImage(

275 masked_image.image,

276 mask=masked_image.mask,

277 variance=masked_image.variance,

278 psf=psf,

279 projection=projection,

280 obs_info=model.obs_info,

281 )._finish_deserialize(model)

282

283 @staticmethod

284 def _get_archive_tree_type[P: pydantic.BaseModel](

285 pointer_type: type[P],

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

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

288 type that uses the given pointer type.

289 """

290 return VisitImageSerializationModel[pointer_type] # type: ignore

291

292 # write_fits and read_fits inherited from MaskedImage.

293

294 @staticmethod

295 def from_legacy(

296 legacy: Any,

297 *,

298 unit: astropy.units.Unit | None = None,

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

300 instrument: str | None = None,

301 visit: int | None = None,

302 component: str | None = None,

303 ) -> VisitImage:

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

305

306 Parameters

307 ----------

308 legacy

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

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

311 unit

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

313 key will be used, if available.

314 plane_map

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

316 description. If `None` (default)

317 `get_legacy_visit_image_mask_planes` is used.

318 instrument

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

320 provided.

321 visit

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

323 """

324 if plane_map is None:

325 plane_map = get_legacy_visit_image_mask_planes()

326 md = legacy.getMetadata()

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

328 instrument = _extract_or_check_header(

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

330 )

331 visit = _extract_or_check_header("LSST BUTLER DATAID VISIT", visit, md, None, int)

332 unit = _extract_or_check_header(

333 "BUNIT", unit, md, None, lambda x: astropy.units.Unit(x, format="fits")

334 )

335 legacy_wcs = legacy.getWcs()

336 if legacy_wcs is None:

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

338 legacy_detector = legacy.getDetector()

339 if legacy_detector is None:

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

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

342 opaque_fits_metadata = FitsOpaqueMetadata()

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

344 with warnings.catch_warnings():

345 # Silence warnings about long keys becoming HIERARCH.

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

347 primary_header.update(md.toOrderedDict())

348 opaque_fits_metadata.extract_legacy_primary_header(primary_header)

349 projection = Projection.from_legacy(

350 legacy_wcs,

351 DetectorFrame(

352 instrument=instrument,

353 visit=visit,

354 detector=legacy_detector.getId(),

355 bbox=detector_bbox,

356 ),

357 )

358 legacy_psf = legacy.getPsf()

359 if legacy_psf is None:

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

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

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

363 result = VisitImage(

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

365 mask=masked_image.mask,

366 variance=masked_image.variance,

367 projection=projection,

368 psf=psf,

369 obs_info=obs_info,

370 )

371

372 result._opaque_metadata = opaque_fits_metadata

373 return result

374

375 @overload # type: ignore[override]

376 @staticmethod

377 def read_legacy( 377 ↛ exitline 377 didn't return from function 'read_legacy' because

378 filename: str,

379 *,

380 component: Literal["bbox"],

381 ) -> Box: ...

382

383 @overload

384 @staticmethod

385 def read_legacy( 385 ↛ exitline 385 didn't return from function 'read_legacy' because

386 filename: str,

387 *,

388 preserve_quantization: bool = False,

389 instrument: str | None = None,

390 visit: int | None = None,

391 component: Literal["image"],

392 ) -> Image: ...

393

394 @overload

395 @staticmethod

396 def read_legacy( 396 ↛ exitline 396 didn't return from function 'read_legacy' because

397 filename: str,

398 *,

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

400 instrument: str | None = None,

401 visit: int | None = None,

402 component: Literal["mask"],

403 ) -> Mask: ...

404

405 @overload

406 @staticmethod

407 def read_legacy( 407 ↛ exitline 407 didn't return from function 'read_legacy' because

408 filename: str,

409 *,

410 preserve_quantization: bool = False,

411 instrument: str | None = None,

412 visit: int | None = None,

413 component: Literal["variance"],

414 ) -> Image: ...

415

416 @overload

417 @staticmethod

418 def read_legacy( 418 ↛ exitline 418 didn't return from function 'read_legacy' because

419 filename: str,

420 *,

421 instrument: str | None = None,

422 visit: int | None = None,

423 component: Literal["projection"],

424 ) -> Projection[DetectorFrame]: ...

425

426 @overload

427 @staticmethod

428 def read_legacy( 428 ↛ exitline 428 didn't return from function 'read_legacy' because

429 filename: str,

430 *,

431 component: Literal["psf"],

432 ) -> PointSpreadFunction: ...

433

434 @overload

435 @staticmethod

436 def read_legacy( 436 ↛ exitline 436 didn't return from function 'read_legacy' because

437 filename: str,

438 *,

439 component: Literal["obs_info"],

440 ) -> ObservationInfo: ...

441

442 @overload

443 @staticmethod

444 def read_legacy( 444 ↛ exitline 444 didn't return from function 'read_legacy' because

445 filename: str,

446 *,

447 preserve_quantization: bool = False,

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

449 instrument: str | None = None,

450 visit: int | None = None,

451 component: None = None,

452 ) -> VisitImage: ...

453

454 @staticmethod

455 def read_legacy( # type: ignore[override]

456 filename: str,

457 *,

458 preserve_quantization: bool = False,

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

460 instrument: str | None = None,

461 visit: int | None = None,

462 component: Literal["bbox", "image", "mask", "variance", "projection", "psf", "obs_info"]

463 | None = None,

464 ) -> Any:

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

466

467 Parameters

468 ----------

469 filename

470 Full name of the file.

471 preserve_quantization

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

473 exactly preserve quantization-compressed pixel values. This causes

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

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

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

477 not transferred to the copy.

478 plane_map

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

480 description. If `None` (default)

481 `get_legacy_visit_image_mask_planes` is used.

482 instrument

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

484 provided.

485 visit

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

487 provided.

488 component

489 A component to read instead of the full image.

490 """

491 from lsst.afw.image import ExposureFitsReader

492

493 reader = ExposureFitsReader(filename)

494 if component == "bbox":

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

496 legacy_detector = reader.readDetector()

497 if legacy_detector is None:

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

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

500 legacy_wcs = None

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

502 legacy_wcs = reader.readWcs()

503 if legacy_wcs is None:

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

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

506 legacy_psf = reader.readPsf()

507 if legacy_psf is None:

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

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

510 if component == "psf":

511 return psf

512 assert component in (None, "image", "mask", "variance", "projection", "obs_info"), (

513 component

514 ) # for MyPy

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

516 primary_header = hdu_list[0].header

517 obs_info = _obs_info_from_md(primary_header)

518 if component == "obs_info":

519 return obs_info

520 instrument = _extract_or_check_header(

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

522 )

523 visit = _extract_or_check_header("LSST BUTLER DATAID VISIT", visit, primary_header, None, int)

524 projection = Projection.from_legacy(

525 legacy_wcs,

526 DetectorFrame(

527 instrument=instrument,

528 visit=visit,

529 detector=legacy_detector.getId(),

530 bbox=detector_bbox,

531 ),

532 )

533 if component == "projection":

534 return projection

535 if plane_map is None:

536 plane_map = get_legacy_visit_image_mask_planes()

537 assert component != "psf", component # for MyPy

538 from_masked_image = MaskedImage._read_legacy_hdus(

539 hdu_list,

540 filename,

541 preserve_quantization=preserve_quantization,

542 plane_map=plane_map,

543 component=component,

544 )

545 if component is not None:

546 # This is the image, mask, or variance; attach the projection

547 # and return

548 return from_masked_image.view(projection=projection)

549 result = VisitImage(

550 from_masked_image.image,

551 mask=from_masked_image.mask,

552 variance=from_masked_image.variance,

553 projection=projection,

554 psf=psf,

555 obs_info=obs_info,

556 )

557 result._opaque_metadata = from_masked_image._opaque_metadata

558 return result

559

560

561class VisitImageSerializationModel[P: pydantic.BaseModel](MaskedImageSerializationModel[P]):

562 """A Pydantic model used to represent a serialized `VisitImage`."""

563

564 # Inherited attributes are duplicated because that improves the docs

565 # (some limitation in the sphinx/pydantic integration), and these are

566 # important docs.

567

568 image: ImageSerializationModel[P] = pydantic.Field(description="The main data image.")

569 mask: MaskSerializationModel[P] = pydantic.Field(

570 description="Bitmask that annotates the main image's pixels."

571 )

572 variance: ImageSerializationModel[P] = pydantic.Field(

573 description="Per-pixel variance estimates for the main image."

574 )

575 projection: ProjectionSerializationModel[P] = pydantic.Field(

576 description="Projection that maps the pixel grid to the sky.",

577 )

578 psf: PiffSerializationModel | PSFExSerializationModel | GaussianPSFSerializationModel | Any = (

579 pydantic.Field(union_mode="left_to_right", description="PSF model for the image.")

580 )

581 obs_info: ObservationInfo = pydantic.Field(

582 description="Standardized description of visit metadata",

583 )

584

585 def deserialize_psf(self, archive: InputArchive[Any]) -> PointSpreadFunction | ArchiveReadError:

586 """Finish deserializing the PSF model, or *return* any exception

587 raised in the attempt.

588 """

589 try:

590 match self.psf:

591 case PiffSerializationModel():

592 return PiffWrapper.deserialize(self.psf, archive)

593 case PSFExSerializationModel():

594 return PSFExWrapper.deserialize(self.psf, archive)

595 case GaussianPSFSerializationModel():

596 return GaussianPointSpreadFunction.deserialize(self.psf, archive)

597 case _:

598 raise ArchiveReadError("PSF model type not recognized.")

599 except ArchiveReadError as err:

600 return err

601

602

603def _extract_or_check_value[T](

604 key: str,

605 given_value: T | None,

606 *sources: tuple[str, T | None],

607) -> T:

608 # Compare given value against multiple sources. If given value is not

609 # supplied return the first non-None value in the reference sources.

610 if given_value is not None:

611 for source_name, source_value in sources:

612 if source_value is not None and source_value != given_value:

613 raise ValueError(

614 f"Given value {given_value!r} does not match {source_value!r} from {source_name}."

615 )

616 if source_value is not None:

617 # Only check the first non-None source rather than checking

618 # all supplied values.

619 break

620 return given_value

621

622 for _, source_value in sources:

623 if source_value is not None:

624 return source_value

625

626 raise ValueError(f"No value found for {key}.")

627

628

629def _extract_or_check_header[T](

630 key: str, given_value: T | None, header: Any, obs_info_value: T | None, coerce: Callable[[Any], T]

631) -> T:

632 hdr_value: T | None = None

633 if (hdr_raw_value := header.get(key)) is not None:

634 hdr_value = coerce(hdr_raw_value)

635 return _extract_or_check_value(

636 key, given_value, ("ObservationInfo", obs_info_value), (f"header key {key}", hdr_value)

637 )

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

232 statements