Coverage for tests/test_visit_image.py: 69%
696 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-23 03:09 -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_repeated_metadata_keys_legacy_round_trip(
455 visit_image_components: dict[str, Any],
456 reset_afw_mask_planes: None, # noqa: F811
457) -> None:
458 """Verify that a header key set more than once keeps all of its values
459 when a VisitImage is converted to a legacy Exposure and back.
460 """
461 from lsst.afw.detection import GaussianPsf
463 opaque_metadata = FitsOpaqueMetadata()
464 header = astropy.io.fits.Header()
465 header.append(("PLATFORM", "lsstcam"), end=True)
466 # SubtractBackgroundTask writes one BGMEAN card per background fit.
467 header.append(("BGMEAN", 1.5), end=True)
468 header.append(("BGMEAN", 2.5), end=True)
469 opaque_metadata.extract_legacy_primary_header(header)
470 visit_image = VisitImage(
471 visit_image_components["image"],
472 variance=visit_image_components["variance"],
473 # A legacy-backed PSF, so that to_legacy can attach it to the
474 # Exposure and from_legacy can read it back.
475 psf=PointSpreadFunction.from_legacy(GaussianPsf(33, 33, 2.5), bounds=Box.factory[0:1024, 0:1024]),
476 mask_schema=visit_image_components["mask_schema"],
477 sky_projection=visit_image_components["sky_projection"],
478 detector=visit_image_components["detector"],
479 obs_info=visit_image_components["obs_info"],
480 band="r",
481 )
482 visit_image._opaque_metadata = opaque_metadata
484 legacy_exposure = visit_image.to_legacy()
485 legacy_metadata = legacy_exposure.getMetadata()
486 assert legacy_metadata.getArray("BGMEAN") == [1.5, 2.5]
487 assert legacy_metadata["PLATFORM"] == "lsstcam"
488 # Add a random key directly
489 legacy_metadata.add("DUMMYVAR", 0.5)
490 legacy_metadata.add("DUMMYVAR", 1.5)
492 round_tripped = VisitImage.from_legacy(
493 legacy_exposure,
494 instrument=visit_image_components["obs_info"].instrument,
495 visit=visit_image_components["sky_projection"].pixel_frame.visit,
496 )
497 round_tripped_header = round_tripped._opaque_metadata.headers[ExtensionKey()]
498 assert "DUMMYVAR" in round_tripped_header
499 assert "BGMEAN" in round_tripped_header
500 assert round_tripped_header["PLATFORM"] == "lsstcam"
501 assert [card.value for card in round_tripped_header.cards if card.keyword == "BGMEAN"] == [1.5, 2.5]
502 assert [card.value for card in round_tripped_header.cards if card.keyword == "DUMMYVAR"] == [0.5, 1.5]
505@skip_no_h5py
506def test_round_trip_ndf(visit_image_components: dict[str, Any]) -> None:
507 """Verify NDF round-trip produces a VisitImage equal to the original."""
508 visit_image = make_visit_image(visit_image_components)
509 with RoundtripNdf(visit_image, "VisitImage") as roundtrip:
510 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False)
513@skip_no_h5py
514def test_fits_ndf_consistency(visit_image_components: dict[str, Any]) -> None:
515 """Verify FITS and NDF backends produce equal VisitImages on round-trip."""
516 visit_image = make_visit_image(visit_image_components)
517 with RoundtripFits(visit_image) as fits_rt, RoundtripNdf(visit_image) as ndf_rt:
518 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False)
519 assert_visit_images_equal(visit_image, ndf_rt.result, expect_view=False)
520 assert_visit_images_equal(fits_rt.result, ndf_rt.result, expect_view=False)
523def test_fits_json_consistency(visit_image_components: dict[str, Any]) -> None:
524 """Verify FITS and JSON backends produce equal VisitImages."""
525 visit_image = make_visit_image(visit_image_components)
526 with (
527 RoundtripFits(visit_image) as fits_rt,
528 RoundtripJson(visit_image) as json_rt,
529 ):
530 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False)
531 assert_visit_images_equal(visit_image, json_rt.result, expect_view=False)
532 assert_visit_images_equal(fits_rt.result, json_rt.result, expect_view=False)
535def test_read_write(visit_image_components: dict[str, Any]) -> None:
536 """Verify a VisitImage round-trips through FITS with correct compression.
538 Checks compression headers, subimage reads, equality, and opaque
539 metadata. Contains only butler-free assertions; component reads live
540 in `test_read_write_components`.
541 """
542 visit_image = make_visit_image(visit_image_components)
543 with RoundtripFits(visit_image, "VisitImage") as roundtrip:
544 # Check that we're still using the right compression, and that we
545 # wrote WCSs.
546 fits = roundtrip.inspect()
547 assert fits[1].header["ZCMPTYPE"] == "GZIP_2"
548 assert fits[1].header["CTYPE1"] == "RA---TAN"
549 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
550 assert fits[2].header["CTYPE1"] == "RA---TAN"
551 assert fits[3].header["ZCMPTYPE"] == "GZIP_2"
552 assert fits[3].header["CTYPE1"] == "RA---TAN"
553 # Check a subimage read (no component arg — does not trigger a skip).
554 subbox = Box.factory[8:13, 9:30]
555 subimage = roundtrip.get(bbox=subbox)
556 assert_masked_images_equal(subimage, visit_image[subbox], expect_view=False)
558 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False)
559 # Check that the round-tripped headers are the same (up to card order).
560 assert len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]) == 1
561 assert dict(visit_image._opaque_metadata.headers[ExtensionKey()]) == dict(
562 roundtrip.result._opaque_metadata.headers[ExtensionKey()]
563 )
564 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")]
565 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")]
566 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")]
567 # Spot-check the concrete background contents (names, field types,
568 # subtracted entry) against the known fixture, so the equality check
569 # above is not vacuously satisfied by empty background maps.
570 assert isinstance(roundtrip.result.backgrounds, BackgroundMap)
571 assert roundtrip.result.backgrounds.keys() == {"standard"}
572 assert isinstance(roundtrip.result.backgrounds["standard"].field, ChebyshevField)
573 assert roundtrip.result.backgrounds.subtracted.name == "standard"
574 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image."
577def test_read_write_components(visit_image_components: dict[str, Any]) -> None:
578 """Verify component reads and storage-class overrides round-trip correctly.
580 Requires a butler; skips when `lsst.daf.butler` is absent.
581 Butler-free assertions live in `test_read_write`.
582 """
583 c = visit_image_components
584 visit_image = make_visit_image(c)
585 with RoundtripFits(visit_image, "VisitImage") as roundtrip:
586 subbox = Box.factory[8:13, 9:30]
587 subimage = roundtrip.get(bbox=subbox)
589 # Get an explicit masked image to compare with the subimage.
590 subimage_masked = roundtrip.get("masked_image", bbox=subbox)
591 assert_masked_images_equal(subimage_masked, subimage, expect_view=False)
593 # Get the same masked image in a multi-component get and ensure
594 # it is the same thing.
595 components = roundtrip.get("components", components=["masked_image", "psf"], bbox=subbox)
596 assert set(components) == {"masked_image", "psf"}
597 assert_masked_images_equal(components["masked_image"], subimage_masked, expect_view=False)
599 assert roundtrip.get("bbox") == visit_image.bbox
601 obs_info = roundtrip.get("obs_info")
602 assert isinstance(obs_info, ObservationInfo)
603 assert obs_info == visit_image.obs_info
605 summary_stats = roundtrip.get("summary_stats")
606 assert isinstance(summary_stats, ObservationSummaryStats)
607 assert summary_stats == visit_image.summary_stats
609 psf = roundtrip.get("psf")
610 assert isinstance(psf, GaussianPointSpreadFunction)
611 assert psf.kernel_bbox == c["gaussian_psf"].kernel_bbox
613 backgrounds = roundtrip.get("backgrounds")
614 assert isinstance(backgrounds, BackgroundMap)
615 assert backgrounds.keys() == {"standard"}
616 assert isinstance(backgrounds["standard"].field, ChebyshevField)
617 assert backgrounds.subtracted.name == "standard"
618 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image."
620 # Test some components get edge cases.
621 components = roundtrip.get("components", components="image")
622 assert isinstance(components["image"], Image)
624 components = roundtrip.get("components")
625 assert set(components) == {
626 "image",
627 "variance",
628 "psf",
629 "bbox",
630 "mask",
631 "obs_info",
632 "backgrounds",
633 "detector",
634 "aperture_corrections",
635 "sky_projection",
636 "summary_stats",
637 "photometric_scaling",
638 }
640 # Butler morphs RuntimeError to ValueError.
641 with pytest.raises(ValueError, match="component nonexistent"):
642 roundtrip.get("components", components=["image", "nonexistent"])
644 with pytest.raises(ValueError, match="should not be specified"):
645 roundtrip.get("components", components=["image", "components"])
647 with pytest.raises(ValueError, match="empty request"):
648 roundtrip.get("components", components=[])
650 with pytest.raises(ValueError, match="not known to any of the requested components"):
651 # PSF does not know how to use bbox so this fails.
652 roundtrip.get("components", components="psf", bbox=subbox)
655def test_sum_background_round_trip_fits(visit_image_components: dict[str, Any]) -> None:
656 """Verify FITS backend keeps two same-named SumField operands as distinct
657 EXTVERs.
658 """
659 visit_image = make_visit_image(visit_image_components)
660 visit = _make_sum_background_visit_image(visit_image_components, visit_image)
661 with RoundtripFits(visit) as roundtrip:
662 _check_sum_background_round_trip(roundtrip.result, visit)
665@skip_no_h5py
666def test_sum_background_round_trip_ndf(visit_image_components: dict[str, Any]) -> None:
667 """Verify NDF backend disambiguates the repeated ``data`` leaf, just as
668 the FITS backend does.
669 """
670 visit_image = make_visit_image(visit_image_components)
671 visit = _make_sum_background_visit_image(visit_image_components, visit_image)
672 with RoundtripNdf(visit) as roundtrip:
673 _check_sum_background_round_trip(roundtrip.result, visit)
676@pytest.mark.parametrize(
677 ("scaling_unit", "operation"),
678 [(u.electron / u.nJy, "multiply"), (u.nJy / u.electron, "divide")],
679 ids=["multiply", "divide"],
680)
681def test_convert_unit_subimage(
682 visit_image_components: dict[str, Any],
683 scaling_unit: u.UnitBase,
684 operation: Literal["multiply", "divide"],
685) -> None:
686 """Verify that converting the units of a subimage applies only the portion
687 of the photometric scaling that overlaps the subimage.
689 A photometric scaling keeps the bounds it was modeled over when the image
690 is subset, so both branches of the conversion must render it over the
691 subimage's bbox rather than over its own bounds.
692 """
693 visit_image = make_visit_image(visit_image_components)
694 scaling = ChebyshevField(
695 visit_image.bbox,
696 np.array([[4.0, 0.5, 0.125], [0.25, 0.0625, 0.0], [0.03125, 0.0, 0.0]]),
697 unit=scaling_unit,
698 )
699 visit_image.photometric_scaling = scaling
700 # Trim a different number of pixels from each side so a scaling rendered
701 # over the wrong box cannot match by chance.
702 subbox = Box.factory[
703 visit_image.bbox.y.start + 7 : visit_image.bbox.y.stop - 13,
704 visit_image.bbox.x.start + 3 : visit_image.bbox.x.stop - 21,
705 ]
706 subimage = visit_image[subbox]
707 assert subimage.photometric_scaling.bounds.bbox == visit_image.bbox
709 converted = subimage.convert_unit(u.electron)
711 assert converted.unit == u.electron
712 assert converted.image.bbox == subbox
713 scaling_array = scaling.render(subbox, dtype=subimage.image.array.dtype).array
714 if operation == "divide":
715 scaling_array = 1.0 / scaling_array
716 assert_values_equal(converted.image.array, subimage.image.array * scaling_array, rtol=1e-5)
717 assert_values_equal(converted.variance.array, subimage.variance.array * scaling_array**2, rtol=1e-5)
720@dataclasses.dataclass
721class _LegacyTestData:
722 filename: str
723 plane_map: dict[str, MaskPlane] = dataclasses.field(default_factory=get_legacy_visit_image_mask_planes)
724 unit: u.Unit = u.nJy
725 storage_class: str = "VisitImage"
726 read_cls: type[VisitImage] = VisitImage
727 legacy_exposure: LegacyExposure = dataclasses.field(init=False)
729 @classmethod
730 def get(
731 cls, which: Literal["visit_image", "preliminary_visit_image", "difference_image"]
732 ) -> _LegacyTestData:
733 if EXTERNAL_DATA_DIR is None: 733 ↛ 735line 733 didn't jump to line 735 because the condition on line 733 was always true
734 pytest.skip("TESTDATA_IMAGES is not set up.")
735 result = cls(
736 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", f"{which}.fits"),
737 )
738 match which:
739 case "preliminary_visit_image":
740 result.unit = u.electron
741 case "difference_image":
742 result.storage_class = "DifferenceImage"
743 result.read_cls = DifferenceImage
744 result.plane_map = get_legacy_difference_image_mask_planes()
745 case "visit_image":
746 pass
747 try:
748 from lsst.afw.image import ExposureFitsReader
750 result.legacy_exposure = ExposureFitsReader(result.filename).read()
751 except ImportError:
752 pytest.skip("lsst.afw.image is not available; cannot read legacy exposures")
753 result.visit_image = result.read_cls.read_legacy(
754 result.filename, preserve_quantization=True, plane_map=result.plane_map
755 )
756 return result
759@pytest.fixture(params=["visit_image", "preliminary_visit_image", "difference_image"])
760def legacy_test_data(request: pytest.FixtureRequest, reset_afw_mask_planes: None) -> _LegacyTestData: # noqa: F811
761 """Return legacy test data.
763 Tests that depend on this parameterized fixture run on all of the legacy
764 test images.
765 """
766 return _LegacyTestData.get(request.param)
769@pytest.fixture(params=["visit_image", "difference_image"])
770def legacy_test_data_calibrated(
771 request: pytest.FixtureRequest,
772 reset_afw_mask_planes: None, # noqa: F811
773) -> _LegacyTestData:
774 """Return legacy test data for calibrated images only.
776 Tests that depend on this parameterized fixture do not run on
777 preliminary_visit_image, since that has 'electron' pixel units
778 """
779 return _LegacyTestData.get(request.param)
782def _check_legacy_obs_info(obs_info: ObservationInfo | None) -> None:
783 """Assert obs_info carries expected LSSTCam/DP2 field values."""
784 assert isinstance(obs_info, ObservationInfo)
785 assert obs_info.instrument == "LSSTCam"
786 assert obs_info.detector_num == 85, obs_info
787 assert obs_info.detector_unique_name == "R21_S11", obs_info
788 assert obs_info.physical_filter == "r_57", obs_info
791def test_legacy_errors(legacy_test_data: _LegacyTestData) -> None:
792 """Verify that from_legacy and read_legacy raise ValueError on
793 conflicting arguments.
794 """
795 with pytest.raises(ValueError, match="does not match"):
796 VisitImage.from_legacy(legacy_test_data.legacy_exposure, instrument="HSC")
797 with pytest.raises(ValueError, match="does not match"):
798 VisitImage.from_legacy(legacy_test_data.legacy_exposure, visit=123456)
799 with pytest.raises(ValueError, match="BUNIT value .* disagrees with given unit"):
800 VisitImage.from_legacy(legacy_test_data.legacy_exposure, unit=u.mJy)
801 visit = VisitImage.from_legacy(
802 legacy_test_data.legacy_exposure,
803 instrument="LSSTCam",
804 unit=legacy_test_data.unit,
805 visit=2025052000177,
806 )
807 assert visit.unit == legacy_test_data.unit
809 with pytest.raises(ValueError, match="does not match"):
810 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, instrument="HSC")
811 with pytest.raises(ValueError, match="does not match"):
812 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, visit=123456)
815def test_component_reads(legacy_test_data: _LegacyTestData) -> None:
816 """Verify that individual components can be read from a legacy FITS
817 file.
818 """
819 visit = VisitImage.read_legacy(legacy_test_data.filename)
820 proj = VisitImage.read_legacy(legacy_test_data.filename, component="sky_projection")
821 assert_sky_projections_equal(proj, visit.sky_projection, expect_identity=False)
822 image = VisitImage.read_legacy(legacy_test_data.filename, component="image")
823 assert image == visit.image
824 assert_sky_projections_equal(proj, image.sky_projection, expect_identity=False)
825 variance = VisitImage.read_legacy(legacy_test_data.filename, component="variance")
826 assert variance == visit.variance
827 assert_sky_projections_equal(proj, variance.sky_projection, expect_identity=False)
828 mask = VisitImage.read_legacy(legacy_test_data.filename, component="mask")
829 assert mask == visit.mask
830 assert_sky_projections_equal(proj, mask.sky_projection, expect_identity=False)
831 psf = VisitImage.read_legacy(legacy_test_data.filename, component="psf")
832 assert isinstance(psf, PointSpreadFunction)
833 obs_info = VisitImage.read_legacy(legacy_test_data.filename, component="obs_info")
834 _check_legacy_obs_info(obs_info)
835 summary_stats = VisitImage.read_legacy(legacy_test_data.filename, component="summary_stats")
836 assert isinstance(summary_stats, ObservationSummaryStats)
837 assert summary_stats.nPsfStar == legacy_test_data.legacy_exposure.info.getSummaryStats().nPsfStar
838 compare_aperture_corrections_to_legacy(
839 VisitImage.read_legacy(legacy_test_data.filename, component="aperture_corrections"),
840 legacy_test_data.legacy_exposure.info.getApCorrMap(),
841 visit.bbox,
842 )
843 detector = VisitImage.read_legacy(legacy_test_data.filename, component="detector")
844 compare_detector_to_legacy(
845 detector, legacy_test_data.legacy_exposure.getDetector(), is_raw_assembled=True
846 )
847 photometric_scaling = VisitImage.read_legacy(legacy_test_data.filename, component="photometric_scaling")
848 compare_photo_calib_to_legacy(
849 photometric_scaling,
850 legacy_test_data.legacy_exposure.getPhotoCalib(),
851 subimage_bbox=visit.bbox,
852 )
855def test_legacy_obs_info(legacy_test_data: _LegacyTestData) -> None:
856 """Verify that ObservationInfo is constructed correctly from a legacy
857 exposure.
858 """
859 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map)
860 assert legacy.obs_info is not None
861 assert legacy.obs_info == legacy_test_data.visit_image.obs_info
862 assert legacy.obs_info is not None # for mypy
863 assert legacy.obs_info.instrument == "LSSTCam"
864 assert legacy.obs_info.detector_num == 85, legacy.obs_info
865 assert legacy.obs_info.detector_unique_name == "R21_S11", legacy.obs_info
866 assert legacy.obs_info.physical_filter == "r_57", legacy.obs_info
869def test_aperture_corrections_to_legacy(legacy_test_data: _LegacyTestData) -> None:
870 """Verify that aperture corrections round-trip through a legacy
871 ApCorrMap.
872 """
873 ap_corrections = legacy_test_data.visit_image.aperture_corrections
874 legacy_ap_corr_map = aperture_corrections_to_legacy(ap_corrections)
875 compare_aperture_corrections_to_legacy(
876 ap_corrections,
877 legacy_ap_corr_map,
878 legacy_test_data.visit_image.bbox,
879 )
882def _check_legacy_headers(visit_image: VisitImage) -> None:
883 """Assert that primary and extension headers are stripped correctly."""
884 header = visit_image._opaque_metadata.headers[ExtensionKey()]
885 assert "EXPTIME" in header
886 assert header["PLATFORM"] == "lsstcam"
887 assert "LSST BUTLER ID" not in header
888 assert "AR HDU" not in header
889 assert "A_ORDER" not in header
890 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("IMAGE"), astropy.io.fits.Header())
891 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("MASK"), astropy.io.fits.Header())
892 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("VARIANCE"), astropy.io.fits.Header())
895def test_read_legacy_headers(legacy_test_data: _LegacyTestData) -> None:
896 """Verify that headers were stripped and interpreted correctly in
897 read_legacy.
898 """
899 assert legacy_test_data.visit_image.unit == legacy_test_data.unit
900 _check_legacy_headers(legacy_test_data.visit_image)
903def test_from_legacy_headers(legacy_test_data: _LegacyTestData) -> None:
904 """Verify that from_legacy handles primary and extension headers
905 correctly.
906 """
907 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map)
908 assert legacy.unit == legacy_test_data.unit
909 _check_legacy_headers(legacy)
912def test_rewrite(legacy_test_data: _LegacyTestData) -> None:
913 """Verify that a legacy VisitImage can be rewritten and round-trips both
914 pixel values and all components.
915 """
916 with RoundtripFits(legacy_test_data.visit_image, legacy_test_data.storage_class) as roundtrip:
917 fits = roundtrip.inspect()
918 assert fits[1].header["ZCMPTYPE"] == "RICE_1"
919 assert fits[1].header["CTYPE1"] == "RA---TAN-SIP"
920 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
921 assert fits[2].header["CTYPE1"] == "RA---TAN-SIP"
922 assert fits[3].header["ZCMPTYPE"] == "RICE_1"
923 assert fits[3].header["CTYPE1"] == "RA---TAN-SIP"
924 subbox = Box.factory[8:13, 9:30]
925 subimage = roundtrip.get(bbox=subbox)
926 assert_masked_images_equal(subimage, legacy_test_data.visit_image[subbox], expect_view=False)
927 alternates: dict[str, Any] = {}
928 assert roundtrip.get("bbox") == legacy_test_data.visit_image.bbox
929 alternates = {
930 k: roundtrip.get(k)
931 for k in [
932 "sky_projection",
933 "image",
934 "mask",
935 "variance",
936 "psf",
937 "obs_info",
938 "summary_stats",
939 "aperture_corrections",
940 "detector",
941 "photometric_scaling",
942 ]
943 }
944 legacy_exposure = roundtrip.get(storageClass="Exposure")
945 assert isinstance(legacy_exposure, LegacyExposure)
946 compare_visit_image_to_legacy(
947 legacy_test_data.visit_image,
948 legacy_exposure,
949 expect_view=False,
950 plane_map=legacy_test_data.plane_map,
951 **MINIMAL_VISIT_DATA_ID,
952 )
953 if legacy_test_data.visit_image.unit == u.nJy:
954 assert legacy_exposure.getPhotoCalib()._isConstant
955 assert legacy_exposure.getPhotoCalib().getCalibrationMean() == 1.0
956 else:
957 compare_photo_calib_to_legacy(
958 legacy_test_data.visit_image.photometric_scaling,
959 legacy_exposure.getPhotoCalib(),
960 subimage_bbox=subbox,
961 )
962 assert legacy_exposure.info.getId() == legacy_test_data.legacy_exposure.info.getId()
963 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
964 assert isinstance(visit_info, LegacyVisitInfo)
965 assert visit_info.getInstrumentLabel() == "LSSTCam"
967 assert_visit_images_equal(roundtrip.result, legacy_test_data.visit_image, expect_view=False)
968 assert dict(legacy_test_data.visit_image._opaque_metadata.headers[ExtensionKey()]) == dict(
969 roundtrip.result._opaque_metadata.headers[ExtensionKey()]
970 )
971 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")]
972 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")]
973 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")]
974 assert roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"] == "lsstcam"
975 compare_visit_image_to_legacy(
976 roundtrip.result,
977 legacy_test_data.legacy_exposure,
978 expect_view=False,
979 plane_map=legacy_test_data.plane_map,
980 **MINIMAL_VISIT_DATA_ID,
981 alternates=alternates,
982 )
983 compare_visit_image_to_legacy(
984 legacy_test_data.read_cls.from_legacy(
985 legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map
986 ),
987 legacy_test_data.legacy_exposure,
988 expect_view=True,
989 plane_map=legacy_test_data.plane_map,
990 **MINIMAL_VISIT_DATA_ID,
991 )
994def test_butler_converters(legacy_test_data: _LegacyTestData) -> None:
995 """Verify that a VisitImage can be read from a Butler dataset written as
996 an Exposure.
997 """
998 try:
999 from lsst.daf.butler import FileDataset
1000 except ImportError:
1001 pytest.skip("lsst.daf.butler could not be imported.")
1003 with TemporaryButler(legacy="ExposureF") as helper:
1004 helper.butler.ingest(
1005 FileDataset(path=legacy_test_data.filename, refs=[helper.legacy]), transfer="symlink"
1006 )
1007 visit_image_ref = helper.legacy.overrideStorageClass(legacy_test_data.storage_class)
1008 with warnings.catch_warnings():
1009 warnings.filterwarnings("ignore", message=".*filter label mismatch.*", category=UserWarning)
1010 visit_image = helper.butler.get(visit_image_ref)
1011 assert visit_image._opaque_metadata.precompressed.keys() == set()
1012 visit_image = helper.butler.get(visit_image_ref, parameters={"preserve_quantization": True})
1013 assert visit_image._opaque_metadata.precompressed.keys() == {"IMAGE", "VARIANCE"}
1014 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
1015 assert bbox == visit_image.bbox
1016 alternates = {
1017 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
1018 for k in ["image", "mask", "variance", "bbox", "psf", "detector"]
1019 }
1020 compare_visit_image_to_legacy(
1021 visit_image,
1022 legacy_test_data.legacy_exposure,
1023 expect_view=False,
1024 plane_map=legacy_test_data.plane_map,
1025 alternates=alternates,
1026 **MINIMAL_VISIT_DATA_ID,
1027 )
1028 helper.butler.pruneDatasets([helper.legacy], purge=True, unstore=True, disassociate=True)
1029 visit_image.metadata["MixedCaseKey"] = 52
1030 helper.butler.put(visit_image, visit_image_ref)
1031 with warnings.catch_warnings():
1032 warnings.filterwarnings("ignore", message=".*filter label mismatch.*", category=UserWarning)
1033 legacy_exposure = helper.butler.get(helper.legacy)
1034 compare_visit_image_to_legacy(
1035 visit_image,
1036 legacy_exposure,
1037 expect_view=False,
1038 plane_map=legacy_test_data.plane_map,
1039 alternates=alternates,
1040 **MINIMAL_VISIT_DATA_ID,
1041 )
1042 visit_image_2 = helper.butler.get(visit_image_ref)
1043 compare_visit_image_to_legacy(
1044 visit_image_2,
1045 legacy_exposure,
1046 expect_view=False,
1047 plane_map=legacy_test_data.plane_map,
1048 alternates=alternates,
1049 **MINIMAL_VISIT_DATA_ID,
1050 )
1051 assert visit_image_2.metadata["MixedCaseKey"] == 52
1054def test_convert_unit(legacy_test_data_calibrated: _LegacyTestData) -> None:
1055 """Verify convert_unit round-trips between nJy, mJy, and electron via
1056 photometric_scaling.
1057 """
1058 from lsst.afw.table import ExposureCatalog
1060 legacy_test_data = legacy_test_data_calibrated
1061 original = legacy_test_data.visit_image.copy()
1062 with pytest.raises(u.UnitConversionError):
1063 original.convert_unit(u.electron)
1064 visit_image_nJy = original.convert_unit(u.nJy, copy=False)
1065 assert np.may_share_memory(visit_image_nJy.image.array, original.image.array)
1066 assert np.may_share_memory(visit_image_nJy.variance.array, original.variance.array)
1067 with pytest.raises(u.UnitConversionError):
1068 original.convert_unit(u.mJy, copy=False)
1069 visit_image_mJy = original.convert_unit(u.mJy, copy="as-needed")
1070 assert visit_image_mJy.unit == u.mJy
1071 assert_values_equal(visit_image_mJy.image.array, original.image.array * 1e-6, rtol=1e-5)
1072 assert np.may_share_memory(visit_image_nJy.mask.array, original.mask.array)
1073 assert_values_equal(visit_image_mJy.variance.array, original.variance.array * 1e-12, rtol=1e-5)
1074 legacy_exposure_mJy = visit_image_mJy.to_legacy()
1075 assert_values_equal(legacy_exposure_mJy.getPhotoCalib().getCalibrationMean(), 1e6, rtol=1e-14)
1076 legacy_masked_image_nJy = legacy_exposure_mJy.getPhotoCalib().calibrateImage(
1077 legacy_exposure_mJy.maskedImage
1078 )
1079 assert_values_equal(visit_image_nJy.image.array, legacy_masked_image_nJy.image.array, rtol=1e-5)
1080 assert_values_equal(visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array, rtol=1e-5)
1081 assert np.may_share_memory(visit_image_mJy.mask.array, original.mask.array)
1082 assert visit_image_mJy.sky_projection is original.sky_projection
1083 assert visit_image_mJy.obs_info is original.obs_info
1084 assert visit_image_mJy.summary_stats is original.summary_stats
1085 assert visit_image_mJy.psf is original.psf
1086 assert visit_image_mJy.detector is original.detector
1087 assert visit_image_mJy.bounds is original.bounds
1088 assert visit_image_mJy.aperture_corrections is original.aperture_corrections
1089 assert visit_image_mJy.photometric_scaling is original.photometric_scaling
1090 visit_summary = ExposureCatalog.readFits(
1091 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_summary.fits")
1092 )
1093 legacy_photo_calib = visit_summary.find(DP2_VISIT_DETECTOR_DATA_ID["detector"]).getPhotoCalib()
1094 visit_image_nJy.photometric_scaling = field_from_legacy_photo_calib(
1095 legacy_photo_calib, bounds=original.detector.bbox, instrumental_unit=u.electron
1096 )
1097 compare_photo_calib_to_legacy(
1098 visit_image_nJy.photometric_scaling,
1099 legacy_test_data.legacy_exposure.getPhotoCalib(),
1100 applied_legacy_photo_calib=legacy_photo_calib,
1101 subimage_bbox=visit_image_nJy.bbox,
1102 )
1103 with pytest.raises(u.UnitConversionError):
1104 visit_image_nJy.convert_unit(u.mm)
1105 with pytest.raises(u.UnitConversionError):
1106 visit_image_nJy.convert_unit(u.electron, copy=False)
1107 legacy_masked_image_e = legacy_photo_calib.uncalibrateImage(legacy_test_data.legacy_exposure.maskedImage)
1108 visit_image_e = visit_image_nJy.convert_unit(u.electron)
1109 assert_values_equal(visit_image_e.image.array, legacy_masked_image_e.image.array, rtol=1e-5)
1110 assert_values_equal(visit_image_e.variance.array, legacy_masked_image_e.variance.array, rtol=1e-5)
1111 assert not np.may_share_memory(visit_image_e.mask.array, visit_image_nJy.mask.array)
1112 visit_image_mJy.photometric_scaling = visit_image_nJy.photometric_scaling
1113 visit_image_e = visit_image_mJy.convert_unit(u.electron)
1114 assert_values_equal(visit_image_e.image.array, legacy_masked_image_e.image.array, rtol=1e-5)
1115 assert_values_equal(visit_image_e.variance.array, legacy_masked_image_e.variance.array, rtol=1e-5)
1116 visit_image_nJy_2 = visit_image_e.convert_unit(u.nJy)
1117 assert_values_equal(visit_image_nJy_2.image.array, visit_image_nJy.image.array, rtol=1e-5)
1118 assert_values_equal(visit_image_nJy_2.variance.array, original.variance.array, rtol=1e-5)
1119 visit_image_e.photometric_scaling = visit_image_nJy.photometric_scaling * (1e-6 * u.mJy / u.nJy)
1120 visit_image_nJy_3 = visit_image_e.convert_unit(u.nJy)
1121 assert_values_equal(visit_image_nJy_3.image.array, visit_image_nJy.image.array, rtol=1e-5)
1122 assert_values_equal(visit_image_nJy_3.variance.array, original.variance.array, rtol=1e-5)
1123 legacy_exposure_e = visit_image_e.to_legacy()
1124 assert_values_equal(
1125 legacy_exposure_e.getPhotoCalib().getCalibrationMean(),
1126 legacy_photo_calib.getCalibrationMean(),
1127 rtol=1e-5,
1128 )
1129 legacy_masked_image_nJy = legacy_exposure_e.getPhotoCalib().calibrateImage(legacy_exposure_e.maskedImage)
1130 assert_values_equal(visit_image_nJy.image.array, legacy_masked_image_nJy.image.array, rtol=1e-5)
1131 assert_values_equal(visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array, rtol=1e-5)
1134def test_background_map_describe() -> None:
1135 """An empty BackgroundMap._describe reports no backgrounds inline."""
1136 bg_map = BackgroundMap()
1137 assert isinstance(bg_map, DescribableMixin)
1138 report = bg_map._describe()
1139 assert isinstance(report, Report)
1140 assert report.type_name == "BackgroundMap"
1141 assert report.inline
1142 assert report.summary == "no backgrounds"
1143 assert report.children == {}
1146def test_background_map_with_entries_describe() -> None:
1147 """BackgroundMap._describe recurses into each background model."""
1148 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1149 bg_map = BackgroundMap(
1150 [Background("sky", cheby, "Sky model."), Background("fringe", cheby)],
1151 subtracted="sky",
1152 )
1153 report = bg_map._describe()
1154 assert report.type_name == "BackgroundMap"
1155 assert report.inline
1156 # The inline summary leads with the subtracted background and marks it,
1157 # since that is all a composite holding this map will show.
1158 # Only the subtracted background is marked, not the whole line.
1159 assert report.summary == "**sky [SUBTRACTED]**; also: fringe"
1160 assert report.emphasis_markup
1161 # str reads the summary as text, so the delimiters never reach a user.
1162 assert report.to_str() == "sky [SUBTRACTED]; also: fringe"
1163 # Standalone, each background is a child carrying its own model's report.
1164 assert set(report.children) == {"sky", "fringe"}
1165 sky = report.children["sky"]
1166 assert sky.type_name == "ChebyshevField"
1167 # The marker is on the heading, where it is read before the fields.
1168 assert sky.title == "ChebyshevField **[SUBTRACTED]**"
1169 assert sky.emphasis_markup
1170 sky_fields = {f.label: f.value for f in sky.fields}
1171 assert sky_fields["description"] == "Sky model."
1172 # The model's own fields survive alongside the background's attributes.
1173 assert "bounds" in sky_fields
1174 # Only the subtracted one is marked, and an absent description is omitted.
1175 fringe = report.children["fringe"]
1176 assert fringe.title is None
1177 assert fringe.emphasis_markup is False
1178 assert "description" not in {f.label for f in fringe.fields}
1181def test_background_map_describe_none_subtracted() -> None:
1182 """A map with no subtracted background says so rather than being silent."""
1183 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1184 bg_map = BackgroundMap([Background("sky", cheby), Background("fringe", cheby)])
1185 report = bg_map._describe()
1186 assert report.summary == "sky, fringe (none subtracted)"
1187 assert report.emphasis_markup is False
1188 assert not any(child.emphasis_markup for child in report.children.values())
1191def test_background_map_describe_leads_with_subtracted() -> None:
1192 """The subtracted background leads the summary wherever it sits in the
1193 map, and only its child is marked.
1194 """
1195 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1196 bg_map = BackgroundMap(
1197 [Background("sky", cheby), Background("skyCorr", cheby)],
1198 subtracted="skyCorr",
1199 )
1200 report = bg_map._describe()
1201 assert report.to_str() == "skyCorr [SUBTRACTED]; also: sky"
1202 assert report.children["sky"].emphasis_markup is False
1203 assert report.children["skyCorr"].emphasis_markup
1206def test_background_map_describe_brief_skips_children() -> None:
1207 """A brief background map report keeps the summary but not the models."""
1208 cheby = ChebyshevField(Box.factory[0:100, 0:200], np.array([[1.0]]))
1209 bg_map = BackgroundMap([Background("sky", cheby)], subtracted="sky")
1210 report = bg_map._describe(DescribeOptions(brief=True))
1211 assert report.to_str() == "sky [SUBTRACTED]"
1212 assert report.children == {}
1215def test_visit_image_repr_str_with_unreadable_psf() -> None:
1216 """Repr and str succeed even when the PSF stored an ArchiveReadError.
1218 An unreadable component is a supported state; repr and str read only the
1219 cheap fields and summary, so they do not build the child tree at all.
1220 """
1221 path = current_fixture_path(FIXTURE_DIR, "visit_image")
1222 visit_image = read_archive(path)
1223 visit_image._psf = ArchiveReadError("psf unreadable")
1224 assert repr(visit_image).startswith("VisitImage(")
1225 assert str(visit_image).startswith("VisitImage(")
1228def test_visit_image_describe_names_an_unreadable_psf() -> None:
1229 """A full report says the PSF could not be read, and describes the rest.
1231 A PSF model can need a package the reader does not have installed, so an
1232 unreadable PSF must not cost the report of everything else.
1233 """
1234 path = current_fixture_path(FIXTURE_DIR, "visit_image")
1235 visit_image = read_archive(path)
1236 visit_image._psf = ArchiveReadError("Failed to import piff.")
1237 report = visit_image.describe(detail=True)
1238 psf = report.children["psf"]
1239 assert psf.inline
1240 assert psf.to_str() == "unreadable (Failed to import piff.)"
1241 # Every other component is described as usual.
1242 assert "sky_projection" in report.children
1243 assert "detector" in report.children
1244 # Both renderers run over the report that contains it.
1245 assert "unreadable" in report._repr_html_()
1246 report.__rich__()
1249def test_visit_image_slice_preserves_unreadable_psf() -> None:
1250 """Slicing propagates a deferred PSF read failure without raising it."""
1251 path = current_fixture_path(FIXTURE_DIR, "visit_image")
1252 visit_image = read_archive(path)
1253 error = ArchiveReadError("psf unreadable")
1254 visit_image._psf = error
1256 sliced = visit_image[...]
1258 assert sliced._psf is error
1259 with pytest.raises(ArchiveReadError, match="psf unreadable"):
1260 sliced.psf
1263def test_observation_summary_stats_describe() -> None:
1264 """ObservationSummaryStats._describe reports the statistics that are set.
1266 Unset ones are omitted rather than shown as NaN.
1267 """
1268 stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4, ra=180.0, dec=-30.0)
1269 assert isinstance(stats, DescribableMixin)
1270 report = stats._describe()
1271 assert isinstance(report, Report)
1272 assert report.type_name == "ObservationSummaryStats"
1273 # Scalars are packed into one group, ordered by name.
1274 (group,) = report.value_groups
1275 assert group.role is FieldRole.DERIVED
1276 assert [name for name, _ in group.values] == sorted(name for name, _ in group.values)
1277 values = dict(group.values)
1278 assert values["psfSigma"] == 2.5
1279 assert values["zeroPoint"] == 31.4
1280 assert values["ra"] == 180.0
1281 assert values["dec"] == -30.0
1282 # Unset statistics are omitted rather than reported as NaN, and the
1283 # serialization plumbing this class inherits never appears at all.
1284 assert "expTime" not in values
1285 assert "skyBg" not in values
1286 assert not {"metadata", "butler_info", "schema_version"} & set(values)
1289def test_observation_summary_stats_describe_omits_empty_sequences() -> None:
1290 """Sequence statistics appear only when they carry a value."""
1291 empty = ObservationSummaryStats(psfSigma=2.5)
1292 # raCorners defaults to all-NaN, which carries no more information than an
1293 # empty sequence does.
1294 assert all(math.isnan(v) for v in empty.raCorners)
1295 assert not any(f.label == "raCorners" for f in empty._describe().fields)
1297 filled = ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
1298 corners = next(f for f in filled._describe().fields if f.label == "raCorners")
1299 assert corners.value == (5.2, 5.4, 5.4, 5.2)
1300 assert corners.role is FieldRole.DERIVED
1303def test_observation_summary_stats_describe_brief_counts() -> None:
1304 """A brief report gives the number set rather than listing them."""
1305 stats = ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
1306 brief = stats._describe(DescribeOptions(brief=True))
1307 assert brief.type_name == "ObservationSummaryStats"
1308 assert brief.value_groups == []
1309 (field,) = brief.fields
1310 assert field.label == "statistics set"
1311 assert field.role is FieldRole.DERIVED
1313 # The count must agree with what the full report actually shows: two
1314 # scalars set by default plus psfSigma, and raCorners as a sequence.
1315 full = stats._describe()
1316 listed = len(full.fields) + sum(len(g.values) for g in full.value_groups)
1317 assert field.value.startswith(f"{listed} of ")
1318 assert brief.to_str() == full.to_str()
1319 assert f"{listed} of " in brief.to_str()
1322def test_observation_summary_stats_pydantic_repr() -> None:
1323 """ObservationSummaryStats uses pydantic's repr, not the mixin's."""
1324 stats = ObservationSummaryStats(psfSigma=2.5)
1325 r = repr(stats)
1326 assert r.startswith("ObservationSummaryStats(")
1327 assert "psfSigma=2.5" in r
1330def test_observation_summary_stats_str_is_the_report_summary() -> None:
1331 """The str output reports the count, not pydantic's field-by-field dump.
1333 The repr keeps the exhaustive form, since that is the one that
1334 round-trips.
1335 """
1336 stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
1337 assert str(stats) == stats.describe().to_str()
1338 # The two set here, plus the integer counters, which default to a genuine
1339 # zero rather than to NaN.
1340 assert str(stats) == "ObservationSummaryStats(4 of 66 statistics set)"
1341 assert stats.nPsfStar == 0
1342 assert stats.nShapeletsStar == 0
1343 # The unset statistics reach repr but not str.
1344 assert "nan" not in str(stats)
1345 assert "nan" in repr(stats)
1348def test_archive_tree_repr_omits_schema_bookkeeping() -> None:
1349 """Schema version fields mirror class constants, so repr leaves them out.
1351 They are never passed on construction and say nothing the type does not.
1352 """
1353 stats = ObservationSummaryStats(psfSigma=2.5, indirect=[1])
1354 for name in ("schema_version", "min_read_version", "indirect"):
1355 assert name not in repr(stats), name
1356 # They are still real fields, and hiding them from repr does not hide them
1357 # from serialization.
1358 assert stats.schema_version == ObservationSummaryStats.SCHEMA_VERSION
1359 dumped = stats.model_dump()
1360 for name in ("schema_version", "min_read_version", "indirect"):
1361 assert name in dumped, name
1364def _match_detector_to_image(visit_image: VisitImage) -> VisitImage:
1365 """Return the visit image with its detector shrunk to the image bbox.
1367 Parameters
1368 ----------
1369 visit_image
1370 Image to adjust.
1372 Returns
1373 -------
1374 visit_image : `VisitImage`
1375 Image whose `~VisitImage.bbox` equals its detector's, as a full-frame
1376 image read from a butler has.
1377 """
1378 detector = visit_image.detector
1379 visit_image._detector = Detector(
1380 detector._attributes.model_copy(update={"bbox": visit_image.bbox}),
1381 detector.amplifiers,
1382 detector._frames,
1383 detector.visit,
1384 )
1385 return visit_image
1388def test_visit_image_describe_hides_summary_stat_corners_for_a_full_frame(
1389 visit_image_components: dict[str, Any],
1390) -> None:
1391 """The sky corners appear once, in the projection's readable table."""
1392 visit = _match_detector_to_image(make_visit_image(visit_image_components))
1393 visit.summary_stats.raCorners = (1.0, 2.0, 3.0, 4.0)
1394 visit.summary_stats.decCorners = (-1.0, -2.0, -3.0, -4.0)
1395 report = visit.describe()
1396 stats = report.children["summary_stats"]
1397 assert not {"raCorners", "decCorners"} & {field.label for field in stats.fields}
1398 # The projection is where they are read instead, and it still has them.
1399 assert any(table.title == "Corners" for table in report.children["sky_projection"].tables)
1402def test_visit_image_cutout_keeps_summary_stat_corners() -> None:
1403 """A cutout keeps the statistics' corners, which still describe the whole
1404 image the statistics were measured on.
1405 """
1406 # The DP2 variant is a cutout of a real image, so it carries the
1407 # statistics measured on the whole of it.
1408 path = current_fixture_path(FIXTURE_DIR, "visit_image", variant="dp2")
1409 visit_image = read_archive(path)
1410 assert visit_image.bbox != visit_image.detector.bbox
1411 stats = visit_image.describe().children["summary_stats"]
1412 assert {"raCorners", "decCorners"} <= {field.label for field in stats.fields}