Coverage for tests/test_visit_image.py: 57%
495 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 15:24 -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 typing import Any, Literal
20import astropy.io.fits
21import astropy.units as u
22import astropy.wcs
23import numpy as np
24import pytest
25from astro_metadata_translator import ObservationInfo
27from lsst.images import (
28 BackgroundMap,
29 Box,
30 DetectorFrame,
31 DifferenceImage,
32 Image,
33 MaskPlane,
34 MaskSchema,
35 ObservationSummaryStats,
36 Polygon,
37 SkyProjectionAstropyView,
38 TractFrame,
39 VisitImage,
40 get_legacy_difference_image_mask_planes,
41 get_legacy_visit_image_mask_planes,
42)
43from lsst.images.aperture_corrections import ApertureCorrectionMap, aperture_corrections_to_legacy
44from lsst.images.cameras import Detector
45from lsst.images.fields import ChebyshevField, SplineField, SumField, field_from_legacy_photo_calib
46from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata
47from lsst.images.psfs import GaussianPointSpreadFunction, PointSpreadFunction
48from lsst.images.serialization import read_archive
49from lsst.images.tests import (
50 DP2_VISIT_DETECTOR_DATA_ID,
51 RoundtripFits,
52 RoundtripJson,
53 RoundtripNdf,
54 TemporaryButler,
55 assert_close,
56 assert_masked_images_equal,
57 assert_sky_projections_equal,
58 assert_visit_images_equal,
59 compare_aperture_corrections_to_legacy,
60 compare_detector_to_legacy,
61 compare_photo_calib_to_legacy,
62 compare_visit_image_to_legacy,
63 make_random_sky_projection,
64)
66try:
67 import h5py # noqa: F401
69 HAVE_H5PY = True
70except ImportError:
71 HAVE_H5PY = False
73try:
74 from lsst.afw.image import Exposure as LegacyExposure
75 from lsst.afw.image import VisitInfo as LegacyVisitInfo
76except ImportError:
77 type LegacyExposure = Any # type: ignore[no-redef]
78 type LegacyVisitInfo = Any # type: ignore[no-redef]
80EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
81LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
83skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed")
86@pytest.fixture(scope="session")
87def visit_image_components() -> dict[str, Any]:
88 """Return a dictionary of VisitImage components."""
89 rng = np.random.default_rng(500)
90 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096])
91 mask_schema = MaskSchema([MaskPlane("M1", "D1")])
92 obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1")
93 summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
94 gaussian_psf = GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13])
95 aperture_corrections: ApertureCorrectionMap = {
96 "flux1": ChebyshevField(det_frame.bbox, np.array([0.75])),
97 "flux2": ChebyshevField(det_frame.bbox, np.array([0.625])),
98 }
99 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector)
100 image = Image(42, shape=(1024, 1024), unit=u.nJy)
101 variance = Image(5.0, shape=(1024, 1024), unit=u.nJy * u.nJy)
102 # polygon is the lower triangle of the image.
103 polygon = Polygon(x_vertices=[-0.5, 1023.5, -0.5], y_vertices=[-0.5, -0.5, 1023.5])
104 sky_projection = make_random_sky_projection(rng, det_frame, det_frame.bbox)
105 return {
106 "mask_schema": mask_schema,
107 "obs_info": obs_info,
108 "summary_stats": summary_stats,
109 "gaussian_psf": gaussian_psf,
110 "aperture_corrections": aperture_corrections,
111 "detector": detector,
112 "image": image,
113 "variance": variance,
114 "polygon": polygon,
115 "sky_projection": sky_projection,
116 }
119def make_visit_image(components: dict[str, Any]) -> VisitImage:
120 """Construct a new VisitImage with most components populated."""
121 det_frame = components["sky_projection"].pixel_frame
122 opaque = FitsOpaqueMetadata()
123 hdr = astropy.io.fits.Header()
124 with warnings.catch_warnings():
125 # Silence warnings about long keys becoming HIERARCH.
126 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning)
127 hdr.update({"PLATFORM": "lsstcam", "LSST BUTLER ID": "123456789"})
128 opaque.extract_legacy_primary_header(hdr)
129 # API signature suggests sky_projection and obs_info can be None but
130 # they are required (unless you pass them in via the image plane).
131 vi = VisitImage(
132 components["image"],
133 variance=components["variance"],
134 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
135 mask_schema=components["mask_schema"],
136 sky_projection=components["sky_projection"],
137 obs_info=components["obs_info"],
138 summary_stats=components["summary_stats"],
139 detector=components["detector"],
140 bounds=components["polygon"],
141 aperture_corrections=components["aperture_corrections"],
142 band="r",
143 )
144 vi.backgrounds.add(
145 "standard",
146 ChebyshevField(det_frame.bbox, np.array([[2.0]])),
147 description="Background subtracted from the image.",
148 is_subtracted=True,
149 )
150 vi._opaque_metadata = opaque
151 return vi
154def make_simplest_visit_image(components: dict[str, Any]) -> VisitImage:
155 """Construct a VisitImage with the minimal set of components populated."""
156 return VisitImage(
157 components["image"],
158 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
159 mask_schema=components["mask_schema"],
160 sky_projection=components["sky_projection"],
161 detector=components["detector"],
162 obs_info=components["obs_info"],
163 band="r",
164 )
167def _make_sum_background_visit_image(components: dict[str, Any], visit_image: VisitImage) -> VisitImage:
168 """Return a VisitImage whose subtracted background is a SumField."""
169 rng = np.random.default_rng(42)
170 bbox = visit_image.image.sky_projection.pixel_frame.bbox
171 bin_y = bbox.y.linspace(6)
172 bin_x = bbox.x.linspace(7)
173 spline_a = SplineField(
174 bbox,
175 rng.standard_normal(size=(bin_y.size, bin_x.size)),
176 y=bin_y,
177 x=bin_x,
178 )
179 spline_b = SplineField(
180 bbox,
181 rng.standard_normal(size=(bin_y.size, bin_x.size)),
182 y=bin_y,
183 x=bin_x,
184 )
185 sum_field = SumField([spline_a, spline_b])
186 bg_map = BackgroundMap()
187 bg_map.add(
188 "stacked",
189 sum_field,
190 description="Two-operand SumField subtracted background.",
191 is_subtracted=True,
192 )
193 return VisitImage(
194 components["image"],
195 variance=components["variance"],
196 psf=components["gaussian_psf"],
197 mask_schema=components["mask_schema"],
198 sky_projection=components["sky_projection"],
199 obs_info=components["obs_info"],
200 summary_stats=components["summary_stats"],
201 detector=components["detector"],
202 band="r",
203 backgrounds=bg_map,
204 )
207def _check_sum_background_round_trip(result: VisitImage, original: VisitImage) -> None:
208 """Assert that a round-tripped SumField background matches the original."""
209 subtracted = result.backgrounds.subtracted
210 assert subtracted is not None
211 assert isinstance(subtracted.field, SumField)
212 original_subtracted = original.backgrounds.subtracted
213 assert original_subtracted is not None
214 original_field = original_subtracted.field
215 assert isinstance(original_field, SumField)
216 round_field = subtracted.field
217 assert isinstance(round_field, SumField)
218 assert len(round_field.operands) == len(original_field.operands)
219 for round_op, orig_op in zip(round_field.operands, original_field.operands, strict=True):
220 assert round_op == orig_op
223def test_basics(visit_image_components: dict[str, Any]) -> None:
224 """Verify VisitImage constructor patterns and required-argument checks."""
225 c = visit_image_components
226 # Test default fill of variance.
227 visit = make_simplest_visit_image(c)
228 assert visit.variance.array[0, 0] == 1.0
229 assert visit[...] is not visit
230 assert str(visit) == "VisitImage(Image([y=0:1024, x=0:1024], int64), ['M1'])"
231 assert repr(visit) == (
232 "VisitImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
233 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
234 " dtype=dtype('uint8')))"
235 )
237 astropy_wcs = visit.astropy_wcs
238 assert isinstance(astropy_wcs, SkyProjectionAstropyView)
239 approx_wcs = visit.fits_wcs
240 assert isinstance(approx_wcs, astropy.wcs.WCS)
242 with pytest.raises(TypeError):
243 # Requires a PSF.
244 VisitImage(
245 c["image"],
246 mask_schema=c["mask_schema"],
247 sky_projection=c["sky_projection"],
248 obs_info=c["obs_info"],
249 detector=c["detector"],
250 band="r",
251 )
253 with pytest.raises(TypeError):
254 # Requires ObservationInfo.
255 VisitImage(
256 c["image"],
257 psf=c["gaussian_psf"],
258 mask_schema=c["mask_schema"],
259 sky_projection=c["sky_projection"],
260 detector=c["detector"],
261 band="r",
262 )
264 with pytest.raises(TypeError):
265 # Requires a sky_projection.
266 VisitImage(
267 c["image"],
268 psf=c["gaussian_psf"],
269 mask_schema=c["mask_schema"],
270 obs_info=c["obs_info"],
271 detector=c["detector"],
272 band="r",
273 )
275 with pytest.raises(TypeError):
276 # Requires a detector.
277 VisitImage(
278 c["image"],
279 psf=c["gaussian_psf"],
280 mask_schema=c["mask_schema"],
281 sky_projection=c["sky_projection"],
282 obs_info=c["obs_info"],
283 band="r",
284 )
286 with pytest.raises(TypeError):
287 # Requires some form of mask.
288 VisitImage(
289 c["image"],
290 psf=c["gaussian_psf"],
291 sky_projection=c["sky_projection"],
292 obs_info=c["obs_info"],
293 detector=c["detector"],
294 band="r",
295 )
297 with pytest.raises(TypeError):
298 VisitImage(
299 Image(42, shape=(5, 5)),
300 psf=c["gaussian_psf"],
301 mask_schema=c["mask_schema"],
302 sky_projection=c["sky_projection"],
303 obs_info=c["obs_info"],
304 detector=c["detector"],
305 band="r",
306 )
308 # Requires a DetectorFrame.
309 rng = np.random.default_rng(501)
310 tract_frame = TractFrame(skymap="Skymap", tract=1, bbox=Box.factory[1:10, 1:10])
311 tract_proj = make_random_sky_projection(rng, tract_frame, Box.factory[1:4096, 1:4096])
312 with pytest.raises(TypeError):
313 VisitImage(
314 c["image"],
315 sky_projection=tract_proj,
316 psf=c["gaussian_psf"],
317 mask_schema=c["mask_schema"],
318 obs_info=c["obs_info"],
319 detector=c["detector"],
320 band="r",
321 )
323 # Variance unit mismatch.
324 with pytest.raises(ValueError):
325 VisitImage(
326 c["image"],
327 variance=c["image"],
328 psf=c["gaussian_psf"],
329 mask_schema=c["mask_schema"],
330 sky_projection=c["sky_projection"],
331 obs_info=c["obs_info"],
332 detector=c["detector"],
333 band="r",
334 )
337def test_copy_and_slice(visit_image_components: dict[str, Any]) -> None:
338 """Verify that copy deep-copies arrays and components while slice shares
339 them.
340 """
341 c = visit_image_components
342 visit_image = make_visit_image(c)
343 copy = visit_image.copy()
344 copy.image.array[0, 0] = 30.0
345 assert visit_image.image.array[0, 0] == 42.0
346 assert copy.image.array[0, 0] == 30.0
347 subvisit = visit_image[Box.factory[0:5, 0:5]]
348 # Check summary stats.
349 assert copy.summary_stats == visit_image.summary_stats
350 assert copy.summary_stats is not visit_image.summary_stats
351 assert subvisit.summary_stats == visit_image.summary_stats
352 assert subvisit.summary_stats is visit_image.summary_stats
353 # Check aperture corrections.
354 assert copy.aperture_corrections.keys() == visit_image.aperture_corrections.keys()
355 assert copy.aperture_corrections is not visit_image.aperture_corrections
356 assert subvisit.aperture_corrections.keys() == visit_image.aperture_corrections.keys()
357 assert subvisit.aperture_corrections is visit_image.aperture_corrections
358 # Check backgrounds.
359 assert copy.backgrounds.keys() == visit_image.backgrounds.keys()
360 assert copy.backgrounds is not visit_image.backgrounds
361 assert subvisit.backgrounds.keys() == visit_image.backgrounds.keys()
362 assert subvisit.backgrounds is visit_image.backgrounds
363 # Check bounds.
364 assert copy.bounds is c["polygon"]
365 assert subvisit.bounds == subvisit.bbox # original polygon wholly encloses subvisit.bbox
368def test_obs_info(visit_image_components: dict[str, Any]) -> None:
369 """Verify that ObservationInfo is present and carries the expected
370 instrument.
371 """
372 visit_image = make_visit_image(visit_image_components)
373 assert visit_image.obs_info is not None
374 assert visit_image.obs_info.instrument == "LSSTCam"
377def test_summary_stats(visit_image_components: dict[str, Any]) -> None:
378 """Verify ObservationSummaryStats equality and inequality comparisons."""
379 summary_stats = visit_image_components["summary_stats"]
380 assert summary_stats == ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
381 assert summary_stats != ObservationSummaryStats(psfSigma=2.5)
382 assert summary_stats != ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
385def test_summary_stats_to_legacy() -> None:
386 """Verify ObservationSummaryStats round-trips through the legacy
387 ExposureSummaryStats even when this package defines fields that the
388 installed afw does not.
389 """
390 try:
391 from lsst.afw.image import ExposureSummaryStats
392 except ImportError:
393 pytest.skip("lsst.afw.image is not available")
395 summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
396 legacy = summary_stats.to_legacy()
397 assert isinstance(legacy, ExposureSummaryStats)
398 assert legacy.psfSigma == 2.5
399 assert legacy.zeroPoint == 31.4
400 # Empty (NaN) fields unknown to the legacy struct are dropped, so the
401 # round trip reproduces the original.
402 assert ObservationSummaryStats.from_legacy(legacy) == summary_stats
404 # A real value in a field unknown to the legacy struct cannot be
405 # represented and must raise rather than be silently dropped.
406 legacy_fields = {field.name for field in dataclasses.fields(ExposureSummaryStats)}
407 extra_float_fields = [
408 name
409 for name, info in ObservationSummaryStats.model_fields.items()
410 if name not in legacy_fields and info.annotation is float
411 ]
412 if extra_float_fields: 412 ↛ exitline 412 didn't return from function 'test_summary_stats_to_legacy' because the condition on line 412 was always true
413 with pytest.raises(ValueError):
414 ObservationSummaryStats(**{extra_float_fields[0]: 1.0}).to_legacy()
417def test_summary_stats_from_legacy_unknown_field() -> None:
418 """Verify from_legacy drops empty unknown fields but errors on set ones."""
420 @dataclasses.dataclass
421 class FakeLegacy:
422 psfSigma: float = 2.5
423 notARealField: float = math.nan
425 # An unknown field that is empty is dropped.
426 stats = ObservationSummaryStats.from_legacy(FakeLegacy())
427 assert stats.psfSigma == 2.5
429 # An unknown field that holds a real value cannot be represented.
430 with pytest.raises(ValueError):
431 ObservationSummaryStats.from_legacy(FakeLegacy(notARealField=1.0))
434@skip_no_h5py
435def test_round_trip_ndf(visit_image_components: dict[str, Any]) -> None:
436 """Verify NDF round-trip produces a VisitImage equal to the original."""
437 visit_image = make_visit_image(visit_image_components)
438 with RoundtripNdf(visit_image, "VisitImage") as roundtrip:
439 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False)
442@skip_no_h5py
443def test_fits_ndf_consistency(visit_image_components: dict[str, Any]) -> None:
444 """Verify FITS and NDF backends produce equal VisitImages on round-trip."""
445 visit_image = make_visit_image(visit_image_components)
446 with RoundtripFits(visit_image) as fits_rt, RoundtripNdf(visit_image) as ndf_rt:
447 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False)
448 assert_visit_images_equal(visit_image, ndf_rt.result, expect_view=False)
449 assert_visit_images_equal(fits_rt.result, ndf_rt.result, expect_view=False)
452def test_fits_json_consistency(visit_image_components: dict[str, Any]) -> None:
453 """Verify FITS and JSON backends produce equal VisitImages."""
454 visit_image = make_visit_image(visit_image_components)
455 with (
456 RoundtripFits(visit_image) as fits_rt,
457 RoundtripJson(visit_image) as json_rt,
458 ):
459 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False)
460 assert_visit_images_equal(visit_image, json_rt.result, expect_view=False)
461 assert_visit_images_equal(fits_rt.result, json_rt.result, expect_view=False)
464def test_read_write(visit_image_components: dict[str, Any]) -> None:
465 """Verify a VisitImage round-trips through FITS with correct compression.
467 Checks compression headers, subimage reads, equality, and opaque
468 metadata. Contains only butler-free assertions; component reads live
469 in `test_read_write_components`.
470 """
471 visit_image = make_visit_image(visit_image_components)
472 with RoundtripFits(visit_image, "VisitImage") as roundtrip:
473 # Check that we're still using the right compression, and that we
474 # wrote WCSs.
475 fits = roundtrip.inspect()
476 assert fits[1].header["ZCMPTYPE"] == "GZIP_2"
477 assert fits[1].header["CTYPE1"] == "RA---TAN"
478 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
479 assert fits[2].header["CTYPE1"] == "RA---TAN"
480 assert fits[3].header["ZCMPTYPE"] == "GZIP_2"
481 assert fits[3].header["CTYPE1"] == "RA---TAN"
482 # Check a subimage read (no component arg — does not trigger a skip).
483 subbox = Box.factory[8:13, 9:30]
484 subimage = roundtrip.get(bbox=subbox)
485 assert_masked_images_equal(subimage, visit_image[subbox], expect_view=False)
487 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False)
488 # Check that the round-tripped headers are the same (up to card order).
489 assert len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]) == 1
490 assert dict(visit_image._opaque_metadata.headers[ExtensionKey()]) == dict(
491 roundtrip.result._opaque_metadata.headers[ExtensionKey()]
492 )
493 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")]
494 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")]
495 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")]
496 # Spot-check the concrete background contents (names, field types,
497 # subtracted entry) against the known fixture, so the equality check
498 # above is not vacuously satisfied by empty background maps.
499 assert isinstance(roundtrip.result.backgrounds, BackgroundMap)
500 assert roundtrip.result.backgrounds.keys() == {"standard"}
501 assert isinstance(roundtrip.result.backgrounds["standard"].field, ChebyshevField)
502 assert roundtrip.result.backgrounds.subtracted.name == "standard"
503 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image."
506def test_read_write_components(visit_image_components: dict[str, Any]) -> None:
507 """Verify component reads and storage-class overrides round-trip correctly.
509 Requires a butler; skips when `lsst.daf.butler` is absent.
510 Butler-free assertions live in `test_read_write`.
511 """
512 c = visit_image_components
513 visit_image = make_visit_image(c)
514 with RoundtripFits(visit_image, "VisitImage") as roundtrip:
515 subbox = Box.factory[8:13, 9:30]
516 subimage = roundtrip.get(bbox=subbox)
518 # Get an explicit masked image to compare with the subimage.
519 subimage_masked = roundtrip.get("masked_image", bbox=subbox)
520 assert_masked_images_equal(subimage_masked, subimage, expect_view=False)
522 # Get the same masked image in a multi-component get and ensure
523 # it is the same thing.
524 components = roundtrip.get("components", components=["masked_image", "psf"], bbox=subbox)
525 assert set(components) == {"masked_image", "psf"}
526 assert_masked_images_equal(components["masked_image"], subimage_masked, expect_view=False)
528 assert roundtrip.get("bbox") == visit_image.bbox
530 obs_info = roundtrip.get("obs_info")
531 assert isinstance(obs_info, ObservationInfo)
532 assert obs_info == visit_image.obs_info
534 summary_stats = roundtrip.get("summary_stats")
535 assert isinstance(summary_stats, ObservationSummaryStats)
536 assert summary_stats == visit_image.summary_stats
538 psf = roundtrip.get("psf")
539 assert isinstance(psf, GaussianPointSpreadFunction)
540 assert psf.kernel_bbox == c["gaussian_psf"].kernel_bbox
542 backgrounds = roundtrip.get("backgrounds")
543 assert isinstance(backgrounds, BackgroundMap)
544 assert backgrounds.keys() == {"standard"}
545 assert isinstance(backgrounds["standard"].field, ChebyshevField)
546 assert backgrounds.subtracted.name == "standard"
547 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image."
549 # Test some components get edge cases.
550 components = roundtrip.get("components", components="image")
551 assert isinstance(components["image"], Image)
553 components = roundtrip.get("components")
554 assert set(components) == {
555 "image",
556 "variance",
557 "psf",
558 "bbox",
559 "mask",
560 "obs_info",
561 "backgrounds",
562 "detector",
563 "aperture_corrections",
564 "sky_projection",
565 "summary_stats",
566 "photometric_scaling",
567 }
569 # Butler morphs RuntimeError to ValueError.
570 with pytest.raises(ValueError):
571 roundtrip.get("components", components=["image", "nonexistent"])
573 with pytest.raises(ValueError):
574 roundtrip.get("components", components=["image", "components"])
576 with pytest.raises(ValueError):
577 roundtrip.get("components", components=[])
579 with pytest.raises(ValueError):
580 # PSF does not know how to use bbox so this fails.
581 roundtrip.get("components", components="psf", bbox=subbox)
584def test_sum_background_round_trip_fits(visit_image_components: dict[str, Any]) -> None:
585 """Verify FITS backend keeps two same-named SumField operands as distinct
586 EXTVERs.
587 """
588 visit_image = make_visit_image(visit_image_components)
589 visit = _make_sum_background_visit_image(visit_image_components, visit_image)
590 with RoundtripFits(visit) as roundtrip:
591 _check_sum_background_round_trip(roundtrip.result, visit)
594@skip_no_h5py
595def test_sum_background_round_trip_ndf(visit_image_components: dict[str, Any]) -> None:
596 """Verify NDF backend disambiguates the repeated ``data`` leaf, just as
597 the FITS backend does.
598 """
599 visit_image = make_visit_image(visit_image_components)
600 visit = _make_sum_background_visit_image(visit_image_components, visit_image)
601 with RoundtripNdf(visit) as roundtrip:
602 _check_sum_background_round_trip(roundtrip.result, visit)
605@dataclasses.dataclass
606class _LegacyTestData:
607 filename: str
608 plane_map: dict[str, MaskPlane] = dataclasses.field(default_factory=get_legacy_visit_image_mask_planes)
609 unit: u.Unit = u.nJy
610 storage_class: str = "VisitImage"
611 read_cls: type[VisitImage] = VisitImage
612 legacy_exposure: LegacyExposure = dataclasses.field(init=False)
614 @classmethod
615 def get(
616 cls, which: Literal["visit_image", "preliminary_visit_image", "difference_image"]
617 ) -> _LegacyTestData:
618 if EXTERNAL_DATA_DIR is None: 618 ↛ 620line 618 didn't jump to line 620 because the condition on line 618 was always true
619 pytest.skip("TESTDATA_IMAGES is not set up.")
620 result = cls(
621 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", f"{which}.fits"),
622 )
623 match which:
624 case "preliminary_visit_image":
625 result.unit = u.electron
626 case "difference_image":
627 result.storage_class = "DifferenceImage"
628 result.read_cls = DifferenceImage
629 result.plane_map = get_legacy_difference_image_mask_planes()
630 case "visit_image":
631 pass
632 try:
633 from lsst.afw.image import ExposureFitsReader
635 result.legacy_exposure = ExposureFitsReader(result.filename).read()
636 except ImportError:
637 pytest.skip("lsst.afw.image is not available; cannot read legacy exposures")
638 result.visit_image = result.read_cls.read_legacy(
639 result.filename, preserve_quantization=True, plane_map=result.plane_map
640 )
641 return result
644@pytest.fixture(scope="session", params=["visit_image", "preliminary_visit_image", "difference_image"])
645def legacy_test_data(request: pytest.FixtureRequest) -> _LegacyTestData:
646 """Return legacy test data.
648 Tests that depend on this parameterized fixture run on all of the legacy
649 test images.
650 """
651 return _LegacyTestData.get(request.param)
654@pytest.fixture(scope="session", params=["visit_image", "difference_image"])
655def legacy_test_data_calibrated(request: pytest.FixtureRequest) -> _LegacyTestData:
656 """Return legacy test data for calibrated images only.
658 Tests that depend on this parameterized fixture do not run on
659 preliminary_visit_image, since that has 'electron' pixel units
660 """
661 return _LegacyTestData.get(request.param)
664def _check_legacy_obs_info(obs_info: ObservationInfo | None) -> None:
665 """Assert obs_info carries expected LSSTCam/DP2 field values."""
666 assert isinstance(obs_info, ObservationInfo)
667 assert obs_info.instrument == "LSSTCam"
668 assert obs_info.detector_num == 85, obs_info
669 assert obs_info.detector_unique_name == "R21_S11", obs_info
670 assert obs_info.physical_filter == "r_57", obs_info
673def test_legacy_errors(legacy_test_data: _LegacyTestData) -> None:
674 """Verify that from_legacy and read_legacy raise ValueError on
675 conflicting arguments.
676 """
677 with pytest.raises(ValueError):
678 VisitImage.from_legacy(legacy_test_data.legacy_exposure, instrument="HSC")
679 with pytest.raises(ValueError):
680 VisitImage.from_legacy(legacy_test_data.legacy_exposure, visit=123456)
681 with pytest.raises(ValueError):
682 VisitImage.from_legacy(legacy_test_data.legacy_exposure, unit=u.mJy)
683 visit = VisitImage.from_legacy(
684 legacy_test_data.legacy_exposure,
685 instrument="LSSTCam",
686 unit=legacy_test_data.unit,
687 visit=2025052000177,
688 )
689 assert visit.unit == legacy_test_data.unit
691 with pytest.raises(ValueError):
692 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, instrument="HSC")
693 with pytest.raises(ValueError):
694 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, visit=123456)
697def test_component_reads(legacy_test_data: _LegacyTestData) -> None:
698 """Verify that individual components can be read from a legacy FITS
699 file.
700 """
701 visit = VisitImage.read_legacy(legacy_test_data.filename)
702 proj = VisitImage.read_legacy(legacy_test_data.filename, component="sky_projection")
703 assert_sky_projections_equal(proj, visit.sky_projection, expect_identity=False)
704 image = VisitImage.read_legacy(legacy_test_data.filename, component="image")
705 assert image == visit.image
706 assert_sky_projections_equal(proj, image.sky_projection, expect_identity=False)
707 variance = VisitImage.read_legacy(legacy_test_data.filename, component="variance")
708 assert variance == visit.variance
709 assert_sky_projections_equal(proj, variance.sky_projection, expect_identity=False)
710 mask = VisitImage.read_legacy(legacy_test_data.filename, component="mask")
711 assert mask == visit.mask
712 assert_sky_projections_equal(proj, mask.sky_projection, expect_identity=False)
713 psf = VisitImage.read_legacy(legacy_test_data.filename, component="psf")
714 assert isinstance(psf, PointSpreadFunction)
715 obs_info = VisitImage.read_legacy(legacy_test_data.filename, component="obs_info")
716 _check_legacy_obs_info(obs_info)
717 summary_stats = VisitImage.read_legacy(legacy_test_data.filename, component="summary_stats")
718 assert isinstance(summary_stats, ObservationSummaryStats)
719 assert summary_stats.nPsfStar == legacy_test_data.legacy_exposure.info.getSummaryStats().nPsfStar
720 compare_aperture_corrections_to_legacy(
721 VisitImage.read_legacy(legacy_test_data.filename, component="aperture_corrections"),
722 legacy_test_data.legacy_exposure.info.getApCorrMap(),
723 visit.bbox,
724 )
725 detector = VisitImage.read_legacy(legacy_test_data.filename, component="detector")
726 compare_detector_to_legacy(
727 detector, legacy_test_data.legacy_exposure.getDetector(), is_raw_assembled=True
728 )
729 photometric_scaling = VisitImage.read_legacy(legacy_test_data.filename, component="photometric_scaling")
730 compare_photo_calib_to_legacy(
731 photometric_scaling,
732 legacy_test_data.legacy_exposure.getPhotoCalib(),
733 subimage_bbox=visit.bbox,
734 )
737def test_legacy_obs_info(legacy_test_data: _LegacyTestData) -> None:
738 """Verify that ObservationInfo is constructed correctly from a legacy
739 exposure.
740 """
741 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map)
742 assert legacy.obs_info is not None
743 assert legacy.obs_info == legacy_test_data.visit_image.obs_info
744 assert legacy.obs_info is not None # for mypy
745 assert legacy.obs_info.instrument == "LSSTCam"
746 assert legacy.obs_info.detector_num == 85, legacy.obs_info
747 assert legacy.obs_info.detector_unique_name == "R21_S11", legacy.obs_info
748 assert legacy.obs_info.physical_filter == "r_57", legacy.obs_info
751def test_aperture_corrections_to_legacy(legacy_test_data: _LegacyTestData) -> None:
752 """Verify that aperture corrections round-trip through a legacy
753 ApCorrMap.
754 """
755 ap_corrections = legacy_test_data.visit_image.aperture_corrections
756 legacy_ap_corr_map = aperture_corrections_to_legacy(ap_corrections)
757 compare_aperture_corrections_to_legacy(
758 ap_corrections,
759 legacy_ap_corr_map,
760 legacy_test_data.visit_image.bbox,
761 )
764def _check_legacy_headers(visit_image: VisitImage) -> None:
765 """Assert that primary and extension headers are stripped correctly."""
766 header = visit_image._opaque_metadata.headers[ExtensionKey()]
767 assert "EXPTIME" in header
768 assert header["PLATFORM"] == "lsstcam"
769 assert "LSST BUTLER ID" not in header
770 assert "AR HDU" not in header
771 assert "A_ORDER" not in header
772 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("IMAGE"), astropy.io.fits.Header())
773 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("MASK"), astropy.io.fits.Header())
774 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("VARIANCE"), astropy.io.fits.Header())
777def test_read_legacy_headers(legacy_test_data: _LegacyTestData) -> None:
778 """Verify that headers were stripped and interpreted correctly in
779 read_legacy.
780 """
781 assert legacy_test_data.visit_image.unit == legacy_test_data.unit
782 _check_legacy_headers(legacy_test_data.visit_image)
785def test_from_legacy_headers(legacy_test_data: _LegacyTestData) -> None:
786 """Verify that from_legacy handles primary and extension headers
787 correctly.
788 """
789 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map)
790 assert legacy.unit == legacy_test_data.unit
791 _check_legacy_headers(legacy)
794def test_rewrite(legacy_test_data: _LegacyTestData) -> None:
795 """Verify that a legacy VisitImage can be rewritten and round-trips both
796 pixel values and all components.
797 """
798 with RoundtripFits(legacy_test_data.visit_image, legacy_test_data.storage_class) as roundtrip:
799 fits = roundtrip.inspect()
800 assert fits[1].header["ZCMPTYPE"] == "RICE_1"
801 assert fits[1].header["CTYPE1"] == "RA---TAN-SIP"
802 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
803 assert fits[2].header["CTYPE1"] == "RA---TAN-SIP"
804 assert fits[3].header["ZCMPTYPE"] == "RICE_1"
805 assert fits[3].header["CTYPE1"] == "RA---TAN-SIP"
806 subbox = Box.factory[8:13, 9:30]
807 subimage = roundtrip.get(bbox=subbox)
808 assert_masked_images_equal(subimage, legacy_test_data.visit_image[subbox], expect_view=False)
809 alternates: dict[str, Any] = {}
810 assert roundtrip.get("bbox") == legacy_test_data.visit_image.bbox
811 alternates = {
812 k: roundtrip.get(k)
813 for k in [
814 "sky_projection",
815 "image",
816 "mask",
817 "variance",
818 "psf",
819 "obs_info",
820 "summary_stats",
821 "aperture_corrections",
822 "detector",
823 "photometric_scaling",
824 ]
825 }
826 legacy_exposure = roundtrip.get(storageClass="Exposure")
827 assert isinstance(legacy_exposure, LegacyExposure)
828 compare_visit_image_to_legacy(
829 legacy_test_data.visit_image,
830 legacy_exposure,
831 expect_view=False,
832 plane_map=legacy_test_data.plane_map,
833 **DP2_VISIT_DETECTOR_DATA_ID,
834 )
835 if legacy_test_data.visit_image.unit == u.nJy:
836 assert legacy_exposure.getPhotoCalib()._isConstant
837 assert legacy_exposure.getPhotoCalib().getCalibrationMean() == 1.0
838 else:
839 compare_photo_calib_to_legacy(
840 legacy_test_data.visit_image.photometric_scaling,
841 legacy_exposure.getPhotoCalib(),
842 subimage_bbox=subbox,
843 )
844 assert legacy_exposure.info.getId() == legacy_test_data.legacy_exposure.info.getId()
845 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
846 assert isinstance(visit_info, LegacyVisitInfo)
847 assert visit_info.getInstrumentLabel() == "LSSTCam"
849 assert_visit_images_equal(roundtrip.result, legacy_test_data.visit_image, expect_view=False)
850 assert dict(legacy_test_data.visit_image._opaque_metadata.headers[ExtensionKey()]) == dict(
851 roundtrip.result._opaque_metadata.headers[ExtensionKey()]
852 )
853 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")]
854 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")]
855 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")]
856 assert roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"] == "lsstcam"
857 compare_visit_image_to_legacy(
858 roundtrip.result,
859 legacy_test_data.legacy_exposure,
860 expect_view=False,
861 plane_map=legacy_test_data.plane_map,
862 **DP2_VISIT_DETECTOR_DATA_ID,
863 alternates=alternates,
864 )
865 compare_visit_image_to_legacy(
866 legacy_test_data.read_cls.from_legacy(
867 legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map
868 ),
869 legacy_test_data.legacy_exposure,
870 expect_view=True,
871 plane_map=legacy_test_data.plane_map,
872 **DP2_VISIT_DETECTOR_DATA_ID,
873 )
876def test_butler_converters(legacy_test_data: _LegacyTestData) -> None:
877 """Verify that a VisitImage can be read from a Butler dataset written as
878 an Exposure.
879 """
880 try:
881 from lsst.daf.butler import FileDataset
882 except ImportError:
883 pytest.skip("lsst.daf.butler could not be imported.")
885 with TemporaryButler(legacy="ExposureF") as helper:
886 helper.butler.ingest(
887 FileDataset(path=legacy_test_data.filename, refs=[helper.legacy]), transfer="symlink"
888 )
889 visit_image_ref = helper.legacy.overrideStorageClass(legacy_test_data.storage_class)
890 with warnings.catch_warnings():
891 warnings.filterwarnings("ignore", message=".*filter label mismatch.*", category=UserWarning)
892 visit_image = helper.butler.get(visit_image_ref)
893 assert visit_image._opaque_metadata.precompressed.keys() == set()
894 visit_image = helper.butler.get(visit_image_ref, parameters={"preserve_quantization": True})
895 assert visit_image._opaque_metadata.precompressed.keys() == {"IMAGE", "VARIANCE"}
896 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
897 assert bbox == visit_image.bbox
898 alternates = {
899 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
900 for k in ["image", "mask", "variance", "bbox", "psf", "detector"]
901 }
902 compare_visit_image_to_legacy(
903 visit_image,
904 legacy_test_data.legacy_exposure,
905 expect_view=False,
906 plane_map=legacy_test_data.plane_map,
907 alternates=alternates,
908 **DP2_VISIT_DETECTOR_DATA_ID,
909 )
910 helper.butler.pruneDatasets([helper.legacy], purge=True, unstore=True, disassociate=True)
911 visit_image.metadata["MixedCaseKey"] = 52
912 helper.butler.put(visit_image, visit_image_ref)
913 with warnings.catch_warnings():
914 warnings.filterwarnings("ignore", message=".*filter label mismatch.*", category=UserWarning)
915 legacy_exposure = helper.butler.get(helper.legacy)
916 compare_visit_image_to_legacy(
917 visit_image,
918 legacy_exposure,
919 expect_view=False,
920 plane_map=legacy_test_data.plane_map,
921 alternates=alternates,
922 **DP2_VISIT_DETECTOR_DATA_ID,
923 )
924 visit_image_2 = helper.butler.get(visit_image_ref)
925 compare_visit_image_to_legacy(
926 visit_image_2,
927 legacy_exposure,
928 expect_view=False,
929 plane_map=legacy_test_data.plane_map,
930 alternates=alternates,
931 **DP2_VISIT_DETECTOR_DATA_ID,
932 )
933 assert visit_image_2.metadata["MixedCaseKey"] == 52
936def test_convert_unit(legacy_test_data_calibrated: _LegacyTestData) -> None:
937 """Verify convert_unit round-trips between nJy, mJy, and electron via
938 photometric_scaling.
939 """
940 from lsst.afw.table import ExposureCatalog
942 legacy_test_data = legacy_test_data_calibrated
943 original = legacy_test_data.visit_image.copy()
944 with pytest.raises(u.UnitConversionError):
945 original.convert_unit(u.electron)
946 visit_image_nJy = original.convert_unit(u.nJy, copy=False)
947 assert np.may_share_memory(visit_image_nJy.image.array, original.image.array)
948 assert np.may_share_memory(visit_image_nJy.variance.array, original.variance.array)
949 with pytest.raises(u.UnitConversionError):
950 original.convert_unit(u.mJy, copy=False)
951 visit_image_mJy = original.convert_unit(u.mJy, copy="as-needed")
952 assert visit_image_mJy.unit == u.mJy
953 assert_close(visit_image_mJy.image.array, original.image.array * 1e-6)
954 assert np.may_share_memory(visit_image_nJy.mask.array, original.mask.array)
955 assert_close(visit_image_mJy.variance.array, original.variance.array * 1e-12)
956 legacy_exposure_mJy = visit_image_mJy.to_legacy()
957 assert_close(legacy_exposure_mJy.getPhotoCalib().getCalibrationMean(), 1e6)
958 legacy_masked_image_nJy = legacy_exposure_mJy.getPhotoCalib().calibrateImage(
959 legacy_exposure_mJy.maskedImage
960 )
961 assert_close(visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
962 assert_close(visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
963 assert np.may_share_memory(visit_image_mJy.mask.array, original.mask.array)
964 assert visit_image_mJy.sky_projection is original.sky_projection
965 assert visit_image_mJy.obs_info is original.obs_info
966 assert visit_image_mJy.summary_stats is original.summary_stats
967 assert visit_image_mJy.psf is original.psf
968 assert visit_image_mJy.detector is original.detector
969 assert visit_image_mJy.bounds is original.bounds
970 assert visit_image_mJy.aperture_corrections is original.aperture_corrections
971 assert visit_image_mJy.photometric_scaling is original.photometric_scaling
972 visit_summary = ExposureCatalog.readFits(
973 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_summary.fits")
974 )
975 legacy_photo_calib = visit_summary.find(DP2_VISIT_DETECTOR_DATA_ID["detector"]).getPhotoCalib()
976 visit_image_nJy.photometric_scaling = field_from_legacy_photo_calib(
977 legacy_photo_calib, bounds=original.detector.bbox, instrumental_unit=u.electron
978 )
979 compare_photo_calib_to_legacy(
980 visit_image_nJy.photometric_scaling,
981 legacy_test_data.legacy_exposure.getPhotoCalib(),
982 applied_legacy_photo_calib=legacy_photo_calib,
983 subimage_bbox=visit_image_nJy.bbox,
984 )
985 with pytest.raises(u.UnitConversionError):
986 visit_image_nJy.convert_unit(u.mm)
987 with pytest.raises(u.UnitConversionError):
988 visit_image_nJy.convert_unit(u.electron, copy=False)
989 legacy_masked_image_e = legacy_photo_calib.uncalibrateImage(legacy_test_data.legacy_exposure.maskedImage)
990 visit_image_e = visit_image_nJy.convert_unit(u.electron)
991 assert_close(visit_image_e.image.array, legacy_masked_image_e.image.array)
992 assert_close(visit_image_e.variance.array, legacy_masked_image_e.variance.array)
993 assert not np.may_share_memory(visit_image_e.mask.array, visit_image_nJy.mask.array)
994 visit_image_mJy.photometric_scaling = visit_image_nJy.photometric_scaling
995 visit_image_e = visit_image_mJy.convert_unit(u.electron)
996 assert_close(visit_image_e.image.array, legacy_masked_image_e.image.array)
997 assert_close(visit_image_e.variance.array, legacy_masked_image_e.variance.array)
998 visit_image_nJy_2 = visit_image_e.convert_unit(u.nJy)
999 assert_close(visit_image_nJy_2.image.array, visit_image_nJy.image.array)
1000 assert_close(visit_image_nJy_2.variance.array, original.variance.array)
1001 visit_image_e.photometric_scaling = visit_image_nJy.photometric_scaling * (1e-6 * u.mJy / u.nJy)
1002 visit_image_nJy_3 = visit_image_e.convert_unit(u.nJy)
1003 assert_close(visit_image_nJy_3.image.array, visit_image_nJy.image.array)
1004 assert_close(visit_image_nJy_3.variance.array, original.variance.array)
1005 legacy_exposure_e = visit_image_e.to_legacy()
1006 assert_close(
1007 legacy_exposure_e.getPhotoCalib().getCalibrationMean(),
1008 legacy_photo_calib.getCalibrationMean(),
1009 )
1010 legacy_masked_image_nJy = legacy_exposure_e.getPhotoCalib().calibrateImage(legacy_exposure_e.maskedImage)
1011 assert_close(visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
1012 assert_close(visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)