Coverage for tests/test_masked_image.py: 95%
250 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-30 11:45 +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.
12from __future__ import annotations
14import dataclasses
15import os
16from pathlib import Path
17from typing import Any
19import astropy.io.fits
20import astropy.units as u
21import numpy as np
22import pytest
23from astropy.coordinates import Angle, SkyCoord
25from lsst.images import (
26 Box,
27 GeneralFrame,
28 Image,
29 MaskedImage,
30 MaskPlane,
31 MaskSchema,
32 NotContainedError,
33 SkyProjection,
34 get_legacy_visit_image_mask_planes,
35)
36from lsst.images.fits import FitsCompressionOptions, FitsOpaqueMetadata
37from lsst.images.tests import (
38 RoundtripFits,
39 RoundtripJson,
40 RoundtripNdf,
41 assert_masked_images_equal,
42 compare_masked_image_to_legacy,
43 reset_afw_mask_planes, # noqa: F401
44)
46try:
47 import h5py # noqa: F401
49 HAVE_H5PY = True
50except ImportError:
51 HAVE_H5PY = False
53try:
54 from lsst.afw.image import MaskedImageReader as LegacyMaskedImageReader
56except ImportError:
57 type LegacyMaskedImageReader = Any # type: ignore[no-redef]
59EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
61skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed")
64@dataclasses.dataclass
65class _LegacyTestData:
66 masked_image: MaskedImage
67 reader: LegacyMaskedImageReader
68 plane_map: dict[str, MaskPlane]
71@pytest.fixture
72def legacy_test_data(reset_afw_mask_planes: None) -> _LegacyTestData: # noqa: F811
73 """Return a Mask read directly from the legacy test dataset and a legacy
74 reader for that image.
76 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.image is unavailable.
77 """
78 # reset_afw_mask_planes will have already skipped if afw is not available.
79 from lsst.afw.image import MaskedImageFitsReader
81 if EXTERNAL_DATA_DIR is None: 81 ↛ 83line 81 didn't jump to line 83 because the condition on line 81 was always true
82 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
83 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits")
84 plane_map = get_legacy_visit_image_mask_planes()
85 masked_image = MaskedImage.read_legacy(filename, plane_map=plane_map)
86 reader = MaskedImageFitsReader(filename)
87 return _LegacyTestData(masked_image=masked_image, reader=reader, plane_map=plane_map)
90def _make_wcs() -> astropy.wcs.WCS:
91 """Build a gnomonic FITS WCS with 0.1 arcsec pixels at (12, 13) deg.
93 The reference pixel is at 0-based pixel (x=5, y=6).
94 """
95 wcs = astropy.wcs.WCS(naxis=2)
96 # FITS CRPIX is 1-based, so CRPIX (6, 7) is 0-based pixel (x=5, y=6).
97 wcs.wcs.crpix = [6.0, 7.0]
98 wcs.wcs.crval = [12.0, 13.0]
99 scale = 0.1 / 3600.0
100 wcs.wcs.cd = [[-scale, 0.0], [0.0, scale]]
101 wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
102 return wcs
105def make_masked_image() -> MaskedImage:
106 """Return a freshly-constructed MaskedImage with BAD and HUNGRY mask
107 planes set.
108 """
109 rng = np.random.default_rng(500)
110 masked_image = MaskedImage(
111 Image(rng.normal(100.0, 8.0, size=(200, 251)), dtype=np.float64, unit=u.nJy, yx0=(5, 8)),
112 mask_schema=MaskSchema(
113 [
114 MaskPlane("BAD", "Pixel is very bad, possibly downright evil."),
115 MaskPlane("HUNGRY", "Pixel hasn't had enough to eat today."),
116 ]
117 ),
118 metadata={"fifty": "5 * 10"},
119 sky_projection=SkyProjection.from_fits_wcs(_make_wcs(), GeneralFrame(unit=u.pix)),
120 )
121 masked_image.mask.array |= np.multiply.outer(
122 masked_image.image.array < 102.0,
123 masked_image.mask.schema.bitmask("BAD"),
124 )
125 masked_image.mask.array |= np.multiply.outer(
126 masked_image.image.array > 98.0,
127 masked_image.mask.schema.bitmask("HUNGRY"),
128 )
129 masked_image.variance.array = rng.normal(64.0, 0.5, size=masked_image.bbox.shape)
130 return masked_image
133def test_masked_image_repr_str_pinned() -> None:
134 """Pin the exact str and repr output of a MaskedImage."""
135 mi = make_masked_image()
136 assert str(mi) == "MaskedImage(Image([y=5:205, x=8:259], float64), ['BAD', 'HUNGRY'])"
137 assert repr(mi) == (
138 "MaskedImage(Image(..., bbox=Box(y=Interval(start=5, stop=205), x=Interval(start=8, stop=259)), "
139 "dtype=dtype('float64')), mask_schema=MaskSchema([MaskPlane(name='BAD', description='Pixel is "
140 "very bad, possibly downright evil.'), MaskPlane(name='HUNGRY', description=\"Pixel hasn't had "
141 "enough to eat today.\")], dtype=dtype('uint8')))"
142 )
145def test_construction() -> None:
146 """Verify the MaskedImage constructed by make_masked_image has the
147 expected attributes.
148 """
149 mi = make_masked_image()
150 assert mi.bbox == Box.factory[5:205, 8:259]
151 assert mi.mask.bbox == mi.bbox
152 assert mi.variance.bbox == mi.bbox
153 assert mi.image.array.shape == mi.bbox.shape
154 assert mi.mask.array.shape == mi.bbox.shape + (1,)
155 assert mi.variance.array.shape == mi.bbox.shape
156 assert mi.unit == u.nJy
157 assert mi.variance.unit == u.nJy**2
158 assert mi.metadata == {"fifty": "5 * 10"}
159 # The checks below are subject to the vagaries of the RNG, but we want
160 # the seed to be such that they all pass, or other tests will be weaker.
161 assert np.sum(mi.mask.array == mi.mask.schema.bitmask("BAD")) > 0
162 assert np.sum(mi.mask.array == mi.mask.schema.bitmask("HUNGRY")) > 0
163 assert np.sum(mi.mask.array == mi.mask.schema.bitmask("BAD", "HUNGRY")) > 0
165 assert mi[...] is not mi
166 assert str(mi) == "MaskedImage(Image([y=5:205, x=8:259], float64), ['BAD', 'HUNGRY'])"
167 assert (
168 repr(mi)
169 == "MaskedImage(Image(..., bbox=Box(y=Interval(start=5, stop=205), x=Interval(start=8, stop=259)), "
170 "dtype=dtype('float64')), mask_schema=MaskSchema([MaskPlane(name='BAD', description='Pixel is "
171 "very bad, possibly downright evil.'), MaskPlane(name='HUNGRY', description=\"Pixel hasn't had "
172 "enough to eat today.\")], dtype=dtype('uint8')))"
173 )
174 copy = mi.copy()
175 original = mi.image.array[0, 0]
176 copy.image.array[0, 0] = 38.0
177 assert mi.image.array[0, 0] == original
178 assert copy.image.array[0, 0] == 38.0
180 # Test error conditions.
181 with pytest.raises(ValueError, match="bboxes do not agree"):
182 # Disagreement over mask bbox.
183 MaskedImage(Image(42.0, shape=(5, 6)), mask=mi.mask)
184 with pytest.raises(TypeError):
185 # No mask definition.
186 MaskedImage(mi.image, variance=mi.variance)
187 with pytest.raises(TypeError):
188 # Can not provide mask and mask schema.
189 MaskedImage(
190 Image(42.0, shape=(5, 5)),
191 mask=mi.mask,
192 mask_schema=mi.mask.schema,
193 )
194 with pytest.raises(ValueError, match="bboxes do not agree"):
195 # image and variance bbox disagreement.
196 MaskedImage(
197 Image(42.0, shape=(5, 5)),
198 mask_schema=mi.mask.schema,
199 variance=mi.variance,
200 )
201 with pytest.raises(ValueError, match="Image has no units but variance does"):
202 # no image unit but there is variance unit.
203 MaskedImage(
204 Image(42.0, shape=(5, 5)),
205 mask_schema=mi.mask.schema,
206 variance=Image(1.0, shape=(5, 5), unit=u.nJy),
207 )
208 with pytest.raises(ValueError, match="should be the square of the image unit"):
209 # image and variance units disagree.
210 MaskedImage(
211 Image(42.0, shape=(5, 5), unit=u.nJy),
212 mask_schema=mi.mask.schema,
213 variance=Image(1.0, shape=(5, 5), unit=u.nJy),
214 )
217def test_subset() -> None:
218 """Verify assignment of a subset into a MaskedImage copy."""
219 mi = make_masked_image()
220 copy = mi.copy()
221 subset = copy.local[0:10, 20:30].copy()
222 subset.image[...] = Image(42.0, shape=(10, 10), unit=u.nJy)
223 copy[subset.bbox] = subset
224 assert copy.image.array[0, 20] == 42.0
225 assert copy.image.array[0, 0] == mi.image.array[0, 0]
228def test_mask_setter() -> None:
229 """Verify the mask plane can be replaced with one grown by add_plane."""
230 mi = make_masked_image()
231 bad = mi.mask.get("BAD")
232 mi.mask = mi.mask.add_plane("OUTSIDE_STENCIL", "Pixel lies outside the stencil.")
233 assert "OUTSIDE_STENCIL" in mi.mask.schema.names
234 assert mi.mask.bbox == mi.image.bbox
235 np.testing.assert_array_equal(mi.mask.get("BAD"), bad)
236 assert not mi.mask.get("OUTSIDE_STENCIL").any()
237 # A mask whose bounding box disagrees with the image is rejected.
238 with pytest.raises(ValueError, match="bboxes do not agree"):
239 mi.mask = mi.mask[Box.factory[10:20, 12:22]]
242def test_fits_roundtrip() -> None:
243 """Verify MaskedImage round-trips correctly through FITS, including
244 subimage reads.
245 """
246 mi = make_masked_image()
247 subbox = Box.factory[11:20, 25:30]
248 subslices = (slice(6, 15), slice(17, 22))
249 np.testing.assert_array_equal(mi.image.array[subslices], mi.image[subbox].array)
250 with RoundtripFits(mi, "MaskedImageV2") as roundtrip:
251 subimage = roundtrip.get(bbox=subbox)
252 # Check that we used lossless compression (the default).
253 fits = roundtrip.inspect()
254 assert fits[1].header["ZCMPTYPE"] == "GZIP_2"
255 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
256 assert fits[3].header["ZCMPTYPE"] == "GZIP_2"
257 assert_masked_images_equal(roundtrip.result, mi, expect_view=False)
258 assert_masked_images_equal(subimage, roundtrip.result[subbox], expect_view=False)
261def test_fits_roundtrip_legacy_read(reset_afw_mask_planes: None) -> None: # noqa: F811
262 """Verify a round-tripped MaskedImageV2 can be read back as a legacy afw
263 MaskedImage.
264 """
265 try:
266 import lsst.afw.image
267 except ImportError:
268 pytest.skip("afw could not be imported")
269 mi = make_masked_image()
270 with RoundtripFits(mi, "MaskedImageV2") as roundtrip:
271 legacy_masked_image = roundtrip.get(storageClass="MaskedImage")
272 assert isinstance(legacy_masked_image, lsst.afw.image.MaskedImage)
273 compare_masked_image_to_legacy(mi, legacy_masked_image, expect_view=False)
276def test_fits_roundtrip_lossy(tmp_path: Path) -> None:
277 """Verify MaskedImage round-trips correctly through FITS with lossy
278 compression.
279 """
280 mi = make_masked_image()
281 subbox = Box.factory[11:20, 25:30]
282 subslices = (slice(6, 15), slice(17, 22))
283 np.testing.assert_array_equal(mi.image.array[subslices], mi.image[subbox].array)
284 path = tmp_path / "lossy.fits"
285 mi.write(
286 path,
287 compression_options={
288 "image": FitsCompressionOptions.LOSSY,
289 "variance": FitsCompressionOptions.LOSSY,
290 },
291 compression_seed=50,
292 )
293 roundtripped = MaskedImage.read(path)
294 subimage = MaskedImage.read(path, bbox=subbox)
295 with astropy.io.fits.open(path, disable_image_compression=True) as fits:
296 assert fits[1].header["ZCMPTYPE"] == "RICE_1"
297 assert fits[2].header["ZCMPTYPE"] == "GZIP_2"
298 assert fits[3].header["ZCMPTYPE"] == "RICE_1"
299 assert_masked_images_equal(roundtripped, mi, expect_view=False, rtol=0.01)
300 assert_masked_images_equal(subimage, roundtripped[subbox], expect_view=False)
303def test_fits_uncompressed_read_is_native_byte_order(tmp_path: Path) -> None:
304 """Verify that reading uncompressed (big-endian on disk) FITS HDUs
305 yields arrays in native byte order.
306 """
307 mi = make_masked_image()
308 path = tmp_path / "uncompressed.fits"
309 mi.write(path, compression_options={"image": None, "mask": None, "variance": None})
310 with astropy.io.fits.open(path) as hdu_list:
311 assert not hdu_list[1].data.dtype.isnative
312 full = MaskedImage.read(path)
313 subimage = MaskedImage.read(path, bbox=Box.factory[11:20, 25:30])
314 for result in (full, subimage):
315 assert result.image.array.dtype.isnative
316 assert result.mask.array.dtype.isnative
317 assert result.variance.array.dtype.isnative
318 assert_masked_images_equal(full, mi, expect_view=False)
321def test_legacy_uncompressed_read_is_native_byte_order(
322 tmp_path: Path,
323 reset_afw_mask_planes: None, # noqa: F811
324) -> None:
325 """Verify that reading an uncompressed legacy FITS file yields arrays
326 in native byte order that can be converted back to afw.
327 """
328 mi = MaskedImage(
329 Image(np.arange(20, dtype=np.float32).reshape(4, 5)),
330 mask_schema=MaskSchema([MaskPlane("BAD", "Pixel is bad.")]),
331 )
332 path = tmp_path / "legacy_uncompressed.fits"
333 mi.to_legacy().writeFits(str(path))
334 with astropy.io.fits.open(path) as hdu_list:
335 assert not hdu_list[1].data.dtype.isnative
336 masked_result = MaskedImage.read_legacy(path)
337 image_result = Image.read_legacy(path)
338 for image in (masked_result.image, masked_result.variance, image_result):
339 assert image.array.dtype.isnative
340 np.testing.assert_array_equal(masked_result.image.array, mi.image.array)
341 np.testing.assert_array_equal(image_result.array, mi.image.array)
342 masked_result.to_legacy()
343 image_result.to_legacy()
346@skip_no_h5py
347def test_round_trip_ndf_compatible_mask() -> None:
348 """Verify NDF round-trip for a MaskedImage with ≤8 mask planes."""
349 mi = make_masked_image()
350 with RoundtripNdf(mi, "MaskedImageV2") as roundtrip:
351 assert_masked_images_equal(roundtrip.result, mi, expect_view=False)
354@skip_no_h5py
355def test_round_trip_ndf_incompatible_mask() -> None:
356 """Verify NDF round-trip for a MaskedImage with more than 8 mask planes."""
357 rng = np.random.default_rng(7)
358 planes = [MaskPlane(f"P{i}", f"plane {i}") for i in range(12)]
359 wide = MaskedImage(
360 Image(
361 rng.normal(100.0, 8.0, size=(50, 60)),
362 dtype=np.float64,
363 unit=u.nJy,
364 yx0=(0, 0),
365 ),
366 mask_schema=MaskSchema(planes),
367 )
368 wide.variance.array = rng.normal(64.0, 0.5, size=wide.bbox.shape)
369 with RoundtripNdf(wide, "MaskedImageV2") as roundtrip:
370 assert_masked_images_equal(roundtrip.result, wide, expect_view=False)
373@skip_no_h5py
374def test_round_trip_ndf_many_plane_mask() -> None:
375 """Verify NDF round-trip for a mask that needs more than one int32
376 chunk.
377 """
378 rng = np.random.default_rng(11)
379 planes = [MaskPlane(f"P{i}", f"plane {i}") for i in range(40)]
380 wide = MaskedImage(
381 Image(
382 rng.normal(100.0, 8.0, size=(10, 12)),
383 dtype=np.float64,
384 unit=u.nJy,
385 yx0=(0, 0),
386 ),
387 mask_schema=MaskSchema(planes),
388 )
389 wide.mask.set("P0", wide.image.array > 100.0)
390 wide.mask.set("P17", wide.image.array < 95.0)
391 wide.mask.set("P39", wide.image.array > 110.0)
392 wide.variance.array = rng.normal(64.0, 0.5, size=wide.bbox.shape)
393 with RoundtripNdf(wide, "MaskedImageV2") as roundtrip:
394 assert_masked_images_equal(roundtrip.result, wide, expect_view=False)
397@skip_no_h5py
398def test_fits_ndf_consistency() -> None:
399 """Verify FITS and NDF backends produce equal MaskedImages on
400 round-trip.
401 """
402 mi = make_masked_image()
403 with (
404 RoundtripFits(mi) as fits_rt,
405 RoundtripNdf(mi) as ndf_rt,
406 ):
407 assert_masked_images_equal(mi, fits_rt.result, expect_view=False)
408 assert_masked_images_equal(mi, ndf_rt.result, expect_view=False)
409 assert_masked_images_equal(fits_rt.result, ndf_rt.result, expect_view=False)
412def test_fits_json_consistency() -> None:
413 """Verify FITS and JSON backends produce equal MaskedImages on
414 round-trip.
415 """
416 mi = make_masked_image()
417 with (
418 RoundtripFits(mi) as fits_rt,
419 RoundtripJson(mi) as json_rt,
420 ):
421 assert_masked_images_equal(mi, fits_rt.result, expect_view=False)
422 assert_masked_images_equal(mi, json_rt.result, expect_view=False)
423 assert_masked_images_equal(fits_rt.result, json_rt.result, expect_view=False)
426def test_legacy(legacy_test_data: _LegacyTestData) -> None:
427 """Test MaskedImage.read_legacy, MaskedImage.to_legacy, and
428 MaskedImage.from_legacy.
429 """
430 legacy_masked_image = legacy_test_data.reader.read()
431 compare_masked_image_to_legacy(
432 legacy_test_data.masked_image,
433 legacy_masked_image,
434 plane_map=legacy_test_data.plane_map,
435 expect_view=False,
436 )
437 compare_masked_image_to_legacy(
438 legacy_test_data.masked_image,
439 legacy_test_data.masked_image.to_legacy(plane_map=legacy_test_data.plane_map),
440 plane_map=legacy_test_data.plane_map,
441 expect_view=True,
442 )
443 compare_masked_image_to_legacy(
444 MaskedImage.from_legacy(legacy_masked_image, plane_map=legacy_test_data.plane_map),
445 legacy_masked_image,
446 expect_view=True,
447 plane_map=legacy_test_data.plane_map,
448 )
451def test_repeated_legacy_metadata_keys(reset_afw_mask_planes: None) -> None: # noqa: F811
452 """Test that a key present on more than one FITS card keeps all of its
453 values in the legacy metadata.
454 """
455 from lsst.daf.base import PropertyList
457 opaque_metadata = FitsOpaqueMetadata()
458 header = astropy.io.fits.Header()
459 # SubtractBackgroundTask writes one BGMEAN card per background fit.
460 header.append(("BGMEAN", 1.5), end=True)
461 header.append(("BGMEAN", 2.5), end=True)
462 header.append(("PLATFORM", "lsstcam"), end=True)
463 opaque_metadata.extract_legacy_primary_header(header)
464 masked_image = make_masked_image()
465 masked_image._opaque_metadata = opaque_metadata
466 legacy_metadata = PropertyList()
467 masked_image._fill_legacy_metadata(legacy_metadata)
468 assert legacy_metadata.getArray("BGMEAN") == [1.5, 2.5]
469 assert legacy_metadata["PLATFORM"] == "lsstcam"
472def test_sky_circle_bbox() -> None:
473 """Test that we can extract a bounding box from a sky circle."""
474 mi = make_masked_image()
476 # This position is on the reference pixel (x=5, y=6), which is just
477 # outside the image (the bbox starts at x=8, y=5), so the box must be
478 # clipped on the low-x and low-y sides. 0.1 arcsec pixels.
479 bbox = mi.bbox_from_sky_circle(
480 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec), clip=True
481 )
482 # The circle has a ~10 pixel radius (1 arcsec at 0.1 arcsec per pixel),
483 # spanning x [-5, 15] and y [-4, 16] before clipping to the image bounds.
484 assert bbox == Box.factory[5:17, 8:16]
486 with pytest.raises(NotContainedError):
487 # Partially off the edge but clipping not requested.
488 mi.bbox_from_sky_circle(
489 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
490 )
492 # Fully inside the image. The image is only ~200 pixels across or
493 # ~20 arcsec.
494 bbox = mi.bbox_from_sky_circle(
495 SkyCoord(ra=12.0 * u.deg - 5.0 * u.arcsec, dec=13.0 * u.deg + 5.0 * u.arcsec, frame="icrs"),
496 Angle(1.0 * u.arcsec),
497 )
498 # The center is offset from the reference pixel by +50 pixels in y and
499 # by +48.7 pixels in x (the 5 arcsec RA offset scales by cos(dec)),
500 # placing it at (x=53.7, y=56) with a ~10 pixel radius.
501 assert bbox == Box.factory[46:67, 44:65]
503 # Fully off the image.
504 with pytest.raises(NotContainedError):
505 mi.bbox_from_sky_circle(
506 SkyCoord(ra=13.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
507 )
509 # Fully off the image with clipping requested.
510 with pytest.raises(NotContainedError):
511 mi.bbox_from_sky_circle(
512 SkyCoord(ra=13.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec), clip=True
513 )
515 # Non-scalar center and radius are rejected.
516 with pytest.raises(ValueError, match="scalar SkyCoord"):
517 mi.bbox_from_sky_circle(
518 SkyCoord(ra=[12.0, 12.1] * u.deg, dec=[13.0, 13.1] * u.deg, frame="icrs"),
519 Angle(1.0 * u.arcsec),
520 )
521 with pytest.raises(ValueError, match="scalar Angle"):
522 mi.bbox_from_sky_circle(
523 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle([1.0, 2.0] * u.arcsec)
524 )
526 # An image without a sky projection cannot calculate a bounding box.
527 no_wcs = Image(0.0, shape=(10, 10), dtype=np.float64)
528 with pytest.raises(ValueError, match="sky projection"):
529 no_wcs.bbox_from_sky_circle(
530 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
531 )