Coverage for tests / test_visit_image.py: 14%
260 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:07 +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 os
15import unittest
16import warnings
17from typing import Any
19import astropy.io.fits
20import astropy.units as u
21import astropy.wcs
22import numpy as np
23from astro_metadata_translator import ObservationInfo
25from lsst.images import (
26 Box,
27 DetectorFrame,
28 Image,
29 MaskPlane,
30 MaskSchema,
31 ObservationSummaryStats,
32 ProjectionAstropyView,
33 TractFrame,
34 VisitImage,
35 get_legacy_visit_image_mask_planes,
36)
37from lsst.images.aperture_corrections import ApertureCorrectionMap, aperture_corrections_to_legacy
38from lsst.images.fields import ChebyshevField
39from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata
40from lsst.images.psfs import GaussianPointSpreadFunction, PointSpreadFunction
41from lsst.images.tests import (
42 DP2_VISIT_DETECTOR_DATA_ID,
43 RoundtripFits,
44 TemporaryButler,
45 assert_masked_images_equal,
46 assert_projections_equal,
47 compare_aperture_corrections_to_legacy,
48 compare_visit_image_to_legacy,
49 make_random_projection,
50)
52DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
55class VisitImageTestCase(unittest.TestCase):
56 """Basic Tests for VisitImage."""
58 @classmethod
59 def setUpClass(cls) -> None:
60 cls.rng = np.random.default_rng(500)
61 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096])
62 cls.projection = make_random_projection(cls.rng, det_frame, Box.factory[1:4096, 1:4096])
63 cls.mask_schema = MaskSchema([MaskPlane("M1", "D1")])
64 cls.obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4)
65 cls.summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
66 cls.gaussian_psf = GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13])
67 cls.aperture_corrections: ApertureCorrectionMap = {
68 "flux1": ChebyshevField(det_frame.bbox, np.array([0.75])),
69 "flux2": ChebyshevField(det_frame.bbox, np.array([0.625])),
70 }
72 opaque = FitsOpaqueMetadata()
73 hdr = astropy.io.fits.Header()
74 with warnings.catch_warnings():
75 # Silence warnings about long keys becoming HIERARCH.
76 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning)
77 hdr.update({"PLATFORM": "lsstcam", "LSST BUTLER ID": "123456789"})
78 opaque.extract_legacy_primary_header(hdr)
80 cls.image = Image(42, shape=(1024, 1024), unit=u.nJy)
81 cls.variance = Image(5.0, shape=(1024, 1024), unit=u.nJy * u.nJy)
82 # API signature suggests projection and obs_info can be None but they
83 # are required (unless you pass them in via the image plane).
84 cls.visit_image = VisitImage(
85 cls.image,
86 variance=cls.variance,
87 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
88 mask_schema=cls.mask_schema,
89 projection=cls.projection,
90 obs_info=cls.obs_info,
91 summary_stats=cls.summary_stats,
92 aperture_corrections=cls.aperture_corrections,
93 )
94 cls.visit_image._opaque_metadata = opaque
95 cls.simplest_visit_image = VisitImage(
96 cls.image,
97 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
98 mask_schema=cls.mask_schema,
99 projection=cls.projection,
100 obs_info=cls.obs_info,
101 )
103 def test_basics(self) -> None:
104 """Test basic constructor patterns."""
105 # Test default fill of variance.
106 visit = self.simplest_visit_image
107 self.assertEqual(visit.variance.array[0, 0], 1.0)
108 self.assertIs(visit[...], visit)
109 self.assertEqual(str(visit), "VisitImage(Image([y=0:1024, x=0:1024], int64), ['M1'])")
110 self.assertEqual(
111 repr(visit),
112 "VisitImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
113 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
114 " dtype=dtype('uint8')))",
115 )
117 astropy_wcs = visit.astropy_wcs
118 self.assertIsInstance(astropy_wcs, ProjectionAstropyView)
119 approx_wcs = visit.fits_wcs
120 self.assertIsInstance(approx_wcs, astropy.wcs.WCS)
122 with self.assertRaises(TypeError):
123 # Requires a PSF.
124 VisitImage(
125 self.image,
126 mask_schema=self.mask_schema,
127 projection=self.projection,
128 obs_info=self.obs_info,
129 )
131 with self.assertRaises(TypeError):
132 # Requires ObservationInfo.
133 VisitImage(
134 self.image,
135 psf=self.gaussian_psf,
136 mask_schema=self.mask_schema,
137 projection=self.projection,
138 )
140 with self.assertRaises(TypeError):
141 # Requires a projection.
142 VisitImage(
143 self.image,
144 psf=self.gaussian_psf,
145 mask_schema=self.mask_schema,
146 obs_info=self.obs_info,
147 )
149 with self.assertRaises(TypeError):
150 # Requires some form of mask.
151 VisitImage(
152 self.image,
153 psf=self.gaussian_psf,
154 projection=self.projection,
155 obs_info=self.obs_info,
156 )
158 with self.assertRaises(TypeError):
159 VisitImage(
160 Image(42, shape=(5, 5)),
161 psf=self.gaussian_psf,
162 mask_schema=self.mask_schema,
163 projection=self.projection,
164 obs_info=self.obs_info,
165 )
167 # Requires a DetectorFrame.
168 tract_frame = TractFrame(skymap="Skymap", tract=1, bbox=Box.factory[1:10, 1:10])
169 tract_proj = make_random_projection(self.rng, tract_frame, Box.factory[1:4096, 1:4096])
170 with self.assertRaises(TypeError):
171 VisitImage(
172 self.image,
173 projection=tract_proj,
174 psf=self.gaussian_psf,
175 mask_schema=self.mask_schema,
176 obs_info=self.obs_info,
177 )
179 # Variance unit mismatch.
180 with self.assertRaises(ValueError):
181 VisitImage(
182 self.image,
183 variance=self.image,
184 psf=self.gaussian_psf,
185 mask_schema=self.mask_schema,
186 projection=self.projection,
187 obs_info=self.obs_info,
188 )
190 def test_copy_and_slice(self) -> None:
191 """Test that arrays and components are copied (when not immutable) by
192 'copy' and referenced by 'slice'.
193 """
194 visit = self.visit_image
195 copy = visit.copy()
196 copy.image.array[0, 0] = 30.0
197 self.assertEqual(visit.image.array[0, 0], 42.0)
198 self.assertEqual(copy.image.array[0, 0], 30.0)
199 subvisit = visit[Box.factory[0:5, 0:5]]
200 # Check summary stats.
201 self.assertEqual(copy.summary_stats, visit.summary_stats)
202 self.assertIsNot(copy.summary_stats, visit.summary_stats)
203 self.assertEqual(subvisit.summary_stats, visit.summary_stats)
204 self.assertIs(subvisit.summary_stats, visit.summary_stats)
205 # Check aperture corrections.
206 self.assertEqual(copy.aperture_corrections.keys(), visit.aperture_corrections.keys())
207 self.assertIsNot(copy.aperture_corrections, visit.aperture_corrections)
208 self.assertEqual(subvisit.aperture_corrections.keys(), visit.aperture_corrections.keys())
209 self.assertIs(subvisit.aperture_corrections, visit.aperture_corrections)
211 def test_obs_info(self) -> None:
212 """Check that ObservationInfo has been constructed."""
213 visit = self.visit_image
214 self.assertIsNotNone(visit.obs_info)
215 self.maxDiff = None
216 assert visit.obs_info is not None # for mypy.
217 self.assertEqual(visit.obs_info.instrument, "LSSTCam")
219 def test_summary_stats(self) -> None:
220 """Test the comparisons and attributes of ObservationSummaryStats."""
221 self.assertEqual(self.summary_stats, ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4))
222 self.assertNotEqual(self.summary_stats, ObservationSummaryStats(psfSigma=2.5))
223 self.assertNotEqual(
224 self.summary_stats, ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
225 )
227 def test_read_write(self) -> None:
228 """Test that a visit can round trip through a FITS file."""
229 with RoundtripFits(self, self.visit_image, "VisitImage") as roundtrip:
230 # Check that we're still using the right compression, and that we
231 # wrote WCSs.
232 fits = roundtrip.inspect()
233 self.assertEqual(fits[1].header["ZCMPTYPE"], "GZIP_2")
234 self.assertEqual(fits[1].header["CTYPE1"], "RA---TAN")
235 self.assertEqual(fits[2].header["ZCMPTYPE"], "GZIP_2")
236 self.assertEqual(fits[2].header["CTYPE1"], "RA---TAN")
237 self.assertEqual(fits[3].header["ZCMPTYPE"], "GZIP_2")
238 self.assertEqual(fits[3].header["CTYPE1"], "RA---TAN")
239 # Check a subimage read.
240 subbox = Box.factory[8:13, 9:30]
241 subimage = roundtrip.get(bbox=subbox)
242 assert_masked_images_equal(self, subimage, self.visit_image[subbox], expect_view=False)
243 with self.subTest():
244 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
245 with self.subTest():
246 obs_info = roundtrip.get("obs_info")
247 self.assertIsInstance(obs_info, ObservationInfo)
248 self.assertEqual(obs_info, self.visit_image.obs_info)
249 with self.subTest():
250 summary_stats = roundtrip.get("summary_stats")
251 self.assertIsInstance(summary_stats, ObservationSummaryStats)
252 self.assertEqual(summary_stats, self.visit_image.summary_stats)
253 with self.subTest():
254 psf = roundtrip.get("psf")
255 self.assertIsInstance(psf, GaussianPointSpreadFunction)
256 self.assertEqual(psf.kernel_bbox, self.gaussian_psf.kernel_bbox)
258 assert_masked_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
259 # Check that the round-tripped headers are the same (up to card order).
260 self.assertEqual(len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]), 1)
261 self.assertEqual(
262 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
263 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
264 )
265 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
266 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
267 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
268 self.assertEqual(roundtrip.result.obs_info, self.visit_image.obs_info)
269 self.assertIsNotNone(roundtrip.result.summary_stats)
270 self.assertEqual(
271 roundtrip.result.summary_stats.psfSigma,
272 self.visit_image.summary_stats.psfSigma,
273 )
274 self.assertEqual(
275 roundtrip.result.summary_stats.zeroPoint,
276 self.visit_image.summary_stats.zeroPoint,
277 )
280@unittest.skipUnless(DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
281class VisitImageLegacyTestCase(unittest.TestCase):
282 """Tests for the VisitImage class and the basics of the archive system.
284 Requires legacy code.
285 """
287 @classmethod
288 def setUpClass(cls) -> None:
289 assert DATA_DIR is not None, "Guaranteed by decorator."
290 cls.filename = os.path.join(DATA_DIR, "dp2", "legacy", "visit_image.fits")
291 try:
292 from lsst.afw.image import ExposureFitsReader
294 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
295 except ImportError:
296 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
297 cls.plane_map = plane_map = get_legacy_visit_image_mask_planes()
298 cls.visit_image = VisitImage.read_legacy(
299 cls.filename, preserve_quantization=True, plane_map=plane_map
300 )
302 def test_legacy_errors(self) -> None:
303 """Legacy read failure modes."""
304 with self.assertRaises(ValueError):
305 VisitImage.from_legacy(self.legacy_exposure, instrument="HSC")
306 with self.assertRaises(ValueError):
307 VisitImage.from_legacy(self.legacy_exposure, visit=123456)
308 with self.assertRaises(ValueError):
309 VisitImage.from_legacy(self.legacy_exposure, unit=u.mJy)
310 visit = VisitImage.from_legacy(
311 self.legacy_exposure, instrument="LSSTCam", unit=u.nJy, visit=2025052000177
312 )
313 self.assertEqual(visit.unit, u.nJy)
315 with self.assertRaises(ValueError):
316 VisitImage.read_legacy(self.filename, instrument="HSC")
317 with self.assertRaises(ValueError):
318 VisitImage.read_legacy(self.filename, visit=123456)
320 def test_component_reads(self) -> None:
321 """Test reads of components from legacy file."""
322 visit = VisitImage.read_legacy(self.filename)
323 proj = VisitImage.read_legacy(self.filename, component="projection")
324 assert_projections_equal(self, proj, visit.projection, expect_identity=False)
325 image = VisitImage.read_legacy(self.filename, component="image")
326 self.assertEqual(image, visit.image)
327 self.check_legacy_obs_info(image.obs_info)
328 assert_projections_equal(self, proj, image.projection, expect_identity=False)
329 variance = VisitImage.read_legacy(self.filename, component="variance")
330 self.assertEqual(variance, visit.variance)
331 assert_projections_equal(self, proj, variance.projection, expect_identity=False)
332 self.check_legacy_obs_info(variance.obs_info)
333 mask = VisitImage.read_legacy(self.filename, component="mask")
334 self.assertEqual(mask, visit.mask)
335 assert_projections_equal(self, proj, mask.projection, expect_identity=False)
336 self.check_legacy_obs_info(mask.obs_info)
337 psf = VisitImage.read_legacy(self.filename, component="psf")
338 self.assertIsInstance(psf, PointSpreadFunction)
339 obs_info = VisitImage.read_legacy(self.filename, component="obs_info")
340 self.check_legacy_obs_info(obs_info)
341 summary_stats = VisitImage.read_legacy(self.filename, component="summary_stats")
342 self.assertIsInstance(summary_stats, ObservationSummaryStats)
343 self.assertEqual(summary_stats.nPsfStar, 93)
344 compare_aperture_corrections_to_legacy(
345 self,
346 VisitImage.read_legacy(self.filename, component="aperture_corrections"),
347 self.legacy_exposure.info.getApCorrMap(),
348 visit.bbox,
349 )
351 def check_legacy_obs_info(self, obs_info: ObservationInfo | None) -> None:
352 """Check that an `ObservationInfo` instance is not `None`, and that it
353 matches the one in the legacy test data file.
354 """
355 self.assertIsInstance(obs_info, ObservationInfo)
356 self.assertEqual(obs_info.instrument, "LSSTCam")
357 self.assertEqual(obs_info.detector_num, 85, obs_info)
358 self.assertEqual(obs_info.detector_unique_name, "R21_S11", obs_info)
359 self.assertEqual(obs_info.physical_filter, "r_57", obs_info)
361 def test_obs_info(self) -> None:
362 """Check that ObservationInfo has been constructed."""
363 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
364 self.assertIsNotNone(legacy.obs_info)
365 self.maxDiff = None
366 self.assertEqual(legacy.obs_info, self.visit_image.obs_info)
367 assert legacy.obs_info is not None # for mypy.
368 self.assertEqual(legacy.obs_info.instrument, "LSSTCam")
369 self.assertEqual(legacy.obs_info.detector_num, 85, legacy.obs_info)
370 self.assertEqual(legacy.obs_info.detector_unique_name, "R21_S11", legacy.obs_info)
371 self.assertEqual(legacy.obs_info.physical_filter, "r_57", legacy.obs_info)
373 def test_aperture_corrections_to_legacy(self) -> None:
374 """Test that we can convert an aperture correction map back to a
375 legacy `lsst.afw.image.ApCorrMap`.
376 """
377 legacy_ap_corr_map = aperture_corrections_to_legacy(self.visit_image.aperture_corrections)
378 compare_aperture_corrections_to_legacy(
379 self, self.visit_image.aperture_corrections, legacy_ap_corr_map, self.visit_image.bbox
380 )
382 def test_read_legacy_headers(self) -> None:
383 """Test that headers were correctly stripped and interpreted in
384 `VisitImage.read_legacy`.
385 """
386 # Check that we read the units from BUNIT.
387 self.assertEqual(self.visit_image.unit, astropy.units.nJy)
388 # Check that the primary header has the keys we want, and none of the
389 # keys we don't want.
390 header = self.visit_image._opaque_metadata.headers[ExtensionKey()]
391 self.assertIn("EXPTIME", header)
392 self.assertEqual(header["PLATFORM"], "lsstcam")
393 self.assertNotIn("LSST BUTLER ID", header)
394 self.assertNotIn("AR HDU", header)
395 self.assertNotIn("A_ORDER", header)
396 # Check that the extension HDUs do not have any custom headers.
397 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
398 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
399 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
401 def test_from_legacy_headers(self) -> None:
402 """Test that from_legacy handles headers properly."""
403 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
404 header = legacy._opaque_metadata.headers[ExtensionKey()]
405 self.assertIn("EXPTIME", header)
406 self.assertEqual(header["PLATFORM"], "lsstcam")
407 self.assertNotIn("LSST BUTLER ID", header)
408 self.assertNotIn("AR HDU", header)
409 self.assertNotIn("A_ORDER", header)
410 # Check that the extension HDUs do not have any custom headers.
411 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
412 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
413 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
415 def test_rewrite(self) -> None:
416 """Test that we can rewrite the visit image and preserve both
417 lossy-compressed pixel values and components exactly.
418 """
419 with RoundtripFits(self, self.visit_image, "VisitImage") as roundtrip:
420 # Check that we're still using the right compression, and that we
421 # wrote WCSs.
422 fits = roundtrip.inspect()
423 self.assertEqual(fits[1].header["ZCMPTYPE"], "RICE_1")
424 self.assertEqual(fits[1].header["CTYPE1"], "RA---TAN-SIP")
425 self.assertEqual(fits[2].header["ZCMPTYPE"], "GZIP_2")
426 self.assertEqual(fits[2].header["CTYPE1"], "RA---TAN-SIP")
427 self.assertEqual(fits[3].header["ZCMPTYPE"], "RICE_1")
428 self.assertEqual(fits[3].header["CTYPE1"], "RA---TAN-SIP")
429 # Check a subimage read.
430 subbox = Box.factory[8:13, 9:30]
431 subimage = roundtrip.get(bbox=subbox)
432 assert_masked_images_equal(self, subimage, self.visit_image[subbox], expect_view=False)
433 alternates: dict[str, Any] = {}
434 with self.subTest():
435 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
436 alternates = {
437 k: roundtrip.get(k)
438 for k in [
439 "projection",
440 "image",
441 "mask",
442 "variance",
443 "psf",
444 "obs_info",
445 "summary_stats",
446 "aperture_corrections",
447 ]
448 }
449 # Try to do a butler get of a component with storage class
450 # override.
451 with self.subTest():
452 if self.legacy_exposure is not None:
453 import lsst.afw.image
455 # We have VisitInfo available.
456 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
457 self.assertIsInstance(visit_info, lsst.afw.image.VisitInfo)
458 self.assertEqual(visit_info.getInstrumentLabel(), "LSSTCam")
459 else:
460 raise unittest.SkipTest("Can not test VisitInfo conversion without afw")
462 assert_masked_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
463 # Check that the round-tripped headers are the same (up to card order).
464 self.assertEqual(
465 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
466 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
467 )
468 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
469 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
470 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
471 self.assertEqual(roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"], "lsstcam")
472 compare_visit_image_to_legacy(
473 self,
474 roundtrip.result,
475 self.legacy_exposure,
476 expect_view=False,
477 plane_map=self.plane_map,
478 **DP2_VISIT_DETECTOR_DATA_ID,
479 alternates=alternates,
480 )
481 # Check converting from the legacy object in-memory.
482 compare_visit_image_to_legacy(
483 self,
484 VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map),
485 self.legacy_exposure,
486 expect_view=True,
487 plane_map=self.plane_map,
488 **DP2_VISIT_DETECTOR_DATA_ID,
489 )
491 def test_butler_converters(self) -> None:
492 """Test that we can read a VisitImage and its components from a butler
493 dataset written as an `lsst.afw.image.Exposure`.
494 """
495 if self.legacy_exposure is None:
496 raise unittest.SkipTest("lsst.afw.image.afw could not be imported.")
497 with TemporaryButler(legacy="ExposureF") as helper:
498 from lsst.daf.butler import FileDataset
500 helper.butler.ingest(FileDataset(path=self.filename, refs=[helper.legacy]), transfer="symlink")
501 visit_image_ref = helper.legacy.overrideStorageClass("VisitImage")
502 visit_image = helper.butler.get(visit_image_ref)
503 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
504 self.assertEqual(bbox, visit_image.bbox)
505 alternates = {
506 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
507 # TODO: including "projection" or "obs_info" here fails because
508 # there's code in daf_butler that expects any component to be
509 # valid for the *internal* storage class, not the requested
510 # one, and that's difficult to fix because it's tied up with
511 # the data ID standardization logic.
512 for k in ["image", "mask", "variance", "psf"]
513 }
514 compare_visit_image_to_legacy(
515 self,
516 visit_image,
517 self.legacy_exposure,
518 expect_view=False,
519 plane_map=self.plane_map,
520 alternates=alternates,
521 **DP2_VISIT_DETECTOR_DATA_ID,
522 )
525if __name__ == "__main__":
526 unittest.main()