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