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