Coverage for tests/test_visit_image.py: 71%
745 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-30 04:21 -0700
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.
12from __future__ import annotations
14import dataclasses
15import math
16import os
17import warnings
18from pathlib import Path
19from typing import Any, Literal
21import astropy.io.fits
22import astropy.units as u
23import astropy.wcs
24import numpy as np
25import pytest
26from astro_metadata_translator import ObservationInfo
28from lsst.images import (
29 Background,
30 BackgroundMap,
31 Box,
32 DetectorFrame,
33 DifferenceImage,
34 Image,
35 MaskPlane,
36 MaskSchema,
37 ObservationSummaryStats,
38 Polygon,
39 SkyProjectionAstropyView,
40 TractFrame,
41 VisitImage,
42 get_legacy_difference_image_mask_planes,
43 get_legacy_visit_image_mask_planes,
44)
45from lsst.images.aperture_corrections import ApertureCorrectionMap, aperture_corrections_to_legacy
46from lsst.images.cameras import Detector
47from lsst.images.describe import DescribableMixin, DescribeOptions, FieldRole, Report
48from lsst.images.fields import ChebyshevField, SplineField, SumField, field_from_legacy_photo_calib
49from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata
50from lsst.images.psfs import GaussianPointSpreadFunction, PointSpreadFunction
51from lsst.images.serialization import ArchiveReadError, read_archive
52from lsst.images.tests import (
53 DP2_VISIT_DETECTOR_DATA_ID,
54 RoundtripFits,
55 RoundtripJson,
56 RoundtripNdf,
57 TemporaryButler,
58 assert_masked_images_equal,
59 assert_sky_projections_equal,
60 assert_values_equal,
61 assert_visit_images_equal,
62 compare_aperture_corrections_to_legacy,
63 compare_detector_to_legacy,
64 compare_photo_calib_to_legacy,
65 compare_visit_image_to_legacy,
66 current_fixture_path,
67 make_random_sky_projection,
68 reset_afw_mask_planes, # noqa: F401
69)
71try:
72 import h5py # noqa: F401
74 HAVE_H5PY = True
75except ImportError:
76 HAVE_H5PY = False
78try:
79 from lsst.afw.image import Exposure as LegacyExposure
80 from lsst.afw.image import VisitInfo as LegacyVisitInfo
81except ImportError:
82 type LegacyExposure = Any # type: ignore[no-redef]
83 type LegacyVisitInfo = Any # type: ignore[no-redef]
85EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
86LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
87FIXTURE_DIR = Path(__file__).parent / "data" / "schemas"
88MINIMAL_VISIT_DATA_ID = {key: DP2_VISIT_DETECTOR_DATA_ID[key] for key in ("instrument", "visit", "detector")}
90skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed")
93@pytest.fixture(scope="session")
94def visit_image_components() -> dict[str, Any]:
95 """Return a dictionary of VisitImage components."""
96 rng = np.random.default_rng(500)
97 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096])
98 mask_schema = MaskSchema([MaskPlane("M1", "D1")])
99 obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1")
100 summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
101 gaussian_psf = GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13])
102 aperture_corrections: ApertureCorrectionMap = {
103 "flux1": ChebyshevField(det_frame.bbox, np.array([0.75])),
104 "flux2": ChebyshevField(det_frame.bbox, np.array([0.625])),
105 }
106 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector)
107 # Real visit images have float pixels, and some operations (e.g. rendering
108 # a photometric scaling) are only defined for floating-point images.
109 image = Image(42.0, shape=(1024, 1024), unit=u.nJy, dtype=np.float32)
110 variance = Image(5.0, shape=(1024, 1024), unit=u.nJy * u.nJy, dtype=np.float32)
111 # polygon is the lower triangle of the image.
112 polygon = Polygon(x_vertices=[-0.5, 1023.5, -0.5], y_vertices=[-0.5, -0.5, 1023.5])
113 sky_projection = make_random_sky_projection(rng, det_frame, det_frame.bbox)
114 return {
115 "mask_schema": mask_schema,
116 "obs_info": obs_info,
117 "summary_stats": summary_stats,
118 "gaussian_psf": gaussian_psf,
119 "aperture_corrections": aperture_corrections,
120 "detector": detector,
121 "image": image,
122 "variance": variance,
123 "polygon": polygon,
124 "sky_projection": sky_projection,
125 }
128def make_visit_image(components: dict[str, Any]) -> VisitImage:
129 """Construct a new VisitImage with most components populated."""
130 det_frame = components["sky_projection"].pixel_frame
131 opaque = FitsOpaqueMetadata()
132 hdr = astropy.io.fits.Header()
133 with warnings.catch_warnings():
134 # Silence warnings about long keys becoming HIERARCH.
135 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning)
136 hdr.update({"PLATFORM": "lsstcam", "LSST BUTLER ID": "123456789"})
137 opaque.extract_legacy_primary_header(hdr)
138 # API signature suggests sky_projection and obs_info can be None but
139 # they are required (unless you pass them in via the image plane).
140 vi = VisitImage(
141 components["image"],
142 variance=components["variance"],
143 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
144 mask_schema=components["mask_schema"],
145 sky_projection=components["sky_projection"],
146 obs_info=components["obs_info"],
147 summary_stats=components["summary_stats"],
148 detector=components["detector"],
149 bounds=components["polygon"],
150 aperture_corrections=components["aperture_corrections"],
151 band="r",
152 )
153 vi.backgrounds.add(
154 "standard",
155 ChebyshevField(det_frame.bbox, np.array([[2.0]])),
156 description="Background subtracted from the image.",
157 is_subtracted=True,
158 )
159 vi._opaque_metadata = opaque
160 return vi
163def make_simplest_visit_image(components: dict[str, Any]) -> VisitImage:
164 """Construct a VisitImage with the minimal set of components populated."""
165 return VisitImage(
166 components["image"],
167 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
168 mask_schema=components["mask_schema"],
169 sky_projection=components["sky_projection"],
170 detector=components["detector"],
171 obs_info=components["obs_info"],
172 band="r",
173 )
176def _make_sum_background_visit_image(components: dict[str, Any], visit_image: VisitImage) -> VisitImage:
177 """Return a VisitImage whose subtracted background is a SumField."""
178 rng = np.random.default_rng(42)
179 bbox = visit_image.image.sky_projection.pixel_frame.bbox
180 bin_y = bbox.y.linspace(6)
181 bin_x = bbox.x.linspace(7)
182 spline_a = SplineField(
183 bbox,
184 rng.standard_normal(size=(bin_y.size, bin_x.size)),
185 y=bin_y,
186 x=bin_x,
187 )
188 spline_b = SplineField(
189 bbox,
190 rng.standard_normal(size=(bin_y.size, bin_x.size)),
191 y=bin_y,
192 x=bin_x,
193 )
194 sum_field = SumField([spline_a, spline_b])
195 bg_map = BackgroundMap()
196 bg_map.add(
197 "stacked",
198 sum_field,
199 description="Two-operand SumField subtracted background.",
200 is_subtracted=True,
201 )
202 return VisitImage(
203 components["image"],
204 variance=components["variance"],
205 psf=components["gaussian_psf"],
206 mask_schema=components["mask_schema"],
207 sky_projection=components["sky_projection"],
208 obs_info=components["obs_info"],
209 summary_stats=components["summary_stats"],
210 detector=components["detector"],
211 band="r",
212 backgrounds=bg_map,
213 )
216def _check_sum_background_round_trip(result: VisitImage, original: VisitImage) -> None:
217 """Assert that a round-tripped SumField background matches the original."""
218 subtracted = result.backgrounds.subtracted
219 assert subtracted is not None
220 assert isinstance(subtracted.field, SumField)
221 original_subtracted = original.backgrounds.subtracted
222 assert original_subtracted is not None
223 original_field = original_subtracted.field
224 assert isinstance(original_field, SumField)
225 round_field = subtracted.field
226 assert isinstance(round_field, SumField)
227 assert len(round_field.operands) == len(original_field.operands)
228 for round_op, orig_op in zip(round_field.operands, original_field.operands, strict=True):
229 assert round_op == orig_op
232def test_visit_image_repr_str_pinned(visit_image_components: dict[str, Any]) -> None:
233 """Pin the exact str and repr output of a VisitImage."""
234 visit = make_simplest_visit_image(visit_image_components)
235 assert str(visit) == "VisitImage(Image([y=0:1024, x=0:1024], float32), ['M1'])"
236 assert repr(visit) == (
237 "VisitImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
238 " dtype=dtype('float32')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
239 " dtype=dtype('uint8')))"
240 )
243def test_basics(visit_image_components: dict[str, Any]) -> None:
244 """Verify VisitImage constructor patterns and required-argument checks."""
245 c = visit_image_components
246 # Test default fill of variance.
247 visit = make_simplest_visit_image(c)
248 assert visit.variance.array[0, 0] == 1.0
249 assert visit[...] is not visit
250 assert str(visit) == "VisitImage(Image([y=0:1024, x=0:1024], float32), ['M1'])"
251 assert repr(visit) == (
252 "VisitImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
253 " dtype=dtype('float32')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
254 " dtype=dtype('uint8')))"
255 )
257 astropy_wcs = visit.astropy_wcs
258 assert isinstance(astropy_wcs, SkyProjectionAstropyView)
259 approx_wcs = visit.fits_wcs
260 assert isinstance(approx_wcs, astropy.wcs.WCS)
262 with pytest.raises(TypeError):
263 # Requires a PSF.
264 VisitImage(
265 c["image"],
266 mask_schema=c["mask_schema"],
267 sky_projection=c["sky_projection"],
268 obs_info=c["obs_info"],
269 detector=c["detector"],
270 band="r",
271 )
273 with pytest.raises(TypeError):
274 # Requires ObservationInfo.
275 VisitImage(
276 c["image"],
277 psf=c["gaussian_psf"],
278 mask_schema=c["mask_schema"],
279 sky_projection=c["sky_projection"],
280 detector=c["detector"],
281 band="r",
282 )
284 with pytest.raises(TypeError):
285 # Requires a sky_projection.
286 VisitImage(
287 c["image"],
288 psf=c["gaussian_psf"],
289 mask_schema=c["mask_schema"],
290 obs_info=c["obs_info"],
291 detector=c["detector"],
292 band="r",
293 )
295 with pytest.raises(TypeError):
296 # Requires a detector.
297 VisitImage(
298 c["image"],
299 psf=c["gaussian_psf"],
300 mask_schema=c["mask_schema"],
301 sky_projection=c["sky_projection"],
302 obs_info=c["obs_info"],
303 band="r",
304 )
306 with pytest.raises(TypeError):
307 # Requires some form of mask.
308 VisitImage(
309 c["image"],
310 psf=c["gaussian_psf"],
311 sky_projection=c["sky_projection"],
312 obs_info=c["obs_info"],
313 detector=c["detector"],
314 band="r",
315 )
317 with pytest.raises(TypeError):
318 VisitImage(
319 Image(42, shape=(5, 5)),
320 psf=c["gaussian_psf"],
321 mask_schema=c["mask_schema"],
322 sky_projection=c["sky_projection"],
323 obs_info=c["obs_info"],
324 detector=c["detector"],
325 band="r",
326 )
328 # Requires a DetectorFrame.
329 rng = np.random.default_rng(501)
330 tract_frame = TractFrame(skymap="Skymap", tract=1, bbox=Box.factory[1:10, 1:10])
331 tract_proj = make_random_sky_projection(rng, tract_frame, Box.factory[1:4096, 1:4096])
332 with pytest.raises(TypeError):
333 VisitImage(
334 c["image"],
335 sky_projection=tract_proj,
336 psf=c["gaussian_psf"],
337 mask_schema=c["mask_schema"],
338 obs_info=c["obs_info"],
339 detector=c["detector"],
340 band="r",
341 )
343 # Variance unit mismatch.
344 with pytest.raises(ValueError, match="should be the square of the image unit"):
345 VisitImage(
346 c["image"],
347 variance=c["image"],
348 psf=c["gaussian_psf"],
349 mask_schema=c["mask_schema"],
350 sky_projection=c["sky_projection"],
351 obs_info=c["obs_info"],
352 detector=c["detector"],
353 band="r",
354 )
357def test_copy_and_slice(visit_image_components: dict[str, Any]) -> None:
358 """Verify that copy deep-copies arrays and components while slice shares
359 them.
360 """
361 c = visit_image_components
362 visit_image = make_visit_image(c)
363 copy = visit_image.copy()
364 copy.image.array[0, 0] = 30.0
365 assert visit_image.image.array[0, 0] == 42.0
366 assert copy.image.array[0, 0] == 30.0
367 subvisit = visit_image[Box.factory[0:5, 0:5]]
368 # Check summary stats.
369 assert copy.summary_stats == visit_image.summary_stats
370 assert copy.summary_stats is not visit_image.summary_stats
371 assert subvisit.summary_stats == visit_image.summary_stats
372 assert subvisit.summary_stats is visit_image.summary_stats
373 # Check aperture corrections.
374 assert copy.aperture_corrections.keys() == visit_image.aperture_corrections.keys()
375 assert copy.aperture_corrections is not visit_image.aperture_corrections
376 assert subvisit.aperture_corrections.keys() == visit_image.aperture_corrections.keys()
377 assert subvisit.aperture_corrections is visit_image.aperture_corrections
378 # Check backgrounds.
379 assert copy.backgrounds.keys() == visit_image.backgrounds.keys()
380 assert copy.backgrounds is not visit_image.backgrounds
381 assert subvisit.backgrounds.keys() == visit_image.backgrounds.keys()
382 assert subvisit.backgrounds is visit_image.backgrounds
383 # Check bounds.
384 assert copy.bounds is c["polygon"]
385 assert subvisit.bounds == subvisit.bbox # original polygon wholly encloses subvisit.bbox
388def test_obs_info(visit_image_components: dict[str, Any]) -> None:
389 """Verify that ObservationInfo is present and carries the expected
390 instrument.
391 """
392 visit_image = make_visit_image(visit_image_components)
393 assert visit_image.obs_info is not None
394 assert visit_image.obs_info.instrument == "LSSTCam"
397def test_summary_stats(visit_image_components: dict[str, Any]) -> None:
398 """Verify ObservationSummaryStats equality and inequality comparisons."""
399 summary_stats = visit_image_components["summary_stats"]
400 assert summary_stats == ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
401 assert summary_stats != ObservationSummaryStats(psfSigma=2.5)
402 assert summary_stats != ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
405def test_summary_stats_to_legacy() -> None:
406 """Verify ObservationSummaryStats round-trips through the legacy
407 ExposureSummaryStats even when this package defines fields that the
408 installed afw does not.
409 """
410 try:
411 from lsst.afw.image import ExposureSummaryStats
412 except ImportError:
413 pytest.skip("lsst.afw.image is not available")
415 summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
416 legacy = summary_stats.to_legacy()
417 assert isinstance(legacy, ExposureSummaryStats)
418 assert legacy.psfSigma == 2.5
419 assert legacy.zeroPoint == 31.4
420 # Empty (NaN) fields unknown to the legacy struct are dropped, so the
421 # round trip reproduces the original.
422 assert ObservationSummaryStats.from_legacy(legacy) == summary_stats
424 # A real value in a field unknown to the legacy struct cannot be
425 # represented and must raise rather than be silently dropped.
426 legacy_fields = {field.name for field in dataclasses.fields(ExposureSummaryStats)}
427 extra_float_fields = [
428 name
429 for name, info in ObservationSummaryStats.model_fields.items()
430 if name not in legacy_fields and info.annotation is float
431 ]
432 if extra_float_fields: 432 ↛ 433line 432 didn't jump to line 433 because the condition on line 432 was never true
433 with pytest.raises(ValueError, match="not supported by this version"):
434 ObservationSummaryStats(**{extra_float_fields[0]: 1.0}).to_legacy()
437def test_summary_stats_from_legacy_unknown_field() -> None:
438 """Verify from_legacy drops empty unknown fields but errors on set ones."""
440 @dataclasses.dataclass
441 class FakeLegacy:
442 psfSigma: float = 2.5
443 notARealField: float = math.nan
445 # An unknown field that is empty is dropped.
446 stats = ObservationSummaryStats.from_legacy(FakeLegacy())
447 assert stats.psfSigma == 2.5
449 # An unknown field that holds a real value cannot be represented.
450 with pytest.raises(ValueError, match="is not known to ObservationSummaryStats"):
451 ObservationSummaryStats.from_legacy(FakeLegacy(notARealField=1.0))
454def test_summary_stats_legacy_butler_read(reset_afw_mask_planes: None) -> None: # noqa: F811
455 """Test that a round-tripped ObservationSummaryStats can be read back as
456 a legacy afw ExposureSummaryStats via the converter registered on the
457 ExposureSummaryStats storage class.
458 """
459 from lsst.afw.image import ExposureSummaryStats as LegacyExposureSummaryStats
461 summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
462 with RoundtripFits(summary_stats, storage_class="ObservationSummaryStats") as roundtrip:
463 legacy_stats = roundtrip.get(storageClass="ExposureSummaryStats")
464 assert isinstance(legacy_stats, LegacyExposureSummaryStats)
465 assert legacy_stats.psfSigma == 2.5
466 assert legacy_stats.zeroPoint == 31.4
467 assert ObservationSummaryStats.from_legacy(legacy_stats) == summary_stats
470def test_detector_legacy_butler_read(reset_afw_mask_planes: None) -> None: # noqa: F811
471 """Test that a round-tripped Detector (written under the DetectorV2
472 storage class) can be read back as a legacy lsst.afw.cameraGeom.Detector
473 via the converter registered on the Detector storage class.
474 """
475 from lsst.afw.cameraGeom import Detector as LegacyDetector
477 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector)
478 with RoundtripFits(detector, storage_class="DetectorV2") as roundtrip:
479 legacy_detector = roundtrip.get(storageClass="Detector")
480 assert isinstance(legacy_detector, LegacyDetector)
481 compare_detector_to_legacy(detector, legacy_detector, is_raw_assembled=True)
484def test_repeated_metadata_keys_legacy_round_trip(
485 visit_image_components: dict[str, Any],
486 reset_afw_mask_planes: None, # noqa: F811
487) -> None:
488 """Verify that a header key set more than once keeps all of its values
489 when a VisitImage is converted to a legacy Exposure and back.
490 """
491 from lsst.afw.detection import GaussianPsf
493 opaque_metadata = FitsOpaqueMetadata()
494 header = astropy.io.fits.Header()
495 header.append(("PLATFORM", "lsstcam"), end=True)
496 # SubtractBackgroundTask writes one BGMEAN card per background fit.
497 header.append(("BGMEAN", 1.5), end=True)
498 header.append(("BGMEAN", 2.5), end=True)
499 opaque_metadata.extract_legacy_primary_header(header)
500 visit_image = VisitImage(
501 visit_image_components["image"],
502 variance=visit_image_components["variance"],
503 # A legacy-backed PSF, so that to_legacy can attach it to the
504 # Exposure and from_legacy can read it back.
505 psf=PointSpreadFunction.from_legacy(GaussianPsf(33, 33, 2.5), bounds=Box.factory[0:1024, 0:1024]),
506 mask_schema=visit_image_components["mask_schema"],
507 sky_projection=visit_image_components["sky_projection"],
508 detector=visit_image_components["detector"],
509 obs_info=visit_image_components["obs_info"],
510 band="r",
511 )
512 visit_image._opaque_metadata = opaque_metadata
514 legacy_exposure = visit_image.to_legacy()
515 legacy_metadata = legacy_exposure.getMetadata()
516 assert legacy_metadata.getArray("BGMEAN") == [1.5, 2.5]
517 assert legacy_metadata["PLATFORM"] == "lsstcam"
518 # Add a random key directly
519 legacy_metadata.add("DUMMYVAR", 0.5)
520 legacy_metadata.add("DUMMYVAR", 1.5)
522 round_tripped = VisitImage.from_legacy(
523 legacy_exposure,
524 instrument=visit_image_components["obs_info"].instrument,
525 visit=visit_image_components["sky_projection"].pixel_frame.visit,
526 )
527 round_tripped_header = round_tripped._opaque_metadata.headers[ExtensionKey()]
528 assert "DUMMYVAR" in round_tripped_header
529 assert "BGMEAN" in round_tripped_header
530 assert round_tripped_header["PLATFORM"] == "lsstcam"
531 assert [card.value for card in round_tripped_header.cards if card.keyword == "BGMEAN"] == [1.5, 2.5]
532 assert [card.value for card in round_tripped_header.cards if card.keyword == "DUMMYVAR"] == [0.5, 1.5]
535def test_external_metadata_legacy_round_trip(
536 visit_image_components: dict[str, Any],
537 reset_afw_mask_planes: None, # noqa: F811
538 tmp_path: Path,
539) -> None:
540 """Verify that native metadata is written to a legacy file as
541 ``LSST IMAGES`` cards and is not visible as external metadata when the
542 file is read back.
543 """
544 from lsst.afw.detection import GaussianPsf
546 opaque_metadata = FitsOpaqueMetadata()
547 header = astropy.io.fits.Header()
548 header.append(("PLATFORM", "lsstcam"), end=True)
549 header.append(("BGMEAN", 1.5), end=True)
550 header.append(("BGMEAN", 2.5), end=True)
551 # An external ID card must neither block nor replace the native "id".
552 header.append(("ID", 99), end=True)
553 opaque_metadata.extract_legacy_primary_header(header)
554 visit_image = VisitImage(
555 visit_image_components["image"],
556 variance=visit_image_components["variance"],
557 psf=PointSpreadFunction.from_legacy(GaussianPsf(33, 33, 2.5), bounds=Box.factory[0:1024, 0:1024]),
558 mask_schema=visit_image_components["mask_schema"],
559 sky_projection=visit_image_components["sky_projection"],
560 detector=visit_image_components["detector"],
561 obs_info=visit_image_components["obs_info"],
562 band="r",
563 metadata={"native_key": 7, "MixedCase": "yes"},
564 )
565 visit_image._opaque_metadata = opaque_metadata
566 path = tmp_path / "legacy.fits"
567 visit_image.to_legacy().writeFits(str(path))
569 with astropy.io.fits.open(path) as hdu_list:
570 primary = hdu_list[0].header
571 native_cards = {}
572 n = 1
573 while f"LSST IMAGES KEY {n}" in primary:
574 native_cards[primary[f"LSST IMAGES KEY {n}"]] = primary[f"LSST IMAGES VALUE {n}"]
575 n += 1
576 assert native_cards == {"native_key": 7, "MixedCase": "yes"}
577 assert [card.value for card in primary.cards if card.keyword == "BGMEAN"] == [1.5, 2.5]
579 result = VisitImage.read_legacy(
580 str(path),
581 instrument=visit_image_components["obs_info"].instrument,
582 visit=visit_image_components["sky_projection"].pixel_frame.visit,
583 )
584 assert result.metadata.native["native_key"] == 7
585 assert result.metadata.native["MixedCase"] == "yes"
586 assert "id" in result.metadata.native
587 assert not [key for key in result.metadata.external if key.startswith("LSST IMAGES")]
588 assert result.metadata.external.get_all("BGMEAN") == (1.5, 2.5)
589 assert result.metadata["platform"] == "lsstcam"
590 assert result.metadata.external["ID"] == 99
592 # Converting the result back must not duplicate the external cards.
593 legacy_metadata = result.to_legacy().getMetadata()
594 assert legacy_metadata.getArray("BGMEAN") == [1.5, 2.5]
595 assert legacy_metadata["LSST IMAGES KEY 1"] == "native_key"
598@skip_no_h5py
599def test_round_trip_ndf(visit_image_components: dict[str, Any]) -> None:
600 """Verify NDF round-trip produces a VisitImage equal to the original."""
601 visit_image = make_visit_image(visit_image_components)
602 with RoundtripNdf(visit_image, "VisitImage") as roundtrip:
603 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False)
606@skip_no_h5py
607def test_fits_ndf_consistency(visit_image_components: dict[str, Any]) -> None:
608 """Verify FITS and NDF backends produce equal VisitImages on round-trip."""
609 visit_image = make_visit_image(visit_image_components)
610 with RoundtripFits(visit_image) as fits_rt, RoundtripNdf(visit_image) as ndf_rt:
611 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False)
612 assert_visit_images_equal(visit_image, ndf_rt.result, expect_view=False)
613 assert_visit_images_equal(fits_rt.result, ndf_rt.result, expect_view=False)
616def test_fits_json_consistency(visit_image_components: dict[str, Any]) -> None:
617 """Verify FITS and JSON backends produce equal VisitImages."""
618 visit_image = make_visit_image(visit_image_components)
619 with (
620 RoundtripFits(visit_image) as fits_rt,
621 RoundtripJson(visit_image) as json_rt,
622 ):
623 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False)
624 assert_visit_images_equal(visit_image, json_rt.result, expect_view=False)
625 assert_visit_images_equal(fits_rt.result, json_rt.result, expect_view=False)
628def test_read_write(visit_image_components: dict[str, Any]) -> None:
629 """Verify a VisitImage round-trips through FITS with correct compression.
631 Checks compression headers, subimage reads, equality, and opaque
632 metadata. Contains only butler-free assertions; component reads live
633 in `test_read_write_components`.
634 """
635 visit_image = make_visit_image(visit_image_components)
636 with RoundtripFits(visit_image, "VisitImage") as roundtrip:
637 # Check that we're still using the right compression, and that we
638 # wrote WCSs.
639 fits = roundtrip.inspect()
640 assert fits[1].header["ZCMPTYPE"] == "GZIP_2"
641 assert fits[1].header["CTYPE1"] == "RA---TAN"
642 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
643 assert fits[2].header["CTYPE1"] == "RA---TAN"
644 assert fits[3].header["ZCMPTYPE"] == "GZIP_2"
645 assert fits[3].header["CTYPE1"] == "RA---TAN"
646 # Check a subimage read (no component arg — does not trigger a skip).
647 subbox = Box.factory[8:13, 9:30]
648 subimage = roundtrip.get(bbox=subbox)
649 assert_masked_images_equal(subimage, visit_image[subbox], expect_view=False)
651 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False)
652 # Check that the round-tripped headers are the same (up to card order).
653 assert len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]) == 1
654 assert dict(visit_image._opaque_metadata.headers[ExtensionKey()]) == dict(
655 roundtrip.result._opaque_metadata.headers[ExtensionKey()]
656 )
657 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")]
658 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")]
659 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")]
660 # Spot-check the concrete background contents (names, field types,
661 # subtracted entry) against the known fixture, so the equality check
662 # above is not vacuously satisfied by empty background maps.
663 assert isinstance(roundtrip.result.backgrounds, BackgroundMap)
664 assert roundtrip.result.backgrounds.keys() == {"standard"}
665 assert isinstance(roundtrip.result.backgrounds["standard"].field, ChebyshevField)
666 assert roundtrip.result.backgrounds.subtracted.name == "standard"
667 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image."
670def test_read_write_components(visit_image_components: dict[str, Any]) -> None:
671 """Verify component reads and storage-class overrides round-trip correctly.
673 Requires a butler; skips when `lsst.daf.butler` is absent.
674 Butler-free assertions live in `test_read_write`.
675 """
676 c = visit_image_components
677 visit_image = make_visit_image(c)
678 with RoundtripFits(visit_image, "VisitImage") as roundtrip:
679 subbox = Box.factory[8:13, 9:30]
680 subimage = roundtrip.get(bbox=subbox)
682 # Get an explicit masked image to compare with the subimage.
683 subimage_masked = roundtrip.get("masked_image", bbox=subbox)
684 assert_masked_images_equal(subimage_masked, subimage, expect_view=False)
686 # Get the same masked image in a multi-component get and ensure
687 # it is the same thing.
688 components = roundtrip.get("components", components=["masked_image", "psf"], bbox=subbox)
689 assert set(components) == {"masked_image", "psf"}
690 assert_masked_images_equal(components["masked_image"], subimage_masked, expect_view=False)
692 assert roundtrip.get("bbox") == visit_image.bbox
694 obs_info = roundtrip.get("obs_info")
695 assert isinstance(obs_info, ObservationInfo)
696 assert obs_info == visit_image.obs_info
698 summary_stats = roundtrip.get("summary_stats")
699 assert isinstance(summary_stats, ObservationSummaryStats)
700 assert summary_stats == visit_image.summary_stats
702 psf = roundtrip.get("psf")
703 assert isinstance(psf, GaussianPointSpreadFunction)
704 assert psf.kernel_bbox == c["gaussian_psf"].kernel_bbox
706 backgrounds = roundtrip.get("backgrounds")
707 assert isinstance(backgrounds, BackgroundMap)
708 assert backgrounds.keys() == {"standard"}
709 assert isinstance(backgrounds["standard"].field, ChebyshevField)
710 assert backgrounds.subtracted.name == "standard"
711 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image."
713 # Test some components get edge cases.
714 components = roundtrip.get("components", components="image")
715 assert isinstance(components["image"], Image)
717 components = roundtrip.get("components")
718 assert set(components) == {
719 "image",
720 "variance",
721 "psf",
722 "bbox",
723 "mask",
724 "obs_info",
725 "backgrounds",
726 "detector",
727 "aperture_corrections",
728 "sky_projection",
729 "summary_stats",
730 "photometric_scaling",
731 }
733 # Butler morphs RuntimeError to ValueError.
734 with pytest.raises(ValueError, match="component nonexistent"):
735 roundtrip.get("components", components=["image", "nonexistent"])
737 with pytest.raises(ValueError, match="should not be specified"):
738 roundtrip.get("components", components=["image", "components"])
740 with pytest.raises(ValueError, match="empty request"):
741 roundtrip.get("components", components=[])
743 with pytest.raises(ValueError, match="not known to any of the requested components"):
744 # PSF does not know how to use bbox so this fails.
745 roundtrip.get("components", components="psf", bbox=subbox)
748def test_sum_background_round_trip_fits(visit_image_components: dict[str, Any]) -> None:
749 """Verify FITS backend keeps two same-named SumField operands as distinct
750 EXTVERs.
751 """
752 visit_image = make_visit_image(visit_image_components)
753 visit = _make_sum_background_visit_image(visit_image_components, visit_image)
754 with RoundtripFits(visit) as roundtrip:
755 _check_sum_background_round_trip(roundtrip.result, visit)
758@skip_no_h5py
759def test_sum_background_round_trip_ndf(visit_image_components: dict[str, Any]) -> None:
760 """Verify NDF backend disambiguates the repeated ``data`` leaf, just as
761 the FITS backend does.
762 """
763 visit_image = make_visit_image(visit_image_components)
764 visit = _make_sum_background_visit_image(visit_image_components, visit_image)
765 with RoundtripNdf(visit) as roundtrip:
766 _check_sum_background_round_trip(roundtrip.result, visit)
769@pytest.mark.parametrize(
770 ("scaling_unit", "operation"),
771 [(u.electron / u.nJy, "multiply"), (u.nJy / u.electron, "divide")],
772 ids=["multiply", "divide"],
773)
774def test_convert_unit_subimage(
775 visit_image_components: dict[str, Any],
776 scaling_unit: u.UnitBase,
777 operation: Literal["multiply", "divide"],
778) -> None:
779 """Verify that converting the units of a subimage applies only the portion
780 of the photometric scaling that overlaps the subimage.
782 A photometric scaling keeps the bounds it was modeled over when the image
783 is subset, so both branches of the conversion must render it over the
784 subimage's bbox rather than over its own bounds.
785 """
786 visit_image = make_visit_image(visit_image_components)
787 scaling = ChebyshevField(
788 visit_image.bbox,
789 np.array([[4.0, 0.5, 0.125], [0.25, 0.0625, 0.0], [0.03125, 0.0, 0.0]]),
790 unit=scaling_unit,
791 )
792 visit_image.photometric_scaling = scaling
793 # Trim a different number of pixels from each side so a scaling rendered
794 # over the wrong box cannot match by chance.
795 subbox = Box.factory[
796 visit_image.bbox.y.start + 7 : visit_image.bbox.y.stop - 13,
797 visit_image.bbox.x.start + 3 : visit_image.bbox.x.stop - 21,
798 ]
799 subimage = visit_image[subbox]
800 assert subimage.photometric_scaling.bounds.bbox == visit_image.bbox
802 converted = subimage.convert_unit(u.electron)
804 assert converted.unit == u.electron
805 assert converted.image.bbox == subbox
806 scaling_array = scaling.render(subbox, dtype=subimage.image.array.dtype).array
807 if operation == "divide":
808 scaling_array = 1.0 / scaling_array
809 assert_values_equal(converted.image.array, subimage.image.array * scaling_array, rtol=1e-5)
810 assert_values_equal(converted.variance.array, subimage.variance.array * scaling_array**2, rtol=1e-5)
813@dataclasses.dataclass
814class _LegacyTestData:
815 filename: str
816 plane_map: dict[str, MaskPlane] = dataclasses.field(default_factory=get_legacy_visit_image_mask_planes)
817 unit: u.Unit = u.nJy
818 storage_class: str = "VisitImage"
819 read_cls: type[VisitImage] = VisitImage
820 legacy_exposure: LegacyExposure = dataclasses.field(init=False)
822 @classmethod
823 def get(
824 cls, which: Literal["visit_image", "preliminary_visit_image", "difference_image"]
825 ) -> _LegacyTestData:
826 if EXTERNAL_DATA_DIR is None: 826 ↛ 828line 826 didn't jump to line 828 because the condition on line 826 was always true
827 pytest.skip("TESTDATA_IMAGES is not set up.")
828 result = cls(
829 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", f"{which}.fits"),
830 )
831 match which:
832 case "preliminary_visit_image":
833 result.unit = u.electron
834 case "difference_image":
835 result.storage_class = "DifferenceImage"
836 result.read_cls = DifferenceImage
837 result.plane_map = get_legacy_difference_image_mask_planes()
838 case "visit_image":
839 pass
840 try:
841 from lsst.afw.image import ExposureFitsReader
843 result.legacy_exposure = ExposureFitsReader(result.filename).read()
844 except ImportError:
845 pytest.skip("lsst.afw.image is not available; cannot read legacy exposures")
846 result.visit_image = result.read_cls.read_legacy(
847 result.filename, preserve_quantization=True, plane_map=result.plane_map
848 )
849 return result
852@pytest.fixture(params=["visit_image", "preliminary_visit_image", "difference_image"])
853def legacy_test_data(request: pytest.FixtureRequest, reset_afw_mask_planes: None) -> _LegacyTestData: # noqa: F811
854 """Return legacy test data.
856 Tests that depend on this parameterized fixture run on all of the legacy
857 test images.
858 """
859 return _LegacyTestData.get(request.param)
862@pytest.fixture(params=["visit_image", "difference_image"])
863def legacy_test_data_calibrated(
864 request: pytest.FixtureRequest,
865 reset_afw_mask_planes: None, # noqa: F811
866) -> _LegacyTestData:
867 """Return legacy test data for calibrated images only.
869 Tests that depend on this parameterized fixture do not run on
870 preliminary_visit_image, since that has 'electron' pixel units
871 """
872 return _LegacyTestData.get(request.param)
875def _check_legacy_obs_info(obs_info: ObservationInfo | None) -> None:
876 """Assert obs_info carries expected LSSTCam/DP2 field values."""
877 assert isinstance(obs_info, ObservationInfo)
878 assert obs_info.instrument == "LSSTCam"
879 assert obs_info.detector_num == 85, obs_info
880 assert obs_info.detector_unique_name == "R21_S11", obs_info
881 assert obs_info.physical_filter == "r_57", obs_info
884def test_legacy_errors(legacy_test_data: _LegacyTestData) -> None:
885 """Verify that from_legacy and read_legacy raise ValueError on
886 conflicting arguments.
887 """
888 with pytest.raises(ValueError, match="does not match"):
889 VisitImage.from_legacy(legacy_test_data.legacy_exposure, instrument="HSC")
890 with pytest.raises(ValueError, match="does not match"):
891 VisitImage.from_legacy(legacy_test_data.legacy_exposure, visit=123456)
892 with pytest.raises(ValueError, match="BUNIT value .* disagrees with given unit"):
893 VisitImage.from_legacy(legacy_test_data.legacy_exposure, unit=u.mJy)
894 visit = VisitImage.from_legacy(
895 legacy_test_data.legacy_exposure,
896 instrument="LSSTCam",
897 unit=legacy_test_data.unit,
898 visit=2025052000177,
899 )
900 assert visit.unit == legacy_test_data.unit
902 with pytest.raises(ValueError, match="does not match"):
903 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, instrument="HSC")
904 with pytest.raises(ValueError, match="does not match"):
905 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, visit=123456)
908def test_component_reads(legacy_test_data: _LegacyTestData) -> None:
909 """Verify that individual components can be read from a legacy FITS
910 file.
911 """
912 visit = VisitImage.read_legacy(legacy_test_data.filename)
913 proj = VisitImage.read_legacy(legacy_test_data.filename, component="sky_projection")
914 assert_sky_projections_equal(proj, visit.sky_projection, expect_identity=False)
915 image = VisitImage.read_legacy(legacy_test_data.filename, component="image")
916 assert image == visit.image
917 assert_sky_projections_equal(proj, image.sky_projection, expect_identity=False)
918 variance = VisitImage.read_legacy(legacy_test_data.filename, component="variance")
919 assert variance == visit.variance
920 assert_sky_projections_equal(proj, variance.sky_projection, expect_identity=False)
921 mask = VisitImage.read_legacy(legacy_test_data.filename, component="mask")
922 assert mask == visit.mask
923 assert_sky_projections_equal(proj, mask.sky_projection, expect_identity=False)
924 psf = VisitImage.read_legacy(legacy_test_data.filename, component="psf")
925 assert isinstance(psf, PointSpreadFunction)
926 obs_info = VisitImage.read_legacy(legacy_test_data.filename, component="obs_info")
927 _check_legacy_obs_info(obs_info)
928 summary_stats = VisitImage.read_legacy(legacy_test_data.filename, component="summary_stats")
929 assert isinstance(summary_stats, ObservationSummaryStats)
930 assert summary_stats.nPsfStar == legacy_test_data.legacy_exposure.info.getSummaryStats().nPsfStar
931 compare_aperture_corrections_to_legacy(
932 VisitImage.read_legacy(legacy_test_data.filename, component="aperture_corrections"),
933 legacy_test_data.legacy_exposure.info.getApCorrMap(),
934 visit.bbox,
935 )
936 detector = VisitImage.read_legacy(legacy_test_data.filename, component="detector")
937 compare_detector_to_legacy(
938 detector, legacy_test_data.legacy_exposure.getDetector(), is_raw_assembled=True
939 )
940 photometric_scaling = VisitImage.read_legacy(legacy_test_data.filename, component="photometric_scaling")
941 compare_photo_calib_to_legacy(
942 photometric_scaling,
943 legacy_test_data.legacy_exposure.getPhotoCalib(),
944 subimage_bbox=visit.bbox,
945 )
948def test_legacy_obs_info(legacy_test_data: _LegacyTestData) -> None:
949 """Verify that ObservationInfo is constructed correctly from a legacy
950 exposure.
951 """
952 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map)
953 assert legacy.obs_info is not None
954 assert legacy.obs_info == legacy_test_data.visit_image.obs_info
955 assert legacy.obs_info is not None # for mypy
956 assert legacy.obs_info.instrument == "LSSTCam"
957 assert legacy.obs_info.detector_num == 85, legacy.obs_info
958 assert legacy.obs_info.detector_unique_name == "R21_S11", legacy.obs_info
959 assert legacy.obs_info.physical_filter == "r_57", legacy.obs_info
962def test_aperture_corrections_to_legacy(legacy_test_data: _LegacyTestData) -> None:
963 """Verify that aperture corrections round-trip through a legacy
964 ApCorrMap.
965 """
966 ap_corrections = legacy_test_data.visit_image.aperture_corrections
967 legacy_ap_corr_map = aperture_corrections_to_legacy(ap_corrections)
968 compare_aperture_corrections_to_legacy(
969 ap_corrections,
970 legacy_ap_corr_map,
971 legacy_test_data.visit_image.bbox,
972 )
975def _check_legacy_headers(visit_image: VisitImage) -> None:
976 """Assert that primary and extension headers are stripped correctly."""
977 header = visit_image._opaque_metadata.headers[ExtensionKey()]
978 assert "EXPTIME" in header
979 assert header["PLATFORM"] == "lsstcam"
980 assert "LSST BUTLER ID" not in header
981 assert "AR HDU" not in header
982 assert "A_ORDER" not in header
983 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("IMAGE"), astropy.io.fits.Header())
984 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("MASK"), astropy.io.fits.Header())
985 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("VARIANCE"), astropy.io.fits.Header())
988def test_read_legacy_headers(legacy_test_data: _LegacyTestData) -> None:
989 """Verify that headers were stripped and interpreted correctly in
990 read_legacy.
991 """
992 assert legacy_test_data.visit_image.unit == legacy_test_data.unit
993 _check_legacy_headers(legacy_test_data.visit_image)
996def test_from_legacy_headers(legacy_test_data: _LegacyTestData) -> None:
997 """Verify that from_legacy handles primary and extension headers
998 correctly.
999 """
1000 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map)
1001 assert legacy.unit == legacy_test_data.unit
1002 _check_legacy_headers(legacy)
1005def test_rewrite(legacy_test_data: _LegacyTestData) -> None:
1006 """Verify that a legacy VisitImage can be rewritten and round-trips both
1007 pixel values and all components.
1008 """
1009 with RoundtripFits(legacy_test_data.visit_image, legacy_test_data.storage_class) as roundtrip:
1010 fits = roundtrip.inspect()
1011 assert fits[1].header["ZCMPTYPE"] == "RICE_1"
1012 assert fits[1].header["CTYPE1"] == "RA---TAN-SIP"
1013 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
1014 assert fits[2].header["CTYPE1"] == "RA---TAN-SIP"
1015 assert fits[3].header["ZCMPTYPE"] == "RICE_1"
1016 assert fits[3].header["CTYPE1"] == "RA---TAN-SIP"
1017 subbox = Box.factory[8:13, 9:30]
1018 subimage = roundtrip.get(bbox=subbox)
1019 assert_masked_images_equal(subimage, legacy_test_data.visit_image[subbox], expect_view=False)
1020 alternates: dict[str, Any] = {}
1021 assert roundtrip.get("bbox") == legacy_test_data.visit_image.bbox
1022 alternates = {
1023 k: roundtrip.get(k)
1024 for k in [
1025 "sky_projection",
1026 "image",
1027 "mask",
1028 "variance",
1029 "psf",
1030 "obs_info",
1031 "summary_stats",
1032 "aperture_corrections",
1033 "detector",
1034 "photometric_scaling",
1035 ]
1036 }
1037 legacy_exposure = roundtrip.get(storageClass="Exposure")
1038 assert isinstance(legacy_exposure, LegacyExposure)
1039 compare_visit_image_to_legacy(
1040 legacy_test_data.visit_image,
1041 legacy_exposure,
1042 expect_view=False,
1043 plane_map=legacy_test_data.plane_map,
1044 **MINIMAL_VISIT_DATA_ID,
1045 )
1046 if legacy_test_data.visit_image.unit == u.nJy:
1047 assert legacy_exposure.getPhotoCalib()._isConstant
1048 assert legacy_exposure.getPhotoCalib().getCalibrationMean() == 1.0
1049 else:
1050 compare_photo_calib_to_legacy(
1051 legacy_test_data.visit_image.photometric_scaling,
1052 legacy_exposure.getPhotoCalib(),
1053 subimage_bbox=subbox,
1054 )
1055 assert legacy_exposure.info.getId() == legacy_test_data.legacy_exposure.info.getId()
1056 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
1057 assert isinstance(visit_info, LegacyVisitInfo)
1058 assert visit_info.getInstrumentLabel() == "LSSTCam"
1060 assert_visit_images_equal(roundtrip.result, legacy_test_data.visit_image, expect_view=False)
1061 assert dict(legacy_test_data.visit_image._opaque_metadata.headers[ExtensionKey()]) == dict(
1062 roundtrip.result._opaque_metadata.headers[ExtensionKey()]
1063 )
1064 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")]
1065 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")]
1066 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")]
1067 assert roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"] == "lsstcam"
1068 compare_visit_image_to_legacy(
1069 roundtrip.result,
1070 legacy_test_data.legacy_exposure,
1071 expect_view=False,
1072 plane_map=legacy_test_data.plane_map,
1073 **MINIMAL_VISIT_DATA_ID,
1074 alternates=alternates,
1075 )
1076 compare_visit_image_to_legacy(
1077 legacy_test_data.read_cls.from_legacy(
1078 legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map
1079 ),
1080 legacy_test_data.legacy_exposure,
1081 expect_view=True,
1082 plane_map=legacy_test_data.plane_map,
1083 **MINIMAL_VISIT_DATA_ID,
1084 )
1087def test_butler_converters(legacy_test_data: _LegacyTestData) -> None:
1088 """Verify that a VisitImage can be read from a Butler dataset written as
1089 an Exposure.
1090 """
1091 try:
1092 from lsst.daf.butler import FileDataset
1093 except ImportError:
1094 pytest.skip("lsst.daf.butler could not be imported.")
1096 with TemporaryButler(legacy="ExposureF") as helper:
1097 helper.butler.ingest(
1098 FileDataset(path=legacy_test_data.filename, refs=[helper.legacy]), transfer="symlink"
1099 )
1100 visit_image_ref = helper.legacy.overrideStorageClass(legacy_test_data.storage_class)
1101 with warnings.catch_warnings():
1102 warnings.filterwarnings("ignore", message=".*filter label mismatch.*", category=UserWarning)
1103 visit_image = helper.butler.get(visit_image_ref)
1104 assert visit_image._opaque_metadata.precompressed.keys() == set()
1105 visit_image = helper.butler.get(visit_image_ref, parameters={"preserve_quantization": True})
1106 assert visit_image._opaque_metadata.precompressed.keys() == {"IMAGE", "VARIANCE"}
1107 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
1108 assert bbox == visit_image.bbox
1109 alternates = {
1110 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
1111 for k in ["image", "mask", "variance", "bbox", "psf", "detector"]
1112 }
1113 compare_visit_image_to_legacy(
1114 visit_image,
1115 legacy_test_data.legacy_exposure,
1116 expect_view=False,
1117 plane_map=legacy_test_data.plane_map,
1118 alternates=alternates,
1119 **MINIMAL_VISIT_DATA_ID,
1120 )
1121 helper.butler.pruneDatasets([helper.legacy], purge=True, unstore=True, disassociate=True)
1122 visit_image.metadata["MixedCaseKey"] = 52
1123 helper.butler.put(visit_image, visit_image_ref)
1124 with warnings.catch_warnings():
1125 warnings.filterwarnings("ignore", message=".*filter label mismatch.*", category=UserWarning)
1126 legacy_exposure = helper.butler.get(helper.legacy)
1127 compare_visit_image_to_legacy(
1128 visit_image,
1129 legacy_exposure,
1130 expect_view=False,
1131 plane_map=legacy_test_data.plane_map,
1132 alternates=alternates,
1133 **MINIMAL_VISIT_DATA_ID,
1134 )
1135 visit_image_2 = helper.butler.get(visit_image_ref)
1136 compare_visit_image_to_legacy(
1137 visit_image_2,
1138 legacy_exposure,
1139 expect_view=False,
1140 plane_map=legacy_test_data.plane_map,
1141 alternates=alternates,
1142 **MINIMAL_VISIT_DATA_ID,
1143 )
1144 assert visit_image_2.metadata["MixedCaseKey"] == 52
1147def test_convert_unit(legacy_test_data_calibrated: _LegacyTestData) -> None:
1148 """Verify convert_unit round-trips between nJy, mJy, and electron via
1149 photometric_scaling.
1150 """
1151 from lsst.afw.table import ExposureCatalog
1153 legacy_test_data = legacy_test_data_calibrated
1154 original = legacy_test_data.visit_image.copy()
1155 with pytest.raises(u.UnitConversionError):
1156 original.convert_unit(u.electron)
1157 visit_image_nJy = original.convert_unit(u.nJy, copy=False)
1158 assert np.may_share_memory(visit_image_nJy.image.array, original.image.array)
1159 assert np.may_share_memory(visit_image_nJy.variance.array, original.variance.array)
1160 with pytest.raises(u.UnitConversionError):
1161 original.convert_unit(u.mJy, copy=False)
1162 visit_image_mJy = original.convert_unit(u.mJy, copy="as-needed")
1163 assert visit_image_mJy.unit == u.mJy
1164 assert_values_equal(visit_image_mJy.image.array, original.image.array * 1e-6, rtol=1e-5)
1165 assert np.may_share_memory(visit_image_nJy.mask.array, original.mask.array)
1166 assert_values_equal(visit_image_mJy.variance.array, original.variance.array * 1e-12, rtol=1e-5)
1167 legacy_exposure_mJy = visit_image_mJy.to_legacy()
1168 assert_values_equal(legacy_exposure_mJy.getPhotoCalib().getCalibrationMean(), 1e6, rtol=1e-14)
1169 legacy_masked_image_nJy = legacy_exposure_mJy.getPhotoCalib().calibrateImage(
1170 legacy_exposure_mJy.maskedImage
1171 )
1172 assert_values_equal(visit_image_nJy.image.array, legacy_masked_image_nJy.image.array, rtol=1e-5)
1173 assert_values_equal(visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array, rtol=1e-5)
1174 assert np.may_share_memory(visit_image_mJy.mask.array, original.mask.array)
1175 assert visit_image_mJy.sky_projection is original.sky_projection
1176 assert visit_image_mJy.obs_info is original.obs_info
1177 assert visit_image_mJy.summary_stats is original.summary_stats
1178 assert visit_image_mJy.psf is original.psf
1179 assert visit_image_mJy.detector is original.detector
1180 assert visit_image_mJy.bounds is original.bounds
1181 assert visit_image_mJy.aperture_corrections is original.aperture_corrections
1182 assert visit_image_mJy.photometric_scaling is original.photometric_scaling
1183 visit_summary = ExposureCatalog.readFits(
1184 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_summary.fits")
1185 )
1186 legacy_photo_calib = visit_summary.find(DP2_VISIT_DETECTOR_DATA_ID["detector"]).getPhotoCalib()
1187 visit_image_nJy.photometric_scaling = field_from_legacy_photo_calib(
1188 legacy_photo_calib, bounds=original.detector.bbox, instrumental_unit=u.electron
1189 )
1190 compare_photo_calib_to_legacy(
1191 visit_image_nJy.photometric_scaling,
1192 legacy_test_data.legacy_exposure.getPhotoCalib(),
1193 applied_legacy_photo_calib=legacy_photo_calib,
1194 subimage_bbox=visit_image_nJy.bbox,
1195 )
1196 with pytest.raises(u.UnitConversionError):
1197 visit_image_nJy.convert_unit(u.mm)
1198 with pytest.raises(u.UnitConversionError):
1199 visit_image_nJy.convert_unit(u.electron, copy=False)
1200 legacy_masked_image_e = legacy_photo_calib.uncalibrateImage(legacy_test_data.legacy_exposure.maskedImage)
1201 visit_image_e = visit_image_nJy.convert_unit(u.electron)
1202 assert_values_equal(visit_image_e.image.array, legacy_masked_image_e.image.array, rtol=1e-5)
1203 assert_values_equal(visit_image_e.variance.array, legacy_masked_image_e.variance.array, rtol=1e-5)
1204 assert not np.may_share_memory(visit_image_e.mask.array, visit_image_nJy.mask.array)
1205 visit_image_mJy.photometric_scaling = visit_image_nJy.photometric_scaling
1206 visit_image_e = visit_image_mJy.convert_unit(u.electron)
1207 assert_values_equal(visit_image_e.image.array, legacy_masked_image_e.image.array, rtol=1e-5)
1208 assert_values_equal(visit_image_e.variance.array, legacy_masked_image_e.variance.array, rtol=1e-5)
1209 visit_image_nJy_2 = visit_image_e.convert_unit(u.nJy)
1210 assert_values_equal(visit_image_nJy_2.image.array, visit_image_nJy.image.array, rtol=1e-5)
1211 assert_values_equal(visit_image_nJy_2.variance.array, original.variance.array, rtol=1e-5)
1212 visit_image_e.photometric_scaling = visit_image_nJy.photometric_scaling * (1e-6 * u.mJy / u.nJy)
1213 visit_image_nJy_3 = visit_image_e.convert_unit(u.nJy)
1214 assert_values_equal(visit_image_nJy_3.image.array, visit_image_nJy.image.array, rtol=1e-5)
1215 assert_values_equal(visit_image_nJy_3.variance.array, original.variance.array, rtol=1e-5)
1216 legacy_exposure_e = visit_image_e.to_legacy()
1217 assert_values_equal(
1218 legacy_exposure_e.getPhotoCalib().getCalibrationMean(),
1219 legacy_photo_calib.getCalibrationMean(),
1220 rtol=1e-5,
1221 )
1222 legacy_masked_image_nJy = legacy_exposure_e.getPhotoCalib().calibrateImage(legacy_exposure_e.maskedImage)
1223 assert_values_equal(visit_image_nJy.image.array, legacy_masked_image_nJy.image.array, rtol=1e-5)
1224 assert_values_equal(visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array, rtol=1e-5)
1227def test_background_map_describe() -> None:
1228 """An empty BackgroundMap._describe reports no backgrounds inline."""
1229 bg_map = BackgroundMap()
1230 assert isinstance(bg_map, DescribableMixin)
1231 report = bg_map._describe()
1232 assert isinstance(report, Report)
1233 assert report.type_name == "BackgroundMap"
1234 assert report.inline
1235 assert report.summary == "no backgrounds"
1236 assert report.children == {}
1239def test_background_map_with_entries_describe() -> None:
1240 """BackgroundMap._describe recurses into each background model."""
1241 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1242 bg_map = BackgroundMap(
1243 [Background("sky", cheby, "Sky model."), Background("fringe", cheby)],
1244 subtracted="sky",
1245 )
1246 report = bg_map._describe()
1247 assert report.type_name == "BackgroundMap"
1248 assert report.inline
1249 # The inline summary leads with the subtracted background and marks it,
1250 # since that is all a composite holding this map will show.
1251 # Only the subtracted background is marked, not the whole line.
1252 assert report.summary == "**sky [SUBTRACTED]**; also: fringe"
1253 assert report.emphasis_markup
1254 # str reads the summary as text, so the delimiters never reach a user.
1255 assert report.to_str() == "sky [SUBTRACTED]; also: fringe"
1256 # Standalone, each background is a child carrying its own model's report.
1257 assert set(report.children) == {"sky", "fringe"}
1258 sky = report.children["sky"]
1259 assert sky.type_name == "ChebyshevField"
1260 # The marker is on the heading, where it is read before the fields.
1261 assert sky.title == "ChebyshevField **[SUBTRACTED]**"
1262 assert sky.emphasis_markup
1263 sky_fields = {f.label: f.value for f in sky.fields}
1264 assert sky_fields["description"] == "Sky model."
1265 # The model's own fields survive alongside the background's attributes.
1266 assert "bounds" in sky_fields
1267 # Only the subtracted one is marked, and an absent description is omitted.
1268 fringe = report.children["fringe"]
1269 assert fringe.title is None
1270 assert fringe.emphasis_markup is False
1271 assert "description" not in {f.label for f in fringe.fields}
1274def test_background_map_describe_none_subtracted() -> None:
1275 """A map with no subtracted background says so rather than being silent."""
1276 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1277 bg_map = BackgroundMap([Background("sky", cheby), Background("fringe", cheby)])
1278 report = bg_map._describe()
1279 assert report.summary == "sky, fringe (none subtracted)"
1280 assert report.emphasis_markup is False
1281 assert not any(child.emphasis_markup for child in report.children.values())
1284def test_background_map_describe_leads_with_subtracted() -> None:
1285 """The subtracted background leads the summary wherever it sits in the
1286 map, and only its child is marked.
1287 """
1288 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1289 bg_map = BackgroundMap(
1290 [Background("sky", cheby), Background("skyCorr", cheby)],
1291 subtracted="skyCorr",
1292 )
1293 report = bg_map._describe()
1294 assert report.to_str() == "skyCorr [SUBTRACTED]; also: sky"
1295 assert report.children["sky"].emphasis_markup is False
1296 assert report.children["skyCorr"].emphasis_markup
1299def test_background_map_describe_brief_skips_children() -> None:
1300 """A brief background map report keeps the summary but not the models."""
1301 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1302 bg_map = BackgroundMap([Background("sky", cheby)], subtracted="sky")
1303 report = bg_map._describe(DescribeOptions(brief=True))
1304 assert report.to_str() == "sky [SUBTRACTED]"
1305 assert report.children == {}
1308def test_visit_image_repr_str_with_unreadable_psf() -> None:
1309 """Repr and str succeed even when the PSF stored an ArchiveReadError.
1311 An unreadable component is a supported state; repr and str read only the
1312 cheap fields and summary, so they do not build the child tree at all.
1313 """
1314 path = current_fixture_path(FIXTURE_DIR, "visit_image")
1315 visit_image = read_archive(path)
1316 visit_image._psf = ArchiveReadError("psf unreadable")
1317 assert repr(visit_image).startswith("VisitImage(")
1318 assert str(visit_image).startswith("VisitImage(")
1321def test_visit_image_describe_names_an_unreadable_psf() -> None:
1322 """A full report says the PSF could not be read, and describes the rest.
1324 A PSF model can need a package the reader does not have installed, so an
1325 unreadable PSF must not cost the report of everything else.
1326 """
1327 path = current_fixture_path(FIXTURE_DIR, "visit_image")
1328 visit_image = read_archive(path)
1329 visit_image._psf = ArchiveReadError("Failed to import piff.")
1330 report = visit_image.describe(detail=True)
1331 psf = report.children["psf"]
1332 assert psf.inline
1333 assert psf.to_str() == "unreadable (Failed to import piff.)"
1334 # Every other component is described as usual.
1335 assert "sky_projection" in report.children
1336 assert "detector" in report.children
1337 # Both renderers run over the report that contains it.
1338 assert "unreadable" in report._repr_html_()
1339 report.__rich__()
1342def test_visit_image_slice_preserves_unreadable_psf() -> None:
1343 """Slicing propagates a deferred PSF read failure without raising it."""
1344 path = current_fixture_path(FIXTURE_DIR, "visit_image")
1345 visit_image = read_archive(path)
1346 error = ArchiveReadError("psf unreadable")
1347 visit_image._psf = error
1349 sliced = visit_image[...]
1351 assert sliced._psf is error
1352 with pytest.raises(ArchiveReadError, match="psf unreadable"):
1353 sliced.psf
1356def test_observation_summary_stats_describe() -> None:
1357 """ObservationSummaryStats._describe reports the statistics that are set.
1359 Unset ones are omitted rather than shown as NaN.
1360 """
1361 stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4, ra=180.0, dec=-30.0)
1362 assert isinstance(stats, DescribableMixin)
1363 report = stats._describe()
1364 assert isinstance(report, Report)
1365 assert report.type_name == "ObservationSummaryStats"
1366 # Scalars are packed into one group, ordered by name.
1367 (group,) = report.value_groups
1368 assert group.role is FieldRole.DERIVED
1369 assert [name for name, _ in group.values] == sorted(name for name, _ in group.values)
1370 values = dict(group.values)
1371 assert values["psfSigma"] == 2.5
1372 assert values["zeroPoint"] == 31.4
1373 assert values["ra"] == 180.0
1374 assert values["dec"] == -30.0
1375 # Unset statistics are omitted rather than reported as NaN, and the
1376 # serialization plumbing this class inherits never appears at all.
1377 assert "expTime" not in values
1378 assert "skyBg" not in values
1379 assert not {"metadata", "butler_info", "schema_version"} & set(values)
1382def test_observation_summary_stats_describe_omits_empty_sequences() -> None:
1383 """Sequence statistics appear only when they carry a value."""
1384 empty = ObservationSummaryStats(psfSigma=2.5)
1385 # raCorners defaults to all-NaN, which carries no more information than an
1386 # empty sequence does.
1387 assert all(math.isnan(v) for v in empty.raCorners)
1388 assert not any(f.label == "raCorners" for f in empty._describe().fields)
1390 filled = ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
1391 corners = next(f for f in filled._describe().fields if f.label == "raCorners")
1392 assert corners.value == (5.2, 5.4, 5.4, 5.2)
1393 assert corners.role is FieldRole.DERIVED
1396def test_observation_summary_stats_describe_brief_counts() -> None:
1397 """A brief report gives the number set rather than listing them."""
1398 stats = ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
1399 brief = stats._describe(DescribeOptions(brief=True))
1400 assert brief.type_name == "ObservationSummaryStats"
1401 assert brief.value_groups == []
1402 (field,) = brief.fields
1403 assert field.label == "statistics set"
1404 assert field.role is FieldRole.DERIVED
1406 # The count must agree with what the full report actually shows: two
1407 # scalars set by default plus psfSigma, and raCorners as a sequence.
1408 full = stats._describe()
1409 listed = len(full.fields) + sum(len(g.values) for g in full.value_groups)
1410 assert field.value.startswith(f"{listed} of ")
1411 assert brief.to_str() == full.to_str()
1412 assert f"{listed} of " in brief.to_str()
1415def test_observation_summary_stats_pydantic_repr() -> None:
1416 """ObservationSummaryStats uses pydantic's repr, not the mixin's."""
1417 stats = ObservationSummaryStats(psfSigma=2.5)
1418 r = repr(stats)
1419 assert r.startswith("ObservationSummaryStats(")
1420 assert "psfSigma=2.5" in r
1423def test_observation_summary_stats_str_is_the_report_summary() -> None:
1424 """The str output reports the count, not pydantic's field-by-field dump.
1426 The repr keeps the exhaustive form, since that is the one that
1427 round-trips.
1428 """
1429 stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
1430 assert str(stats) == stats.describe().to_str()
1431 # The two set here, plus the integer counters, which default to a genuine
1432 # zero rather than to NaN.
1433 assert str(stats) == "ObservationSummaryStats(4 of 66 statistics set)"
1434 assert stats.nPsfStar == 0
1435 assert stats.nShapeletsStar == 0
1436 # The unset statistics reach repr but not str.
1437 assert "nan" not in str(stats)
1438 assert "nan" in repr(stats)
1441def test_archive_tree_repr_omits_schema_bookkeeping() -> None:
1442 """Schema version fields mirror class constants, so repr leaves them out.
1444 They are never passed on construction and say nothing the type does not.
1445 """
1446 stats = ObservationSummaryStats(psfSigma=2.5, indirect=[1])
1447 for name in ("schema_version", "min_read_version", "indirect"):
1448 assert name not in repr(stats), name
1449 # They are still real fields, and hiding them from repr does not hide them
1450 # from serialization.
1451 assert stats.schema_version == ObservationSummaryStats.SCHEMA_VERSION
1452 dumped = stats.model_dump()
1453 for name in ("schema_version", "min_read_version", "indirect"):
1454 assert name in dumped, name
1457def _match_detector_to_image(visit_image: VisitImage) -> VisitImage:
1458 """Return the visit image with its detector shrunk to the image bbox.
1460 Parameters
1461 ----------
1462 visit_image
1463 Image to adjust.
1465 Returns
1466 -------
1467 visit_image : `VisitImage`
1468 Image whose `~VisitImage.bbox` equals its detector's, as a full-frame
1469 image read from a butler has.
1470 """
1471 detector = visit_image.detector
1472 visit_image._detector = Detector(
1473 detector._attributes.model_copy(update={"bbox": visit_image.bbox}),
1474 detector.amplifiers,
1475 detector._frames,
1476 detector.visit,
1477 )
1478 return visit_image
1481def test_visit_image_describe_hides_summary_stat_corners_for_a_full_frame(
1482 visit_image_components: dict[str, Any],
1483) -> None:
1484 """The sky corners appear once, in the projection's readable table."""
1485 visit = _match_detector_to_image(make_visit_image(visit_image_components))
1486 visit.summary_stats.raCorners = (1.0, 2.0, 3.0, 4.0)
1487 visit.summary_stats.decCorners = (-1.0, -2.0, -3.0, -4.0)
1488 report = visit.describe()
1489 stats = report.children["summary_stats"]
1490 assert not {"raCorners", "decCorners"} & {field.label for field in stats.fields}
1491 # The projection is where they are read instead, and it still has them.
1492 assert any(table.title == "Corners" for table in report.children["sky_projection"].tables)
1495def test_visit_image_cutout_keeps_summary_stat_corners() -> None:
1496 """A cutout keeps the statistics' corners, which still describe the whole
1497 image the statistics were measured on.
1498 """
1499 # The DP2 variant is a cutout of a real image, so it carries the
1500 # statistics measured on the whole of it.
1501 path = current_fixture_path(FIXTURE_DIR, "visit_image", variant="dp2")
1502 visit_image = read_archive(path)
1503 assert visit_image.bbox != visit_image.detector.bbox
1504 stats = visit_image.describe().children["summary_stats"]
1505 assert {"raCorners", "decCorners"} <= {field.label for field in stats.fields}