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