Coverage for python/lsst/images/tests/_checks.py: 40%
441 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-29 22:40 +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_cell_coadds_equal",
17 "assert_close",
18 "assert_equal_allow_nan",
19 "assert_images_equal",
20 "assert_masked_images_equal",
21 "assert_masks_equal",
22 "assert_psfs_equal",
23 "assert_sky_projections_equal",
24 "assert_visit_images_equal",
25 "check_archive_tree_class_invariants",
26 "check_astropy_wcs_interface",
27 "check_projection",
28 "check_transform",
29 "compare_amplifier_to_legacy",
30 "compare_aperture_corrections_to_legacy",
31 "compare_cell_coadd_to_legacy",
32 "compare_detector_to_legacy",
33 "compare_field_to_legacy",
34 "compare_image_to_legacy",
35 "compare_mask_to_legacy",
36 "compare_masked_image_to_legacy",
37 "compare_observation_summary_stats_to_legacy",
38 "compare_photo_calib_to_legacy",
39 "compare_psf_to_legacy",
40 "compare_sky_projection_to_legacy_wcs",
41 "compare_visit_image_to_legacy",
42 "iter_concrete_archive_tree_subclasses",
43 "legacy_coords_to_astropy",
44 "legacy_points_to_xy_array",
45)
47import dataclasses
48import math
49import unittest
50from collections.abc import Iterator, Mapping
51from typing import TYPE_CHECKING, Any, Literal, cast
53import astropy.units as u
54import astropy.wcs.wcsapi
55import numpy as np
56from astropy.coordinates import SkyCoord
58from .._geom import XY, YX, Box
59from .._image import Image
60from .._mask import Mask, MaskPlane
61from .._masked_image import MaskedImage
62from .._observation_summary_stats import ObservationSummaryStats
63from .._transforms import DetectorFrame, Frame, SkyFrame, SkyProjection, TractFrame, Transform
64from .._visit_image import VisitImage
65from ..cameras import Amplifier, Detector, DetectorType, ReadoutCorner
66from ..cells import CellCoadd, CellIJ, CoaddProvenance
67from ..fields import BaseField, ChebyshevField
68from ..psfs import PointSpreadFunction
69from ..serialization import ArchiveTree
71if TYPE_CHECKING:
72 try:
73 from lsst.cell_coadds import MultipleCellCoadd
74 except ImportError:
75 type MultipleCellCoadd = Any # type: ignore[no-redef]
76 try:
77 from lsst.afw.image import PhotoCalib as LegacyPhotoCalib
78 except ImportError:
79 type LegacyPhotoCalib = Any # type: ignore[no-redef]
82def assert_close(
83 tc: unittest.TestCase,
84 a: np.ndarray | u.Quantity | float,
85 b: np.ndarray | u.Quantity | float,
86 **kwargs: Any,
87) -> None:
88 """Test that two arrays, floats, or quantities are almost equal.
90 Parameters
91 ----------
92 tc
93 Test case object with assert methods to use.
94 a
95 Array, float, or quantity to compare.
96 b
97 Array, float, or quantity to compare.
98 **kwargs
99 Forwarded to `astropy.units.allclose`.
100 """
101 tc.assertTrue(u.allclose(a, b, **kwargs), msg=f"{a} != {b}")
104def assert_equal_allow_nan(tc: unittest.TestCase, a: float, b: float) -> None:
105 """Test that two floating point values are equal, with nan == nan.
107 Parameters
108 ----------
109 tc
110 Test case providing the assertion methods to use.
111 a
112 First value to compare.
113 b
114 Second value to compare.
115 """
116 try:
117 tc.assertEqual(a, b)
118 except AssertionError:
119 if not (math.isnan(a) and math.isnan(b)):
120 raise
123def assert_images_equal(
124 tc: unittest.TestCase,
125 a: Image,
126 b: Image,
127 *,
128 rtol: float = 0.0,
129 atol: float = 0.0,
130 expect_view: bool | Literal["array"] | None = None,
131) -> None:
132 """Assert that two images are equal or nearly equal.
134 Parameters
135 ----------
136 tc
137 Test case providing the assertion methods to use.
138 a
139 First image to compare.
140 b
141 Second image to compare.
142 rtol
143 Relative tolerance for the pixel comparison.
144 atol
145 Absolute tolerance for the pixel comparison.
146 expect_view
147 If not `None`, also assert whether ``b`` shares memory with ``a``
148 (i.e. is a view); ``"array"`` checks only the pixel arrays.
149 """
150 tc.assertEqual(a.bbox, b.bbox)
151 tc.assertEqual(a.unit, b.unit)
152 assert_sky_projections_equal(tc, a.sky_projection, b.sky_projection)
153 if expect_view is not None:
154 tc.assertEqual(np.may_share_memory(a.array, b.array), bool(expect_view))
155 if expect_view == "array":
156 tc.assertEqual(a.metadata, b.metadata)
157 else:
158 tc.assertEqual(a.metadata is b.metadata, expect_view)
159 if not expect_view:
160 assert_close(tc, a.array, b.array, atol=atol, rtol=rtol)
161 tc.assertEqual(a.metadata, b.metadata)
164def assert_masks_equal(tc: unittest.TestCase, a: Mask, b: Mask) -> None:
165 """Assert that two masks are equal or nearly equal.
167 Parameters
168 ----------
169 tc
170 Test case providing the assertion methods to use.
171 a
172 First mask to compare.
173 b
174 Second mask to compare.
175 """
176 tc.assertEqual(a.bbox, b.bbox)
177 tc.assertEqual(a.schema, b.schema)
178 tc.assertEqual(a.metadata, b.metadata)
179 assert_sky_projections_equal(tc, a.sky_projection, b.sky_projection)
180 np.testing.assert_array_equal(a.array, b.array)
183def assert_masked_images_equal(
184 tc: unittest.TestCase,
185 a: MaskedImage,
186 b: MaskedImage,
187 *,
188 rtol: float = 0.0,
189 atol: float = 0.0,
190 expect_view: bool | None = None,
191) -> None:
192 """Assert that two masked images are equal or nearly equal.
194 Parameters
195 ----------
196 tc
197 Test case providing the assertion methods to use.
198 a
199 First masked image to compare.
200 b
201 Second masked image to compare.
202 rtol
203 Relative tolerance for the pixel comparison.
204 atol
205 Absolute tolerance for the pixel comparison.
206 expect_view
207 If not `None`, also assert whether ``b`` shares memory with ``a``
208 (i.e. is a view).
209 """
210 tc.assertEqual(a.metadata, b.metadata)
211 assert_sky_projections_equal(tc, a.sky_projection, b.sky_projection)
212 assert_images_equal(tc, a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view)
213 assert_masks_equal(tc, a.mask, b.mask)
214 assert_images_equal(tc, a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view)
217def assert_psfs_equal(
218 tc: unittest.TestCase,
219 psf1: PointSpreadFunction,
220 psf2: PointSpreadFunction,
221 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
222) -> int:
223 """Compare two PSF objets.
225 Parameters
226 ----------
227 tc
228 Test case object with assert methods to use.
229 psf1
230 Point-spread function to test.
231 psf2
232 The other point-spread function to test.
233 points
234 Points to evaluate the PSFs at. If not provided, the intersection of
235 the PSF bounding boxes are used to generate points on a grid.
237 Returns
238 -------
239 `int`
240 The number of points actually tested.
241 """
242 if points is None: 242 ↛ 245line 242 didn't jump to line 245 because the condition on line 242 was always true
243 points = psf1.bounds.bbox.intersection(psf2.bounds.bbox).meshgrid(3).map(np.ravel)
245 tc.assertEqual(psf1.kernel_bbox, psf2.kernel_bbox)
247 n_points_tested: int = 0
248 for x, y in zip(points.x, points.y):
249 # The two PSFs must agree on which points fall inside their input
250 # domain. Querying ``.contains`` directly (rather than relying on
251 # ``compute_kernel_image`` to raise) makes this test tolerant of
252 # implementations that do not raise on out-of-domain points -- in
253 # particular ``CellPointSpreadFunction``, where evaluating in a
254 # missing cell does not always raise ``BoundsError``.
255 contains1 = psf1.bounds.contains(x=x, y=y)
256 contains2 = psf2.bounds.contains(x=x, y=y)
257 tc.assertEqual(
258 contains1,
259 contains2,
260 f"PSFs disagree on whether ({x}, {y}) is in-bounds: psf1={contains1}, psf2={contains2}",
261 )
262 if not contains1: 262 ↛ 263line 262 didn't jump to line 263 because the condition on line 262 was never true
263 continue
264 tc.assertEqual(psf1.compute_kernel_image(x=x, y=y), psf2.compute_kernel_image(x=x, y=y))
265 tc.assertEqual(psf1.compute_stellar_bbox(x=x, y=y), psf2.compute_stellar_bbox(x=x, y=y))
266 tc.assertEqual(psf1.compute_stellar_image(x=x, y=y), psf2.compute_stellar_image(x=x, y=y))
267 n_points_tested += 1
268 return n_points_tested
271def assert_visit_images_equal(
272 tc: unittest.TestCase,
273 a: VisitImage,
274 b: VisitImage,
275 *,
276 expect_view: bool | None = None,
277) -> None:
278 """Assert that two `.VisitImage` instances carry the same persistent state.
280 Extends `assert_masked_images_equal` with the VisitImage-specific
281 attributes (PSF, filter, observation info, detector, aperture
282 corrections, photometric scaling, backgrounds, polygon bounds,
283 summary stats) so a round-trip check on a `.VisitImage` does not
284 silently miss differences in any of them.
286 Parameters
287 ----------
288 tc
289 Test case providing the assertion methods to use.
290 a
291 First visit image to compare.
292 b
293 Second visit image to compare.
294 expect_view
295 If not `None`, also assert whether ``b`` shares memory with ``a``
296 (i.e. is a view).
297 """
298 assert_masked_images_equal(tc, a, b, expect_view=expect_view)
299 tc.assertEqual(a.summary_stats, b.summary_stats)
300 tc.assertEqual(a.physical_filter, b.physical_filter)
301 tc.assertEqual(a.band, b.band)
302 tc.assertEqual(a.obs_info, b.obs_info)
303 tc.assertEqual(a.detector, b.detector)
304 tc.assertEqual(dict(a.aperture_corrections), dict(b.aperture_corrections))
305 tc.assertEqual(a.photometric_scaling, b.photometric_scaling)
306 tc.assertEqual(dict(a.backgrounds), dict(b.backgrounds))
307 tc.assertEqual(a.backgrounds.subtracted, b.backgrounds.subtracted)
308 tc.assertEqual(a.bounds, b.bounds)
309 assert_psfs_equal(tc, a.psf, b.psf)
312def assert_cell_coadds_equal(
313 tc: unittest.TestCase,
314 a: CellCoadd,
315 b: CellCoadd,
316 *,
317 expect_view: bool | None = None,
318) -> None:
319 """Assert that two `.CellCoadd` instances carry the same persistent state.
321 Extends the masked-image-style equality check with the
322 CellCoadd-specific attributes (PSF, cell grid, missing cells,
323 backgrounds, patch/tract, band) so a round-trip check does not
324 silently miss differences in any of them.
326 Parameters
327 ----------
328 tc
329 Test case providing the assertion methods to use.
330 a
331 First cell coadd to compare.
332 b
333 Second cell coadd to compare.
334 expect_view
335 If not `None`, also assert whether ``b`` shares memory with ``a``
336 (i.e. is a view).
337 """
338 assert_masked_images_equal(tc, a, b, expect_view=expect_view)
339 tc.assertEqual(a.band, b.band)
340 tc.assertEqual(a.patch, b.patch)
341 tc.assertEqual(a.tract, b.tract)
342 tc.assertEqual(a.grid, b.grid)
343 tc.assertEqual(a.bounds.missing, b.bounds.missing)
344 tc.assertEqual(dict(a.backgrounds), dict(b.backgrounds))
345 tc.assertEqual(a.backgrounds.subtracted, b.backgrounds.subtracted)
346 assert_psfs_equal(tc, a.psf, b.psf)
349def compare_image_to_legacy(
350 tc: unittest.TestCase, image: Image, legacy_image: Any, expect_view: bool | None = None
351) -> None:
352 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object.
354 Parameters
355 ----------
356 tc
357 Test case providing the assertion methods to use.
358 image
359 Image to compare.
360 legacy_image
361 Legacy `lsst.afw.image.Image` to compare against.
362 expect_view
363 If not `None`, also assert whether ``image`` shares memory with
364 ``legacy_image`` (i.e. is a view).
365 """
366 tc.assertEqual(image.bbox, Box.from_legacy(legacy_image.getBBox()))
367 if expect_view is not None: 367 ↛ 369line 367 didn't jump to line 369 because the condition on line 367 was always true
368 tc.assertEqual(np.may_share_memory(image.array, legacy_image.array), expect_view)
369 if not expect_view: 369 ↛ exitline 369 didn't return from function 'compare_image_to_legacy' because the condition on line 369 was always true
370 np.testing.assert_array_equal(image.array, legacy_image.array)
373def compare_mask_to_legacy(
374 tc: unittest.TestCase, mask: Mask, legacy_mask: Any, plane_map: Mapping[str, MaskPlane] | None = None
375) -> None:
376 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object.
378 Parameters
379 ----------
380 tc
381 Test case providing the assertion methods to use.
382 mask
383 Mask to compare.
384 legacy_mask
385 Legacy `lsst.afw.image.Mask` to compare against.
386 plane_map
387 Mapping from legacy plane name to the new mask plane; defaults to
388 the planes in ``mask.schema``.
389 """
390 tc.assertEqual(mask.bbox, Box.from_legacy(legacy_mask.getBBox()))
391 if plane_map is None: 391 ↛ 393line 391 didn't jump to line 393 because the condition on line 391 was always true
392 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
393 for old_name, new_plane in plane_map.items():
394 np.testing.assert_array_equal(
395 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
396 mask.get(new_plane.name),
397 )
400def compare_masked_image_to_legacy(
401 tc: unittest.TestCase,
402 masked_image: MaskedImage,
403 legacy_masked_image: Any,
404 *,
405 plane_map: Mapping[str, MaskPlane] | None = None,
406 expect_view: bool | None = None,
407 alternates: Mapping[str, Any] | None = None,
408) -> None:
409 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
410 object.
412 Parameters
413 ----------
414 tc
415 Test case to use for asserts.
416 masked_image
417 New image to test.
418 legacy_masked_image
419 Legacy image to test against.
420 plane_map
421 Mapping between new and legacy mask planes.
422 expect_view
423 Whether to test that the image and variance arrays do or do not share
424 memory.
425 alternates
426 A mapping of other versions of one or more (new) components to also
427 check against the legacy versions of those components.
428 """
429 compare_image_to_legacy(tc, masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view)
430 compare_mask_to_legacy(tc, masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map)
431 compare_image_to_legacy(
432 tc, masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view
433 )
434 if alternates: 434 ↛ 435line 434 didn't jump to line 435 because the condition on line 434 was never true
435 if image := alternates.get("image"):
436 compare_image_to_legacy(tc, image, legacy_masked_image.getImage(), expect_view=expect_view)
437 if mask := alternates.get("mask"):
438 compare_mask_to_legacy(tc, mask, legacy_masked_image.getMask(), plane_map=plane_map)
439 if variance := alternates.get("variance"):
440 compare_image_to_legacy(tc, variance, legacy_masked_image.getVariance(), expect_view=expect_view)
443def compare_visit_image_to_legacy(
444 tc: unittest.TestCase,
445 visit_image: VisitImage,
446 legacy_exposure: Any,
447 *,
448 plane_map: Mapping[str, MaskPlane] | None = None,
449 expect_view: bool | None = None,
450 instrument: str,
451 visit: int,
452 detector: int,
453 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
454 alternates: Mapping[str, Any] | None = None,
455) -> None:
456 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
457 object.
459 Parameters
460 ----------
461 tc
462 Test case to use for asserts.
463 visit_image
464 New image to test.
465 legacy_exposure
466 Legacy image to test against.
467 plane_map
468 Mapping between new and legacy mask planes.
469 expect_view
470 Whether to test that the image and variance arrays do or do not share
471 memory.
472 instrument
473 Expected instrument name.
474 visit
475 Expected visit ID.
476 detector
477 Expected detector ID.
478 applied_legacy_photo_calib
479 Legacy `lsst.afw.image.PhotoCalib` already applied to
480 ``legacy_exposure``, used when comparing photometric scaling.
481 alternates
482 A mapping of other versions of one or more (new) components to also
483 check against the legacy versions of those components.
484 """
485 compare_masked_image_to_legacy(
486 tc,
487 visit_image,
488 legacy_exposure.getMaskedImage(),
489 plane_map=plane_map,
490 expect_view=expect_view,
491 alternates=alternates,
492 )
493 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
494 compare_sky_projection_to_legacy_wcs(
495 tc,
496 visit_image.sky_projection,
497 legacy_exposure.getWcs(),
498 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
499 visit_image.bbox,
500 )
501 tc.assertIs(visit_image.sky_projection, visit_image.mask.sky_projection)
502 tc.assertIs(visit_image.sky_projection, visit_image.variance.sky_projection)
503 compare_psf_to_legacy(tc, visit_image.psf, legacy_exposure.getPsf())
504 compare_observation_summary_stats_to_legacy(
505 tc, visit_image.summary_stats, legacy_exposure.info.getSummaryStats()
506 )
507 compare_detector_to_legacy(tc, visit_image.detector, legacy_exposure.getDetector(), is_raw_assembled=True)
508 # Make a tiny box for Field comparisons that need to make arrays; that can
509 # get expensive otherwisre.
510 tiny_bbox = detector_bbox.local[2:4, 3:6]
511 compare_aperture_corrections_to_legacy(
512 tc, visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
513 )
514 compare_photo_calib_to_legacy(
515 tc,
516 visit_image.photometric_scaling,
517 legacy_exposure.info.getPhotoCalib(),
518 applied_legacy_photo_calib=applied_legacy_photo_calib,
519 subimage_bbox=tiny_bbox,
520 )
521 if alternates:
522 if (bbox := alternates.get("bbox")) is not None:
523 tc.assertEqual(bbox, visit_image.bbox)
524 if sky_projection := alternates.get("sky_projection"):
525 compare_sky_projection_to_legacy_wcs(
526 tc,
527 sky_projection,
528 legacy_exposure.getWcs(),
529 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
530 visit_image.bbox,
531 )
532 if psf := alternates.get("psf"):
533 compare_psf_to_legacy(tc, psf, legacy_exposure.getPsf())
534 if summary_stats := alternates.get("summary_stats"):
535 compare_observation_summary_stats_to_legacy(
536 tc, summary_stats, legacy_exposure.info.getSummaryStats()
537 )
538 if detector_obj := alternates.get("detector"):
539 compare_detector_to_legacy(tc, detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True)
540 if obs_info := alternates.get("obs_info"):
541 visitInfo = legacy_exposure.visitInfo
542 tc.assertEqual(obs_info.instrument, visitInfo.getInstrumentLabel())
543 if aperture_corrections := alternates.get("aperture_corrections"):
544 compare_aperture_corrections_to_legacy(
545 tc, aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
546 )
547 if (photometric_scaling := alternates.get("photometic_scaling", ...)) is not ...:
548 compare_photo_calib_to_legacy(
549 tc,
550 photometric_scaling,
551 legacy_exposure.info.getPhotoCalib(),
552 applied_legacy_photo_calib=applied_legacy_photo_calib,
553 subimage_bbox=tiny_bbox,
554 )
557def compare_photo_calib_to_legacy(
558 tc: unittest.TestCase,
559 photometric_scaling: BaseField | None,
560 legacy_photo_calib: LegacyPhotoCalib,
561 *,
562 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
563 subimage_bbox: Box,
564) -> None:
565 if legacy_photo_calib._isConstant:
566 if legacy_photo_calib.getCalibrationMean() == 1.0:
567 if applied_legacy_photo_calib is None:
568 tc.assertIsNone(photometric_scaling)
569 return
570 else:
571 legacy_photo_calib = applied_legacy_photo_calib
572 if legacy_photo_calib._isConstant:
573 assert isinstance(photometric_scaling, ChebyshevField)
574 assert_close(
575 tc, photometric_scaling.coefficients, np.array([[legacy_photo_calib.getCalibrationMean()]])
576 )
577 else:
578 assert photometric_scaling is not None
579 compare_field_to_legacy(
580 tc,
581 photometric_scaling / legacy_photo_calib.getCalibrationMean(),
582 legacy_photo_calib.computeScaledCalibration(),
583 subimage_bbox,
584 )
587def compare_cell_coadd_to_legacy(
588 tc: unittest.TestCase,
589 cell_coadd: CellCoadd,
590 legacy_cell_coadd: MultipleCellCoadd,
591 *,
592 tract_bbox: Box,
593 plane_map: Mapping[str, MaskPlane] | None = None,
594 alternates: Mapping[str, Any] | None = None,
595 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
596) -> None:
597 """Compare a `.cells.CellCoadd` object to a legacy
598 `lsst.cell_coadds.MultipleCellCoadd` object.
600 Parameters
601 ----------
602 tc
603 Test case to use for asserts.
604 cell_coadd
605 New coadd to test.
606 legacy_cell_coadd
607 Legacy coadd to test against.
608 tract_bbox
609 Bounding box of the full tract.
610 plane_map
611 Mapping between new and legacy mask planes.
612 alternates
613 A mapping of other versions of one or more (new) components to also
614 check against the legacy versions of those components.
615 psf_points
616 Points to use to compare the PSFs.
617 """
618 legacy_stitched = legacy_cell_coadd.stitch(cell_coadd.bbox.to_legacy())
619 compare_image_to_legacy(tc, cell_coadd.image, legacy_stitched.image, expect_view=False)
620 compare_mask_to_legacy(tc, cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map)
621 compare_image_to_legacy(tc, cell_coadd.variance, legacy_stitched.variance, expect_view=False)
622 if legacy_stitched.mask_fractions is not None:
623 compare_image_to_legacy(
624 tc, cell_coadd.mask_fractions["rejected"], legacy_stitched.mask_fractions, expect_view=False
625 )
626 for n in range(legacy_stitched.n_noise_realizations):
627 compare_image_to_legacy(
628 tc, cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False
629 )
630 tc.assertEqual(cell_coadd.skymap, legacy_stitched.identifiers.skymap)
631 tc.assertEqual(cell_coadd.tract, legacy_stitched.identifiers.tract)
632 tc.assertEqual(cell_coadd.patch.index.x, legacy_stitched.identifiers.patch.x)
633 tc.assertEqual(cell_coadd.patch.index.y, legacy_stitched.identifiers.patch.y)
634 tc.assertEqual(cell_coadd.band, legacy_stitched.identifiers.band)
635 tc.assertTrue(tract_bbox.contains(cell_coadd.patch.outer_bbox))
636 tc.assertTrue(cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox))
637 tc.assertTrue(cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox))
638 tc.assertEqual(cell_coadd.unit, u.Unit(legacy_cell_coadd.common.units.value))
639 tc.assertTrue(cell_coadd.bounds.bbox.contains(cell_coadd.bbox))
640 tc.assertTrue(cell_coadd.grid.bbox.contains(cell_coadd.bbox))
641 compare_sky_projection_to_legacy_wcs(
642 tc,
643 cell_coadd.sky_projection,
644 legacy_cell_coadd.wcs,
645 TractFrame(
646 skymap=legacy_cell_coadd.identifiers.skymap,
647 tract=legacy_cell_coadd.identifiers.tract,
648 bbox=tract_bbox,
649 ),
650 cell_coadd.bbox,
651 is_fits=True,
652 )
653 tc.assertIs(cell_coadd.sky_projection, cell_coadd.mask.sky_projection)
654 tc.assertIs(cell_coadd.sky_projection, cell_coadd.variance.sky_projection)
655 compare_psf_to_legacy(
656 tc, cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points
657 )
658 compare_aperture_corrections_to_legacy(
659 tc, cell_coadd.aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
660 )
661 compare_cell_coadd_provenance_to_legacy(tc, cell_coadd.provenance, legacy_cell_coadd)
662 if alternates:
663 if sky_projection := alternates.get("sky_projection"):
664 compare_sky_projection_to_legacy_wcs(
665 tc,
666 sky_projection,
667 legacy_stitched.wcs,
668 TractFrame(
669 skymap=legacy_cell_coadd.identifiers.skymap,
670 tract=legacy_cell_coadd.identifiers.tract,
671 bbox=tract_bbox,
672 ),
673 cell_coadd.bbox,
674 is_fits=True,
675 )
676 if psf := alternates.get("psf"):
677 compare_psf_to_legacy(tc, psf, legacy_stitched.psf, points=psf_points)
678 if aperture_corrections := alternates.get("aperture_corrections"):
679 compare_aperture_corrections_to_legacy(
680 tc, aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
681 )
682 if provenance := alternates.get("provenance"):
683 compare_cell_coadd_provenance_to_legacy(tc, provenance, legacy_cell_coadd)
686def compare_cell_coadd_provenance_to_legacy(
687 tc: unittest.TestCase, provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd
688) -> None:
689 """Compare a `.cells.CoaddProvenance` object to a legacy
690 `lsst.cell_coadds.MultipleCellCoadd` object.
692 Parameters
693 ----------
694 tc
695 Test case to use for asserts.
696 provenance
697 New provenance object to test.
698 legacy_cell_coadd
699 Legacy coadd to test against.
700 """
701 from lsst.cell_coadds import ObservationIdentifiers
703 for legacy_cell in legacy_cell_coadd.cells.values():
704 cell_index = CellIJ.from_legacy(legacy_cell.identifiers.cell)
705 prov = provenance[cell_index]
706 legacy_table = astropy.table.Table(
707 rows=[
708 [
709 ids.instrument,
710 ids.visit,
711 ids.detector,
712 ids.day_obs,
713 ids.physical_filter,
714 legacy_input.overlaps_center,
715 legacy_input.overlap_fraction,
716 legacy_input.weight,
717 legacy_input.psf_shape.getIxx(),
718 legacy_input.psf_shape.getIyy(),
719 legacy_input.psf_shape.getIxy(),
720 legacy_input.psf_shape_flag,
721 ]
722 for ids, legacy_input in legacy_cell.inputs.items()
723 ],
724 dtype=[
725 np.object_,
726 np.uint64,
727 np.uint16,
728 np.uint32,
729 np.object_,
730 np.bool_,
731 np.float64,
732 np.float64,
733 np.float64,
734 np.float64,
735 np.float64,
736 np.bool_,
737 ],
738 names=[
739 "instrument",
740 "visit",
741 "detector",
742 "day_obs",
743 "physical_filter",
744 "overlaps_center",
745 "overlap_fraction",
746 "weight",
747 "psf_shape_xx",
748 "psf_shape_yy",
749 "psf_shape_xy",
750 "psf_shape_flag",
751 ],
752 )
753 # For a single cell all 'inputs' are also 'contributions'.
754 tc.assertEqual(len(legacy_cell.inputs), len(prov.inputs))
755 tc.assertEqual(len(legacy_cell.inputs), len(prov.contributions))
756 prov.inputs.sort(["instrument", "visit", "detector"])
757 prov.contributions.sort(["instrument", "visit", "detector"])
758 legacy_table.sort(["instrument", "visit", "detector"])
759 np.testing.assert_array_equal(prov.inputs["instrument"], prov.contributions["instrument"])
760 np.testing.assert_array_equal(prov.inputs["visit"], prov.contributions["visit"])
761 np.testing.assert_array_equal(prov.inputs["detector"], prov.contributions["detector"])
762 np.testing.assert_array_equal(prov.inputs["instrument"], legacy_table["instrument"])
763 np.testing.assert_array_equal(prov.inputs["visit"], legacy_table["visit"])
764 np.testing.assert_array_equal(prov.inputs["detector"], legacy_table["detector"])
765 np.testing.assert_array_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"])
766 np.testing.assert_array_equal(prov.inputs["day_obs"], legacy_table["day_obs"])
767 np.testing.assert_array_equal(prov.contributions["overlaps_center"], legacy_table["overlaps_center"])
768 np.testing.assert_array_equal(
769 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"]
770 )
771 np.testing.assert_array_equal(prov.contributions["weight"], legacy_table["weight"])
772 np.testing.assert_array_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"])
773 np.testing.assert_array_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"])
774 np.testing.assert_array_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"])
775 np.testing.assert_array_equal(prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"])
776 for row in prov.inputs:
777 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"})
778 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key]
779 tc.assertEqual(legacy_polygon, row["polygon"].to_legacy())
782def compare_psf_to_legacy(
783 tc: unittest.TestCase,
784 psf: PointSpreadFunction,
785 legacy_psf: Any,
786 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
787 expect_legacy_raise_on_out_of_bounds: bool = False,
788) -> int:
789 """Compare a PSF model object to its legacy interface.
791 Parameters
792 ----------
793 tc
794 Test case object with assert methods to use.
795 psf
796 Point-spread function to test.
797 legacy_psf
798 Legacy `lsst.afw.detection.Psf` instance to compare with.
799 points
800 Points to evaluate the PSFs at. If not provided, the intersection of
801 the PSF bounding boxes are used to generate points on a grid.
802 expect_legacy_raise_on_out_of_bounds
803 If `True`, expect ``legacy_psf`` to raise
804 `lsst.afw.detection.InvalidPsfError` when evaluated at a position
805 considered out-of-bounds by ``psf``.
807 Returns
808 -------
809 `int`
810 The number of points actually tested.
811 """
812 from lsst.afw.detection import InvalidPsfError
814 if points is None:
815 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel)
816 legacy_points = arrays_to_legacy_points(points.x, points.y)
817 n_points_tested: int = 0
818 for p in legacy_points:
819 if not psf.bounds.contains(x=p.x, y=p.y):
820 if expect_legacy_raise_on_out_of_bounds:
821 with tc.assertRaises(InvalidPsfError):
822 legacy_psf.computeKernelImage(p)
823 continue
824 tc.assertEqual(psf.kernel_bbox, Box.from_legacy(legacy_psf.computeKernelBBox(p)))
825 tc.assertEqual(
826 psf.compute_kernel_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeKernelImage(p))
827 )
828 tc.assertEqual(
829 psf.compute_stellar_bbox(x=p.x, y=p.y), Box.from_legacy(legacy_psf.computeImageBBox(p))
830 )
831 tc.assertEqual(psf.compute_stellar_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeImage(p)))
832 n_points_tested += 1
833 return n_points_tested
836def compare_field_to_legacy(
837 tc: unittest.TestCase,
838 field: BaseField,
839 legacy_field: Any,
840 subimage_bbox: Box,
841) -> None:
842 """Test a Field object by comparing it to an equivalent
843 `lsst.afw.math.BoundedField`.
845 Parameters
846 ----------
847 tc
848 Test case object with assert methods to use.
849 field
850 Field to test.
851 legacy_field : ``lsst.afw.math.BoundedField``
852 Equivalent legacy bounded field.
853 subimage_bbox
854 Bounding box for full-image tests.
855 """
856 from lsst.afw.math import BoundedField as LegacyBoundedField
858 tc.assertEqual(field.bounds.bbox, Box.from_legacy(legacy_field.getBBox()))
859 # Pixel coordinates to test the numpy array interface with.
860 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
861 if not isinstance(field.bounds, Box): 861 ↛ 862line 861 didn't jump to line 862 because the condition on line 861 was never true
862 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y)
863 pixel_xy = pixel_xy.map(lambda v: v[mask])
864 try:
865 assert_close(
866 tc,
867 field(x=pixel_xy.x, y=pixel_xy.y),
868 legacy_field.evaluate(pixel_xy.x, pixel_xy.y),
869 equal_nan=True,
870 )
871 except AssertionError as err:
872 err.add_note(f"evaluated at {pixel_xy}")
873 raise
874 if not isinstance(legacy_field, LegacyBoundedField): 874 ↛ 877line 874 didn't jump to line 877 because the condition on line 874 was never true
875 # Legacy StitchedApertureCorrection objects are not true BoundedFields
876 # and don't have addToImage.
877 return
878 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
879 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
880 assert_images_equal(
881 tc, field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13
882 )
885def compare_aperture_corrections_to_legacy(
886 tc: unittest.TestCase,
887 aperture_corrections: Mapping[str, BaseField],
888 legacy_ap_corr_map: Any,
889 subimage_bbox: Box,
890) -> None:
891 """Test an aperture correction `dict` by comparing it to an equivalent
892 `lsst.afw.image.ApCorrMap`.
894 Parameters
895 ----------
896 tc
897 Test case object with assert methods to use.
898 aperture_corrections
899 Dictionary to test.
900 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
901 Equivalent legacy aperture correction map.
902 subimage_bbox
903 Bounding box for full-image tests.
904 """
905 tc.assertEqual(aperture_corrections.keys(), set(legacy_ap_corr_map.keys()))
906 for name, field in aperture_corrections.items():
907 compare_field_to_legacy(tc, field, legacy_ap_corr_map[name], subimage_bbox)
910def compare_observation_summary_stats_to_legacy(
911 tc: unittest.TestCase,
912 summary_stats: ObservationSummaryStats,
913 legacy_summary_stats: Any,
914) -> None:
915 """Test an ObservationSummaryStats object by comparing it to an equivalent
916 `lsst.afw.image.ExposureSummaryStats`.
918 Parameters
919 ----------
920 tc
921 Test case object with assert methods to use.
922 summary_stats
923 Struct to test.
924 legacy_summary_stats : ``lsst.afw.image.ExposureSummaryStats``
925 Equivalent legacy struct.
926 """
927 for field in dataclasses.fields(legacy_summary_stats):
928 a = getattr(legacy_summary_stats, field.name)
929 b = getattr(summary_stats, field.name)
930 if isinstance(b, tuple):
931 for ai, bi in zip(a, b):
932 tc.assertTrue(ai == bi or (math.isnan(ai) and math.isnan(bi)), f"{field.name}: {a} != {b}")
933 else:
934 tc.assertTrue(a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}")
937def compare_sky_projection_to_legacy_wcs[F: Frame](
938 tc: unittest.TestCase,
939 sky_projection: SkyProjection[F],
940 legacy_wcs: Any,
941 pixel_frame: F,
942 subimage_bbox: Box,
943 is_fits: bool = False,
944) -> None:
945 """Test a Projection object by comparing it to an equivalent
946 `lsst.afw.geom.SkyWcs`.
948 Parameters
949 ----------
950 tc
951 Test case object with assert methods to use.
952 sky_projection
953 Projection to test.
954 legacy_wcs : ``lsst.afw.geom.SkyWcs``
955 Equivalent legacy WCS.
956 pixel_frame
957 Expected pixel frame for the sky_projection.
958 subimage_bbox
959 Bounding box of points to generate for tests.
960 is_fits
961 Whether this sky_projection is expected to be exactly representable as
962 a FITS WCS. If `False` it is assumed to have a FITS approximation
963 attached instead.
964 """
965 # Pixel coordinates to test on over the subimage region of interest:
966 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
967 # Array indices of those pixel values (subtract off bbox starts):
968 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
969 sky_coords = legacy_coords_to_astropy(
970 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
971 )
972 # Test transforming with the Projection itself, which also tests its
973 # nested Transform and an Astropy High-Level WCS view with no origin
974 # change.
975 check_projection(tc, sky_projection, pixel_xy, sky_coords, pixel_frame)
976 # Also test the Astropy High-Level WCS view with an origin change to
977 # array indices.
978 check_astropy_wcs_interface(
979 tc, sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
980 )
981 if is_fits:
982 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
983 assert fits_wcs is not None
984 check_astropy_wcs_interface(tc, fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
985 # Use that FITS approximation to check that we can make a
986 # Projection from a FITS WCS, too.
987 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
988 check_projection(
989 tc,
990 fits_projection,
991 subimage_array_xy,
992 sky_coords,
993 pixel_frame,
994 pixel_atol=1e-5,
995 )
996 # We want Projections we create from a FITS WCS to be backed by an
997 # AST FrameSet so we can convert them into legacy
998 # `lsst.afw.geom.SkyWcs` objects if desired.
999 tc.assertIn("Begin FrameSet", fits_projection.show())
1000 else:
1001 tc.assertIsNone(sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False))
1002 # The legacy SkyWcs should instead have a FITS approximation
1003 # attached; run the same tests on that.
1004 assert sky_projection.fits_approximation is not None
1005 compare_sky_projection_to_legacy_wcs(
1006 tc,
1007 sky_projection.fits_approximation,
1008 legacy_wcs.getFitsApproximation(),
1009 pixel_frame,
1010 subimage_bbox,
1011 is_fits=True,
1012 )
1015def check_transform[I: Frame, O: Frame](
1016 tc: unittest.TestCase,
1017 transform: Transform[I, O],
1018 input_xy: XY[np.ndarray],
1019 output_xy: XY[np.ndarray],
1020 in_frame: Frame,
1021 out_frame: Frame,
1022 *,
1023 check_inverted: bool = True,
1024 in_atol: u.Quantity | None = None,
1025 out_atol: u.Quantity | None = None,
1026) -> None:
1027 """Test Transform against known arrays of input and output points.
1029 Parameters
1030 ----------
1031 tc
1032 Test case object with assert methods to use.
1033 transform
1034 Transform to test.
1035 input_xy
1036 Arrays of input points.
1037 output_xy
1038 Arrays of output points.
1039 in_frame
1040 Expected input frame.
1041 out_frame
1042 Expected output frame.
1043 check_inverted
1044 If `True`, recurse (once) to test the inverse transform.
1045 in_atol
1046 Expected absolute precision of input points.
1047 out_atol
1048 Expected absolute precision of output points.
1049 """
1050 tc.assertEqual(transform.in_frame, in_frame)
1051 tc.assertEqual(transform.out_frame, out_frame)
1052 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None
1053 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None
1054 # Test array interfaces.
1055 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
1056 assert_close(tc, test_output_xy.x, output_xy.x, atol=out_atol_v)
1057 assert_close(tc, test_output_xy.y, output_xy.y, atol=out_atol_v)
1058 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
1059 assert_close(tc, test_input_xy.x, input_xy.x, atol=in_atol_v)
1060 assert_close(tc, test_input_xy.y, input_xy.y, atol=in_atol_v)
1061 # Test scalar interfaces with numpy scalars.
1062 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
1063 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v)
1064 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v)
1065 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v)
1066 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v)
1067 # Test quantity array interfaces.
1068 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
1069 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
1070 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
1071 assert_close(tc, test_output_q_xy.x, output_q_xy.x, atol=out_atol)
1072 assert_close(tc, test_output_q_xy.y, output_q_xy.y, atol=out_atol)
1073 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
1074 assert_close(tc, test_input_q_xy.x, input_q_xy.x, atol=in_atol)
1075 assert_close(tc, test_input_q_xy.y, input_q_xy.y, atol=in_atol)
1076 # Test quantity scalar interfaces.
1077 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
1078 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
1079 ):
1080 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol)
1081 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol)
1082 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol)
1083 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol)
1084 if check_inverted:
1085 # Test the inverse transform.
1086 check_transform(
1087 tc,
1088 transform.inverted(),
1089 output_xy,
1090 input_xy,
1091 out_frame,
1092 in_frame,
1093 check_inverted=False,
1094 out_atol=in_atol,
1095 in_atol=out_atol,
1096 )
1099def check_projection[P: Frame](
1100 tc: unittest.TestCase,
1101 sky_projection: SkyProjection[P],
1102 pixel_xy: XY[np.ndarray],
1103 sky_coords: SkyCoord,
1104 pixel_frame: Frame,
1105 *,
1106 pixel_atol: float | None = None,
1107 sky_atol: u.Quantity | None = None,
1108) -> None:
1109 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1110 coordinates.
1112 Parameters
1113 ----------
1114 tc
1115 Test case object with assert methods to use.
1116 sky_projection
1117 Projection to test.
1118 pixel_xy
1119 Arrays of pixel coordinates.
1120 sky_coords
1121 Corresponding sky coordinates.
1122 pixel_frame
1123 Expected pixel frame.
1124 pixel_atol
1125 Expected absolute precision of pixel points.
1126 sky_atol
1127 Expected absolute precision of sky coordinates.
1128 """
1129 tc.assertEqual(sky_projection.pixel_frame, pixel_frame)
1130 tc.assertEqual(sky_projection.sky_frame, SkyFrame.ICRS)
1131 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None
1132 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None
1133 # Test array interfaces.
1134 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1135 assert_close(tc, test_pixel_xy.x, pixel_xy.x, atol=pixel_atol)
1136 assert_close(tc, test_pixel_xy.y, pixel_xy.y, atol=pixel_atol)
1137 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1138 assert_close(tc, test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v)
1139 assert_close(tc, test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v)
1140 # Test scalar interfaces.
1141 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1142 assert_close(tc, sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol)
1143 assert_close(tc, sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol)
1144 assert_close(tc, sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v)
1145 assert_close(
1146 tc, sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v
1147 )
1148 # Test the underlying Transform object.
1149 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1150 check_transform(
1151 tc,
1152 sky_projection.pixel_to_sky_transform,
1153 pixel_xy,
1154 sky_xy,
1155 pixel_frame,
1156 SkyFrame.ICRS,
1157 check_inverted=False,
1158 in_atol=pixel_atol_q,
1159 out_atol=sky_atol,
1160 )
1161 check_transform(
1162 tc,
1163 sky_projection.sky_to_pixel_transform,
1164 sky_xy,
1165 pixel_xy,
1166 SkyFrame.ICRS,
1167 pixel_frame,
1168 check_inverted=False,
1169 in_atol=sky_atol,
1170 out_atol=pixel_atol_q,
1171 )
1172 # Test the Astropy interface adapter.
1173 check_astropy_wcs_interface(
1174 tc, sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1175 )
1178def assert_sky_projections_equal(
1179 tc: unittest.TestCase,
1180 a: SkyProjection[Any] | None,
1181 b: SkyProjection[Any] | None,
1182 expect_identity: bool | None = None,
1183) -> None:
1184 """Test that two `.SkyProjection` instances are equivalent.
1186 Parameters
1187 ----------
1188 tc
1189 Test case providing the assertion methods to use.
1190 a
1191 First sky projection to compare.
1192 b
1193 Second sky projection to compare.
1194 expect_identity
1195 If not `None`, assert whether ``a`` and ``b`` are the same object.
1196 """
1197 if a is None and b is None:
1198 return
1199 assert a is not None and b is not None
1200 match expect_identity:
1201 case True:
1202 tc.assertIs(a, b)
1203 return
1204 case False:
1205 tc.assertIsNot(a, b)
1206 case None if a is b:
1207 return
1208 tc.assertEqual(a, b)
1211def check_astropy_wcs_interface(
1212 tc: unittest.TestCase,
1213 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1214 pixel_xy: XY[np.ndarray],
1215 sky_coords: SkyCoord,
1216 *,
1217 pixel_atol: float | None = None,
1218 sky_atol: u.Quantity | None = None,
1219) -> None:
1220 """Test an Astropy WCS instance against known arrays of pixel and
1221 sky coordinates.
1223 Parameters
1224 ----------
1225 tc
1226 Test case object with assert methods to use.
1227 wcs
1228 Astropy WCS object to test.
1229 pixel_xy
1230 Arrays of pixel coordinates.
1231 sky_coords
1232 Corresponding sky coordinates.
1233 pixel_atol
1234 Expected absolute precision of pixel points.
1235 sky_atol
1236 Expected absolute precision of sky coordinates.
1237 """
1238 test_x, test_y = wcs.world_to_pixel(sky_coords)
1239 assert_close(tc, test_x, pixel_xy.x, atol=pixel_atol)
1240 assert_close(tc, test_y, pixel_xy.y, atol=pixel_atol)
1241 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
1242 assert_close(tc, test_sky_coords.ra, sky_coords.ra, atol=sky_atol)
1243 assert_close(tc, test_sky_coords.dec, sky_coords.dec, atol=sky_atol)
1246def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
1247 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array.
1249 Parameters
1250 ----------
1251 legacy_points
1252 Legacy ``lsst.geom.Point2D`` objects to convert.
1253 """
1254 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
1257def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
1258 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
1259 coordinate object.
1261 Parameters
1262 ----------
1263 legacy_coords
1264 Legacy ``lsst.geom.SpherePoint`` objects to convert.
1265 """
1266 return SkyCoord(
1267 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
1268 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
1269 )
1272def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
1273 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``.
1275 Parameters
1276 ----------
1277 x
1278 X coordinates of the points.
1279 y
1280 Y coordinates of the points.
1281 """
1282 from lsst.geom import Point2D
1284 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
1287def compare_amplifier_to_legacy(
1288 tc: unittest.TestCase,
1289 amplifier: Amplifier,
1290 legacy_amplifier: Any,
1291 *,
1292 is_raw_assembled: bool,
1293 expect_nominal_calibrations: bool = True,
1294) -> None:
1295 """Compare an `~.cameras.Amplifier` to a legacy
1296 `lsst.afw.cameraGeom.Amplifier`.
1298 Parameters
1299 ----------
1300 tc
1301 Test case providing the assertion methods to use.
1302 amplifier
1303 Amplifier to compare.
1304 legacy_amplifier
1305 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against.
1306 is_raw_assembled
1307 Whether the raw geometry is expected to be the assembled-raw
1308 geometry (`True`) or the unassembled-raw geometry (`False`).
1309 expect_nominal_calibrations
1310 Whether the amplifier is expected to carry nominal calibrations.
1311 """
1312 tc.assertEqual(legacy_amplifier.getName(), amplifier.name)
1313 tc.assertEqual(Box.from_legacy(legacy_amplifier.getBBox()), amplifier.bbox)
1314 if is_raw_assembled:
1315 raw_geom = amplifier.assembled_raw_geometry
1316 else:
1317 raw_geom = amplifier.unassembled_raw_geometry
1318 assert raw_geom is not None
1319 tc.assertEqual(ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()), raw_geom.readout_corner)
1320 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawBBox()), raw_geom.bbox)
1321 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawDataBBox()), raw_geom.data_bbox)
1322 tc.assertEqual(legacy_amplifier.getRawFlipX(), raw_geom.flip_x)
1323 tc.assertEqual(legacy_amplifier.getRawFlipY(), raw_geom.flip_y)
1324 tc.assertEqual(legacy_amplifier.getRawXYOffset().getX(), raw_geom.x_offset)
1325 tc.assertEqual(legacy_amplifier.getRawXYOffset().getY(), raw_geom.y_offset)
1326 tc.assertEqual(
1327 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()), raw_geom.horizontal_overscan_bbox
1328 )
1329 tc.assertEqual(
1330 Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()), raw_geom.vertical_overscan_bbox
1331 )
1332 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawPrescanBBox()), raw_geom.horizontal_prescan_bbox)
1333 if expect_nominal_calibrations:
1334 assert amplifier.nominal_calibrations is not None
1335 assert_equal_allow_nan(tc, legacy_amplifier.getGain(), amplifier.nominal_calibrations.gain)
1336 assert_equal_allow_nan(tc, legacy_amplifier.getReadNoise(), amplifier.nominal_calibrations.read_noise)
1337 assert_equal_allow_nan(
1338 tc, legacy_amplifier.getSaturation(), amplifier.nominal_calibrations.saturation
1339 )
1340 assert_equal_allow_nan(
1341 tc, legacy_amplifier.getSuspectLevel(), amplifier.nominal_calibrations.suspect_level
1342 )
1343 np.testing.assert_array_equal(
1344 legacy_amplifier.getLinearityCoeffs(), amplifier.nominal_calibrations.linearity_coefficients
1345 )
1346 tc.assertEqual(legacy_amplifier.getLinearityType(), amplifier.nominal_calibrations.linearity_type)
1349def compare_detector_to_legacy(
1350 tc: unittest.TestCase,
1351 detector: Detector,
1352 legacy_detector: Any,
1353 *,
1354 is_raw_assembled: bool,
1355 expect_nominal_calibrations: bool = True,
1356) -> None:
1357 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`.
1359 Parameters
1360 ----------
1361 tc
1362 Test case providing the assertion methods to use.
1363 detector
1364 Detector to compare.
1365 legacy_detector
1366 Legacy `lsst.afw.cameraGeom.Detector` to compare against.
1367 is_raw_assembled
1368 Whether the raw geometry is expected to be the assembled-raw
1369 geometry (`True`) or the unassembled-raw geometry (`False`).
1370 expect_nominal_calibrations
1371 Whether the detector's amplifiers are expected to carry nominal
1372 calibrations.
1373 """
1374 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
1376 tc.assertEqual(legacy_detector.getName(), detector.name)
1377 tc.assertEqual(legacy_detector.getId(), detector.id)
1378 tc.assertEqual(DetectorType.from_legacy(legacy_detector.getType()), detector.type)
1379 tc.assertEqual(Box.from_legacy(legacy_detector.getBBox()), detector.bbox)
1380 tc.assertEqual(legacy_detector.getSerial(), detector.serial)
1381 legacy_orientation = legacy_detector.getOrientation()
1382 tc.assertEqual(legacy_orientation.getFpPosition3().getX(), detector.orientation.focal_plane_x)
1383 tc.assertEqual(legacy_orientation.getFpPosition3().getY(), detector.orientation.focal_plane_y)
1384 tc.assertEqual(legacy_orientation.getFpPosition3().getZ(), detector.orientation.focal_plane_z)
1385 tc.assertEqual(legacy_orientation.getReferencePoint().getX(), detector.orientation.pixel_reference_x)
1386 tc.assertEqual(legacy_orientation.getReferencePoint().getY(), detector.orientation.pixel_reference_y)
1387 tc.assertEqual(legacy_orientation.getYaw().asRadians(), detector.orientation.yaw.to_value(u.rad))
1388 tc.assertEqual(legacy_orientation.getPitch().asRadians(), detector.orientation.pitch.to_value(u.rad))
1389 tc.assertEqual(legacy_orientation.getRoll().asRadians(), detector.orientation.roll.to_value(u.rad))
1390 tc.assertEqual(legacy_detector.getPixelSize().getX(), detector.pixel_size)
1391 tc.assertEqual(legacy_detector.getPhysicalType(), detector.physical_type)
1392 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
1393 compare_amplifier_to_legacy(
1394 tc,
1395 amplifier,
1396 legacy_amplifier,
1397 is_raw_assembled=is_raw_assembled,
1398 expect_nominal_calibrations=expect_nominal_calibrations,
1399 )
1400 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
1401 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
1402 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
1403 check_transform(
1404 tc,
1405 detector.to_focal_plane,
1406 pixel_xy,
1407 legacy_points_to_xy_array(fp_legacy_points),
1408 detector.frame,
1409 detector.to_focal_plane.out_frame,
1410 in_atol=1e-9 * u.pix,
1411 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
1412 )
1413 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
1414 check_transform(
1415 tc,
1416 detector.to_field_angle,
1417 pixel_xy,
1418 legacy_points_to_xy_array(fa_legacy_points),
1419 detector.frame,
1420 detector.to_field_angle.out_frame,
1421 in_atol=1e-9 * u.pix,
1422 out_atol=1e-7 * u.arcsec,
1423 )
1426def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
1427 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
1429 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
1430 classes. Importing this module already imports every ``lsst.images``
1431 module that defines a subclass, so the tree is fully populated by the time
1432 this is called.
1434 This discovery is deliberately separate from
1435 `check_archive_tree_class_invariants` so that the per-class check stays
1436 usable on a single class even if this metaprogramming is removed later.
1437 """
1438 seen: set[type] = set()
1439 stack: list[type] = [ArchiveTree]
1440 while stack:
1441 kls = stack.pop()
1442 for sub in kls.__subclasses__():
1443 if sub in seen: 1443 ↛ 1444line 1443 didn't jump to line 1444 because the condition on line 1443 was never true
1444 continue
1445 seen.add(sub)
1446 stack.append(sub)
1447 if not getattr(sub, "__abstractmethods__", None): 1447 ↛ 1442line 1447 didn't jump to line 1442 because the condition on line 1447 was always true
1448 yield sub
1451def check_archive_tree_class_invariants(tc: unittest.TestCase, cls: type[ArchiveTree]) -> None:
1452 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
1453 well-formed schema-version constants and an in-memory type.
1455 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
1456 ``PUBLIC_TYPE`` are present and well-typed, that the version is
1457 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
1458 schema major.
1460 Parameters
1461 ----------
1462 tc
1463 Test case object with assert methods to use.
1464 cls
1465 The concrete `.serialization.ArchiveTree` subclass to check.
1466 """
1467 tc.assertTrue(hasattr(cls, "SCHEMA_NAME"), f"{cls.__name__} lacks SCHEMA_NAME")
1468 tc.assertTrue(hasattr(cls, "SCHEMA_VERSION"), f"{cls.__name__} lacks SCHEMA_VERSION")
1469 tc.assertTrue(hasattr(cls, "MIN_READ_VERSION"), f"{cls.__name__} lacks MIN_READ_VERSION")
1470 tc.assertTrue(hasattr(cls, "PUBLIC_TYPE"), f"{cls.__name__} lacks PUBLIC_TYPE")
1471 tc.assertIsInstance(cls.SCHEMA_NAME, str)
1472 tc.assertGreater(len(cls.SCHEMA_NAME), 0)
1473 tc.assertRegex(cls.SCHEMA_VERSION, r"^\d+\.\d+\.\d+$")
1474 tc.assertIsInstance(cls.MIN_READ_VERSION, int)
1475 tc.assertGreaterEqual(cls.MIN_READ_VERSION, 1)
1476 tc.assertIsInstance(cls.PUBLIC_TYPE, type)
1477 major = int(cls.SCHEMA_VERSION.split(".")[0])
1478 tc.assertLessEqual(cls.MIN_READ_VERSION, major)