Coverage for tests/test_visit_image.py: 16%
419 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-09 02:02 -0700
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 ProjectionAstropyView,
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, 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_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_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.projection = make_random_projection(cls.rng, det_frame, Box.factory[1:4096, 1:4096])
105 # API signature suggests projection and obs_info can be None but they
106 # 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 projection=cls.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 projection=cls.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, ProjectionAstropyView)
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 projection=self.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 projection=self.projection,
174 detector=self.detector,
175 band="r",
176 )
178 with self.assertRaises(TypeError):
179 # Requires a 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 projection=self.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 projection=self.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 projection=self.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_projection(self.rng, tract_frame, Box.factory[1:4096, 1:4096])
225 with self.assertRaises(TypeError):
226 VisitImage(
227 self.image,
228 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 projection=self.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 )
381class VisitImageLegacyTestMixin:
382 """Tests for the VisitImage class and the basics of the archive, to be
383 specialized for a particular test image.
385 `setUp` or `setUpClass` must be implemented to set the attributes declared
386 in the class.
387 """
389 filename: str
390 legacy_exposure: Any
391 plane_map: dict[str, MaskPlane]
392 visit_image: VisitImage
393 unit: u.UnitBase
394 storage_class: ClassVar[str] = "VisitImage"
396 def test_legacy_errors(self) -> None:
397 """Legacy read failure modes."""
398 with self.assertRaises(ValueError):
399 VisitImage.from_legacy(self.legacy_exposure, instrument="HSC")
400 with self.assertRaises(ValueError):
401 VisitImage.from_legacy(self.legacy_exposure, visit=123456)
402 with self.assertRaises(ValueError):
403 VisitImage.from_legacy(self.legacy_exposure, unit=u.mJy)
404 visit = VisitImage.from_legacy(
405 self.legacy_exposure, instrument="LSSTCam", unit=self.unit, visit=2025052000177
406 )
407 self.assertEqual(visit.unit, self.unit)
409 with self.assertRaises(ValueError):
410 VisitImage.read_legacy(self.filename, instrument="HSC")
411 with self.assertRaises(ValueError):
412 VisitImage.read_legacy(self.filename, visit=123456)
414 def test_component_reads(self) -> None:
415 """Test reads of components from legacy file."""
416 visit = VisitImage.read_legacy(self.filename)
417 proj = VisitImage.read_legacy(self.filename, component="projection")
418 assert_projections_equal(self, proj, visit.projection, expect_identity=False)
419 image = VisitImage.read_legacy(self.filename, component="image")
420 self.assertEqual(image, visit.image)
421 assert_projections_equal(self, proj, image.projection, expect_identity=False)
422 variance = VisitImage.read_legacy(self.filename, component="variance")
423 self.assertEqual(variance, visit.variance)
424 assert_projections_equal(self, proj, variance.projection, expect_identity=False)
425 mask = VisitImage.read_legacy(self.filename, component="mask")
426 self.assertEqual(mask, visit.mask)
427 assert_projections_equal(self, proj, mask.projection, expect_identity=False)
428 psf = VisitImage.read_legacy(self.filename, component="psf")
429 self.assertIsInstance(psf, PointSpreadFunction)
430 obs_info = VisitImage.read_legacy(self.filename, component="obs_info")
431 self.check_legacy_obs_info(obs_info)
432 summary_stats = VisitImage.read_legacy(self.filename, component="summary_stats")
433 self.assertIsInstance(summary_stats, ObservationSummaryStats)
434 self.assertEqual(summary_stats.nPsfStar, self.legacy_exposure.info.getSummaryStats().nPsfStar)
435 compare_aperture_corrections_to_legacy(
436 self,
437 VisitImage.read_legacy(self.filename, component="aperture_corrections"),
438 self.legacy_exposure.info.getApCorrMap(),
439 visit.bbox,
440 )
441 detector = VisitImage.read_legacy(self.filename, component="detector")
442 compare_detector_to_legacy(self, detector, self.legacy_exposure.getDetector(), is_raw_assembled=True)
443 photometric_scaling = VisitImage.read_legacy(self.filename, component="photometric_scaling")
444 compare_photo_calib_to_legacy(
445 self,
446 photometric_scaling,
447 self.legacy_exposure.getPhotoCalib(),
448 subimage_bbox=visit.bbox,
449 )
451 def check_legacy_obs_info(self, obs_info: ObservationInfo | None) -> None:
452 """Check that an `ObservationInfo` instance is not `None`, and that it
453 matches the one in the legacy test data file.
454 """
455 self.assertIsInstance(obs_info, ObservationInfo)
456 self.assertEqual(obs_info.instrument, "LSSTCam")
457 self.assertEqual(obs_info.detector_num, 85, obs_info)
458 self.assertEqual(obs_info.detector_unique_name, "R21_S11", obs_info)
459 self.assertEqual(obs_info.physical_filter, "r_57", obs_info)
461 def test_obs_info(self) -> None:
462 """Check that ObservationInfo has been constructed."""
463 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
464 self.assertIsNotNone(legacy.obs_info)
465 self.maxDiff = None
466 self.assertEqual(legacy.obs_info, self.visit_image.obs_info)
467 assert legacy.obs_info is not None # for mypy.
468 self.assertEqual(legacy.obs_info.instrument, "LSSTCam")
469 self.assertEqual(legacy.obs_info.detector_num, 85, legacy.obs_info)
470 self.assertEqual(legacy.obs_info.detector_unique_name, "R21_S11", legacy.obs_info)
471 self.assertEqual(legacy.obs_info.physical_filter, "r_57", legacy.obs_info)
473 def test_aperture_corrections_to_legacy(self) -> None:
474 """Test that we can convert an aperture correction map back to a
475 legacy `lsst.afw.image.ApCorrMap`.
476 """
477 legacy_ap_corr_map = aperture_corrections_to_legacy(self.visit_image.aperture_corrections)
478 compare_aperture_corrections_to_legacy(
479 self, self.visit_image.aperture_corrections, legacy_ap_corr_map, self.visit_image.bbox
480 )
482 def test_read_legacy_headers(self) -> None:
483 """Test that headers were correctly stripped and interpreted in
484 `VisitImage.read_legacy`.
485 """
486 # Check that we read the units from BUNIT.
487 self.assertEqual(self.visit_image.unit, self.unit)
488 # Check that the primary header has the keys we want, and none of the
489 # keys we don't want.
490 header = self.visit_image._opaque_metadata.headers[ExtensionKey()]
491 self.assertIn("EXPTIME", header)
492 self.assertEqual(header["PLATFORM"], "lsstcam")
493 self.assertNotIn("LSST BUTLER ID", header)
494 self.assertNotIn("AR HDU", header)
495 self.assertNotIn("A_ORDER", header)
496 # Check that the extension HDUs do not have any custom headers.
497 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
498 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
499 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
501 def test_from_legacy_headers(self) -> None:
502 """Test that from_legacy handles headers properly."""
503 legacy = VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map)
504 header = legacy._opaque_metadata.headers[ExtensionKey()]
505 self.assertIn("EXPTIME", header)
506 self.assertEqual(header["PLATFORM"], "lsstcam")
507 self.assertNotIn("LSST BUTLER ID", header)
508 self.assertNotIn("AR HDU", header)
509 self.assertNotIn("A_ORDER", header)
510 # Check that the extension HDUs do not have any custom headers.
511 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("IMAGE")])
512 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("MASK")])
513 self.assertFalse(self.visit_image._opaque_metadata.headers[ExtensionKey("VARIANCE")])
515 def test_rewrite(self) -> None:
516 """Test that we can rewrite the visit image and preserve both
517 lossy-compressed pixel values and components exactly.
518 """
519 import lsst.afw.image
521 with RoundtripFits(self, self.visit_image, self.storage_class) as roundtrip:
522 # Check that we're still using the right compression, and that we
523 # wrote WCSs.
524 fits = roundtrip.inspect()
525 self.assertEqual(fits[1].header["ZCMPTYPE"], "RICE_1")
526 self.assertEqual(fits[1].header["CTYPE1"], "RA---TAN-SIP")
527 self.assertEqual(fits[2].header["ZCMPTYPE"], "GZIP_2")
528 self.assertEqual(fits[2].header["CTYPE1"], "RA---TAN-SIP")
529 self.assertEqual(fits[3].header["ZCMPTYPE"], "RICE_1")
530 self.assertEqual(fits[3].header["CTYPE1"], "RA---TAN-SIP")
531 # Check a subimage read.
532 subbox = Box.factory[8:13, 9:30]
533 subimage = roundtrip.get(bbox=subbox)
534 assert_masked_images_equal(self, subimage, self.visit_image[subbox], expect_view=False)
535 alternates: dict[str, Any] = {}
536 with self.subTest():
537 self.assertEqual(roundtrip.get("bbox"), self.visit_image.bbox)
538 alternates = {
539 k: roundtrip.get(k)
540 for k in [
541 "projection",
542 "image",
543 "mask",
544 "variance",
545 "psf",
546 "obs_info",
547 "summary_stats",
548 "aperture_corrections",
549 "detector",
550 "photometric_scaling",
551 ]
552 }
553 # Test reading back in as an Exposure.
554 with self.subTest():
555 legacy_exposure = roundtrip.get(storageClass="Exposure")
556 self.assertIsInstance(legacy_exposure, lsst.afw.image.Exposure)
557 # This covers most of the compnents, which have clean 1-1
558 # mappings from legacy to new:
559 compare_visit_image_to_legacy(
560 self,
561 self.visit_image,
562 legacy_exposure,
563 expect_view=False,
564 plane_map=self.plane_map,
565 **DP2_VISIT_DETECTOR_DATA_ID,
566 )
567 # A few components are different enough to merit extra
568 # attention:
569 if self.visit_image.unit == u.nJy:
570 self.assertTrue(legacy_exposure.getPhotoCalib()._isConstant)
571 self.assertEqual(legacy_exposure.getPhotoCalib().getCalibrationMean(), 1.0)
572 else:
573 compare_photo_calib_to_legacy(
574 self,
575 self.visit_image.photometric_scaling,
576 legacy_exposure.getPhotoCalib(),
577 subimage_bbox=subbox,
578 )
579 self.assertEqual(legacy_exposure.info.getId(), self.legacy_exposure.info.getId())
580 # Try to do a butler get of a component with storage class
581 # override.
582 with self.subTest():
583 # We have VisitInfo available.
584 visit_info = roundtrip.get("obs_info", storageClass="VisitInfo")
585 self.assertIsInstance(visit_info, lsst.afw.image.VisitInfo)
586 self.assertEqual(visit_info.getInstrumentLabel(), "LSSTCam")
588 assert_visit_images_equal(self, roundtrip.result, self.visit_image, expect_view=False)
589 # Check that the round-tripped headers are the same (up to card order).
590 self.assertEqual(
591 dict(self.visit_image._opaque_metadata.headers[ExtensionKey()]),
592 dict(roundtrip.result._opaque_metadata.headers[ExtensionKey()]),
593 )
594 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")])
595 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")])
596 self.assertFalse(roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")])
597 self.assertEqual(roundtrip.result._opaque_metadata.headers[ExtensionKey()]["PLATFORM"], "lsstcam")
598 compare_visit_image_to_legacy(
599 self,
600 roundtrip.result,
601 self.legacy_exposure,
602 expect_view=False,
603 plane_map=self.plane_map,
604 **DP2_VISIT_DETECTOR_DATA_ID,
605 alternates=alternates,
606 )
607 # Check converting from the legacy object in-memory.
608 compare_visit_image_to_legacy(
609 self,
610 VisitImage.from_legacy(self.legacy_exposure, plane_map=self.plane_map),
611 self.legacy_exposure,
612 expect_view=True,
613 plane_map=self.plane_map,
614 **DP2_VISIT_DETECTOR_DATA_ID,
615 )
617 def test_butler_converters(self) -> None:
618 """Test that we can read a VisitImage and its components from a butler
619 dataset written as an `lsst.afw.image.Exposure`.
620 """
621 if self.legacy_exposure is None:
622 raise unittest.SkipTest("lsst.afw.image.afw could not be imported.")
623 with TemporaryButler(legacy="ExposureF") as helper:
624 from lsst.daf.butler import FileDataset
626 helper.butler.ingest(FileDataset(path=self.filename, refs=[helper.legacy]), transfer="symlink")
627 visit_image_ref = helper.legacy.overrideStorageClass(self.storage_class)
628 with warnings.catch_warnings():
629 # Silence warnings about data ID and filter label disagreeing.
630 warnings.simplefilter("ignore", category=UserWarning)
631 visit_image = helper.butler.get(visit_image_ref)
632 # We didn't ask for the quantization to be preserved, so it
633 # shouldn't be.
634 self.assertEqual(visit_image._opaque_metadata.precompressed.keys(), set())
635 # This time preserve the quantization
636 visit_image = helper.butler.get(visit_image_ref, parameters={"preserve_quantization": True})
637 self.assertEqual(visit_image._opaque_metadata.precompressed.keys(), {"IMAGE", "VARIANCE"})
638 bbox = helper.butler.get(visit_image_ref.makeComponentRef("bbox"))
639 self.assertEqual(bbox, visit_image.bbox)
640 alternates = {
641 k: helper.butler.get(visit_image_ref.makeComponentRef(k))
642 # TODO: including "projection" or "obs_info" here fails because
643 # there's code in daf_butler that expects any component to be
644 # valid for the *internal* storage class, not the requested
645 # one, and that's difficult to fix because it's tied up with
646 # the data ID standardization logic.
647 for k in ["image", "mask", "variance", "bbox", "psf", "detector"]
648 }
649 compare_visit_image_to_legacy(
650 self,
651 visit_image,
652 self.legacy_exposure,
653 expect_view=False,
654 plane_map=self.plane_map,
655 alternates=alternates,
656 **DP2_VISIT_DETECTOR_DATA_ID,
657 )
658 # Add some metadata to the new VisitImage and then do a converting
659 # `put` that should write to the old format (we have to delete the
660 # old one first, which just deletes a symlink).
661 helper.butler.pruneDatasets([helper.legacy], purge=True, unstore=True, disassociate=True)
662 visit_image.metadata["MixedCaseKey"] = 52
663 helper.butler.put(visit_image, visit_image_ref)
664 # Check that we can read *that* back in as a legacy exposure.
665 legacy_exposure = helper.butler.get(helper.legacy)
666 compare_visit_image_to_legacy(
667 self,
668 visit_image,
669 legacy_exposure,
670 expect_view=False,
671 plane_map=self.plane_map,
672 alternates=alternates,
673 **DP2_VISIT_DETECTOR_DATA_ID,
674 )
675 # Check that we can read it back in as a VisitImage, and that the
676 # new metadata is preserved.
677 visit_image_2 = helper.butler.get(visit_image_ref)
678 compare_visit_image_to_legacy(
679 self,
680 visit_image_2,
681 legacy_exposure,
682 expect_view=False,
683 plane_map=self.plane_map,
684 alternates=alternates,
685 **DP2_VISIT_DETECTOR_DATA_ID,
686 )
687 self.assertEqual(visit_image_2.metadata["MixedCaseKey"], 52)
690@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
691class VisitImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
692 """Tests for the VisitImage class using a DRP-final visit_image dataset.
694 Requires legacy code.
695 """
697 @classmethod
698 def setUpClass(cls) -> None:
699 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
700 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits")
701 try:
702 from lsst.afw.image import ExposureFitsReader
704 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
705 except ImportError:
706 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
707 cls.plane_map = get_legacy_visit_image_mask_planes()
708 cls.visit_image = VisitImage.read_legacy(
709 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
710 )
711 cls.unit = u.nJy
713 def test_convert_unit(self) -> None:
714 """Test using the ``photometric_scaling`` to swap between
715 calibrated and instrumental units.
717 This includes tests of the `VisitImage.to_legacy` logic for units that
718 don't map directly to `lsst.afw.image.PhotoCalib` conventions.
719 """
720 from lsst.afw.table import ExposureCatalog
722 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
723 # Make a copy of class state so we can modify it without breaking
724 # other tests.
725 original = self.visit_image.copy()
726 # We should not be able to convert to instrumental units when there is
727 # no photometric scaling.
728 with self.assertRaises(u.UnitConversionError):
729 original.convert_unit(u.electron)
730 # Converting to the current unit should be a no-op that does not need
731 # to copy.
732 visit_image_nJy = original.convert_unit(u.nJy, copy=False)
733 self.assertTrue(np.may_share_memory(visit_image_nJy.image.array, original.image.array))
734 self.assertTrue(np.may_share_memory(visit_image_nJy.variance.array, original.variance.array))
735 # Even without a photometric_scaling attached, we should be able to
736 # convert to a compatible unit, but only if we allow copies.
737 with self.assertRaises(u.UnitConversionError):
738 original.convert_unit(u.mJy, copy=False)
739 visit_image_mJy = original.convert_unit(u.mJy, copy="as-needed")
740 self.assertEqual(visit_image_mJy.unit, u.mJy)
741 assert_close(self, visit_image_mJy.image.array, original.image.array * 1e-6)
742 self.assertTrue(np.may_share_memory(visit_image_nJy.mask.array, original.mask.array))
743 assert_close(self, visit_image_mJy.variance.array, original.variance.array * 1e-12)
744 # Converting a mJy image to legacy should make a PhotoCalib that maps
745 # mJy to nJy.
746 legacy_exposure_mJy = visit_image_mJy.to_legacy()
747 assert_close(self, legacy_exposure_mJy.getPhotoCalib().getCalibrationMean(), 1e6)
748 legacy_masked_image_nJy = legacy_exposure_mJy.getPhotoCalib().calibrateImage(
749 legacy_exposure_mJy.maskedImage
750 )
751 assert_close(self, visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
752 assert_close(self, visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
753 # Test that we haven't dropped any component objects along the way,
754 # and that they're all still the same objects or thin views.
755 self.assertTrue(np.may_share_memory(visit_image_mJy.mask.array, original.mask.array))
756 self.assertIs(visit_image_mJy.projection, original.projection)
757 self.assertIs(visit_image_mJy.obs_info, original.obs_info)
758 self.assertIs(visit_image_mJy.summary_stats, original.summary_stats)
759 self.assertIs(visit_image_mJy.psf, original.psf)
760 self.assertIs(visit_image_mJy.detector, original.detector)
761 self.assertIs(visit_image_mJy.bounds, original.bounds)
762 self.assertIs(visit_image_mJy.aperture_corrections, original.aperture_corrections)
763 self.assertIs(visit_image_mJy.photometric_scaling, original.photometric_scaling)
764 # Attach the final PhotoCalib (this isn't stored with the legacy file
765 # because that is the mapping to nJy, which is trivial).
766 visit_summary = ExposureCatalog.readFits(
767 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_summary.fits")
768 )
769 legacy_photo_calib = visit_summary.find(DP2_VISIT_DETECTOR_DATA_ID["detector"]).getPhotoCalib()
770 visit_image_nJy.photometric_scaling = field_from_legacy_photo_calib(
771 legacy_photo_calib, bounds=original.detector.bbox, instrumental_unit=u.electron
772 )
773 compare_photo_calib_to_legacy(
774 self,
775 visit_image_nJy.photometric_scaling,
776 self.legacy_exposure.getPhotoCalib(),
777 applied_legacy_photo_calib=legacy_photo_calib,
778 subimage_bbox=visit_image_nJy.bbox,
779 )
780 # We still can't convert to completely unrelated units.
781 with self.assertRaises(u.UnitConversionError):
782 visit_image_nJy.convert_unit(u.mm)
783 # Uncalibrating via the photometric_scaling matches what legacy code
784 # does, and by default it copies everything.
785 with self.assertRaises(u.UnitConversionError):
786 visit_image_nJy.convert_unit(u.electron, copy=False)
787 legacy_masked_image_e = legacy_photo_calib.uncalibrateImage(self.legacy_exposure.maskedImage)
788 visit_image_e = visit_image_nJy.convert_unit(u.electron)
789 assert_close(self, visit_image_e.image.array, legacy_masked_image_e.image.array)
790 assert_close(self, visit_image_e.variance.array, legacy_masked_image_e.variance.array)
791 self.assertFalse(np.may_share_memory(visit_image_e.mask.array, visit_image_nJy.mask.array))
792 # We can also uncalibrate if we start with an image that has units
793 # that are compatible with the photometric_scaling but not identical
794 # to it.
795 visit_image_mJy.photometric_scaling = visit_image_nJy.photometric_scaling
796 visit_image_e = visit_image_mJy.convert_unit(u.electron)
797 assert_close(self, visit_image_e.image.array, legacy_masked_image_e.image.array)
798 assert_close(self, visit_image_e.variance.array, legacy_masked_image_e.variance.array)
799 # We can re-apply the scaling to go back to calibrated units.
800 visit_image_nJy_2 = visit_image_e.convert_unit(u.nJy)
801 assert_close(self, visit_image_nJy_2.image.array, visit_image_nJy.image.array)
802 assert_close(self, visit_image_nJy_2.variance.array, original.variance.array)
803 # Try calibrating an image with a scaling that has units other than
804 # nJy in the numerator.
805 visit_image_e.photometric_scaling = visit_image_nJy.photometric_scaling * (1e-6 * u.mJy / u.nJy)
806 visit_image_nJy_3 = visit_image_e.convert_unit(u.nJy)
807 assert_close(self, visit_image_nJy_3.image.array, visit_image_nJy.image.array)
808 assert_close(self, visit_image_nJy_3.variance.array, original.variance.array)
809 # Try converting that uncalibrated image to legacy; the extra mJy/nJy
810 # factor should get included in the PhotoCalib to recover the original
811 # PhotoCalib.
812 legacy_exposure_e = visit_image_e.to_legacy()
813 assert_close(
814 self,
815 legacy_exposure_e.getPhotoCalib().getCalibrationMean(),
816 legacy_photo_calib.getCalibrationMean(),
817 )
818 legacy_masked_image_nJy = legacy_exposure_e.getPhotoCalib().calibrateImage(
819 legacy_exposure_e.maskedImage
820 )
821 assert_close(self, visit_image_nJy.image.array, legacy_masked_image_nJy.image.array)
822 assert_close(self, visit_image_nJy.variance.array, legacy_masked_image_nJy.variance.array)
825@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
826class PreliminaryVisitImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
827 """Tests for the VisitImage class using a DRP preliminary_visit_image
828 dataset.
830 Requires legacy code.
831 """
833 @classmethod
834 def setUpClass(cls) -> None:
835 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
836 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "preliminary_visit_image.fits")
837 try:
838 from lsst.afw.image import ExposureFitsReader
840 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
841 except ImportError:
842 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
843 cls.plane_map = get_legacy_visit_image_mask_planes()
844 cls.visit_image = VisitImage.read_legacy(
845 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
846 )
847 cls.unit = u.electron
850@unittest.skipUnless(EXTERNAL_DATA_DIR is not None, "TESTDATA_IMAGES_DIR is not in the environment.")
851class DifferenceImageLegacyTestCase(unittest.TestCase, VisitImageLegacyTestMixin):
852 """Tests for the DifferenceImage class using a DRP difference_image
853 dataset.
855 Because DifferenceImage is a trivial subclass of VisitImage (it may be
856 extended in the future), we run the VisitImage tests to make sure nothing
857 has gone wrong in anything that wasn't trivially inherited.
859 Requires legacy code.
860 """
862 storage_class: ClassVar[str] = "DifferenceImage"
864 @classmethod
865 def setUpClass(cls) -> None:
866 assert EXTERNAL_DATA_DIR is not None, "Guaranteed by decorator."
867 cls.filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "difference_image.fits")
868 try:
869 from lsst.afw.image import ExposureFitsReader
871 cls.legacy_exposure = ExposureFitsReader(cls.filename).read()
872 except ImportError:
873 raise unittest.SkipTest("afw not available; cannot read legacy visit images") from None
874 cls.plane_map = get_legacy_difference_image_mask_planes()
875 cls.visit_image = DifferenceImage.read_legacy(
876 cls.filename, preserve_quantization=True, plane_map=cls.plane_map
877 )
878 cls.unit = u.nJy
881if __name__ == "__main__":
882 unittest.main()