Coverage for python / lsst / images / tests / _checks.py: 15%
179 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 02:28 -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
14__all__ = (
15 "arrays_to_legacy_points",
16 "assert_close",
17 "assert_images_equal",
18 "assert_masked_images_equal",
19 "assert_masks_equal",
20 "assert_projections_equal",
21 "check_astropy_wcs_interface",
22 "check_projection",
23 "check_transform",
24 "compare_image_to_legacy",
25 "compare_mask_to_legacy",
26 "compare_masked_image_to_legacy",
27 "compare_projection_to_legacy_wcs",
28 "compare_psf_to_legacy",
29 "compare_visit_image_to_legacy",
30 "legacy_coords_to_astropy",
31 "legacy_points_to_xy_array",
32)
34import unittest
35from collections.abc import Mapping
36from typing import Any, Literal, cast
38import astropy.units as u
39import astropy.wcs.wcsapi
40import numpy as np
41from astropy.coordinates import SkyCoord
43from .._geom import XY, Box
44from .._image import Image
45from .._mask import Mask, MaskPlane
46from .._masked_image import MaskedImage
47from .._transforms import DetectorFrame, Frame, Projection, SkyFrame, Transform
48from .._visit_image import VisitImage
49from ..psfs import PointSpreadFunction
52def assert_close(
53 tc: unittest.TestCase,
54 a: np.ndarray | u.Quantity | float,
55 b: np.ndarray | u.Quantity | float,
56 **kwargs: Any,
57) -> None:
58 """Test that two arrays, floats, or quantities are almost equal.
60 Parameters
61 ----------
62 tc
63 Test case object with assert methods to use.
64 a
65 Array, float, or quantity to compare.
66 b
67 Array, float, or quantity to compare.
68 **kwargs
69 Forwarded to `astropy.units.allclose`.
70 """
71 tc.assertTrue(u.allclose(a, b, **kwargs), msg=f"{a} != {b}")
74def assert_images_equal(
75 tc: unittest.TestCase,
76 a: Image,
77 b: Image,
78 *,
79 rtol: float = 0.0,
80 atol: float = 0.0,
81 expect_view: bool | Literal["array"] | None = None,
82) -> None:
83 """Assert that two images are equal or nearly equal."""
84 tc.assertEqual(a.bbox, b.bbox)
85 tc.assertEqual(a.unit, b.unit)
86 assert_projections_equal(tc, a.projection, b.projection)
87 if expect_view is not None:
88 tc.assertEqual(np.may_share_memory(a.array, b.array), bool(expect_view))
89 if expect_view == "array":
90 tc.assertEqual(a.metadata, b.metadata)
91 else:
92 tc.assertEqual(a.metadata is b.metadata, expect_view)
93 if not expect_view:
94 assert_close(tc, a.array, b.array, atol=atol, rtol=rtol)
95 tc.assertEqual(a.metadata, b.metadata)
98def assert_masks_equal(tc: unittest.TestCase, a: Mask, b: Mask) -> None:
99 """Assert that two masks are equal or nearly equal."""
100 tc.assertEqual(a.bbox, b.bbox)
101 tc.assertEqual(a.schema, b.schema)
102 tc.assertEqual(a.metadata, b.metadata)
103 assert_projections_equal(tc, a.projection, b.projection)
104 np.testing.assert_array_equal(a.array, b.array)
107def assert_masked_images_equal(
108 tc: unittest.TestCase,
109 a: MaskedImage,
110 b: MaskedImage,
111 *,
112 rtol: float = 0.0,
113 atol: float = 0.0,
114 expect_view: bool | None = None,
115) -> None:
116 """Assert that two masked images are equal or nearly equal."""
117 tc.assertEqual(a.metadata, b.metadata)
118 assert_projections_equal(tc, a.projection, b.projection)
119 assert_images_equal(tc, a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view)
120 assert_masks_equal(tc, a.mask, b.mask)
121 assert_images_equal(tc, a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view)
124def compare_image_to_legacy(
125 tc: unittest.TestCase, image: Image, legacy_image: Any, expect_view: bool | None = None
126) -> None:
127 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object."""
128 tc.assertEqual(image.bbox, Box.from_legacy(legacy_image.getBBox()))
129 if expect_view is not None:
130 tc.assertEqual(np.may_share_memory(image.array, legacy_image.array), expect_view)
131 if not expect_view:
132 np.testing.assert_array_equal(image.array, legacy_image.array)
135def compare_mask_to_legacy(
136 tc: unittest.TestCase, mask: Mask, legacy_mask: Any, plane_map: Mapping[str, MaskPlane] | None = None
137) -> None:
138 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object."""
139 tc.assertEqual(mask.bbox, Box.from_legacy(legacy_mask.getBBox()))
140 if plane_map is None:
141 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
142 for old_name, new_plane in plane_map.items():
143 np.testing.assert_array_equal(
144 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
145 mask.get(new_plane.name),
146 )
149def compare_masked_image_to_legacy(
150 tc: unittest.TestCase,
151 masked_image: MaskedImage,
152 legacy_masked_image: Any,
153 *,
154 plane_map: Mapping[str, MaskPlane] | None = None,
155 expect_view: bool | None = None,
156 alternates: Mapping[str, Any] | None = None,
157) -> None:
158 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
159 object.
161 Parameters
162 ----------
163 tc
164 Test case to use for asserts.
165 masked_image
166 New image to test.
167 legacy_masked_image
168 Legacy image to test against.
169 plane_map
170 Mapping between new and legacy mask planes.
171 expect_view
172 Whether to test that the image and variance arrays do or do not share
173 memory.
174 alternates
175 A mapping of other versions of one or more (new) components to also
176 check against the legacy versions of those components.
177 """
178 compare_image_to_legacy(tc, masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view)
179 compare_mask_to_legacy(tc, masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map)
180 compare_image_to_legacy(
181 tc, masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view
182 )
183 if alternates:
184 if image := alternates.get("image"):
185 compare_image_to_legacy(tc, image, legacy_masked_image.getImage(), expect_view=expect_view)
186 if mask := alternates.get("mask"):
187 compare_mask_to_legacy(tc, mask, legacy_masked_image.getMask(), plane_map=plane_map)
188 if variance := alternates.get("variance"):
189 compare_image_to_legacy(tc, variance, legacy_masked_image.getVariance(), expect_view=expect_view)
192def compare_visit_image_to_legacy(
193 tc: unittest.TestCase,
194 visit_image: VisitImage,
195 legacy_exposure: Any,
196 *,
197 plane_map: Mapping[str, MaskPlane] | None = None,
198 expect_view: bool | None = None,
199 instrument: str,
200 visit: int,
201 detector: int,
202 alternates: Mapping[str, Any] | None = None,
203) -> None:
204 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
205 object.
207 Parameters
208 ----------
209 tc
210 Test case to use for asserts.
211 visit_image
212 New image to test.
213 legacy_exposure
214 Legacy image to test against.
215 plane_map
216 Mapping between new and legacy mask planes.
217 expect_view
218 Whether to test that the image and variance arrays do or do not share
219 memory.
220 instrument
221 Expected instrument name.
222 visit
223 Expected visit ID.
224 detetector
225 Expected detector ID.
226 alternates
227 A mapping of other versions of one or more (new) components to also
228 check against the legacy versions of those components.
229 """
230 compare_masked_image_to_legacy(
231 tc,
232 visit_image,
233 legacy_exposure.getMaskedImage(),
234 plane_map=plane_map,
235 expect_view=expect_view,
236 alternates=alternates,
237 )
238 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
239 compare_projection_to_legacy_wcs(
240 tc,
241 visit_image.projection,
242 legacy_exposure.getWcs(),
243 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
244 visit_image.bbox,
245 )
246 tc.assertIs(visit_image.projection, visit_image.mask.projection)
247 tc.assertIs(visit_image.projection, visit_image.variance.projection)
248 compare_psf_to_legacy(tc, visit_image.psf, legacy_exposure.getPsf())
249 if alternates:
250 if projection := alternates.get("projection"):
251 compare_projection_to_legacy_wcs(
252 tc,
253 projection,
254 legacy_exposure.getWcs(),
255 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
256 visit_image.bbox,
257 )
258 if psf := alternates.get("psf"):
259 compare_psf_to_legacy(tc, psf, legacy_exposure.getPsf())
262def compare_psf_to_legacy(tc: unittest.TestCase, psf: PointSpreadFunction, legacy_psf: Any) -> None:
263 """Compare a PSF model object to its legacy interface.
265 Parameters
266 ----------
267 tc
268 Test case object with assert methods to use.
269 """
270 from lsst.geom import Point2D
272 for p in [Point2D(50.0, 60.0), Point2D(801.2, 322.8), Point2D(33.5, 22.1)]:
273 tc.assertEqual(psf.kernel_bbox, Box.from_legacy(legacy_psf.computeKernelBBox(p)))
274 tc.assertEqual(
275 psf.compute_kernel_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeKernelImage(p))
276 )
277 tc.assertEqual(
278 psf.compute_stellar_bbox(x=p.x, y=p.y), Box.from_legacy(legacy_psf.computeImageBBox(p))
279 )
280 tc.assertEqual(psf.compute_stellar_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeImage(p)))
283def compare_projection_to_legacy_wcs[F: Frame](
284 tc: unittest.TestCase,
285 projection: Projection[F],
286 legacy_wcs: Any,
287 pixel_frame: F,
288 subimage_bbox: Box,
289 is_fits: bool = False,
290) -> None:
291 """Test a Projection object by comparing it to an equivalent
292 `lsst.afw.geom.SkyWcs`.
294 Parameters
295 ----------
296 tc
297 Test case object with assert methods to use.
298 projection
299 Projection to test.
300 legacy_wcs : ``lsst.afw.geom.SkyWcs``
301 Equivalent legacy WCS.
302 pixel_frame
303 Expected pixel frame for the projection.
304 subimage_bbox
305 Bounding box of points to generate for tests.
306 is_fits
307 Whether this projection is expected to be exactly representable as a
308 FIT WCS. If `False` it is assumed to have a FITS approximation
309 attached instead.
310 """
311 # Pixel coordinates to test on over the subimage region of interest:
312 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
313 # Array indices of those pixel values (subtract off bbox starts):
314 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
315 sky_coords = legacy_coords_to_astropy(
316 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
317 )
318 # Test transforming with the Projection itself, which also tests its
319 # nested Transform and an Astropy High-Level WCS view with no origin
320 # change.
321 check_projection(tc, projection, pixel_xy, sky_coords, pixel_frame)
322 # Also test the Astropy High-Level WCS view with an origin change to
323 # array indices.
324 check_astropy_wcs_interface(
325 tc, projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
326 )
327 if is_fits:
328 fits_wcs = projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
329 assert fits_wcs is not None
330 check_astropy_wcs_interface(tc, fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
331 # Use that FITS approximation to check that we can make a
332 # Projection from a FITS WCS, too.
333 fits_projection = Projection.from_fits_wcs(fits_wcs, pixel_frame)
334 check_projection(
335 tc,
336 fits_projection,
337 subimage_array_xy,
338 sky_coords,
339 pixel_frame,
340 pixel_atol=1e-5,
341 )
342 else:
343 tc.assertIsNone(projection.as_fits_wcs(subimage_bbox, allow_approximation=False))
344 # The legacy SkyWcs should instead have a FITS approximation
345 # attached; run the same tests on that.
346 assert projection.fits_approximation is not None
347 compare_projection_to_legacy_wcs(
348 tc,
349 projection.fits_approximation,
350 legacy_wcs.getFitsApproximation(),
351 pixel_frame,
352 subimage_bbox,
353 is_fits=True,
354 )
357def check_transform[I: Frame, O: Frame](
358 tc: unittest.TestCase,
359 transform: Transform[I, O],
360 input_xy: XY[np.ndarray],
361 output_xy: XY[np.ndarray],
362 in_frame: Frame,
363 out_frame: Frame,
364 *,
365 check_inverted: bool = True,
366 in_atol: u.Quantity | None = None,
367 out_atol: u.Quantity | None = None,
368) -> None:
369 """Test Transform against known arrays of input and output points.
371 Parameters
372 ----------
373 tc
374 Test case object with assert methods to use.
375 transform
376 Transform to test.
377 input_xy
378 Arrays of input points.
379 outout_xy
380 Arrays of output points.
381 in_frame
382 Expected input frame.
383 out_frame
384 Expected output frame.
385 check_inverted
386 If `True`, recurse (once) to test the inverse transform.
387 in_atol
388 Expected absolute precision of input points.
389 out_atol
390 Expected absolute precision of output points.
391 """
392 tc.assertEqual(transform.in_frame, in_frame)
393 tc.assertEqual(transform.out_frame, out_frame)
394 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None
395 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None
396 # Test array interfaces.
397 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
398 assert_close(tc, test_output_xy.x, output_xy.x, atol=out_atol_v)
399 assert_close(tc, test_output_xy.y, output_xy.y, atol=out_atol_v)
400 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
401 assert_close(tc, test_input_xy.x, input_xy.x, atol=in_atol_v)
402 assert_close(tc, test_input_xy.y, input_xy.y, atol=in_atol_v)
403 # Test scalar interfaces with numpy scalars.
404 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
405 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v)
406 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v)
407 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v)
408 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v)
409 # Test quantity array interfaces.
410 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
411 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
412 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
413 assert_close(tc, test_output_q_xy.x, output_q_xy.x, atol=out_atol)
414 assert_close(tc, test_output_q_xy.y, output_q_xy.y, atol=out_atol)
415 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
416 assert_close(tc, test_input_q_xy.x, input_q_xy.x, atol=in_atol)
417 assert_close(tc, test_input_q_xy.y, input_q_xy.y, atol=in_atol)
418 # Test quantity scalar interfaces.
419 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
420 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
421 ):
422 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol)
423 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol)
424 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol)
425 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol)
426 if check_inverted:
427 # Test the inverse transform.
428 check_transform(
429 tc,
430 transform.inverted(),
431 output_xy,
432 input_xy,
433 out_frame,
434 in_frame,
435 check_inverted=False,
436 out_atol=in_atol,
437 in_atol=out_atol,
438 )
441def check_projection[P: Frame](
442 tc: unittest.TestCase,
443 projection: Projection[P],
444 pixel_xy: XY[np.ndarray],
445 sky_coords: SkyCoord,
446 pixel_frame: Frame,
447 *,
448 pixel_atol: float | None = None,
449 sky_atol: u.Quantity | None = None,
450) -> None:
451 """Test a `.Projection` instance against known arrays of pixel and sky
452 coordinates.
454 Parameters
455 ----------
456 tc
457 Test case object with assert methods to use.
458 projection
459 Projection to test.
460 pixel_xy
461 Arrays of pixel coordinates.
462 sky_coords
463 Corresponding sky coordinates.
464 pixel_frame
465 Expected pixel frame.
466 pixel_atol
467 Expected absolute precision of pixel points.
468 sky_atol
469 Expected absolute precision of sky coordinates.
470 """
471 tc.assertEqual(projection.pixel_frame, pixel_frame)
472 tc.assertEqual(projection.sky_frame, SkyFrame.ICRS)
473 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None
474 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None
475 # Test array interfaces.
476 test_pixel_xy = cast(XY[np.ndarray], projection.sky_to_pixel(sky_coords))
477 assert_close(tc, test_pixel_xy.x, pixel_xy.x, atol=pixel_atol)
478 assert_close(tc, test_pixel_xy.y, pixel_xy.y, atol=pixel_atol)
479 test_sky_astropy = projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
480 assert_close(tc, test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v)
481 assert_close(tc, test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v)
482 # Test scalar interfaces.
483 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
484 assert_close(tc, projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol)
485 assert_close(tc, projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol)
486 assert_close(tc, projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v)
487 assert_close(tc, projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v)
488 # Test the underlying Transform object.
489 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
490 check_transform(
491 tc,
492 projection.pixel_to_sky_transform,
493 pixel_xy,
494 sky_xy,
495 pixel_frame,
496 SkyFrame.ICRS,
497 check_inverted=False,
498 in_atol=pixel_atol_q,
499 out_atol=sky_atol,
500 )
501 check_transform(
502 tc,
503 projection.sky_to_pixel_transform,
504 sky_xy,
505 pixel_xy,
506 SkyFrame.ICRS,
507 pixel_frame,
508 check_inverted=False,
509 in_atol=sky_atol,
510 out_atol=pixel_atol_q,
511 )
512 # Test the Astropy interface adapter.
513 check_astropy_wcs_interface(
514 tc, projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
515 )
518def assert_projections_equal(
519 tc: unittest.TestCase,
520 a: Projection[Any] | None,
521 b: Projection[Any] | None,
522 expect_identity: bool | None = None,
523) -> None:
524 """Test that two `.Projection` instances are equivalent."""
525 if a is None and b is None:
526 return
527 assert a is not None and b is not None
528 match expect_identity:
529 case True:
530 tc.assertIs(a, b)
531 return
532 case False:
533 tc.assertIsNot(a, b)
534 case None if a is b:
535 return
536 tc.assertEqual(a.pixel_frame, b.pixel_frame)
537 tc.assertEqual(a.show(simplified=True), b.show(simplified=True))
538 assert_projections_equal(
539 tc, a.fits_approximation, cast(Projection[Any], b.fits_approximation), expect_identity=False
540 )
543def check_astropy_wcs_interface(
544 tc: unittest.TestCase,
545 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
546 pixel_xy: XY[np.ndarray],
547 sky_coords: SkyCoord,
548 *,
549 pixel_atol: float | None = None,
550 sky_atol: u.Quantity | None = None,
551) -> None:
552 """Test an Astropy WCS instance against known arrays of pixel and
553 sky coordinates.
555 Parameters
556 ----------
557 tc
558 Test case object with assert methods to use.
559 wcs
560 Astropy WCS object to test.
561 pixel_xy
562 Arrays of pixel coordinates.
563 sky_coords
564 Corresponding sky coordinates.
565 pixel_atol
566 Expected absolute precision of pixel points.
567 sky_atol
568 Expected absolute precision of sky coordinates.
569 """
570 test_x, test_y = wcs.world_to_pixel(sky_coords)
571 assert_close(tc, test_x, pixel_xy.x, atol=pixel_atol)
572 assert_close(tc, test_y, pixel_xy.y, atol=pixel_atol)
573 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
574 assert_close(tc, test_sky_coords.ra, sky_coords.ra, atol=sky_atol)
575 assert_close(tc, test_sky_coords.dec, sky_coords.dec, atol=sky_atol)
578def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
579 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array."""
580 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
583def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
584 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
585 coordinate object.
586 """
587 return SkyCoord(
588 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
589 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
590 )
593def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
594 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``."""
595 from lsst.geom import Point2D
597 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]