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