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