Coverage for tests/test_visit_image.py: 53%
471 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-01 09:14 +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, ClassVar
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 BackgroundMap,
27 Box,
28 DetectorFrame,
29 DifferenceImage,
30 Image,
31 MaskPlane,
32 MaskSchema,
33 ObservationSummaryStats,
34 Polygon,
35 SkyProjectionAstropyView,
36 TractFrame,
37 VisitImage,
38 get_legacy_difference_image_mask_planes,
39 get_legacy_visit_image_mask_planes,
40)
41from lsst.images.aperture_corrections import ApertureCorrectionMap, aperture_corrections_to_legacy
42from lsst.images.cameras import Detector
43from lsst.images.fields import ChebyshevField, SplineField, SumField, field_from_legacy_photo_calib
44from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata
45from lsst.images.psfs import GaussianPointSpreadFunction, PointSpreadFunction
46from lsst.images.serialization import read
47from lsst.images.tests import (
48 DP2_VISIT_DETECTOR_DATA_ID,
49 RoundtripFits,
50 RoundtripJson,
51 RoundtripNdf,
52 TemporaryButler,
53 assert_close,
54 assert_masked_images_equal,
55 assert_sky_projections_equal,
56 assert_visit_images_equal,
57 compare_aperture_corrections_to_legacy,
58 compare_detector_to_legacy,
59 compare_photo_calib_to_legacy,
60 compare_visit_image_to_legacy,
61 make_random_sky_projection,
62)
64try:
65 import h5py # noqa: F401
67 HAVE_H5PY = True
68except ImportError:
69 HAVE_H5PY = False
71EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
72LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data")
75class VisitImageTestCase(unittest.TestCase):
76 """Basic Tests for VisitImage."""
78 @classmethod
79 def setUpClass(cls) -> None:
80 cls.rng = np.random.default_rng(500)
81 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096])
82 cls.mask_schema = MaskSchema([MaskPlane("M1", "D1")])
83 cls.obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1")
84 cls.summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4)
85 cls.gaussian_psf = GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13])
86 cls.aperture_corrections: ApertureCorrectionMap = {
87 "flux1": ChebyshevField(det_frame.bbox, np.array([0.75])),
88 "flux2": ChebyshevField(det_frame.bbox, np.array([0.625])),
89 }
90 cls.detector = read(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector)
92 opaque = FitsOpaqueMetadata()
93 hdr = astropy.io.fits.Header()
94 with warnings.catch_warnings():
95 # Silence warnings about long keys becoming HIERARCH.
96 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning)
97 hdr.update({"PLATFORM": "lsstcam", "LSST BUTLER ID": "123456789"})
98 opaque.extract_legacy_primary_header(hdr)
100 cls.image = Image(42, shape=(1024, 1024), unit=u.nJy)
101 cls.variance = Image(5.0, shape=(1024, 1024), unit=u.nJy * u.nJy)
102 # polygon is the lower triangle of the image.
103 cls.polygon = Polygon(x_vertices=[-0.5, 1023.5, -0.5], y_vertices=[-0.5, -0.5, 1023.5])
104 cls.sky_projection = make_random_sky_projection(cls.rng, det_frame, Box.factory[1:4096, 1:4096])
105 # API signature suggests sky_projection and obs_info can be None but
106 # they are required (unless you pass them in via the image plane).
107 cls.visit_image = VisitImage(
108 cls.image,
109 variance=cls.variance,
110 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
111 mask_schema=cls.mask_schema,
112 sky_projection=cls.sky_projection,
113 obs_info=cls.obs_info,
114 summary_stats=cls.summary_stats,
115 detector=cls.detector,
116 bounds=cls.polygon,
117 aperture_corrections=cls.aperture_corrections,
118 band="r",
119 )
120 cls.visit_image.backgrounds.add(
121 "standard",
122 ChebyshevField(det_frame.bbox, np.array([[2.0]])),
123 description="Background subtracted from the image.",
124 is_subtracted=True,
125 )
126 cls.visit_image._opaque_metadata = opaque
127 cls.simplest_visit_image = VisitImage(
128 cls.image,
129 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]),
130 mask_schema=cls.mask_schema,
131 sky_projection=cls.sky_projection,
132 detector=cls.detector,
133 obs_info=cls.obs_info,
134 band="r",
135 )
137 def test_basics(self) -> None:
138 """Test basic constructor patterns."""
139 # Test default fill of variance.
140 visit = self.simplest_visit_image
141 self.assertEqual(visit.variance.array[0, 0], 1.0)
142 self.assertIs(visit[...], visit)
143 self.assertEqual(str(visit), "VisitImage(Image([y=0:1024, x=0:1024], int64), ['M1'])")
144 self.assertEqual(
145 repr(visit),
146 "VisitImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024)),"
147 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')],"
148 " dtype=dtype('uint8')))",
149 )
151 astropy_wcs = visit.astropy_wcs
152 self.assertIsInstance(astropy_wcs, SkyProjectionAstropyView)
153 approx_wcs = visit.fits_wcs
154 self.assertIsInstance(approx_wcs, astropy.wcs.WCS)
156 with self.assertRaises(TypeError):
157 # Requires a PSF.
158 VisitImage(
159 self.image,
160 mask_schema=self.mask_schema,
161 sky_projection=self.sky_projection,
162 obs_info=self.obs_info,
163 detector=self.detector,
164 band="r",
165 )
167 with self.assertRaises(TypeError):
168 # Requires ObservationInfo.
169 VisitImage(
170 self.image,
171 psf=self.gaussian_psf,
172 mask_schema=self.mask_schema,
173 sky_projection=self.sky_projection,
174 detector=self.detector,
175 band="r",
176 )
178 with self.assertRaises(TypeError):
179 # Requires a sky_projection.
180 VisitImage(
181 self.image,
182 psf=self.gaussian_psf,
183 mask_schema=self.mask_schema,
184 obs_info=self.obs_info,
185 detector=self.detector,
186 band="r",
187 )
189 with self.assertRaises(TypeError):
190 # Requires a detector.
191 VisitImage(
192 self.image,
193 psf=self.gaussian_psf,
194 mask_schema=self.mask_schema,
195 sky_projection=self.sky_projection,
196 obs_info=self.obs_info,
197 band="r",
198 )
200 with self.assertRaises(TypeError):
201 # Requires some form of mask.
202 VisitImage(
203 self.image,
204 psf=self.gaussian_psf,
205 sky_projection=self.sky_projection,
206 obs_info=self.obs_info,
207 detector=self.detector,
208 band="r",
209 )
211 with self.assertRaises(TypeError):
212 VisitImage(
213 Image(42, shape=(5, 5)),
214 psf=self.gaussian_psf,
215 mask_schema=self.mask_schema,
216 sky_projection=self.sky_projection,
217 obs_info=self.obs_info,
218 detector=self.detector,
219 band="r",
220 )
222 # Requires a DetectorFrame.
223 tract_frame = TractFrame(skymap="Skymap", tract=1, bbox=Box.factory[1:10, 1:10])
224 tract_proj = make_random_sky_projection(self.rng, tract_frame, Box.factory[1:4096, 1:4096])
225 with self.assertRaises(TypeError):
226 VisitImage(
227 self.image,
228 sky_projection=tract_proj,
229 psf=self.gaussian_psf,
230 mask_schema=self.mask_schema,
231 obs_info=self.obs_info,
232 detector=self.detector,
233 band="r",
234 )
236 # Variance unit mismatch.
237 with self.assertRaises(ValueError):
238 VisitImage(
239 self.image,
240 variance=self.image,
241 psf=self.gaussian_psf,
242 mask_schema=self.mask_schema,
243 sky_projection=self.sky_projection,
244 obs_info=self.obs_info,
245 detector=self.detector,
246 band="r",
247 )
249 def test_copy_and_slice(self) -> None:
250 """Test that arrays and components are copied (when not immutable) by
251 'copy' and referenced by 'slice'.
252 """
253 visit = self.visit_image
254 copy = visit.copy()
255 copy.image.array[0, 0] = 30.0
256 self.assertEqual(visit.image.array[0, 0], 42.0)
257 self.assertEqual(copy.image.array[0, 0], 30.0)
258 subvisit = visit[Box.factory[0:5, 0:5]]
259 # Check summary stats.
260 self.assertEqual(copy.summary_stats, visit.summary_stats)
261 self.assertIsNot(copy.summary_stats, visit.summary_stats)
262 self.assertEqual(subvisit.summary_stats, visit.summary_stats)
263 self.assertIs(subvisit.summary_stats, visit.summary_stats)
264 # Check aperture corrections.
265 self.assertEqual(copy.aperture_corrections.keys(), visit.aperture_corrections.keys())
266 self.assertIsNot(copy.aperture_corrections, visit.aperture_corrections)
267 self.assertEqual(subvisit.aperture_corrections.keys(), visit.aperture_corrections.keys())
268 self.assertIs(subvisit.aperture_corrections, visit.aperture_corrections)
269 # Check backgrounds.
270 self.assertEqual(copy.backgrounds.keys(), visit.backgrounds.keys())
271 self.assertIsNot(copy.backgrounds, visit.backgrounds)
272 self.assertEqual(subvisit.backgrounds.keys(), visit.backgrounds.keys())
273 self.assertIs(subvisit.backgrounds, visit.backgrounds)
274 # Check bounds.
275 self.assertIs(copy.bounds, self.polygon)
276 self.assertEqual(subvisit.bounds, subvisit.bbox) # original polygon wholly encloses subvisit.bbox
278 def test_obs_info(self) -> None:
279 """Check that ObservationInfo has been constructed."""
280 visit = self.visit_image
281 self.assertIsNotNone(visit.obs_info)
282 self.maxDiff = None
283 assert visit.obs_info is not None # for mypy.
284 self.assertEqual(visit.obs_info.instrument, "LSSTCam")
286 def test_summary_stats(self) -> None:
287 """Test the comparisons and attributes of ObservationSummaryStats."""
288 self.assertEqual(self.summary_stats, ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4))
289 self.assertNotEqual(self.summary_stats, ObservationSummaryStats(psfSigma=2.5))
290 self.assertNotEqual(
291 self.summary_stats, ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2))
292 )
294 @unittest.skipUnless(HAVE_H5PY, "h5py is not installed")
295 def test_round_trip_ndf(self) -> None:
296 """NDF round-trip for VisitImage."""
297 with RoundtripNdf(self, self.visit_image, "VisitImage") as roundtrip:
298 assert_visit_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
300 @unittest.skipUnless(HAVE_H5PY, "h5py is not installed")
301 def test_fits_ndf_consistency(self) -> None:
302 """FITS and NDF backends produce equal VisitImages on round-trip."""
303 with RoundtripFits(self, self.visit_image) as fits_rt, RoundtripNdf(self, self.visit_image) as ndf_rt:
304 assert_visit_images_equal(self, self.visit_image, fits_rt.result, expect_view=False)
305 assert_visit_images_equal(self, self.visit_image, ndf_rt.result, expect_view=False)
306 assert_visit_images_equal(self, fits_rt.result, ndf_rt.result, expect_view=False)
308 def test_fits_json_consistency(self) -> None:
309 """FITS and JSON backends produce equal VisitImages on round-trip."""
310 with (
311 RoundtripFits(self, self.visit_image) as fits_rt,
312 RoundtripJson(self, self.visit_image) as json_rt,
313 ):
314 assert_visit_images_equal(self, self.visit_image, fits_rt.result, expect_view=False)
315 assert_visit_images_equal(self, self.visit_image, json_rt.result, expect_view=False)
316 assert_visit_images_equal(self, fits_rt.result, json_rt.result, expect_view=False)
318 def test_read_write(self) -> None:
319 """Test that a visit can round trip through a FITS file."""
320 with RoundtripFits(self, self.visit_image, "VisitImage") as roundtrip:
321 # Check that we're still using the right compression, and that we
322 # wrote WCSs.
323 fits = roundtrip.inspect()
324 self.assertEqual(fits[1].header["ZCMPTYPE"], "GZIP_2")
325 self.assertEqual(fits[1].header["CTYPE1"], "RA---TAN")
326 self.assertEqual(fits[2].header["ZCMPTYPE"], "GZIP_2")
327 self.assertEqual(fits[2].header["CTYPE1"], "RA---TAN")
328 self.assertEqual(fits[3].header["ZCMPTYPE"], "GZIP_2")
329 self.assertEqual(fits[3].header["CTYPE1"], "RA---TAN")
330 # Check a subimage read.
331 subbox = Box.factory[8:13, 9:30]
332 subimage = roundtrip.get(bbox=subbox)
333 assert_masked_images_equal(self, subimage, self.visit_image[subbox], expect_view=False)
335 # Get an explicit masked image to compare with the subimage
336 with self.subTest():
337 subimage_masked = roundtrip.get("masked_image", bbox=subbox)
338 assert_masked_images_equal(self, subimage_masked, subimage, expect_view=False)
340 # Get the same masked image in a multi-component get and ensure
341 # it is the same thing.
342 components = roundtrip.get("components", components=["masked_image", "psf"], bbox=subbox)
343 self.assertEqual(set(components), {"masked_image", "psf"})
344 assert_masked_images_equal(
345 self, components["masked_image"], subimage_masked, expect_view=False
346 )
348 with self.subTest():
349 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
350 with self.subTest():
351 obs_info = roundtrip.get("obs_info")
352 self.assertIsInstance(obs_info, ObservationInfo)
353 self.assertEqual(obs_info, self.visit_image.obs_info)
354 with self.subTest():
355 summary_stats = roundtrip.get("summary_stats")
356 self.assertIsInstance(summary_stats, ObservationSummaryStats)
357 self.assertEqual(summary_stats, self.visit_image.summary_stats)
358 with self.subTest():
359 psf = roundtrip.get("psf")
360 self.assertIsInstance(psf, GaussianPointSpreadFunction)
361 self.assertEqual(psf.kernel_bbox, self.gaussian_psf.kernel_bbox)
362 with self.subTest():
363 backgrounds = roundtrip.get("backgrounds")
364 self.assertIsInstance(backgrounds, BackgroundMap)
365 self.assertEqual(backgrounds.keys(), {"standard"})
366 self.assertIsInstance(backgrounds["standard"].field, ChebyshevField)
367 self.assertEqual(backgrounds.subtracted.name, "standard")
368 self.assertEqual(
369 roundtrip.result.backgrounds.subtracted.description,
370 "Background subtracted from the image.",
371 )
373 with self.subTest(components="components edge cases"):
374 # Test some components get edge cases.
375 components = roundtrip.get("components", components="image")
376 self.assertIsInstance(components["image"], Image)
378 components = roundtrip.get("components")
379 self.assertEqual(
380 set(components),
381 {
382 "image",
383 "variance",
384 "psf",
385 "bbox",
386 "mask",
387 "obs_info",
388 "backgrounds",
389 "detector",
390 "aperture_corrections",
391 "sky_projection",
392 "summary_stats",
393 "photometric_scaling",
394 },
395 )
397 # Butler morphs RuntimeError to ValueError.
398 with self.assertRaises(ValueError):
399 roundtrip.get("components", components=["image", "nonexistent"])
401 with self.assertRaises(ValueError):
402 roundtrip.get("components", components=["image", "components"])
404 with self.assertRaises(ValueError):
405 roundtrip.get("components", components=[])
407 with self.assertRaises(ValueError):
408 # PSF does not know how to use bbox so this fails.
409 roundtrip.get("components", components="psf", bbox=subbox)
411 assert_visit_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
412 # Check that the round-tripped headers are the same (up to card order).
413 self.assertEqual(len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]), 1)
414 self.assertEqual(
415 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
416 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
417 )
418 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
419 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
420 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
421 # Spot-check the concrete background contents (names, field types,
422 # subtracted entry) against the known fixture, so the equality check
423 # above is not vacuously satisfied by empty background maps.
424 self.assertIsInstance(roundtrip.result.backgrounds, BackgroundMap)
425 self.assertEqual(roundtrip.result.backgrounds.keys(), {"standard"})
426 self.assertIsInstance(roundtrip.result.backgrounds["standard"].field, ChebyshevField)
427 self.assertEqual(roundtrip.result.backgrounds.subtracted.name, "standard")
428 self.assertEqual(
429 roundtrip.result.backgrounds.subtracted.description, "Background subtracted from the image."
430 )
432 def _make_sum_background_visit_image(self) -> VisitImage:
433 """Return a VisitImage whose subtracted background is a SumField.
435 Each operand of the SumField calls ``add_array(name="data")`` from
436 the same nested archive, exercising the per-name disambiguation
437 the output archives perform via `_register_name`.
438 """
439 det_frame = self.visit_image.image.sky_projection.pixel_frame
440 bbox = det_frame.bbox
441 bin_y = bbox.y.linspace(6)
442 bin_x = bbox.x.linspace(7)
443 spline_a = SplineField(
444 bbox,
445 self.rng.standard_normal(size=(bin_y.size, bin_x.size)),
446 y=bin_y,
447 x=bin_x,
448 )
449 spline_b = SplineField(
450 bbox,
451 self.rng.standard_normal(size=(bin_y.size, bin_x.size)),
452 y=bin_y,
453 x=bin_x,
454 )
455 sum_field = SumField([spline_a, spline_b])
456 bg_map = BackgroundMap()
457 bg_map.add(
458 "stacked",
459 sum_field,
460 description="Two-operand SumField subtracted background.",
461 is_subtracted=True,
462 )
463 return VisitImage(
464 self.image,
465 variance=self.variance,
466 psf=self.gaussian_psf,
467 mask_schema=self.mask_schema,
468 sky_projection=self.sky_projection,
469 obs_info=self.obs_info,
470 summary_stats=self.summary_stats,
471 detector=self.detector,
472 band="r",
473 backgrounds=bg_map,
474 )
476 def test_sum_background_round_trip_fits(self) -> None:
477 """Two operands of a SumField background each call ``add_array``
478 with the same name; the FITS backend must keep them as distinct
479 EXTVERs rather than overwriting.
480 """
481 visit = self._make_sum_background_visit_image()
482 with RoundtripFits(self, visit) as roundtrip:
483 self._check_sum_background_round_trip(roundtrip.result, visit)
485 @unittest.skipUnless(HAVE_H5PY, "h5py is not installed")
486 def test_sum_background_round_trip_ndf(self) -> None:
487 """NDF must disambiguate the repeated ``data`` leaf the same way."""
488 visit = self._make_sum_background_visit_image()
489 with RoundtripNdf(self, visit) as roundtrip:
490 self._check_sum_background_round_trip(roundtrip.result, visit)
492 def _check_sum_background_round_trip(self, result: VisitImage, original: VisitImage) -> None:
493 subtracted = result.backgrounds.subtracted
494 assert subtracted is not None
495 self.assertIsInstance(subtracted.field, SumField)
496 original_subtracted = original.backgrounds.subtracted
497 assert original_subtracted is not None
498 original_field = original_subtracted.field
499 assert isinstance(original_field, SumField)
500 round_field = subtracted.field
501 assert isinstance(round_field, SumField)
502 self.assertEqual(len(round_field.operands), len(original_field.operands))
503 for round_op, orig_op in zip(round_field.operands, original_field.operands, strict=True):
504 self.assertEqual(round_op, orig_op)
507class VisitImageLegacyTestMixin:
508 """Tests for the VisitImage class and the basics of the archive, to be
509 specialized for a particular test image.
511 `setUp` or `setUpClass` must be implemented to set the attributes declared
512 in the class.
513 """
515 filename: str
516 legacy_exposure: Any
517 plane_map: dict[str, MaskPlane]
518 visit_image: VisitImage
519 unit: u.UnitBase
520 storage_class: ClassVar[str] = "VisitImage"
522 def test_legacy_errors(self) -> None:
523 """Legacy read failure modes."""
524 with self.assertRaises(ValueError):
525 VisitImage.from_legacy(self.legacy_exposure, instrument="HSC")
526 with self.assertRaises(ValueError):
527 VisitImage.from_legacy(self.legacy_exposure, visit=123456)
528 with self.assertRaises(ValueError):
529 VisitImage.from_legacy(self.legacy_exposure, unit=u.mJy)
530 visit = VisitImage.from_legacy(
531 self.legacy_exposure, instrument="LSSTCam", unit=self.unit, visit=2025052000177
532 )
533 self.assertEqual(visit.unit, self.unit)
535 with self.assertRaises(ValueError):
536 VisitImage.read_legacy(self.filename, instrument="HSC")
537 with self.assertRaises(ValueError):
538 VisitImage.read_legacy(self.filename, visit=123456)
540 def test_component_reads(self) -> None:
541 """Test reads of components from legacy file."""
542 visit = VisitImage.read_legacy(self.filename)
543 proj = VisitImage.read_legacy(self.filename, component="sky_projection")
544 assert_sky_projections_equal(self, proj, visit.sky_projection, expect_identity=False)
545 image = VisitImage.read_legacy(self.filename, component="image")
546 self.assertEqual(image, visit.image)
547 assert_sky_projections_equal(self, proj, image.sky_projection, expect_identity=False)
548 variance = VisitImage.read_legacy(self.filename, component="variance")
549 self.assertEqual(variance, visit.variance)
550 assert_sky_projections_equal(self, proj, variance.sky_projection, expect_identity=False)
551 mask = VisitImage.read_legacy(self.filename, component="mask")
552 self.assertEqual(mask, visit.mask)
553 assert_sky_projections_equal(self, proj, mask.sky_projection, expect_identity=False)
554 psf = VisitImage.read_legacy(self.filename, component="psf")
555 self.assertIsInstance(psf, PointSpreadFunction)
556 obs_info = VisitImage.read_legacy(self.filename, component="obs_info")
557 self.check_legacy_obs_info(obs_info)
558 summary_stats = VisitImage.read_legacy(self.filename, component="summary_stats")
559 self.assertIsInstance(summary_stats, ObservationSummaryStats)
560 self.assertEqual(summary_stats.nPsfStar, self.legacy_exposure.info.getSummaryStats().nPsfStar)
561 compare_aperture_corrections_to_legacy(
562 self,
563 VisitImage.read_legacy(self.filename, component="aperture_corrections"),
564 self.legacy_exposure.info.getApCorrMap(),
565 visit.bbox,
566 )
567 detector = VisitImage.read_legacy(self.filename, component="detector")
568 compare_detector_to_legacy(self, detector, self.legacy_exposure.getDetector(), is_raw_assembled=True)
569 photometric_scaling = VisitImage.read_legacy(self.filename, component="photometric_scaling")
570 compare_photo_calib_to_legacy(
571 self,
572 photometric_scaling,
573 self.legacy_exposure.getPhotoCalib(),
574 subimage_bbox=visit.bbox,
575 )
577 def check_legacy_obs_info(self, obs_info: ObservationInfo | None) -> None:
578 """Check that an `ObservationInfo` instance is not `None`, and that it
579 matches the one in the legacy test data file.
580 """
581 self.assertIsInstance(obs_info, ObservationInfo)
582 self.assertEqual(obs_info.instrument, "LSSTCam")
583 self.assertEqual(obs_info.detector_num, 85, obs_info)
584 self.assertEqual(obs_info.detector_unique_name, "R21_S11", obs_info)
585 self.assertEqual(obs_info.physical_filter, "r_57", obs_info)
587 def test_obs_info(self) -> None:
588 """Check that ObservationInfo has been constructed."""
589 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
590 self.assertIsNotNone(legacy.obs_info)
591 self.maxDiff = None
592 self.assertEqual(legacy.obs_info, self.visit_image.obs_info)
593 assert legacy.obs_info is not None # for mypy.
594 self.assertEqual(legacy.obs_info.instrument, "LSSTCam")
595 self.assertEqual(legacy.obs_info.detector_num, 85, legacy.obs_info)
596 self.assertEqual(legacy.obs_info.detector_unique_name, "R21_S11", legacy.obs_info)
597 self.assertEqual(legacy.obs_info.physical_filter, "r_57", legacy.obs_info)
599 def test_aperture_corrections_to_legacy(self) -> None:
600 """Test that we can convert an aperture correction map back to a
601 legacy `lsst.afw.image.ApCorrMap`.
602 """
603 legacy_ap_corr_map = aperture_corrections_to_legacy(self.visit_image.aperture_corrections)
604 compare_aperture_corrections_to_legacy(
605 self, self.visit_image.aperture_corrections, legacy_ap_corr_map, self.visit_image.bbox
606 )
608 def test_read_legacy_headers(self) -> None:
609 """Test that headers were correctly stripped and interpreted in
610 `VisitImage.read_legacy`.
611 """
612 # Check that we read the units from BUNIT.
613 self.assertEqual(self.visit_image.unit, self.unit)
614 # Check that the primary header has the keys we want, and none of the
615 # keys we don't want.
616 header = self.visit_image._opaque_metadata.headers[ExtensionKey()]
617 self.assertIn("EXPTIME", header)
618 self.assertEqual(header["PLATFORM"], "lsstcam")
619 self.assertNotIn("LSST BUTLER ID", header)
620 self.assertNotIn("AR HDU", header)
621 self.assertNotIn("A_ORDER", header)
622 # Check that the extension HDUs do not have any custom headers.
623 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
624 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
625 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
627 def test_from_legacy_headers(self) -> None:
628 """Test that from_legacy handles headers properly."""
629 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
630 header = legacy._opaque_metadata.headers[ExtensionKey()]
631 self.assertIn("EXPTIME", header)
632 self.assertEqual(header["PLATFORM"], "lsstcam")
633 self.assertNotIn("LSST BUTLER ID", header)
634 self.assertNotIn("AR HDU", header)
635 self.assertNotIn("A_ORDER", header)
636 # Check that the extension HDUs do not have any custom headers.
637 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
638 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
639 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
641 def test_rewrite(self) -> None:
642 """Test that we can rewrite the visit image and preserve both
643 lossy-compressed pixel values and components exactly.
644 """
645 import lsst.afw.image
647 with RoundtripFits(self, self.visit_image, self.storage_class) as roundtrip:
648 # Check that we're still using the right compression, and that we
649 # wrote WCSs.
650 fits = roundtrip.inspect()
651 self.assertEqual(fits[1].header["ZCMPTYPE"], "RICE_1")
652 self.assertEqual(fits[1].header["CTYPE1"], "RA---TAN-SIP")
653 self.assertEqual(fits[2].header["ZCMPTYPE"], "GZIP_2")
654 self.assertEqual(fits[2].header["CTYPE1"], "RA---TAN-SIP")
655 self.assertEqual(fits[3].header["ZCMPTYPE"], "RICE_1")
656 self.assertEqual(fits[3].header["CTYPE1"], "RA---TAN-SIP")
657 # Check a subimage read.
658 subbox = Box.factory[8:13, 9:30]
659 subimage = roundtrip.get(bbox=subbox)
660 assert_masked_images_equal(self, subimage, self.visit_image[subbox], expect_view=False)
661 alternates: dict[str, Any] = {}
662 with self.subTest():
663 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
664 alternates = {
665 k: roundtrip.get(k)
666 for k in [
667 "sky_projection",
668 "image",
669 "mask",
670 "variance",
671 "psf",
672 "obs_info",
673 "summary_stats",
674 "aperture_corrections",
675 "detector",
676 "photometric_scaling",
677 ]
678 }
679 # Test reading back in as an Exposure.
680 with self.subTest():
681 legacy_exposure = roundtrip.get(storageClass="Exposure")
682 self.assertIsInstance(legacy_exposure, lsst.afw.image.Exposure)
683 # This covers most of the compnents, which have clean 1-1
684 # mappings from legacy to new:
685 compare_visit_image_to_legacy(
686 self,
687 self.visit_image,
688 legacy_exposure,
689 expect_view=False,
690 plane_map=self.plane_map,
691 **DP2_VISIT_DETECTOR_DATA_ID,
692 )
693 # A few components are different enough to merit extra
694 # attention:
695 if self.visit_image.unit == u.nJy:
696 self.assertTrue(legacy_exposure.getPhotoCalib()._isConstant)
697 self.assertEqual(legacy_exposure.getPhotoCalib().getCalibrationMean(), 1.0)
698 else:
699 compare_photo_calib_to_legacy(
700 self,
701 self.visit_image.photometric_scaling,
702 legacy_exposure.getPhotoCalib(),
703 subimage_bbox=subbox,
704 )
705 self.assertEqual(legacy_exposure.info.getId(), self.legacy_exposure.info.getId())
706 # Try to do a butler get of a component with storage class
707 # override.
708 with self.subTest():
709 # We have VisitInfo available.
710 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
711 self.assertIsInstance(visit_info, lsst.afw.image.VisitInfo)
712 self.assertEqual(visit_info.getInstrumentLabel(), "LSSTCam")
714 assert_visit_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
715 # Check that the round-tripped headers are the same (up to card order).
716 self.assertEqual(
717 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
718 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
719 )
720 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
721 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
722 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
723 self.assertEqual(roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"], "lsstcam")
724 compare_visit_image_to_legacy(
725 self,
726 roundtrip.result,
727 self.legacy_exposure,
728 expect_view=False,
729 plane_map=self.plane_map,
730 **DP2_VISIT_DETECTOR_DATA_ID,
731 alternates=alternates,
732 )
733 # Check converting from the legacy object in-memory.
734 compare_visit_image_to_legacy(
735 self,
736 VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map),
737 self.legacy_exposure,
738 expect_view=True,
739 plane_map=self.plane_map,
740 **DP2_VISIT_DETECTOR_DATA_ID,
741 )
743 def test_butler_converters(self) -> None:
744 """Test that we can read a VisitImage and its components from a butler
745 dataset written as an `lsst.afw.image.Exposure`.
746 """
747 if self.legacy_exposure is None:
748 raise unittest.SkipTest("lsst.afw.image.afw could not be imported.")
749 with TemporaryButler(legacy="ExposureF") as helper:
750 from lsst.daf.butler import FileDataset
752 helper.butler.ingest(FileDataset(path=self.filename, refs=[helper.legacy]), transfer="symlink")
753 visit_image_ref = helper.legacy.overrideStorageClass(self.storage_class)
754 with warnings.catch_warnings():
755 # Silence warnings about data ID and filter label disagreeing.
756 warnings.simplefilter("ignore", category=UserWarning)
757 visit_image = helper.butler.get(visit_image_ref)
758 # We didn't ask for the quantization to be preserved, so it
759 # shouldn't be.
760 self.assertEqual(visit_image._opaque_metadata.precompressed.keys(), set())
761 # This time preserve the quantization
762 visit_image = helper.butler.get(visit_image_ref, parameters={"preserve_quantization": True})
763 self.assertEqual(visit_image._opaque_metadata.precompressed.keys(), {"IMAGE", "VARIANCE"})
764 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
765 self.assertEqual(bbox, visit_image.bbox)
766 alternates = {
767 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
768 # TODO: including "sky_projection" or "obs_info" here fails
769 # because there's code in daf_butler that expects any component
770 # to be valid for the *internal* storage class, not the
771 # requested one, and that's difficult to fix because it's tied
772 # up with the data ID standardization logic.
773 for k in ["image", "mask", "variance", "bbox", "psf", "detector"]
774 }
775 compare_visit_image_to_legacy(
776 self,
777 visit_image,
778 self.legacy_exposure,
779 expect_view=False,
780 plane_map=self.plane_map,
781 alternates=alternates,
782 **DP2_VISIT_DETECTOR_DATA_ID,
783 )
784 # Add some metadata to the new VisitImage and then do a converting
785 # `put` that should write to the old format (we have to delete the
786 # old one first, which just deletes a symlink).
787 helper.butler.pruneDatasets([helper.legacy], purge=True, unstore=True, disassociate=True)
788 visit_image.metadata["MixedCaseKey"] = 52
789 helper.butler.put(visit_image, visit_image_ref)
790 # Check that we can read *that* back in as a legacy exposure.
791 legacy_exposure = helper.butler.get(helper.legacy)
792 compare_visit_image_to_legacy(
793 self,
794 visit_image,
795 legacy_exposure,
796 expect_view=False,
797 plane_map=self.plane_map,
798 alternates=alternates,
799 **DP2_VISIT_DETECTOR_DATA_ID,
800 )
801 # Check that we can read it back in as a VisitImage, and that the
802 # new metadata is preserved.
803 visit_image_2 = helper.butler.get(visit_image_ref)
804 compare_visit_image_to_legacy(
805 self,
806 visit_image_2,
807 legacy_exposure,
808 expect_view=False,
809 plane_map=self.plane_map,
810 alternates=alternates,
811 **DP2_VISIT_DETECTOR_DATA_ID,
812 )
813 self.assertEqual(visit_image_2.metadata["MixedCaseKey"], 52)
816@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
817class VisitImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
818 """Tests for the VisitImage class using a DRP-final visit_image dataset.
820 Requires legacy code.
821 """
823 @classmethod
824 def setUpClass(cls) -> None:
825 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
826 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits")
827 try:
828 from lsst.afw.image import ExposureFitsReader
830 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
831 except ImportError:
832 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
833 cls.plane_map = get_legacy_visit_image_mask_planes()
834 cls.visit_image = VisitImage.read_legacy(
835 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
836 )
837 cls.unit = u.nJy
839 def test_convert_unit(self) -> None:
840 """Test using the ``photometric_scaling`` to swap between
841 calibrated and instrumental units.
843 This includes tests of the `VisitImage.to_legacy` logic for units that
844 don't map directly to `lsst.afw.image.PhotoCalib` conventions.
845 """
846 from lsst.afw.table import ExposureCatalog
848 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
849 # Make a copy of class state so we can modify it without breaking
850 # other tests.
851 original = self.visit_image.copy()
852 # We should not be able to convert to instrumental units when there is
853 # no photometric scaling.
854 with self.assertRaises(u.UnitConversionError):
855 original.convert_unit(u.electron)
856 # Converting to the current unit should be a no-op that does not need
857 # to copy.
858 visit_image_nJy = original.convert_unit(u.nJy, copy=False)
859 self.assertTrue(np.may_share_memory(visit_image_nJy.image.array, original.image.array))
860 self.assertTrue(np.may_share_memory(visit_image_nJy.variance.array, original.variance.array))
861 # Even without a photometric_scaling attached, we should be able to
862 # convert to a compatible unit, but only if we allow copies.
863 with self.assertRaises(u.UnitConversionError):
864 original.convert_unit(u.mJy, copy=False)
865 visit_image_mJy = original.convert_unit(u.mJy, copy="as-needed")
866 self.assertEqual(visit_image_mJy.unit, u.mJy)
867 assert_close(self, visit_image_mJy.image.array, original.image.array * 1e-6)
868 self.assertTrue(np.may_share_memory(visit_image_nJy.mask.array, original.mask.array))
869 assert_close(self, visit_image_mJy.variance.array, original.variance.array * 1e-12)
870 # Converting a mJy image to legacy should make a PhotoCalib that maps
871 # mJy to nJy.
872 legacy_exposure_mJy = visit_image_mJy.to_legacy()
873 assert_close(self, legacy_exposure_mJy.getPhotoCalib().getCalibrationMean(), 1e6)
874 legacy_masked_image_nJy = legacy_exposure_mJy.getPhotoCalib().calibrateImage(
875 legacy_exposure_mJy.maskedImage
876 )
877 assert_close(self, visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
878 assert_close(self, visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
879 # Test that we haven't dropped any component objects along the way,
880 # and that they're all still the same objects or thin views.
881 self.assertTrue(np.may_share_memory(visit_image_mJy.mask.array, original.mask.array))
882 self.assertIs(visit_image_mJy.sky_projection, original.sky_projection)
883 self.assertIs(visit_image_mJy.obs_info, original.obs_info)
884 self.assertIs(visit_image_mJy.summary_stats, original.summary_stats)
885 self.assertIs(visit_image_mJy.psf, original.psf)
886 self.assertIs(visit_image_mJy.detector, original.detector)
887 self.assertIs(visit_image_mJy.bounds, original.bounds)
888 self.assertIs(visit_image_mJy.aperture_corrections, original.aperture_corrections)
889 self.assertIs(visit_image_mJy.photometric_scaling, original.photometric_scaling)
890 # Attach the final PhotoCalib (this isn't stored with the legacy file
891 # because that is the mapping to nJy, which is trivial).
892 visit_summary = ExposureCatalog.readFits(
893 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_summary.fits")
894 )
895 legacy_photo_calib = visit_summary.find(DP2_VISIT_DETECTOR_DATA_ID["detector"]).getPhotoCalib()
896 visit_image_nJy.photometric_scaling = field_from_legacy_photo_calib(
897 legacy_photo_calib, bounds=original.detector.bbox, instrumental_unit=u.electron
898 )
899 compare_photo_calib_to_legacy(
900 self,
901 visit_image_nJy.photometric_scaling,
902 self.legacy_exposure.getPhotoCalib(),
903 applied_legacy_photo_calib=legacy_photo_calib,
904 subimage_bbox=visit_image_nJy.bbox,
905 )
906 # We still can't convert to completely unrelated units.
907 with self.assertRaises(u.UnitConversionError):
908 visit_image_nJy.convert_unit(u.mm)
909 # Uncalibrating via the photometric_scaling matches what legacy code
910 # does, and by default it copies everything.
911 with self.assertRaises(u.UnitConversionError):
912 visit_image_nJy.convert_unit(u.electron, copy=False)
913 legacy_masked_image_e = legacy_photo_calib.uncalibrateImage(self.legacy_exposure.maskedImage)
914 visit_image_e = visit_image_nJy.convert_unit(u.electron)
915 assert_close(self, visit_image_e.image.array, legacy_masked_image_e.image.array)
916 assert_close(self, visit_image_e.variance.array, legacy_masked_image_e.variance.array)
917 self.assertFalse(np.may_share_memory(visit_image_e.mask.array, visit_image_nJy.mask.array))
918 # We can also uncalibrate if we start with an image that has units
919 # that are compatible with the photometric_scaling but not identical
920 # to it.
921 visit_image_mJy.photometric_scaling = visit_image_nJy.photometric_scaling
922 visit_image_e = visit_image_mJy.convert_unit(u.electron)
923 assert_close(self, visit_image_e.image.array, legacy_masked_image_e.image.array)
924 assert_close(self, visit_image_e.variance.array, legacy_masked_image_e.variance.array)
925 # We can re-apply the scaling to go back to calibrated units.
926 visit_image_nJy_2 = visit_image_e.convert_unit(u.nJy)
927 assert_close(self, visit_image_nJy_2.image.array, visit_image_nJy.image.array)
928 assert_close(self, visit_image_nJy_2.variance.array, original.variance.array)
929 # Try calibrating an image with a scaling that has units other than
930 # nJy in the numerator.
931 visit_image_e.photometric_scaling = visit_image_nJy.photometric_scaling * (1e-6 * u.mJy / u.nJy)
932 visit_image_nJy_3 = visit_image_e.convert_unit(u.nJy)
933 assert_close(self, visit_image_nJy_3.image.array, visit_image_nJy.image.array)
934 assert_close(self, visit_image_nJy_3.variance.array, original.variance.array)
935 # Try converting that uncalibrated image to legacy; the extra mJy/nJy
936 # factor should get included in the PhotoCalib to recover the original
937 # PhotoCalib.
938 legacy_exposure_e = visit_image_e.to_legacy()
939 assert_close(
940 self,
941 legacy_exposure_e.getPhotoCalib().getCalibrationMean(),
942 legacy_photo_calib.getCalibrationMean(),
943 )
944 legacy_masked_image_nJy = legacy_exposure_e.getPhotoCalib().calibrateImage(
945 legacy_exposure_e.maskedImage
946 )
947 assert_close(self, visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
948 assert_close(self, visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
951@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
952class PreliminaryVisitImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
953 """Tests for the VisitImage class using a DRP preliminary_visit_image
954 dataset.
956 Requires legacy code.
957 """
959 @classmethod
960 def setUpClass(cls) -> None:
961 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
962 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "preliminary_visit_image.fits")
963 try:
964 from lsst.afw.image import ExposureFitsReader
966 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
967 except ImportError:
968 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
969 cls.plane_map = get_legacy_visit_image_mask_planes()
970 cls.visit_image = VisitImage.read_legacy(
971 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
972 )
973 cls.unit = u.electron
976@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
977class DifferenceImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
978 """Tests for the DifferenceImage class using a DRP difference_image
979 dataset.
981 Because DifferenceImage is a trivial subclass of VisitImage (it may be
982 extended in the future), we run the VisitImage tests to make sure nothing
983 has gone wrong in anything that wasn't trivially inherited.
985 Requires legacy code.
986 """
988 storage_class: ClassVar[str] = "DifferenceImage"
990 @classmethod
991 def setUpClass(cls) -> None:
992 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
993 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_image.fits")
994 try:
995 from lsst.afw.image import ExposureFitsReader
997 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
998 except ImportError:
999 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
1000 cls.plane_map = get_legacy_difference_image_mask_planes()
1001 cls.visit_image = DifferenceImage.read_legacy(
1002 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
1003 )
1004 cls.unit = u.nJy
1007if __name__ == "__main__":
1008 unittest.main()