Coverage for tests/test_visit_image.py: 50%
452 statements
« prev ^ index » next coverage.py v7.14.2, created at 2026-06-22 08:46 +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):
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):
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):
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)
334 with self.subTest():
335 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
336 with self.subTest():
337 obs_info = roundtrip.get("obs_info")
338 self.assertIsInstance(obs_info, ObservationInfo)
339 self.assertEqual(obs_info, self.visit_image.obs_info)
340 with self.subTest():
341 summary_stats = roundtrip.get("summary_stats")
342 self.assertIsInstance(summary_stats, ObservationSummaryStats)
343 self.assertEqual(summary_stats, self.visit_image.summary_stats)
344 with self.subTest():
345 psf = roundtrip.get("psf")
346 self.assertIsInstance(psf, GaussianPointSpreadFunction)
347 self.assertEqual(psf.kernel_bbox, self.gaussian_psf.kernel_bbox)
348 with self.subTest():
349 backgrounds = roundtrip.get("backgrounds")
350 self.assertIsInstance(backgrounds, BackgroundMap)
351 self.assertEqual(backgrounds.keys(), {"standard"})
352 self.assertIsInstance(backgrounds["standard"].field, ChebyshevField)
353 self.assertEqual(backgrounds.subtracted.name, "standard")
354 self.assertEqual(
355 roundtrip.result.backgrounds.subtracted.description,
356 "Background subtracted from the image.",
357 )
359 assert_visit_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
360 # Check that the round-tripped headers are the same (up to card order).
361 self.assertEqual(len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]), 1)
362 self.assertEqual(
363 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
364 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
365 )
366 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
367 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
368 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
369 # Spot-check the concrete background contents (names, field types,
370 # subtracted entry) against the known fixture, so the equality check
371 # above is not vacuously satisfied by empty background maps.
372 self.assertIsInstance(roundtrip.result.backgrounds, BackgroundMap)
373 self.assertEqual(roundtrip.result.backgrounds.keys(), {"standard"})
374 self.assertIsInstance(roundtrip.result.backgrounds["standard"].field, ChebyshevField)
375 self.assertEqual(roundtrip.result.backgrounds.subtracted.name, "standard")
376 self.assertEqual(
377 roundtrip.result.backgrounds.subtracted.description, "Background subtracted from the image."
378 )
380 def _make_sum_background_visit_image(self) -> VisitImage:
381 """Return a VisitImage whose subtracted background is a SumField.
383 Each operand of the SumField calls ``add_array(name="data")`` from
384 the same nested archive, exercising the per-name disambiguation
385 the output archives perform via `_register_name`.
386 """
387 det_frame = self.visit_image.image.sky_projection.pixel_frame
388 bbox = det_frame.bbox
389 bin_y = bbox.y.linspace(6)
390 bin_x = bbox.x.linspace(7)
391 spline_a = SplineField(
392 bbox,
393 self.rng.standard_normal(size=(bin_y.size, bin_x.size)),
394 y=bin_y,
395 x=bin_x,
396 )
397 spline_b = SplineField(
398 bbox,
399 self.rng.standard_normal(size=(bin_y.size, bin_x.size)),
400 y=bin_y,
401 x=bin_x,
402 )
403 sum_field = SumField([spline_a, spline_b])
404 bg_map = BackgroundMap()
405 bg_map.add(
406 "stacked",
407 sum_field,
408 description="Two-operand SumField subtracted background.",
409 is_subtracted=True,
410 )
411 return VisitImage(
412 self.image,
413 variance=self.variance,
414 psf=self.gaussian_psf,
415 mask_schema=self.mask_schema,
416 sky_projection=self.sky_projection,
417 obs_info=self.obs_info,
418 summary_stats=self.summary_stats,
419 detector=self.detector,
420 band="r",
421 backgrounds=bg_map,
422 )
424 def test_sum_background_round_trip_fits(self) -> None:
425 """Two operands of a SumField background each call ``add_array``
426 with the same name; the FITS backend must keep them as distinct
427 EXTVERs rather than overwriting.
428 """
429 visit = self._make_sum_background_visit_image()
430 with RoundtripFits(self, visit) as roundtrip:
431 self._check_sum_background_round_trip(roundtrip.result, visit)
433 @unittest.skipUnless(HAVE_H5PY, "h5py is not installed")
434 def test_sum_background_round_trip_ndf(self) -> None:
435 """NDF must disambiguate the repeated ``data`` leaf the same way."""
436 visit = self._make_sum_background_visit_image()
437 with RoundtripNdf(self, visit) as roundtrip:
438 self._check_sum_background_round_trip(roundtrip.result, visit)
440 def _check_sum_background_round_trip(self, result: VisitImage, original: VisitImage) -> None:
441 subtracted = result.backgrounds.subtracted
442 assert subtracted is not None
443 self.assertIsInstance(subtracted.field, SumField)
444 original_subtracted = original.backgrounds.subtracted
445 assert original_subtracted is not None
446 original_field = original_subtracted.field
447 assert isinstance(original_field, SumField)
448 round_field = subtracted.field
449 assert isinstance(round_field, SumField)
450 self.assertEqual(len(round_field.operands), len(original_field.operands))
451 for round_op, orig_op in zip(round_field.operands, original_field.operands, strict=True):
452 self.assertEqual(round_op, orig_op)
455class VisitImageLegacyTestMixin:
456 """Tests for the VisitImage class and the basics of the archive, to be
457 specialized for a particular test image.
459 `setUp` or `setUpClass` must be implemented to set the attributes declared
460 in the class.
461 """
463 filename: str
464 legacy_exposure: Any
465 plane_map: dict[str, MaskPlane]
466 visit_image: VisitImage
467 unit: u.UnitBase
468 storage_class: ClassVar[str] = "VisitImage"
470 def test_legacy_errors(self) -> None:
471 """Legacy read failure modes."""
472 with self.assertRaises(ValueError):
473 VisitImage.from_legacy(self.legacy_exposure, instrument="HSC")
474 with self.assertRaises(ValueError):
475 VisitImage.from_legacy(self.legacy_exposure, visit=123456)
476 with self.assertRaises(ValueError):
477 VisitImage.from_legacy(self.legacy_exposure, unit=u.mJy)
478 visit = VisitImage.from_legacy(
479 self.legacy_exposure, instrument="LSSTCam", unit=self.unit, visit=2025052000177
480 )
481 self.assertEqual(visit.unit, self.unit)
483 with self.assertRaises(ValueError):
484 VisitImage.read_legacy(self.filename, instrument="HSC")
485 with self.assertRaises(ValueError):
486 VisitImage.read_legacy(self.filename, visit=123456)
488 def test_component_reads(self) -> None:
489 """Test reads of components from legacy file."""
490 visit = VisitImage.read_legacy(self.filename)
491 proj = VisitImage.read_legacy(self.filename, component="sky_projection")
492 assert_sky_projections_equal(self, proj, visit.sky_projection, expect_identity=False)
493 image = VisitImage.read_legacy(self.filename, component="image")
494 self.assertEqual(image, visit.image)
495 assert_sky_projections_equal(self, proj, image.sky_projection, expect_identity=False)
496 variance = VisitImage.read_legacy(self.filename, component="variance")
497 self.assertEqual(variance, visit.variance)
498 assert_sky_projections_equal(self, proj, variance.sky_projection, expect_identity=False)
499 mask = VisitImage.read_legacy(self.filename, component="mask")
500 self.assertEqual(mask, visit.mask)
501 assert_sky_projections_equal(self, proj, mask.sky_projection, expect_identity=False)
502 psf = VisitImage.read_legacy(self.filename, component="psf")
503 self.assertIsInstance(psf, PointSpreadFunction)
504 obs_info = VisitImage.read_legacy(self.filename, component="obs_info")
505 self.check_legacy_obs_info(obs_info)
506 summary_stats = VisitImage.read_legacy(self.filename, component="summary_stats")
507 self.assertIsInstance(summary_stats, ObservationSummaryStats)
508 self.assertEqual(summary_stats.nPsfStar, self.legacy_exposure.info.getSummaryStats().nPsfStar)
509 compare_aperture_corrections_to_legacy(
510 self,
511 VisitImage.read_legacy(self.filename, component="aperture_corrections"),
512 self.legacy_exposure.info.getApCorrMap(),
513 visit.bbox,
514 )
515 detector = VisitImage.read_legacy(self.filename, component="detector")
516 compare_detector_to_legacy(self, detector, self.legacy_exposure.getDetector(), is_raw_assembled=True)
517 photometric_scaling = VisitImage.read_legacy(self.filename, component="photometric_scaling")
518 compare_photo_calib_to_legacy(
519 self,
520 photometric_scaling,
521 self.legacy_exposure.getPhotoCalib(),
522 subimage_bbox=visit.bbox,
523 )
525 def check_legacy_obs_info(self, obs_info: ObservationInfo | None) -> None:
526 """Check that an `ObservationInfo` instance is not `None`, and that it
527 matches the one in the legacy test data file.
528 """
529 self.assertIsInstance(obs_info, ObservationInfo)
530 self.assertEqual(obs_info.instrument, "LSSTCam")
531 self.assertEqual(obs_info.detector_num, 85, obs_info)
532 self.assertEqual(obs_info.detector_unique_name, "R21_S11", obs_info)
533 self.assertEqual(obs_info.physical_filter, "r_57", obs_info)
535 def test_obs_info(self) -> None:
536 """Check that ObservationInfo has been constructed."""
537 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
538 self.assertIsNotNone(legacy.obs_info)
539 self.maxDiff = None
540 self.assertEqual(legacy.obs_info, self.visit_image.obs_info)
541 assert legacy.obs_info is not None # for mypy.
542 self.assertEqual(legacy.obs_info.instrument, "LSSTCam")
543 self.assertEqual(legacy.obs_info.detector_num, 85, legacy.obs_info)
544 self.assertEqual(legacy.obs_info.detector_unique_name, "R21_S11", legacy.obs_info)
545 self.assertEqual(legacy.obs_info.physical_filter, "r_57", legacy.obs_info)
547 def test_aperture_corrections_to_legacy(self) -> None:
548 """Test that we can convert an aperture correction map back to a
549 legacy `lsst.afw.image.ApCorrMap`.
550 """
551 legacy_ap_corr_map = aperture_corrections_to_legacy(self.visit_image.aperture_corrections)
552 compare_aperture_corrections_to_legacy(
553 self, self.visit_image.aperture_corrections, legacy_ap_corr_map, self.visit_image.bbox
554 )
556 def test_read_legacy_headers(self) -> None:
557 """Test that headers were correctly stripped and interpreted in
558 `VisitImage.read_legacy`.
559 """
560 # Check that we read the units from BUNIT.
561 self.assertEqual(self.visit_image.unit, self.unit)
562 # Check that the primary header has the keys we want, and none of the
563 # keys we don't want.
564 header = self.visit_image._opaque_metadata.headers[ExtensionKey()]
565 self.assertIn("EXPTIME", header)
566 self.assertEqual(header["PLATFORM"], "lsstcam")
567 self.assertNotIn("LSST BUTLER ID", header)
568 self.assertNotIn("AR HDU", header)
569 self.assertNotIn("A_ORDER", header)
570 # Check that the extension HDUs do not have any custom headers.
571 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
572 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
573 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
575 def test_from_legacy_headers(self) -> None:
576 """Test that from_legacy handles headers properly."""
577 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
578 header = legacy._opaque_metadata.headers[ExtensionKey()]
579 self.assertIn("EXPTIME", header)
580 self.assertEqual(header["PLATFORM"], "lsstcam")
581 self.assertNotIn("LSST BUTLER ID", header)
582 self.assertNotIn("AR HDU", header)
583 self.assertNotIn("A_ORDER", header)
584 # Check that the extension HDUs do not have any custom headers.
585 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
586 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
587 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
589 def test_rewrite(self) -> None:
590 """Test that we can rewrite the visit image and preserve both
591 lossy-compressed pixel values and components exactly.
592 """
593 import lsst.afw.image
595 with RoundtripFits(self, self.visit_image, self.storage_class) as roundtrip:
596 # Check that we're still using the right compression, and that we
597 # wrote WCSs.
598 fits = roundtrip.inspect()
599 self.assertEqual(fits[1].header["ZCMPTYPE"], "RICE_1")
600 self.assertEqual(fits[1].header["CTYPE1"], "RA---TAN-SIP")
601 self.assertEqual(fits[2].header["ZCMPTYPE"], "GZIP_2")
602 self.assertEqual(fits[2].header["CTYPE1"], "RA---TAN-SIP")
603 self.assertEqual(fits[3].header["ZCMPTYPE"], "RICE_1")
604 self.assertEqual(fits[3].header["CTYPE1"], "RA---TAN-SIP")
605 # Check a subimage read.
606 subbox = Box.factory[8:13, 9:30]
607 subimage = roundtrip.get(bbox=subbox)
608 assert_masked_images_equal(self, subimage, self.visit_image[subbox], expect_view=False)
609 alternates: dict[str, Any] = {}
610 with self.subTest():
611 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
612 alternates = {
613 k: roundtrip.get(k)
614 for k in [
615 "sky_projection",
616 "image",
617 "mask",
618 "variance",
619 "psf",
620 "obs_info",
621 "summary_stats",
622 "aperture_corrections",
623 "detector",
624 "photometric_scaling",
625 ]
626 }
627 # Test reading back in as an Exposure.
628 with self.subTest():
629 legacy_exposure = roundtrip.get(storageClass="Exposure")
630 self.assertIsInstance(legacy_exposure, lsst.afw.image.Exposure)
631 # This covers most of the compnents, which have clean 1-1
632 # mappings from legacy to new:
633 compare_visit_image_to_legacy(
634 self,
635 self.visit_image,
636 legacy_exposure,
637 expect_view=False,
638 plane_map=self.plane_map,
639 **DP2_VISIT_DETECTOR_DATA_ID,
640 )
641 # A few components are different enough to merit extra
642 # attention:
643 if self.visit_image.unit == u.nJy:
644 self.assertTrue(legacy_exposure.getPhotoCalib()._isConstant)
645 self.assertEqual(legacy_exposure.getPhotoCalib().getCalibrationMean(), 1.0)
646 else:
647 compare_photo_calib_to_legacy(
648 self,
649 self.visit_image.photometric_scaling,
650 legacy_exposure.getPhotoCalib(),
651 subimage_bbox=subbox,
652 )
653 self.assertEqual(legacy_exposure.info.getId(), self.legacy_exposure.info.getId())
654 # Try to do a butler get of a component with storage class
655 # override.
656 with self.subTest():
657 # We have VisitInfo available.
658 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
659 self.assertIsInstance(visit_info, lsst.afw.image.VisitInfo)
660 self.assertEqual(visit_info.getInstrumentLabel(), "LSSTCam")
662 assert_visit_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
663 # Check that the round-tripped headers are the same (up to card order).
664 self.assertEqual(
665 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
666 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
667 )
668 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
669 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
670 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
671 self.assertEqual(roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"], "lsstcam")
672 compare_visit_image_to_legacy(
673 self,
674 roundtrip.result,
675 self.legacy_exposure,
676 expect_view=False,
677 plane_map=self.plane_map,
678 **DP2_VISIT_DETECTOR_DATA_ID,
679 alternates=alternates,
680 )
681 # Check converting from the legacy object in-memory.
682 compare_visit_image_to_legacy(
683 self,
684 VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map),
685 self.legacy_exposure,
686 expect_view=True,
687 plane_map=self.plane_map,
688 **DP2_VISIT_DETECTOR_DATA_ID,
689 )
691 def test_butler_converters(self) -> None:
692 """Test that we can read a VisitImage and its components from a butler
693 dataset written as an `lsst.afw.image.Exposure`.
694 """
695 if self.legacy_exposure is None:
696 raise unittest.SkipTest("lsst.afw.image.afw could not be imported.")
697 with TemporaryButler(legacy="ExposureF") as helper:
698 from lsst.daf.butler import FileDataset
700 helper.butler.ingest(FileDataset(path=self.filename, refs=[helper.legacy]), transfer="symlink")
701 visit_image_ref = helper.legacy.overrideStorageClass(self.storage_class)
702 with warnings.catch_warnings():
703 # Silence warnings about data ID and filter label disagreeing.
704 warnings.simplefilter("ignore", category=UserWarning)
705 visit_image = helper.butler.get(visit_image_ref)
706 # We didn't ask for the quantization to be preserved, so it
707 # shouldn't be.
708 self.assertEqual(visit_image._opaque_metadata.precompressed.keys(), set())
709 # This time preserve the quantization
710 visit_image = helper.butler.get(visit_image_ref, parameters={"preserve_quantization": True})
711 self.assertEqual(visit_image._opaque_metadata.precompressed.keys(), {"IMAGE", "VARIANCE"})
712 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
713 self.assertEqual(bbox, visit_image.bbox)
714 alternates = {
715 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
716 # TODO: including "sky_projection" or "obs_info" here fails
717 # because there's code in daf_butler that expects any component
718 # to be valid for the *internal* storage class, not the
719 # requested one, and that's difficult to fix because it's tied
720 # up with the data ID standardization logic.
721 for k in ["image", "mask", "variance", "bbox", "psf", "detector"]
722 }
723 compare_visit_image_to_legacy(
724 self,
725 visit_image,
726 self.legacy_exposure,
727 expect_view=False,
728 plane_map=self.plane_map,
729 alternates=alternates,
730 **DP2_VISIT_DETECTOR_DATA_ID,
731 )
732 # Add some metadata to the new VisitImage and then do a converting
733 # `put` that should write to the old format (we have to delete the
734 # old one first, which just deletes a symlink).
735 helper.butler.pruneDatasets([helper.legacy], purge=True, unstore=True, disassociate=True)
736 visit_image.metadata["MixedCaseKey"] = 52
737 helper.butler.put(visit_image, visit_image_ref)
738 # Check that we can read *that* back in as a legacy exposure.
739 legacy_exposure = helper.butler.get(helper.legacy)
740 compare_visit_image_to_legacy(
741 self,
742 visit_image,
743 legacy_exposure,
744 expect_view=False,
745 plane_map=self.plane_map,
746 alternates=alternates,
747 **DP2_VISIT_DETECTOR_DATA_ID,
748 )
749 # Check that we can read it back in as a VisitImage, and that the
750 # new metadata is preserved.
751 visit_image_2 = helper.butler.get(visit_image_ref)
752 compare_visit_image_to_legacy(
753 self,
754 visit_image_2,
755 legacy_exposure,
756 expect_view=False,
757 plane_map=self.plane_map,
758 alternates=alternates,
759 **DP2_VISIT_DETECTOR_DATA_ID,
760 )
761 self.assertEqual(visit_image_2.metadata["MixedCaseKey"], 52)
764@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
765class VisitImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
766 """Tests for the VisitImage class using a DRP-final visit_image dataset.
768 Requires legacy code.
769 """
771 @classmethod
772 def setUpClass(cls) -> None:
773 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
774 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits")
775 try:
776 from lsst.afw.image import ExposureFitsReader
778 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
779 except ImportError:
780 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
781 cls.plane_map = get_legacy_visit_image_mask_planes()
782 cls.visit_image = VisitImage.read_legacy(
783 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
784 )
785 cls.unit = u.nJy
787 def test_convert_unit(self) -> None:
788 """Test using the ``photometric_scaling`` to swap between
789 calibrated and instrumental units.
791 This includes tests of the `VisitImage.to_legacy` logic for units that
792 don't map directly to `lsst.afw.image.PhotoCalib` conventions.
793 """
794 from lsst.afw.table import ExposureCatalog
796 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
797 # Make a copy of class state so we can modify it without breaking
798 # other tests.
799 original = self.visit_image.copy()
800 # We should not be able to convert to instrumental units when there is
801 # no photometric scaling.
802 with self.assertRaises(u.UnitConversionError):
803 original.convert_unit(u.electron)
804 # Converting to the current unit should be a no-op that does not need
805 # to copy.
806 visit_image_nJy = original.convert_unit(u.nJy, copy=False)
807 self.assertTrue(np.may_share_memory(visit_image_nJy.image.array, original.image.array))
808 self.assertTrue(np.may_share_memory(visit_image_nJy.variance.array, original.variance.array))
809 # Even without a photometric_scaling attached, we should be able to
810 # convert to a compatible unit, but only if we allow copies.
811 with self.assertRaises(u.UnitConversionError):
812 original.convert_unit(u.mJy, copy=False)
813 visit_image_mJy = original.convert_unit(u.mJy, copy="as-needed")
814 self.assertEqual(visit_image_mJy.unit, u.mJy)
815 assert_close(self, visit_image_mJy.image.array, original.image.array * 1e-6)
816 self.assertTrue(np.may_share_memory(visit_image_nJy.mask.array, original.mask.array))
817 assert_close(self, visit_image_mJy.variance.array, original.variance.array * 1e-12)
818 # Converting a mJy image to legacy should make a PhotoCalib that maps
819 # mJy to nJy.
820 legacy_exposure_mJy = visit_image_mJy.to_legacy()
821 assert_close(self, legacy_exposure_mJy.getPhotoCalib().getCalibrationMean(), 1e6)
822 legacy_masked_image_nJy = legacy_exposure_mJy.getPhotoCalib().calibrateImage(
823 legacy_exposure_mJy.maskedImage
824 )
825 assert_close(self, visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
826 assert_close(self, visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
827 # Test that we haven't dropped any component objects along the way,
828 # and that they're all still the same objects or thin views.
829 self.assertTrue(np.may_share_memory(visit_image_mJy.mask.array, original.mask.array))
830 self.assertIs(visit_image_mJy.sky_projection, original.sky_projection)
831 self.assertIs(visit_image_mJy.obs_info, original.obs_info)
832 self.assertIs(visit_image_mJy.summary_stats, original.summary_stats)
833 self.assertIs(visit_image_mJy.psf, original.psf)
834 self.assertIs(visit_image_mJy.detector, original.detector)
835 self.assertIs(visit_image_mJy.bounds, original.bounds)
836 self.assertIs(visit_image_mJy.aperture_corrections, original.aperture_corrections)
837 self.assertIs(visit_image_mJy.photometric_scaling, original.photometric_scaling)
838 # Attach the final PhotoCalib (this isn't stored with the legacy file
839 # because that is the mapping to nJy, which is trivial).
840 visit_summary = ExposureCatalog.readFits(
841 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_summary.fits")
842 )
843 legacy_photo_calib = visit_summary.find(DP2_VISIT_DETECTOR_DATA_ID["detector"]).getPhotoCalib()
844 visit_image_nJy.photometric_scaling = field_from_legacy_photo_calib(
845 legacy_photo_calib, bounds=original.detector.bbox, instrumental_unit=u.electron
846 )
847 compare_photo_calib_to_legacy(
848 self,
849 visit_image_nJy.photometric_scaling,
850 self.legacy_exposure.getPhotoCalib(),
851 applied_legacy_photo_calib=legacy_photo_calib,
852 subimage_bbox=visit_image_nJy.bbox,
853 )
854 # We still can't convert to completely unrelated units.
855 with self.assertRaises(u.UnitConversionError):
856 visit_image_nJy.convert_unit(u.mm)
857 # Uncalibrating via the photometric_scaling matches what legacy code
858 # does, and by default it copies everything.
859 with self.assertRaises(u.UnitConversionError):
860 visit_image_nJy.convert_unit(u.electron, copy=False)
861 legacy_masked_image_e = legacy_photo_calib.uncalibrateImage(self.legacy_exposure.maskedImage)
862 visit_image_e = visit_image_nJy.convert_unit(u.electron)
863 assert_close(self, visit_image_e.image.array, legacy_masked_image_e.image.array)
864 assert_close(self, visit_image_e.variance.array, legacy_masked_image_e.variance.array)
865 self.assertFalse(np.may_share_memory(visit_image_e.mask.array, visit_image_nJy.mask.array))
866 # We can also uncalibrate if we start with an image that has units
867 # that are compatible with the photometric_scaling but not identical
868 # to it.
869 visit_image_mJy.photometric_scaling = visit_image_nJy.photometric_scaling
870 visit_image_e = visit_image_mJy.convert_unit(u.electron)
871 assert_close(self, visit_image_e.image.array, legacy_masked_image_e.image.array)
872 assert_close(self, visit_image_e.variance.array, legacy_masked_image_e.variance.array)
873 # We can re-apply the scaling to go back to calibrated units.
874 visit_image_nJy_2 = visit_image_e.convert_unit(u.nJy)
875 assert_close(self, visit_image_nJy_2.image.array, visit_image_nJy.image.array)
876 assert_close(self, visit_image_nJy_2.variance.array, original.variance.array)
877 # Try calibrating an image with a scaling that has units other than
878 # nJy in the numerator.
879 visit_image_e.photometric_scaling = visit_image_nJy.photometric_scaling * (1e-6 * u.mJy / u.nJy)
880 visit_image_nJy_3 = visit_image_e.convert_unit(u.nJy)
881 assert_close(self, visit_image_nJy_3.image.array, visit_image_nJy.image.array)
882 assert_close(self, visit_image_nJy_3.variance.array, original.variance.array)
883 # Try converting that uncalibrated image to legacy; the extra mJy/nJy
884 # factor should get included in the PhotoCalib to recover the original
885 # PhotoCalib.
886 legacy_exposure_e = visit_image_e.to_legacy()
887 assert_close(
888 self,
889 legacy_exposure_e.getPhotoCalib().getCalibrationMean(),
890 legacy_photo_calib.getCalibrationMean(),
891 )
892 legacy_masked_image_nJy = legacy_exposure_e.getPhotoCalib().calibrateImage(
893 legacy_exposure_e.maskedImage
894 )
895 assert_close(self, visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
896 assert_close(self, visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
899@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
900class PreliminaryVisitImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
901 """Tests for the VisitImage class using a DRP preliminary_visit_image
902 dataset.
904 Requires legacy code.
905 """
907 @classmethod
908 def setUpClass(cls) -> None:
909 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
910 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "preliminary_visit_image.fits")
911 try:
912 from lsst.afw.image import ExposureFitsReader
914 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
915 except ImportError:
916 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
917 cls.plane_map = get_legacy_visit_image_mask_planes()
918 cls.visit_image = VisitImage.read_legacy(
919 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
920 )
921 cls.unit = u.electron
924@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
925class DifferenceImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
926 """Tests for the DifferenceImage class using a DRP difference_image
927 dataset.
929 Because DifferenceImage is a trivial subclass of VisitImage (it may be
930 extended in the future), we run the VisitImage tests to make sure nothing
931 has gone wrong in anything that wasn't trivially inherited.
933 Requires legacy code.
934 """
936 storage_class: ClassVar[str] = "DifferenceImage"
938 @classmethod
939 def setUpClass(cls) -> None:
940 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
941 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_image.fits")
942 try:
943 from lsst.afw.image import ExposureFitsReader
945 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
946 except ImportError:
947 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
948 cls.plane_map = get_legacy_difference_image_mask_planes()
949 cls.visit_image = DifferenceImage.read_legacy(
950 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
951 )
952 cls.unit = u.nJy
955if __name__ == "__main__":
956 unittest.main()