Coverage for python / lsst / images / tests / _checks.py: 14%
223 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:07 +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
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_aperture_corrections_to_legacy",
25 "compare_field_to_legacy",
26 "compare_image_to_legacy",
27 "compare_mask_to_legacy",
28 "compare_masked_image_to_legacy",
29 "compare_observation_summary_stats_to_legacy",
30 "compare_projection_to_legacy_wcs",
31 "compare_psf_to_legacy",
32 "compare_visit_image_to_legacy",
33 "legacy_coords_to_astropy",
34 "legacy_points_to_xy_array",
35)
37import dataclasses
38import math
39import unittest
40from collections.abc import Mapping
41from typing import Any, Literal, cast
43import astropy.units as u
44import astropy.wcs.wcsapi
45import numpy as np
46from astropy.coordinates import SkyCoord
48from .._geom import XY, Box
49from .._image import Image
50from .._mask import Mask, MaskPlane
51from .._masked_image import MaskedImage
52from .._observation_summary_stats import ObservationSummaryStats
53from .._transforms import DetectorFrame, Frame, Projection, SkyFrame, Transform
54from .._visit_image import VisitImage
55from ..aperture_corrections import ApertureCorrectionMap
56from ..fields import BaseField
57from ..psfs import PointSpreadFunction
60def assert_close(
61 tc: unittest.TestCase,
62 a: np.ndarray | u.Quantity | float,
63 b: np.ndarray | u.Quantity | float,
64 **kwargs: Any,
65) -> None:
66 """Test that two arrays, floats, or quantities are almost equal.
68 Parameters
69 ----------
70 tc
71 Test case object with assert methods to use.
72 a
73 Array, float, or quantity to compare.
74 b
75 Array, float, or quantity to compare.
76 **kwargs
77 Forwarded to `astropy.units.allclose`.
78 """
79 tc.assertTrue(u.allclose(a, b, **kwargs), msg=f"{a} != {b}")
82def assert_images_equal(
83 tc: unittest.TestCase,
84 a: Image,
85 b: Image,
86 *,
87 rtol: float = 0.0,
88 atol: float = 0.0,
89 expect_view: bool | Literal["array"] | None = None,
90) -> None:
91 """Assert that two images are equal or nearly equal."""
92 tc.assertEqual(a.bbox, b.bbox)
93 tc.assertEqual(a.unit, b.unit)
94 assert_projections_equal(tc, a.projection, b.projection)
95 if expect_view is not None:
96 tc.assertEqual(np.may_share_memory(a.array, b.array), bool(expect_view))
97 if expect_view == "array":
98 tc.assertEqual(a.metadata, b.metadata)
99 else:
100 tc.assertEqual(a.metadata is b.metadata, expect_view)
101 if not expect_view:
102 assert_close(tc, a.array, b.array, atol=atol, rtol=rtol)
103 tc.assertEqual(a.metadata, b.metadata)
106def assert_masks_equal(tc: unittest.TestCase, a: Mask, b: Mask) -> None:
107 """Assert that two masks are equal or nearly equal."""
108 tc.assertEqual(a.bbox, b.bbox)
109 tc.assertEqual(a.schema, b.schema)
110 tc.assertEqual(a.metadata, b.metadata)
111 assert_projections_equal(tc, a.projection, b.projection)
112 np.testing.assert_array_equal(a.array, b.array)
115def assert_masked_images_equal(
116 tc: unittest.TestCase,
117 a: MaskedImage,
118 b: MaskedImage,
119 *,
120 rtol: float = 0.0,
121 atol: float = 0.0,
122 expect_view: bool | None = None,
123) -> None:
124 """Assert that two masked images are equal or nearly equal."""
125 tc.assertEqual(a.metadata, b.metadata)
126 assert_projections_equal(tc, a.projection, b.projection)
127 assert_images_equal(tc, a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view)
128 assert_masks_equal(tc, a.mask, b.mask)
129 assert_images_equal(tc, a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view)
132def compare_image_to_legacy(
133 tc: unittest.TestCase, image: Image, legacy_image: Any, expect_view: bool | None = None
134) -> None:
135 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object."""
136 tc.assertEqual(image.bbox, Box.from_legacy(legacy_image.getBBox()))
137 if expect_view is not None:
138 tc.assertEqual(np.may_share_memory(image.array, legacy_image.array), expect_view)
139 if not expect_view:
140 np.testing.assert_array_equal(image.array, legacy_image.array)
143def compare_mask_to_legacy(
144 tc: unittest.TestCase, mask: Mask, legacy_mask: Any, plane_map: Mapping[str, MaskPlane] | None = None
145) -> None:
146 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object."""
147 tc.assertEqual(mask.bbox, Box.from_legacy(legacy_mask.getBBox()))
148 if plane_map is None:
149 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
150 for old_name, new_plane in plane_map.items():
151 np.testing.assert_array_equal(
152 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
153 mask.get(new_plane.name),
154 )
157def compare_masked_image_to_legacy(
158 tc: unittest.TestCase,
159 masked_image: MaskedImage,
160 legacy_masked_image: Any,
161 *,
162 plane_map: Mapping[str, MaskPlane] | None = None,
163 expect_view: bool | None = None,
164 alternates: Mapping[str, Any] | None = None,
165) -> None:
166 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
167 object.
169 Parameters
170 ----------
171 tc
172 Test case to use for asserts.
173 masked_image
174 New image to test.
175 legacy_masked_image
176 Legacy image to test against.
177 plane_map
178 Mapping between new and legacy mask planes.
179 expect_view
180 Whether to test that the image and variance arrays do or do not share
181 memory.
182 alternates
183 A mapping of other versions of one or more (new) components to also
184 check against the legacy versions of those components.
185 """
186 compare_image_to_legacy(tc, masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view)
187 compare_mask_to_legacy(tc, masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map)
188 compare_image_to_legacy(
189 tc, masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view
190 )
191 if alternates:
192 if image := alternates.get("image"):
193 compare_image_to_legacy(tc, image, legacy_masked_image.getImage(), expect_view=expect_view)
194 if mask := alternates.get("mask"):
195 compare_mask_to_legacy(tc, mask, legacy_masked_image.getMask(), plane_map=plane_map)
196 if variance := alternates.get("variance"):
197 compare_image_to_legacy(tc, variance, legacy_masked_image.getVariance(), expect_view=expect_view)
200def compare_visit_image_to_legacy(
201 tc: unittest.TestCase,
202 visit_image: VisitImage,
203 legacy_exposure: Any,
204 *,
205 plane_map: Mapping[str, MaskPlane] | None = None,
206 expect_view: bool | None = None,
207 instrument: str,
208 visit: int,
209 detector: int,
210 alternates: Mapping[str, Any] | None = None,
211) -> None:
212 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
213 object.
215 Parameters
216 ----------
217 tc
218 Test case to use for asserts.
219 visit_image
220 New image to test.
221 legacy_exposure
222 Legacy image to test against.
223 plane_map
224 Mapping between new and legacy mask planes.
225 expect_view
226 Whether to test that the image and variance arrays do or do not share
227 memory.
228 instrument
229 Expected instrument name.
230 visit
231 Expected visit ID.
232 detector
233 Expected detector ID.
234 alternates
235 A mapping of other versions of one or more (new) components to also
236 check against the legacy versions of those components.
237 """
238 compare_masked_image_to_legacy(
239 tc,
240 visit_image,
241 legacy_exposure.getMaskedImage(),
242 plane_map=plane_map,
243 expect_view=expect_view,
244 alternates=alternates,
245 )
246 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
247 compare_projection_to_legacy_wcs(
248 tc,
249 visit_image.projection,
250 legacy_exposure.getWcs(),
251 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
252 visit_image.bbox,
253 )
254 tc.assertIs(visit_image.projection, visit_image.mask.projection)
255 tc.assertIs(visit_image.projection, visit_image.variance.projection)
256 compare_psf_to_legacy(tc, visit_image.psf, legacy_exposure.getPsf())
257 compare_observation_summary_stats_to_legacy(
258 tc, visit_image.summary_stats, legacy_exposure.info.getSummaryStats()
259 )
260 compare_aperture_corrections_to_legacy(
261 tc, visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), detector_bbox
262 )
263 if alternates:
264 if projection := alternates.get("projection"):
265 compare_projection_to_legacy_wcs(
266 tc,
267 projection,
268 legacy_exposure.getWcs(),
269 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
270 visit_image.bbox,
271 )
272 if psf := alternates.get("psf"):
273 compare_psf_to_legacy(tc, psf, legacy_exposure.getPsf())
274 if summary_stats := alternates.get("summary_stats"):
275 compare_observation_summary_stats_to_legacy(
276 tc, summary_stats, legacy_exposure.info.getSummaryStats()
277 )
278 if obs_info := alternates.get("obs_info"):
279 visitInfo = legacy_exposure.visitInfo
280 tc.assertEqual(obs_info.instrument, visitInfo.getInstrumentLabel())
281 if aperture_corrections := alternates.get("aperture_corrections"):
282 compare_aperture_corrections_to_legacy(
283 tc, aperture_corrections, legacy_exposure.info.getApCorrMap(), detector_bbox
284 )
287def compare_psf_to_legacy(
288 tc: unittest.TestCase, psf: PointSpreadFunction, legacy_psf: Any, subimage_bbox: Box | None = None
289) -> int:
290 """Compare a PSF model object to its legacy interface.
292 Parameters
293 ----------
294 tc
295 Test case object with assert methods to use.
296 psf
297 Point-spread function to test.
298 legacy_psf
299 Legacy `lsst.afw.detection.Psf` instance to compare with.
300 subimage_bbox
301 Bounding box to draw test points from. Defaults to ``psf.bounds``.
303 Returns
304 -------
305 `int`
306 The number of points actually tested.
307 """
308 if subimage_bbox is None:
309 if isinstance(psf.bounds, Box):
310 subimage_bbox = psf.bounds
311 else:
312 subimage_bbox = psf.bounds.bbox
314 # Pixel coordinates to test on over the subimage region of interest:
315 pixel_xy = subimage_bbox.meshgrid(n=3).map(np.ravel)
316 legacy_points = arrays_to_legacy_points(pixel_xy.x, pixel_xy.y)
318 n_points_tested: int = 0
319 for p in legacy_points:
320 if not psf.bounds.contains(x=p.x, y=p.y):
321 continue
322 tc.assertEqual(psf.kernel_bbox, Box.from_legacy(legacy_psf.computeKernelBBox(p)))
323 tc.assertEqual(
324 psf.compute_kernel_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeKernelImage(p))
325 )
326 tc.assertEqual(
327 psf.compute_stellar_bbox(x=p.x, y=p.y), Box.from_legacy(legacy_psf.computeImageBBox(p))
328 )
329 tc.assertEqual(psf.compute_stellar_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeImage(p)))
330 n_points_tested += 1
331 return n_points_tested
334def compare_field_to_legacy(
335 tc: unittest.TestCase,
336 field: BaseField,
337 legacy_field: Any,
338 subimage_bbox: Box,
339) -> None:
340 """Test a Field object by comparing it to an equivalent
341 `lsst.afw.math.BoundedField`.
343 Parameters
344 ----------
345 tc
346 Test case object with assert methods to use.
347 field
348 Field to test.
349 legacy_field : ``lsst.afw.math.BoundedField``
350 Equivalent legacy bounded field.
351 subimage_bbox
352 Bounding box for full-image tests.
353 """
354 tc.assertEqual(field.bounds.bbox, Box.from_legacy(legacy_field.getBBox()))
355 # Pixel coordinates to test the numpy array interface with.
356 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
357 assert_close(tc, field(x=pixel_xy.x, y=pixel_xy.y), legacy_field.evaluate(pixel_xy.x, pixel_xy.y))
358 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
359 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
360 assert_images_equal(tc, field.render(subimage_bbox), Image.from_legacy(legacy_image_1), rtol=1e-13)
363def compare_aperture_corrections_to_legacy(
364 tc: unittest.TestCase,
365 aperture_corrections: ApertureCorrectionMap,
366 legacy_ap_corr_map: Any,
367 subimage_bbox: Box,
368) -> None:
369 """Test an aperture correction `dict` by comparing it to an equivalent
370 `lsst.afw.image.ApCorrMap`.
372 Parameters
373 ----------
374 tc
375 Test case object with assert methods to use.
376 aperture_corrections
377 Dictionary to test.
378 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
379 Equivalent legacy aperture correction map.
380 subimage_bbox
381 Bounding box for full-image tests.
382 """
383 tc.assertEqual(aperture_corrections.keys(), set(legacy_ap_corr_map.keys()))
384 for name, field in aperture_corrections.items():
385 compare_field_to_legacy(tc, field, legacy_ap_corr_map[name], subimage_bbox)
388def compare_observation_summary_stats_to_legacy(
389 tc: unittest.TestCase,
390 summary_stats: ObservationSummaryStats,
391 legacy_summary_stats: Any,
392) -> None:
393 """Test an ObservationSummaryStats object by comparing it to an equivalent
394 `lsst.afw.image.ExposureSummaryStats`.
396 Parameters
397 ----------
398 tc
399 Test case object with assert methods to use.
400 summary_stats
401 Struct to test.
402 legacy : ``lsst.afw.image.ExposureSummaryStats``
403 Equivalent legacy struct.
404 """
405 for field in dataclasses.fields(legacy_summary_stats):
406 a = getattr(legacy_summary_stats, field.name)
407 b = getattr(summary_stats, field.name)
408 if isinstance(b, tuple):
409 for ai, bi in zip(a, b):
410 tc.assertTrue(ai == bi or (math.isnan(ai) and math.isnan(bi)), f"{field.name}: {a} != {b}")
411 else:
412 tc.assertTrue(a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}")
415def compare_projection_to_legacy_wcs[F: Frame](
416 tc: unittest.TestCase,
417 projection: Projection[F],
418 legacy_wcs: Any,
419 pixel_frame: F,
420 subimage_bbox: Box,
421 is_fits: bool = False,
422) -> None:
423 """Test a Projection object by comparing it to an equivalent
424 `lsst.afw.geom.SkyWcs`.
426 Parameters
427 ----------
428 tc
429 Test case object with assert methods to use.
430 projection
431 Projection to test.
432 legacy_wcs : ``lsst.afw.geom.SkyWcs``
433 Equivalent legacy WCS.
434 pixel_frame
435 Expected pixel frame for the projection.
436 subimage_bbox
437 Bounding box of points to generate for tests.
438 is_fits
439 Whether this projection is expected to be exactly representable as a
440 FIT WCS. If `False` it is assumed to have a FITS approximation
441 attached instead.
442 """
443 # Pixel coordinates to test on over the subimage region of interest:
444 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
445 # Array indices of those pixel values (subtract off bbox starts):
446 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
447 sky_coords = legacy_coords_to_astropy(
448 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
449 )
450 # Test transforming with the Projection itself, which also tests its
451 # nested Transform and an Astropy High-Level WCS view with no origin
452 # change.
453 check_projection(tc, projection, pixel_xy, sky_coords, pixel_frame)
454 # Also test the Astropy High-Level WCS view with an origin change to
455 # array indices.
456 check_astropy_wcs_interface(
457 tc, projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
458 )
459 if is_fits:
460 fits_wcs = projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
461 assert fits_wcs is not None
462 check_astropy_wcs_interface(tc, fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
463 # Use that FITS approximation to check that we can make a
464 # Projection from a FITS WCS, too.
465 fits_projection = Projection.from_fits_wcs(fits_wcs, pixel_frame)
466 check_projection(
467 tc,
468 fits_projection,
469 subimage_array_xy,
470 sky_coords,
471 pixel_frame,
472 pixel_atol=1e-5,
473 )
474 # We want Projections we create from a FITS WCS to be backed by an
475 # AST FrameSet so we can convert them into legacy
476 # `lsst.afw.geom.SkyWcs` objects if desired.
477 tc.assertIn("Begin FrameSet", fits_projection.show())
478 else:
479 tc.assertIsNone(projection.as_fits_wcs(subimage_bbox, allow_approximation=False))
480 # The legacy SkyWcs should instead have a FITS approximation
481 # attached; run the same tests on that.
482 assert projection.fits_approximation is not None
483 compare_projection_to_legacy_wcs(
484 tc,
485 projection.fits_approximation,
486 legacy_wcs.getFitsApproximation(),
487 pixel_frame,
488 subimage_bbox,
489 is_fits=True,
490 )
493def check_transform[I: Frame, O: Frame](
494 tc: unittest.TestCase,
495 transform: Transform[I, O],
496 input_xy: XY[np.ndarray],
497 output_xy: XY[np.ndarray],
498 in_frame: Frame,
499 out_frame: Frame,
500 *,
501 check_inverted: bool = True,
502 in_atol: u.Quantity | None = None,
503 out_atol: u.Quantity | None = None,
504) -> None:
505 """Test Transform against known arrays of input and output points.
507 Parameters
508 ----------
509 tc
510 Test case object with assert methods to use.
511 transform
512 Transform to test.
513 input_xy
514 Arrays of input points.
515 output_xy
516 Arrays of output points.
517 in_frame
518 Expected input frame.
519 out_frame
520 Expected output frame.
521 check_inverted
522 If `True`, recurse (once) to test the inverse transform.
523 in_atol
524 Expected absolute precision of input points.
525 out_atol
526 Expected absolute precision of output points.
527 """
528 tc.assertEqual(transform.in_frame, in_frame)
529 tc.assertEqual(transform.out_frame, out_frame)
530 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None
531 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None
532 # Test array interfaces.
533 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
534 assert_close(tc, test_output_xy.x, output_xy.x, atol=out_atol_v)
535 assert_close(tc, test_output_xy.y, output_xy.y, atol=out_atol_v)
536 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
537 assert_close(tc, test_input_xy.x, input_xy.x, atol=in_atol_v)
538 assert_close(tc, test_input_xy.y, input_xy.y, atol=in_atol_v)
539 # Test scalar interfaces with numpy scalars.
540 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
541 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v)
542 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v)
543 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v)
544 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v)
545 # Test quantity array interfaces.
546 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
547 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
548 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
549 assert_close(tc, test_output_q_xy.x, output_q_xy.x, atol=out_atol)
550 assert_close(tc, test_output_q_xy.y, output_q_xy.y, atol=out_atol)
551 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
552 assert_close(tc, test_input_q_xy.x, input_q_xy.x, atol=in_atol)
553 assert_close(tc, test_input_q_xy.y, input_q_xy.y, atol=in_atol)
554 # Test quantity scalar interfaces.
555 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
556 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
557 ):
558 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol)
559 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol)
560 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol)
561 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol)
562 if check_inverted:
563 # Test the inverse transform.
564 check_transform(
565 tc,
566 transform.inverted(),
567 output_xy,
568 input_xy,
569 out_frame,
570 in_frame,
571 check_inverted=False,
572 out_atol=in_atol,
573 in_atol=out_atol,
574 )
577def check_projection[P: Frame](
578 tc: unittest.TestCase,
579 projection: Projection[P],
580 pixel_xy: XY[np.ndarray],
581 sky_coords: SkyCoord,
582 pixel_frame: Frame,
583 *,
584 pixel_atol: float | None = None,
585 sky_atol: u.Quantity | None = None,
586) -> None:
587 """Test a `.Projection` instance against known arrays of pixel and sky
588 coordinates.
590 Parameters
591 ----------
592 tc
593 Test case object with assert methods to use.
594 projection
595 Projection to test.
596 pixel_xy
597 Arrays of pixel coordinates.
598 sky_coords
599 Corresponding sky coordinates.
600 pixel_frame
601 Expected pixel frame.
602 pixel_atol
603 Expected absolute precision of pixel points.
604 sky_atol
605 Expected absolute precision of sky coordinates.
606 """
607 tc.assertEqual(projection.pixel_frame, pixel_frame)
608 tc.assertEqual(projection.sky_frame, SkyFrame.ICRS)
609 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None
610 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None
611 # Test array interfaces.
612 test_pixel_xy = cast(XY[np.ndarray], projection.sky_to_pixel(sky_coords))
613 assert_close(tc, test_pixel_xy.x, pixel_xy.x, atol=pixel_atol)
614 assert_close(tc, test_pixel_xy.y, pixel_xy.y, atol=pixel_atol)
615 test_sky_astropy = projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
616 assert_close(tc, test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v)
617 assert_close(tc, test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v)
618 # Test scalar interfaces.
619 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
620 assert_close(tc, projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol)
621 assert_close(tc, projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol)
622 assert_close(tc, projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v)
623 assert_close(tc, projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v)
624 # Test the underlying Transform object.
625 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
626 check_transform(
627 tc,
628 projection.pixel_to_sky_transform,
629 pixel_xy,
630 sky_xy,
631 pixel_frame,
632 SkyFrame.ICRS,
633 check_inverted=False,
634 in_atol=pixel_atol_q,
635 out_atol=sky_atol,
636 )
637 check_transform(
638 tc,
639 projection.sky_to_pixel_transform,
640 sky_xy,
641 pixel_xy,
642 SkyFrame.ICRS,
643 pixel_frame,
644 check_inverted=False,
645 in_atol=sky_atol,
646 out_atol=pixel_atol_q,
647 )
648 # Test the Astropy interface adapter.
649 check_astropy_wcs_interface(
650 tc, projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
651 )
654def assert_projections_equal(
655 tc: unittest.TestCase,
656 a: Projection[Any] | None,
657 b: Projection[Any] | None,
658 expect_identity: bool | None = None,
659) -> None:
660 """Test that two `.Projection` instances are equivalent."""
661 if a is None and b is None:
662 return
663 assert a is not None and b is not None
664 match expect_identity:
665 case True:
666 tc.assertIs(a, b)
667 return
668 case False:
669 tc.assertIsNot(a, b)
670 case None if a is b:
671 return
672 tc.assertEqual(a.pixel_frame, b.pixel_frame)
673 tc.assertEqual(a.show(simplified=True), b.show(simplified=True))
674 assert_projections_equal(
675 tc, a.fits_approximation, cast(Projection[Any], b.fits_approximation), expect_identity=False
676 )
679def check_astropy_wcs_interface(
680 tc: unittest.TestCase,
681 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
682 pixel_xy: XY[np.ndarray],
683 sky_coords: SkyCoord,
684 *,
685 pixel_atol: float | None = None,
686 sky_atol: u.Quantity | None = None,
687) -> None:
688 """Test an Astropy WCS instance against known arrays of pixel and
689 sky coordinates.
691 Parameters
692 ----------
693 tc
694 Test case object with assert methods to use.
695 wcs
696 Astropy WCS object to test.
697 pixel_xy
698 Arrays of pixel coordinates.
699 sky_coords
700 Corresponding sky coordinates.
701 pixel_atol
702 Expected absolute precision of pixel points.
703 sky_atol
704 Expected absolute precision of sky coordinates.
705 """
706 test_x, test_y = wcs.world_to_pixel(sky_coords)
707 assert_close(tc, test_x, pixel_xy.x, atol=pixel_atol)
708 assert_close(tc, test_y, pixel_xy.y, atol=pixel_atol)
709 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
710 assert_close(tc, test_sky_coords.ra, sky_coords.ra, atol=sky_atol)
711 assert_close(tc, test_sky_coords.dec, sky_coords.dec, atol=sky_atol)
714def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
715 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array."""
716 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
719def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
720 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
721 coordinate object.
722 """
723 return SkyCoord(
724 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
725 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
726 )
729def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
730 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``."""
731 from lsst.geom import Point2D
733 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]