Coverage for tests/test_masked_image.py: 92%
212 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-30 08:38 +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
37from lsst.images.tests import (
38 RoundtripFits,
39 RoundtripJson,
40 RoundtripNdf,
41 assert_masked_images_equal,
42 compare_masked_image_to_legacy,
43)
45try:
46 import h5py # noqa: F401
48 HAVE_H5PY = True
49except ImportError:
50 HAVE_H5PY = False
52try:
53 from lsst.afw.image import MaskedImageReader as LegacyMaskedImageReader
55except ImportError:
56 type LegacyMaskedImageReader = Any # type: ignore[no-redef]
58EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
60skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed")
63@dataclasses.dataclass
64class _LegacyTestData:
65 masked_image: MaskedImage
66 reader: LegacyMaskedImageReader
67 plane_map: dict[str, MaskPlane]
70@pytest.fixture(scope="session")
71def legacy_test_data() -> _LegacyTestData:
72 """Return a Mask read directly from the legacy test dataset and a legacy
73 reader for that image.
75 Skips if TESTDATA_IMAGES_DIR is unset or lsst.afw.image is unavailable.
76 """
77 if EXTERNAL_DATA_DIR is None: 77 ↛ 79line 77 didn't jump to line 79 because the condition on line 77 was always true
78 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
79 try:
80 from lsst.afw.image import MaskedImageFitsReader
81 except ImportError:
82 pytest.skip("'lsst.afw.image' could not be imported.")
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):
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):
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):
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):
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):
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() -> None:
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)
303@skip_no_h5py
304def test_round_trip_ndf_compatible_mask() -> None:
305 """Verify NDF round-trip for a MaskedImage with ≤8 mask planes."""
306 mi = make_masked_image()
307 with RoundtripNdf(mi, "MaskedImageV2") as roundtrip:
308 assert_masked_images_equal(roundtrip.result, mi, expect_view=False)
311@skip_no_h5py
312def test_round_trip_ndf_incompatible_mask() -> None:
313 """Verify NDF round-trip for a MaskedImage with more than 8 mask planes."""
314 rng = np.random.default_rng(7)
315 planes = [MaskPlane(f"P{i}", f"plane {i}") for i in range(12)]
316 wide = MaskedImage(
317 Image(
318 rng.normal(100.0, 8.0, size=(50, 60)),
319 dtype=np.float64,
320 unit=u.nJy,
321 yx0=(0, 0),
322 ),
323 mask_schema=MaskSchema(planes),
324 )
325 wide.variance.array = rng.normal(64.0, 0.5, size=wide.bbox.shape)
326 with RoundtripNdf(wide, "MaskedImageV2") as roundtrip:
327 assert_masked_images_equal(roundtrip.result, wide, expect_view=False)
330@skip_no_h5py
331def test_round_trip_ndf_many_plane_mask() -> None:
332 """Verify NDF round-trip for a mask that needs more than one int32
333 chunk.
334 """
335 rng = np.random.default_rng(11)
336 planes = [MaskPlane(f"P{i}", f"plane {i}") for i in range(40)]
337 wide = MaskedImage(
338 Image(
339 rng.normal(100.0, 8.0, size=(10, 12)),
340 dtype=np.float64,
341 unit=u.nJy,
342 yx0=(0, 0),
343 ),
344 mask_schema=MaskSchema(planes),
345 )
346 wide.mask.set("P0", wide.image.array > 100.0)
347 wide.mask.set("P17", wide.image.array < 95.0)
348 wide.mask.set("P39", wide.image.array > 110.0)
349 wide.variance.array = rng.normal(64.0, 0.5, size=wide.bbox.shape)
350 with RoundtripNdf(wide, "MaskedImageV2") as roundtrip:
351 assert_masked_images_equal(roundtrip.result, wide, expect_view=False)
354@skip_no_h5py
355def test_fits_ndf_consistency() -> None:
356 """Verify FITS and NDF backends produce equal MaskedImages on
357 round-trip.
358 """
359 mi = make_masked_image()
360 with (
361 RoundtripFits(mi) as fits_rt,
362 RoundtripNdf(mi) as ndf_rt,
363 ):
364 assert_masked_images_equal(mi, fits_rt.result, expect_view=False)
365 assert_masked_images_equal(mi, ndf_rt.result, expect_view=False)
366 assert_masked_images_equal(fits_rt.result, ndf_rt.result, expect_view=False)
369def test_fits_json_consistency() -> None:
370 """Verify FITS and JSON backends produce equal MaskedImages on
371 round-trip.
372 """
373 mi = make_masked_image()
374 with (
375 RoundtripFits(mi) as fits_rt,
376 RoundtripJson(mi) as json_rt,
377 ):
378 assert_masked_images_equal(mi, fits_rt.result, expect_view=False)
379 assert_masked_images_equal(mi, json_rt.result, expect_view=False)
380 assert_masked_images_equal(fits_rt.result, json_rt.result, expect_view=False)
383def test_legacy(legacy_test_data: _LegacyTestData) -> None:
384 """Test MaskedImage.read_legacy, MaskedImage.to_legacy, and
385 MaskedImage.from_legacy.
386 """
387 legacy_masked_image = legacy_test_data.reader.read()
388 compare_masked_image_to_legacy(
389 legacy_test_data.masked_image,
390 legacy_masked_image,
391 plane_map=legacy_test_data.plane_map,
392 expect_view=False,
393 )
394 compare_masked_image_to_legacy(
395 legacy_test_data.masked_image,
396 legacy_test_data.masked_image.to_legacy(plane_map=legacy_test_data.plane_map),
397 plane_map=legacy_test_data.plane_map,
398 expect_view=True,
399 )
400 compare_masked_image_to_legacy(
401 MaskedImage.from_legacy(legacy_masked_image, plane_map=legacy_test_data.plane_map),
402 legacy_masked_image,
403 expect_view=True,
404 plane_map=legacy_test_data.plane_map,
405 )
408def test_sky_circle_bbox() -> None:
409 """Test that we can extract a bounding box from a sky circle."""
410 mi = make_masked_image()
412 # This position is on the reference pixel (x=5, y=6), which is just
413 # outside the image (the bbox starts at x=8, y=5), so the box must be
414 # clipped on the low-x and low-y sides. 0.1 arcsec pixels.
415 bbox = mi.bbox_from_sky_circle(
416 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec), clip=True
417 )
418 # The circle has a ~10 pixel radius (1 arcsec at 0.1 arcsec per pixel),
419 # spanning x [-5, 15] and y [-4, 16] before clipping to the image bounds.
420 assert bbox == Box.factory[5:17, 8:16]
422 with pytest.raises(NotContainedError):
423 # Partially off the edge but clipping not requested.
424 mi.bbox_from_sky_circle(
425 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
426 )
428 # Fully inside the image. The image is only ~200 pixels across or
429 # ~20 arcsec.
430 bbox = mi.bbox_from_sky_circle(
431 SkyCoord(ra=12.0 * u.deg - 5.0 * u.arcsec, dec=13.0 * u.deg + 5.0 * u.arcsec, frame="icrs"),
432 Angle(1.0 * u.arcsec),
433 )
434 # The center is offset from the reference pixel by +50 pixels in y and
435 # by +48.7 pixels in x (the 5 arcsec RA offset scales by cos(dec)),
436 # placing it at (x=53.7, y=56) with a ~10 pixel radius.
437 assert bbox == Box.factory[46:67, 44:65]
439 # Fully off the image.
440 with pytest.raises(NotContainedError):
441 mi.bbox_from_sky_circle(
442 SkyCoord(ra=13.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
443 )
445 # Fully off the image with clipping requested.
446 with pytest.raises(NotContainedError):
447 mi.bbox_from_sky_circle(
448 SkyCoord(ra=13.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec), clip=True
449 )
451 # Non-scalar center and radius are rejected.
452 with pytest.raises(ValueError, match="scalar SkyCoord"):
453 mi.bbox_from_sky_circle(
454 SkyCoord(ra=[12.0, 12.1] * u.deg, dec=[13.0, 13.1] * u.deg, frame="icrs"),
455 Angle(1.0 * u.arcsec),
456 )
457 with pytest.raises(ValueError, match="scalar Angle"):
458 mi.bbox_from_sky_circle(
459 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle([1.0, 2.0] * u.arcsec)
460 )
462 # An image without a sky projection cannot calculate a bounding box.
463 no_wcs = Image(0.0, shape=(10, 10), dtype=np.float64)
464 with pytest.raises(ValueError, match="sky projection"):
465 no_wcs.bbox_from_sky_circle(
466 SkyCoord(ra=12.0 * u.deg, dec=13.0 * u.deg, frame="icrs"), Angle(1.0 * u.arcsec)
467 )