Coverage for tests/test_from_hdu_list.py: 99%
247 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-10 09:13 +0000
1# This file is part of lsst-images.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12"""Tests for reconstructing Image/MaskedImage from cut-down HDU lists such as
13those written by ``dax_images_cutout``.
14"""
16from __future__ import annotations
18from pathlib import Path
20import astropy.io.fits
21import astropy.units as u
22import numpy as np
23import pytest
25from lsst.images import Box, GeneralFrame, Image, Mask, MaskedImage, MaskPlane, MaskSchema
26from lsst.images import fits as images_fits
27from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata
28from lsst.images.tests import (
29 assert_images_equal,
30 assert_masked_images_equal,
31 make_random_sky_projection,
32)
35def make_hdu_list(
36 tmp_path: Path,
37 *,
38 projection: bool = False,
39 planes: int = 3,
40 rng: np.random.Generator | None = None,
41) -> tuple[MaskedImage, astropy.io.fits.HDUList]:
42 """Return a (MaskedImage, cut-down HDUList) pair for use in tests."""
43 if rng is None:
44 rng = np.random.default_rng(7)
45 shape = (20, 25)
46 yx0 = (5, 8)
47 bbox = Box.from_shape(shape, start=yx0)
48 proj = make_random_sky_projection(rng, GeneralFrame(unit=u.pix), bbox) if projection else None
49 image = Image(rng.normal(100.0, 8.0, shape).astype("float32"), unit=u.nJy, yx0=yx0, sky_projection=proj)
50 schema = MaskSchema([MaskPlane(f"P{i}", f"description {i}") for i in range(planes)])
51 masked_image = MaskedImage(image, mask_schema=schema)
52 for i in range(planes):
53 masked_image.mask.set(f"P{i}", rng.random(shape) > 0.5)
54 masked_image.variance.array = rng.normal(64.0, 0.5, shape)
56 path = str(tmp_path / "x.fits")
57 masked_image.write(path)
58 with astropy.io.fits.open(path) as hdul:
59 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [
60 astropy.io.fits.PrimaryHDU(header=hdul[0].header.copy())
61 ]
62 for name in ["IMAGE", "MASK", "VARIANCE"]:
63 src = hdul[name]
64 hdus.append(
65 astropy.io.fits.ImageHDU(data=np.asarray(src.data), header=src.header.copy(), name=name)
66 )
67 hdu_list = astropy.io.fits.HDUList(hdus)
68 return masked_image, hdu_list
71def _cutdown(
72 obj: object, names: list[str], tmp_path: Path, **write_kwargs: object
73) -> astropy.io.fits.HDUList:
74 """Return an in-memory cut-down HDUList with the requested extension
75 names.
76 """
77 path = str(tmp_path / "cutdown.fits")
78 obj.write(path, **write_kwargs) # type: ignore[attr-defined]
79 with astropy.io.fits.open(path) as hdul:
80 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [
81 astropy.io.fits.PrimaryHDU(header=hdul[0].header.copy())
82 ]
83 for name in names:
84 src = hdul[name]
85 hdus.append(
86 astropy.io.fits.ImageHDU(data=np.asarray(src.data), header=src.header.copy(), name=name)
87 )
88 return astropy.io.fits.HDUList(hdus)
91def make_legacy_masked_image_hdu_list(
92 planes: dict[str, int],
93 set_pixels: dict[str, tuple[int, int]],
94 *,
95 shape: tuple[int, int] = (6, 7),
96 yx0: tuple[int, int] = (5, 8),
97 rng: np.random.Generator | None = None,
98) -> astropy.io.fits.HDUList:
99 """Return an afw-style legacy cut-down HDUList with MP_ cards in the MASK
100 HDU.
101 """
102 if rng is None: 102 ↛ 104line 102 didn't jump to line 104 because the condition on line 102 was always true
103 rng = np.random.default_rng(7)
104 mask_data = np.zeros(shape, dtype=np.int32)
105 for name, (y, x) in set_pixels.items():
106 mask_data[y, x] |= 1 << planes[name]
107 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [astropy.io.fits.PrimaryHDU()]
108 for name, data in [
109 ("IMAGE", rng.normal(0.0, 1.0, shape).astype("float32")),
110 ("MASK", mask_data),
111 ("VARIANCE", rng.normal(1.0, 0.1, shape).astype("float32")),
112 ]:
113 hdu = astropy.io.fits.ImageHDU(data=data, name=name)
114 hdu.header["LTV1"] = -yx0[1]
115 hdu.header["LTV2"] = -yx0[0]
116 hdus.append(hdu)
117 with images_fits.suppress_fits_card_warnings():
118 for name, bit in planes.items():
119 hdus[2].header[f"MP_{name}"] = bit
120 return astropy.io.fits.HDUList(hdus)
123def _schema_index(mask: Mask, name: str) -> int:
124 """Return the zero-based plane index of the named plane in the mask
125 schema.
126 """
127 return next(n for n, plane in enumerate(mask.schema) if plane is not None and plane.name == name)
130AFW_VISIT_BITS = {
131 "BAD": 0,
132 "SAT": 1,
133 "INTRP": 2,
134 "CR": 3,
135 "EDGE": 4,
136 "DETECTED": 5,
137 "DETECTED_NEGATIVE": 6,
138 "SUSPECT": 7,
139 "NO_DATA": 8,
140}
143def _legacy_mask_hdu(
144 set_pixels: dict[str, tuple[int, int]], *, shape: tuple[int, int] = (6, 7)
145) -> astropy.io.fits.ImageHDU:
146 """Return a standalone afw-style MASK ImageHDU with MP_ cards."""
147 data = np.zeros(shape, dtype=np.int32)
148 for name, (y, x) in set_pixels.items():
149 data[y, x] |= 1 << AFW_VISIT_BITS[name]
150 hdu = astropy.io.fits.ImageHDU(data=data, name="MASK")
151 hdu.header["LTV1"] = -8
152 hdu.header["LTV2"] = -5
153 with images_fits.suppress_fits_card_warnings():
154 for name, bit in AFW_VISIT_BITS.items():
155 hdu.header[f"MP_{name}"] = bit
156 return hdu
159def _legacy_full_hdu_list(
160 set_pixels: dict[str, tuple[int, int]], *, shape: tuple[int, int] = (6, 7)
161) -> astropy.io.fits.HDUList:
162 """Return an afw-style full HDUList (PRIMARY+IMAGE+MASK+VARIANCE) with MP_
163 cards.
164 """
165 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [astropy.io.fits.PrimaryHDU()]
166 for name, data in [
167 ("IMAGE", np.zeros(shape, dtype=np.float32)),
168 ("MASK", None),
169 ("VARIANCE", np.ones(shape, dtype=np.float32)),
170 ]:
171 if name == "MASK":
172 hdus.append(_legacy_mask_hdu(set_pixels, shape=shape))
173 else:
174 hdu = astropy.io.fits.ImageHDU(data=data, name=name)
175 hdu.header["LTV1"] = -8
176 hdu.header["LTV2"] = -5
177 hdus.append(hdu)
178 return astropy.io.fits.HDUList(hdus)
181def test_offset_wcs_round_trip() -> None:
182 """Verify add_offset_wcs and read_offset_wcs are inverses of each other."""
183 header = astropy.io.fits.Header()
184 images_fits.add_offset_wcs(header, x=19190, y=22580, key="A")
185 assert images_fits.read_offset_wcs(header, key="A") == (19190, 22580)
188def test_offset_wcs_absent_returns_none() -> None:
189 """Verify read_offset_wcs returns None when no offset WCS key is
190 present.
191 """
192 assert images_fits.read_offset_wcs(astropy.io.fits.Header(), key="A") is None
195def test_offset_wcs_other_key_ignored() -> None:
196 """Verify read_offset_wcs ignores WCS keys written under a different
197 letter.
198 """
199 header = astropy.io.fits.Header()
200 images_fits.add_offset_wcs(header, x=3, y=4, key="A")
201 assert images_fits.read_offset_wcs(header, key="B") is None
204def test_read_yx0_from_offset_wcs() -> None:
205 """Verify read_yx0 recovers yx0 from an offset WCS header."""
206 header = astropy.io.fits.Header()
207 images_fits.add_offset_wcs(header, x=19190, y=22580)
208 yx0 = images_fits.read_yx0(header)
209 assert (yx0.y, yx0.x) == (22580, 19190)
212def test_read_yx0_from_ltv() -> None:
213 """Verify read_yx0 recovers yx0 from LTV1/LTV2 header cards."""
214 header = astropy.io.fits.Header()
215 header["LTV1"] = -19190
216 header["LTV2"] = -22580
217 yx0 = images_fits.read_yx0(header)
218 assert (yx0.y, yx0.x) == (22580, 19190)
221def test_read_yx0_missing_raises() -> None:
222 """Verify read_yx0 raises ValueError when no offset information is
223 present.
224 """
225 with pytest.raises(ValueError, match="records no LTV1/LTV2 cards"):
226 images_fits.read_yx0(astropy.io.fits.Header())
229def test_masked_image_round_trip(tmp_path: Path) -> None:
230 """Verify a cut-down MaskedImage reconstructs to an equal MaskedImage."""
231 masked_image, cutdown = make_hdu_list(tmp_path)
232 result = MaskedImage.from_hdu_list(cutdown)
233 assert_masked_images_equal(result, masked_image)
236def test_legacy_non_cell_coadd_from_hdu_list() -> None:
237 """Verify a legacy afw-style cut-down HDUList reconstructs a MaskedImage
238 with MP_ planes.
239 """
240 planes = {"BAD": 0, "DETECTED": 5, "INEXACT_PSF": 11, "SENSOR_EDGE": 14}
241 set_pixels = {"DETECTED": (1, 2), "INEXACT_PSF": (3, 4), "SENSOR_EDGE": (5, 6)}
242 hdul = make_legacy_masked_image_hdu_list(planes, set_pixels)
243 result = MaskedImage.from_hdu_list(hdul)
244 assert "SENSOR_EDGE" in result.mask.schema.names
245 assert result.mask.get("SENSOR_EDGE")[5, 6]
246 assert result.mask.get("INEXACT_PSF")[3, 4]
247 assert result.mask.bbox.y.start == 5
248 assert result.mask.bbox.x.start == 8
251def test_masked_image_round_trip_with_projection(tmp_path: Path) -> None:
252 """Verify the sky projection is recovered from the FITS WCS in the cut-down
253 HDUList.
254 """
255 masked_image, cutdown = make_hdu_list(tmp_path, projection=True)
256 result = MaskedImage.from_hdu_list(cutdown)
257 assert result.sky_projection is not None
258 center = (masked_image.bbox.x.size / 2, masked_image.bbox.y.size / 2)
259 expected_wcs = masked_image.fits_wcs
260 actual_wcs = result.fits_wcs
261 assert expected_wcs is not None
262 assert actual_wcs is not None
263 expected = expected_wcs.pixel_to_world(*center)
264 actual = actual_wcs.pixel_to_world(*center)
265 assert expected.separation(actual).arcsec < 1e-3
268def test_mask_planes_repacked_across_byte_boundary(tmp_path: Path) -> None:
269 """Verify nine planes stored in an int32 HDU are repacked into the two-byte
270 uint8 layout.
271 """
272 rng = np.random.default_rng(7)
273 masked_image, _ = make_hdu_list(tmp_path, planes=9, rng=rng)
274 assert masked_image.mask.schema.mask_size == 2
275 cutdown = _cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"], tmp_path)
276 assert np.asarray(cutdown["MASK"].data).dtype.kind == "i"
277 result = MaskedImage.from_hdu_list(cutdown)
278 assert result.mask.schema.mask_size == 2
279 assert_masked_images_equal(result, masked_image)
282def test_image_from_hdu_list_reads_first_two_hdus(tmp_path: Path) -> None:
283 """Verify Image.from_hdu_list reads PRIMARY+IMAGE and ignores later
284 HDUs.
285 """
286 masked_image, cutdown = make_hdu_list(tmp_path)
287 result = Image.from_hdu_list(cutdown)
288 assert_images_equal(result, masked_image.image)
289 two = _cutdown(masked_image, ["IMAGE"], tmp_path)
290 assert_images_equal(Image.from_hdu_list(two), masked_image.image)
293def test_missing_mask_schema_raises(tmp_path: Path) -> None:
294 """Verify a MASK HDU without MSK* cards raises ValueError."""
295 _, cutdown = make_hdu_list(tmp_path)
296 for key in [k for k in cutdown["MASK"].header if k.startswith(("MSKN", "MSKM", "MSKD"))]:
297 del cutdown["MASK"].header[key]
298 with pytest.raises(ValueError, match="cannot reconstruct the mask schema"):
299 MaskedImage.from_hdu_list(cutdown)
302def test_multiple_mask_hdus_raises(tmp_path: Path) -> None:
303 """Verify two MASK HDUs raise ValueError rather than silently dropping
304 one.
305 """
306 _, cutdown = make_hdu_list(tmp_path)
307 extra_mask = astropy.io.fits.ImageHDU(
308 data=np.asarray(cutdown["MASK"].data), header=cutdown["MASK"].header.copy(), name="MASK"
309 )
310 extra_mask.header["EXTVER"] = 2
311 cutdown.append(extra_mask)
312 with pytest.raises(ValueError, match="supports only a single MASK extension"):
313 MaskedImage.from_hdu_list(cutdown)
316def test_primary_header_preserved(tmp_path: Path) -> None:
317 """Verify confusing container cards are dropped and other primary cards
318 survive as opaque metadata.
319 """
320 masked_image, _ = make_hdu_list(tmp_path)
322 def add_card(header: astropy.io.fits.Header) -> None:
323 header["MYCARD"] = "hello"
325 cutdown = _cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"], tmp_path, update_header=add_card)
326 result = MaskedImage.from_hdu_list(cutdown)
327 assert isinstance(result._opaque_metadata, FitsOpaqueMetadata)
328 primary = result._opaque_metadata.headers[ExtensionKey()]
329 assert primary["MYCARD"] == "hello"
330 assert "DATAMODL" not in primary
331 assert "INDXADDR" not in primary
332 assert "JSONADDR" not in primary
335def test_legacy_mp_ltv_path() -> None:
336 """Verify the legacy MP_/LTV branch reads a synthetic afw-style MASK HDU
337 correctly.
338 """
339 data = np.zeros((6, 7), dtype=np.int32)
340 data[1, 2] = 0b01 # BAD
341 data[3, 4] = 0b10 # SAT
342 hdu = astropy.io.fits.ImageHDU(data=data, name="MASK")
343 hdu.header["LTV1"] = -8
344 hdu.header["LTV2"] = -5
345 hdu.header["MP_BAD"] = 0
346 hdu.header["MP_SAT"] = 1
347 plane_map = {"BAD": MaskPlane("BAD", "bad"), "SAT": MaskPlane("SATURATED", "saturated")}
348 mask = Mask._read_legacy_hdu(hdu, FitsOpaqueMetadata(), plane_map=plane_map)
349 assert mask.bbox.y.start == 5
350 assert mask.bbox.x.start == 8
351 assert set(mask.schema.names) == {"BAD", "SATURATED"}
352 assert mask.get("BAD")[1, 2]
353 assert mask.get("SATURATED")[3, 4]
354 assert not mask.get("BAD")[3, 4]
357def test_read_legacy_strip_false_keeps_cards() -> None:
358 """Verify MaskPlane.read_legacy with strip=False reads planes but leaves
359 MP_ cards in the header.
360 """
361 header = astropy.io.fits.Header()
362 header["MP_BAD"] = 0
363 header["MP_SAT"] = 1
364 planes = MaskPlane.read_legacy(header, strip=False)
365 assert planes == {"BAD": 0, "SAT": 1}
366 assert "MP_BAD" in header
367 assert "MP_SAT" in header
370def test_read_legacy_default_strips_cards() -> None:
371 """Verify MaskPlane.read_legacy default strips MP_ cards from the
372 header.
373 """
374 header = astropy.io.fits.Header()
375 header["MP_BAD"] = 0
376 header["MP_SAT"] = 1
377 planes = MaskPlane.read_legacy(header)
378 assert planes == {"BAD": 0, "SAT": 1}
379 assert "MP_BAD" not in header
380 assert "MP_SAT" not in header
383def test_read_legacy_hdu_reindexes_retained_cards() -> None:
384 """Verify _read_legacy_hdu with strip_legacy_planes=False rewrites MP_
385 cards to new bit positions.
386 """
387 hdu = _legacy_mask_hdu({"SUSPECT": (1, 2), "NO_DATA": (3, 4), "BAD": (0, 0)})
388 mask = Mask._read_legacy_hdu(hdu, FitsOpaqueMetadata(), strip_legacy_planes=False)
389 assert mask.get("SUSPECT")[1, 2]
390 assert mask.get("NO_DATA")[3, 4]
391 assert _schema_index(mask, "SUSPECT") == 6
392 assert _schema_index(mask, "NO_DATA") == 7
393 assert hdu.header["MP_SUSPECT"] == 6
394 assert hdu.header["MP_NO_DATA"] == 7
395 assert hdu.header["MP_BAD"] == 0
396 assert hdu.header["MP_DETECTED"] == 5
397 assert "MP_DETECTED_NEGATIVE" not in hdu.header
400def test_read_legacy_hdu_default_strips() -> None:
401 """Verify _read_legacy_hdu default behavior strips all MP_ cards."""
402 hdu = _legacy_mask_hdu({"BAD": (0, 0)})
403 Mask._read_legacy_hdu(hdu, FitsOpaqueMetadata())
404 assert not [k for k in hdu.header if k.startswith("MP_")]
407def test_from_hdu_list_round_trips_reindexed_mp_cards(tmp_path: Path) -> None:
408 """Verify a reconstructed legacy MaskedImage re-serializes with correct MP_
409 bit indices.
410 """
411 hdul = _legacy_full_hdu_list({"SUSPECT": (1, 2), "NO_DATA": (3, 4), "BAD": (0, 0)})
412 masked_image = MaskedImage.from_hdu_list(hdul)
413 path = str(tmp_path / "out.fits")
414 masked_image.write(path)
415 with astropy.io.fits.open(path) as out:
416 header = out["MASK"].header
417 array = np.asarray(out["MASK"].data)
418 mskn = {header[k]: int(k.removeprefix("MSKN")) for k in header if k.startswith("MSKN")}
419 assert header["MP_SUSPECT"] == mskn["SUSPECT"]
420 assert header["MP_NO_DATA"] == mskn["NO_DATA"]
421 assert header["MP_SUSPECT"] == 6
422 assert "MP_DETECTED_NEGATIVE" not in header
423 assert array[1, 2] & (1 << header["MP_SUSPECT"])
424 assert array[3, 4] & (1 << header["MP_NO_DATA"])
425 assert array[0, 0] & (1 << header["MP_BAD"])
428def test_normal_read_strips_mp_cards(tmp_path: Path) -> None:
429 """Verify a normal read + rewrite of a legacy-cutout file drops the MP_
430 cards.
431 """
432 hdul = _legacy_full_hdu_list({"SUSPECT": (1, 2), "NO_DATA": (3, 4), "BAD": (0, 0)})
433 masked_image = MaskedImage.from_hdu_list(hdul)
434 legacy_cutout = str(tmp_path / "legacy_cutout.fits")
435 masked_image.write(legacy_cutout)
436 with astropy.io.fits.open(legacy_cutout) as out:
437 assert [k for k in out["MASK"].header if k.startswith("MP_")]
438 rewritten = str(tmp_path / "rewritten.fits")
439 MaskedImage.read(legacy_cutout).write(rewritten)
440 with astropy.io.fits.open(rewritten) as out:
441 assert not [k for k in out["MASK"].header if k.startswith("MP_")]