Coverage for tests/test_imageCutoutsBackend.py: 23%
99 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-27 01:58 -0700
1# This file is part of dax_images_cutout.
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# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import os.path
23import tempfile
24import unittest
26import astropy.io.fits
27import numpy as np
29import lsst.daf.butler
30import lsst.images
31import lsst.sphgeom
32import lsst.utils.tests
33from lsst.dax.images.cutout import CutoutMode, ImageCutoutFactory, projection_finders, stencils
35try:
36 import lsst.afw.image
37 import lsst.geom
39 HAVE_AFW = True
40except ImportError:
41 HAVE_AFW = False
44@unittest.skipUnless(HAVE_AFW, "lsst.afw not available")
45class TestImageCutoutsBackend(lsst.utils.tests.TestCase):
46 """Tests for ImageCutoutsBackend."""
48 @classmethod
49 def setUpClass(cls):
50 try:
51 cls.data_dir = lsst.utils.getPackageDir("testdata_image_cutouts")
52 except LookupError:
53 raise unittest.SkipTest("testdata_image_cutouts not setup.") from None
55 def setUp(self):
56 collection = "2.2i/runs/test-med-1/w_2022_03/DM-33223/20220118T193330Z"
57 self.butler = lsst.daf.butler.Butler(os.path.join(self.data_dir, "repo"), collections=collection)
59 # Centered on a galaxy
60 point = lsst.sphgeom.LonLat.fromDegrees(56.6400770, -36.4492250)
61 radius = lsst.sphgeom.Angle((10 * lsst.geom.arcseconds).asRadians())
62 self.stencil = stencils.SkyCircle(point, radius)
64 self.projectionFinders = (
65 projection_finders.ReadComponents(),
66 projection_finders.ReadComponentsAstropyFits(),
67 )
69 self.dataId = {"patch": 24, "tract": 3828, "band": "r", "skymap": "DC2"}
70 self.ref = self.butler.find_dataset("deepCoadd_calexp", data_id=self.dataId)
72 def test_extract_ref(self):
73 """Test that extract_ref produces a reasonable cutout."""
74 with tempfile.TemporaryDirectory() as tempdir:
75 for cutout_mode in CutoutMode:
76 # Try each available projection finder at least once.
77 proj_finder = self.projectionFinders[cutout_mode.value % len(self.projectionFinders)]
78 cutoutBackend = ImageCutoutFactory(self.butler, proj_finder, tempdir)
79 result = cutoutBackend.extract_ref(self.stencil, self.ref, cutout_mode=cutout_mode)
80 # The galaxy should be near the center of the image.
81 match result.cutout:
82 case lsst.afw.image.Exposure() | lsst.afw.image.MaskedImage():
83 box = result.cutout.getBBox()
84 array = result.cutout.image.array
85 case lsst.afw.image.Image():
86 array = result.cutout.array
87 box = result.cutout.getBBox()
88 case lsst.images.MaskedImage():
89 box = result.cutout.bbox.to_legacy()
90 array = result.cutout.image.array
91 case lsst.images.Image():
92 array = result.cutout.array
93 box = result.cutout.bbox.to_legacy()
94 case _:
95 raise RuntimeError(f"Unexpected cutout type: {type(result.cutout)}")
96 self.assertEqual(box.width, 101)
97 self.assertEqual(box.height, 101)
98 self.assertFloatsAlmostEqual(array[50, 49], 2.083247661590576)
99 self.assertFloatsAlmostEqual(array[1, 1], -0.14575983583927155)
100 self.assertFloatsAlmostEqual(array[100, 100], 0.08674515783786774)
102 output = cutoutBackend.process_ref(self.stencil, self.ref, cutout_mode=cutout_mode)
103 self.assertTrue(output.exists())
105 output = cutoutBackend.process_uuid(self.stencil, self.ref.id, cutout_mode=cutout_mode)
106 self.assertTrue(output.exists())
108 def test_provenance_in_primary_header(self):
109 """Cutout provenance must land in the primary FITS header for every
110 cutout mode, including the native ``lsst.images`` container modes.
111 """
112 provenance_keys = ("CUTVERS",)
113 # Card comments that fit within the 80-column FITS card limit and so
114 # survive the round trip, covering both the provenance keys (DATE-CUT)
115 # and the stencil keys (ST_TYPE). Some cards can potentially have
116 # truncated comments.
117 provenance_comments = {
118 "BTLRUUID": "Butler UUID of full image",
119 "BTLRNAME": "Butler dataset type",
120 "ST_TYPE": "Type of stencil used to create this cutout",
121 "ST_RA": "[deg] Circle center Right Ascension",
122 "DATE-CUT": "Time of cutout extraction",
123 }
125 with tempfile.TemporaryDirectory() as tempdir:
126 cutoutBackend = ImageCutoutFactory(self.butler, self.projectionFinders[0], tempdir)
127 for cutout_mode in CutoutMode:
128 with self.subTest(cutout_mode=str(cutout_mode)):
129 output = cutoutBackend.process_ref(self.stencil, self.ref, cutout_mode=cutout_mode)
130 with output.open("rb") as fh, astropy.io.fits.open(fh) as hdul:
131 header = hdul[0].header
132 for key in provenance_keys:
133 self.assertIn(
134 key,
135 header,
136 f"{key} missing from primary header for {cutout_mode}",
137 )
138 self.assertEqual(header["BTLRUUID"], self.ref.id.hex)
139 self.assertEqual(header["BTLRNAME"], self.ref.datasetType.name)
140 for key, comment in provenance_comments.items():
141 self.assertEqual(
142 header.comments[key],
143 comment,
144 f"{key} comment wrong in primary header for {cutout_mode}",
145 )
147 # This test file has a BGMEAN header in the original
148 # primary header and we need to ensure that it still
149 # exists in the cutout to indicate that the primary
150 # header was propagated. IMAGE does not have it since
151 # IMAGE only reads the second header without merging
152 if cutout_mode != CutoutMode.IMAGE_ONLY:
153 self.assertIn("BGMEAN", header)
155 def test_astropy_masked_image_planes_readable_by_afw(self):
156 """The ``ASTROPY_MASKED_IMAGE`` backend re-serializes the cutout as an
157 ``lsst.images`` file but keeps the legacy ``MP_`` mask-plane cards
158 (re-indexed to the reshuffled schema). ``lsst.afw.image`` can
159 therefore read the cutout, and each mask plane covers the same pixels
160 as the afw-only ``MASKED_IMAGE`` backend even though the bit numbering
161 may differ between the two.
162 """
163 with tempfile.TemporaryDirectory() as tempdir:
164 cutoutBackend = ImageCutoutFactory(self.butler, self.projectionFinders[0], tempdir)
165 astropy_uri = cutoutBackend.process_ref(
166 self.stencil, self.ref, cutout_mode=CutoutMode.ASTROPY_MASKED_IMAGE
167 )
168 afw_uri = cutoutBackend.process_ref(self.stencil, self.ref, cutout_mode=CutoutMode.MASKED_IMAGE)
169 # Both files are read back through afw; the astropy cutout is an
170 # lsst.images file, so this exercises afw reading the re-indexed
171 # MP_ cards.
172 astropy_mask = lsst.afw.image.MaskedImageF(astropy_uri.ospath).mask
173 afw_mask = lsst.afw.image.MaskedImageF(afw_uri.ospath).mask
175 common_planes = set(astropy_mask.getMaskPlaneDict()) & set(afw_mask.getMaskPlaneDict())
176 self.assertTrue(common_planes, "afw read no shared mask planes from the astropy cutout")
178 any_set = False
179 for plane in sorted(common_planes):
180 astropy_set = (astropy_mask.array & astropy_mask.getPlaneBitMask(plane)) != 0
181 afw_set = (afw_mask.array & afw_mask.getPlaneBitMask(plane)) != 0
182 np.testing.assert_array_equal(
183 astropy_set, afw_set, err_msg=f"mask plane {plane!r} coverage differs"
184 )
185 any_set = any_set or bool(afw_set.any())
186 self.assertTrue(any_set, "expected at least one mask plane to be set in the cutout")
189if __name__ == "__main__": 189 ↛ 190line 189 didn't jump to line 190 because the condition on line 189 was never true
190 unittest.main()