Coverage for tests/test_from_hdu_list.py: 100%
249 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-29 15:33 -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.
12"""Tests for reconstructing Image/MaskedImage from cut-down HDU lists such as
13those written by ``dax_images_cutout``.
14"""
16from __future__ import annotations
18import os
19import tempfile
20import unittest
22import astropy.io.fits
23import astropy.units as u
24import numpy as np
26from lsst.images import Box, GeneralFrame, Image, Mask, MaskedImage, MaskPlane, MaskSchema
27from lsst.images import fits as images_fits
28from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata
29from lsst.images.tests import (
30 assert_images_equal,
31 assert_masked_images_equal,
32 make_random_sky_projection,
33)
36class ReadOffsetWcsTestCase(unittest.TestCase):
37 """Tests for the inverse of `lsst.images.fits.add_offset_wcs`."""
39 def test_round_trip(self) -> None:
40 header = astropy.io.fits.Header()
41 images_fits.add_offset_wcs(header, x=19190, y=22580, key="A")
42 self.assertEqual(images_fits.read_offset_wcs(header, key="A"), (19190, 22580))
44 def test_absent_returns_none(self) -> None:
45 self.assertIsNone(images_fits.read_offset_wcs(astropy.io.fits.Header(), key="A"))
47 def test_other_key_ignored(self) -> None:
48 header = astropy.io.fits.Header()
49 images_fits.add_offset_wcs(header, x=3, y=4, key="A")
50 self.assertIsNone(images_fits.read_offset_wcs(header, key="B"))
53class ReadYx0TestCase(unittest.TestCase):
54 """Tests for recovering yx0 from either offset convention."""
56 def test_from_offset_wcs(self) -> None:
57 header = astropy.io.fits.Header()
58 images_fits.add_offset_wcs(header, x=19190, y=22580)
59 yx0 = images_fits.read_yx0(header)
60 self.assertEqual((yx0.y, yx0.x), (22580, 19190))
62 def test_from_ltv(self) -> None:
63 header = astropy.io.fits.Header()
64 header["LTV1"] = -19190
65 header["LTV2"] = -22580
66 yx0 = images_fits.read_yx0(header)
67 self.assertEqual((yx0.y, yx0.x), (22580, 19190))
69 def test_missing_raises(self) -> None:
70 with self.assertRaises(ValueError):
71 images_fits.read_yx0(astropy.io.fits.Header())
74class FromHduListTestCase(unittest.TestCase):
75 """Tests for reconstructing Image/MaskedImage from cut-down HDU lists."""
77 def setUp(self) -> None:
78 self.maxDiff = None
79 self.rng = np.random.default_rng(7)
81 def _build_masked_image(self, *, projection: bool = False, planes: int = 3) -> MaskedImage:
82 shape = (20, 25)
83 yx0 = (5, 8)
84 bbox = Box.from_shape(shape, start=yx0)
85 proj = make_random_sky_projection(self.rng, GeneralFrame(unit=u.pix), bbox) if projection else None
86 image = Image(
87 self.rng.normal(100.0, 8.0, shape).astype("float32"), unit=u.nJy, yx0=yx0, sky_projection=proj
88 )
89 schema = MaskSchema([MaskPlane(f"P{i}", f"description {i}") for i in range(planes)])
90 masked_image = MaskedImage(image, mask_schema=schema)
91 for i in range(planes):
92 masked_image.mask.set(f"P{i}", self.rng.random(shape) > 0.5)
93 masked_image.variance.array = self.rng.normal(64.0, 0.5, shape)
94 return masked_image
96 def _cutdown(self, obj: object, names: list[str], **write_kwargs: object) -> astropy.io.fits.HDUList:
97 """Serialize ``obj`` and return an in-memory cut-down HDU list: the
98 primary HDU plus uncompressed ImageHDUs for ``names`` (JSON and INDEX
99 dropped), mimicking a ``dax_images_cutout`` file.
100 """
101 with tempfile.TemporaryDirectory() as tmp:
102 path = os.path.join(tmp, "x.fits")
103 obj.write(path, **write_kwargs) # type: ignore[attr-defined]
104 with astropy.io.fits.open(path) as hdul:
105 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [
106 astropy.io.fits.PrimaryHDU(header=hdul[0].header.copy())
107 ]
108 for name in names:
109 src = hdul[name]
110 hdus.append(
111 astropy.io.fits.ImageHDU(
112 data=np.asarray(src.data), header=src.header.copy(), name=name
113 )
114 )
115 return astropy.io.fits.HDUList(hdus)
117 def _legacy_masked_image_hdu_list(
118 self,
119 planes: dict[str, int],
120 set_pixels: dict[str, tuple[int, int]],
121 *,
122 shape: tuple[int, int] = (6, 7),
123 yx0: tuple[int, int] = (5, 8),
124 ) -> astropy.io.fits.HDUList:
125 """Build an afw-style legacy cut-down HDU list (PRIMARY + IMAGE + MASK
126 + VARIANCE) whose MASK HDU carries ``MP_`` cards instead of ``MSKN``,
127 mimicking a ``dax_images_cutout`` file made from an afw-written image.
129 ``planes`` maps legacy plane name to bit index; ``set_pixels`` maps
130 legacy plane name to a single ``(y, x)`` pixel to set for that plane.
131 """
132 mask_data = np.zeros(shape, dtype=np.int32)
133 for name, (y, x) in set_pixels.items():
134 mask_data[y, x] |= 1 << planes[name]
135 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [astropy.io.fits.PrimaryHDU()]
136 for name, data in [
137 ("IMAGE", self.rng.normal(0.0, 1.0, shape).astype("float32")),
138 ("MASK", mask_data),
139 ("VARIANCE", self.rng.normal(1.0, 0.1, shape).astype("float32")),
140 ]:
141 hdu = astropy.io.fits.ImageHDU(data=data, name=name)
142 hdu.header["LTV1"] = -yx0[1]
143 hdu.header["LTV2"] = -yx0[0]
144 hdus.append(hdu)
145 with images_fits.suppress_fits_card_warnings():
146 for name, bit in planes.items():
147 hdus[2].header[f"MP_{name}"] = bit
148 return astropy.io.fits.HDUList(hdus)
150 def test_masked_image_round_trip(self) -> None:
151 """A cut-down MaskedImage reconstructs to an equal MaskedImage."""
152 masked_image = self._build_masked_image()
153 cutdown = self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"])
154 result = MaskedImage.from_hdu_list(cutdown)
155 assert_masked_images_equal(self, result, masked_image)
157 def test_legacy_non_cell_coadd_from_hdu_list(self) -> None:
158 """A cut-down afw-style non-cell coadd (``MP_`` cards, ``SENSOR_EDGE``
159 set) is reconstructed, mapping ``SENSOR_EDGE`` to its own plane.
160 """
161 planes = {"BAD": 0, "DETECTED": 5, "INEXACT_PSF": 11, "SENSOR_EDGE": 14}
162 set_pixels = {"DETECTED": (1, 2), "INEXACT_PSF": (3, 4), "SENSOR_EDGE": (5, 6)}
163 hdul = self._legacy_masked_image_hdu_list(planes, set_pixels)
164 result = MaskedImage.from_hdu_list(hdul)
165 self.assertIn("SENSOR_EDGE", result.mask.schema.names)
166 self.assertTrue(result.mask.get("SENSOR_EDGE")[5, 6])
167 self.assertTrue(result.mask.get("INEXACT_PSF")[3, 4])
168 self.assertEqual(result.mask.bbox.y.start, 5)
169 self.assertEqual(result.mask.bbox.x.start, 8)
171 def test_masked_image_round_trip_with_projection(self) -> None:
172 """The sky projection is recovered from the FITS WCS."""
173 masked_image = self._build_masked_image(projection=True)
174 cutdown = self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"])
175 result = MaskedImage.from_hdu_list(cutdown)
176 self.assertIsNotNone(result.sky_projection)
177 center = (masked_image.bbox.x.size / 2, masked_image.bbox.y.size / 2)
178 expected_wcs = masked_image.fits_wcs
179 actual_wcs = result.fits_wcs
180 assert expected_wcs is not None and actual_wcs is not None
181 expected = expected_wcs.pixel_to_world(*center)
182 actual = actual_wcs.pixel_to_world(*center)
183 self.assertLess(expected.separation(actual).arcsec, 1e-3)
185 def test_mask_planes_repacked_across_byte_boundary(self) -> None:
186 """Nine planes stored as one on-disk int32 HDU are repacked into the
187 two-byte uint8 in-memory layout, preserving every plane.
188 """
189 masked_image = self._build_masked_image(planes=9)
190 self.assertEqual(masked_image.mask.schema.mask_size, 2)
191 cutdown = self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"])
192 self.assertEqual(np.asarray(cutdown["MASK"].data).dtype.kind, "i")
193 result = MaskedImage.from_hdu_list(cutdown)
194 self.assertEqual(result.mask.schema.mask_size, 2)
195 assert_masked_images_equal(self, result, masked_image)
197 def test_image_from_hdu_list_reads_first_two_hdus(self) -> None:
198 """Image.from_hdu_list reads PRIMARY+IMAGE and ignores later HDUs."""
199 masked_image = self._build_masked_image()
200 # A full four-HDU list: only PRIMARY and IMAGE are consulted.
201 result = Image.from_hdu_list(self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"]))
202 assert_images_equal(self, result, masked_image.image)
203 # A bare two-HDU list works too.
204 two = self._cutdown(masked_image, ["IMAGE"])
205 assert_images_equal(self, Image.from_hdu_list(two), masked_image.image)
207 def test_missing_mask_schema_raises(self) -> None:
208 """A MASK HDU without MSK* cards is rejected (for now)."""
209 masked_image = self._build_masked_image()
210 cutdown = self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"])
211 for key in [k for k in cutdown["MASK"].header if k.startswith(("MSKN", "MSKM", "MSKD"))]:
212 del cutdown["MASK"].header[key]
213 with self.assertRaises(ValueError):
214 MaskedImage.from_hdu_list(cutdown)
216 def test_multiple_mask_hdus_raises(self) -> None:
217 """Two MASK HDUs (e.g. EXTVER 1 and 2) are rejected rather than
218 silently dropping the mask information from all but the first.
219 """
220 masked_image = self._build_masked_image()
221 cutdown = self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"])
222 extra_mask = astropy.io.fits.ImageHDU(
223 data=np.asarray(cutdown["MASK"].data), header=cutdown["MASK"].header.copy(), name="MASK"
224 )
225 extra_mask.header["EXTVER"] = 2
226 cutdown.append(extra_mask)
227 with self.assertRaises(ValueError):
228 MaskedImage.from_hdu_list(cutdown)
230 def test_primary_header_preserved(self) -> None:
231 """Confusing container cards are dropped; other primary cards survive
232 as opaque metadata.
233 """
234 masked_image = self._build_masked_image()
236 def add_card(header: astropy.io.fits.Header) -> None:
237 header["MYCARD"] = "hello"
239 cutdown = self._cutdown(masked_image, ["IMAGE", "MASK", "VARIANCE"], update_header=add_card)
240 result = MaskedImage.from_hdu_list(cutdown)
241 assert isinstance(result._opaque_metadata, FitsOpaqueMetadata)
242 primary = result._opaque_metadata.headers[ExtensionKey()]
243 self.assertEqual(primary["MYCARD"], "hello")
244 self.assertNotIn("DATAMODL", primary)
245 self.assertNotIn("INDXADDR", primary)
246 self.assertNotIn("JSONADDR", primary)
249class LegacyMaskBranchTestCase(unittest.TestCase):
250 """Characterize the legacy MP_/LTV branch of the generalized mask reader.
252 The full ``read_legacy`` paths require ``lsst.afw.image`` and external test
253 data, so this exercises the legacy MP_/LTV branch directly with a synthetic
254 afw-style HDU instead.
255 """
257 def test_legacy_mp_ltv_path(self) -> None:
258 data = np.zeros((6, 7), dtype=np.int32)
259 data[1, 2] = 0b01 # BAD
260 data[3, 4] = 0b10 # SAT
261 hdu = astropy.io.fits.ImageHDU(data=data, name="MASK")
262 hdu.header["LTV1"] = -8
263 hdu.header["LTV2"] = -5
264 hdu.header["MP_BAD"] = 0
265 hdu.header["MP_SAT"] = 1
266 plane_map = {"BAD": MaskPlane("BAD", "bad"), "SAT": MaskPlane("SATURATED", "saturated")}
267 mask = Mask._read_legacy_hdu(hdu, FitsOpaqueMetadata(), plane_map=plane_map)
268 self.assertEqual(mask.bbox.y.start, 5)
269 self.assertEqual(mask.bbox.x.start, 8)
270 self.assertEqual(set(mask.schema.names), {"BAD", "SATURATED"})
271 self.assertTrue(mask.get("BAD")[1, 2])
272 self.assertTrue(mask.get("SATURATED")[3, 4])
273 self.assertFalse(mask.get("BAD")[3, 4])
276class LegacyPlaneRetentionTestCase(unittest.TestCase):
277 """The legacy-compatibility reconstruction path keeps the ``MP_``
278 mask-plane cards, re-indexed to the reshuffled schema, instead of stripping
279 them.
280 """
282 # A sample set of ``MP_`` bit assignments used purely as a test fixture.
283 # This is not an externally stable layout: the bits afw actually writes
284 # depend on every Mask created in the process, not just an image's own
285 # planes. ``DETECTED_NEGATIVE`` is not part of the visit-image schema, so
286 # it is dropped on reconstruction and every plane after it (``SUSPECT``,
287 # ``NO_DATA``) shifts down by one.
288 AFW_VISIT_BITS = {
289 "BAD": 0,
290 "SAT": 1,
291 "INTRP": 2,
292 "CR": 3,
293 "EDGE": 4,
294 "DETECTED": 5,
295 "DETECTED_NEGATIVE": 6,
296 "SUSPECT": 7,
297 "NO_DATA": 8,
298 }
300 def test_read_legacy_strip_false_keeps_cards(self) -> None:
301 """``MaskPlane.read_legacy(strip=False)`` reads the plane bits but
302 leaves the ``MP_`` cards in the header.
303 """
304 header = astropy.io.fits.Header()
305 header["MP_BAD"] = 0
306 header["MP_SAT"] = 1
307 planes = MaskPlane.read_legacy(header, strip=False)
308 self.assertEqual(planes, {"BAD": 0, "SAT": 1})
309 self.assertIn("MP_BAD", header)
310 self.assertIn("MP_SAT", header)
312 def test_read_legacy_default_strips_cards(self) -> None:
313 """The default still strips the ``MP_`` cards (behavior for code that
314 declares a new-schema-only mask).
315 """
316 header = astropy.io.fits.Header()
317 header["MP_BAD"] = 0
318 header["MP_SAT"] = 1
319 planes = MaskPlane.read_legacy(header)
320 self.assertEqual(planes, {"BAD": 0, "SAT": 1})
321 self.assertNotIn("MP_BAD", header)
322 self.assertNotIn("MP_SAT", header)
324 def _legacy_mask_hdu(
325 self, set_pixels: dict[str, tuple[int, int]], *, shape: tuple[int, int] = (6, 7)
326 ) -> astropy.io.fits.ImageHDU:
327 """Build a standalone afw-style ``MASK`` HDU."""
328 data = np.zeros(shape, dtype=np.int32)
329 for name, (y, x) in set_pixels.items():
330 data[y, x] |= 1 << self.AFW_VISIT_BITS[name]
331 hdu = astropy.io.fits.ImageHDU(data=data, name="MASK")
332 hdu.header["LTV1"] = -8
333 hdu.header["LTV2"] = -5
334 with images_fits.suppress_fits_card_warnings():
335 for name, bit in self.AFW_VISIT_BITS.items():
336 hdu.header[f"MP_{name}"] = bit
337 return hdu
339 def _schema_index(self, mask: Mask, name: str) -> int:
340 return next(n for n, plane in enumerate(mask.schema) if plane is not None and plane.name == name)
342 def test_read_legacy_hdu_reindexes_retained_cards(self) -> None:
343 """With ``strip_legacy_planes=False`` the surviving ``MP_`` cards are
344 rewritten to the bit positions of the reshuffled schema.
345 """
346 hdu = self._legacy_mask_hdu({"SUSPECT": (1, 2), "NO_DATA": (3, 4), "BAD": (0, 0)})
347 mask = Mask._read_legacy_hdu(hdu, FitsOpaqueMetadata(), strip_legacy_planes=False)
348 # The plane pixels are repacked into the canonical visit-image schema.
349 self.assertTrue(mask.get("SUSPECT")[1, 2])
350 self.assertTrue(mask.get("NO_DATA")[3, 4])
351 # SUSPECT/NO_DATA shifted down by one relative to the input layout.
352 self.assertEqual(self._schema_index(mask, "SUSPECT"), 6)
353 self.assertEqual(self._schema_index(mask, "NO_DATA"), 7)
354 # Retained MP_ cards now record the new positions, not the input's.
355 self.assertEqual(hdu.header["MP_SUSPECT"], 6)
356 self.assertEqual(hdu.header["MP_NO_DATA"], 7)
357 # Planes whose positions did not move are unchanged.
358 self.assertEqual(hdu.header["MP_BAD"], 0)
359 self.assertEqual(hdu.header["MP_DETECTED"], 5)
360 # The plane that is not in the new schema is dropped entirely.
361 self.assertNotIn("MP_DETECTED_NEGATIVE", hdu.header)
363 def test_read_legacy_hdu_default_strips(self) -> None:
364 """The default ``_read_legacy_hdu`` behavior still strips ``MP_``."""
365 hdu = self._legacy_mask_hdu({"BAD": (0, 0)})
366 Mask._read_legacy_hdu(hdu, FitsOpaqueMetadata())
367 self.assertFalse([k for k in hdu.header if k.startswith("MP_")])
369 def _legacy_full_hdu_list(
370 self, set_pixels: dict[str, tuple[int, int]], *, shape: tuple[int, int] = (6, 7)
371 ) -> astropy.io.fits.HDUList:
372 hdus: list[astropy.io.fits.hdu.base._BaseHDU] = [astropy.io.fits.PrimaryHDU()]
373 for name, data in [
374 ("IMAGE", np.zeros(shape, dtype=np.float32)),
375 ("MASK", None),
376 ("VARIANCE", np.ones(shape, dtype=np.float32)),
377 ]:
378 if name == "MASK":
379 hdus.append(self._legacy_mask_hdu(set_pixels, shape=shape))
380 else:
381 hdu = astropy.io.fits.ImageHDU(data=data, name=name)
382 hdu.header["LTV1"] = -8
383 hdu.header["LTV2"] = -5
384 hdus.append(hdu)
385 return astropy.io.fits.HDUList(hdus)
387 def test_from_hdu_list_round_trips_reindexed_mp_cards(self) -> None:
388 """A reconstructed legacy MaskedImage re-serializes with ``MP_`` cards
389 whose bit indices match the written ``MSKN`` layout and select the
390 correct pixels, so a legacy reader stays correct.
391 """
392 hdul = self._legacy_full_hdu_list({"SUSPECT": (1, 2), "NO_DATA": (3, 4), "BAD": (0, 0)})
393 masked_image = MaskedImage.from_hdu_list(hdul)
394 with tempfile.TemporaryDirectory() as tmp:
395 path = os.path.join(tmp, "out.fits")
396 masked_image.write(path)
397 with astropy.io.fits.open(path) as out:
398 header = out["MASK"].header
399 array = np.asarray(out["MASK"].data)
400 # New-schema MSKN index for each plane name.
401 mskn = {header[k]: int(k.removeprefix("MSKN")) for k in header if k.startswith("MSKN")}
402 # Retained MP_ cards agree with the MSKN layout.
403 self.assertEqual(header["MP_SUSPECT"], mskn["SUSPECT"])
404 self.assertEqual(header["MP_NO_DATA"], mskn["NO_DATA"])
405 self.assertEqual(header["MP_SUSPECT"], 6)
406 self.assertNotIn("MP_DETECTED_NEGATIVE", header)
407 # The MP_ bit index selects the right pixel in the on-disk array.
408 self.assertTrue(array[1, 2] & (1 << header["MP_SUSPECT"]))
409 self.assertTrue(array[3, 4] & (1 << header["MP_NO_DATA"]))
410 self.assertTrue(array[0, 0] & (1 << header["MP_BAD"]))
412 def test_normal_read_strips_mp_cards(self) -> None:
413 """The normal ``lsst.images`` reader drops any ``MP_`` cards, so they
414 cannot drift out of sync with the authoritative ``MSKN`` schema or be
415 re-propagated on rewrite. They are written only for the legacy-cutout
416 afw-compatibility scenario.
417 """
418 hdul = self._legacy_full_hdu_list({"SUSPECT": (1, 2), "NO_DATA": (3, 4), "BAD": (0, 0)})
419 masked_image = MaskedImage.from_hdu_list(hdul)
420 with tempfile.TemporaryDirectory() as tmp:
421 legacy_cutout = os.path.join(tmp, "legacy_cutout.fits")
422 masked_image.write(legacy_cutout)
423 # The legacy-cutout file carries MP_ cards for afw.
424 with astropy.io.fits.open(legacy_cutout) as out:
425 self.assertTrue([k for k in out["MASK"].header if k.startswith("MP_")])
426 # A normal read + rewrite must not retain them.
427 rewritten = os.path.join(tmp, "rewritten.fits")
428 MaskedImage.read(legacy_cutout).write(rewritten)
429 with astropy.io.fits.open(rewritten) as out:
430 self.assertFalse([k for k in out["MASK"].header if k.startswith("MP_")])
433if __name__ == "__main__":
434 unittest.main()