Coverage for python/lsst/images/tests/_checks.py: 69%
666 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-07 23:31 +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 "annotate_errors",
16 "arrays_to_legacy_points",
17 "assert_cell_coadds_equal",
18 "assert_images_equal",
19 "assert_masked_images_equal",
20 "assert_masks_equal",
21 "assert_obs_metadata_fields_equal",
22 "assert_psfs_equal",
23 "assert_sky_projections_equal",
24 "assert_values_equal",
25 "assert_visit_images_equal",
26 "check_archive_tree_class_invariants",
27 "check_astropy_wcs_interface",
28 "check_bounds_contains_broadcasting",
29 "check_projection",
30 "check_transform",
31 "compare_amplifier_to_legacy",
32 "compare_aperture_corrections_to_legacy",
33 "compare_cell_coadd_to_legacy",
34 "compare_detector_to_legacy",
35 "compare_field_to_legacy",
36 "compare_image_to_legacy",
37 "compare_mask_to_legacy",
38 "compare_masked_image_to_legacy",
39 "compare_observation_summary_stats_to_legacy",
40 "compare_photo_calib_to_legacy",
41 "compare_psf_to_legacy",
42 "compare_sky_projection_to_legacy_wcs",
43 "compare_visit_image_to_legacy",
44 "iter_concrete_archive_tree_subclasses",
45 "legacy_coords_to_astropy",
46 "legacy_points_to_xy_array",
47)
49import dataclasses
50import re
51from abc import ABC, abstractmethod
52from collections.abc import Generator, Iterator, Mapping
53from contextlib import contextmanager
54from typing import TYPE_CHECKING, Any, Literal, NoReturn, Self, cast
56import astropy.coordinates as coords
57import astropy.time
58import astropy.units as u
59import astropy.wcs.wcsapi
60import numpy as np
61import pytest
62from astropy.coordinates import SkyCoord
64from .._geom import XY, YX, Bounds, Box
65from .._image import Image
66from .._mask import Mask, MaskPlane, get_legacy_visit_image_mask_planes
67from .._masked_image import MaskedImage
68from .._observation_summary_stats import ObservationSummaryStats
69from .._transforms import DetectorFrame, Frame, SkyFrame, SkyProjection, TractFrame, Transform
70from .._visit_image import VisitImage
71from ..cameras import Amplifier, Detector, DetectorType, ReadoutCorner
72from ..cells import CellCoadd, CellIJ, CoaddProvenance
73from ..fields import BaseField, ChebyshevField
74from ..psfs import PointSpreadFunction
75from ..serialization import ArchiveTree
77if TYPE_CHECKING:
78 try:
79 from lsst.cell_coadds import MultipleCellCoadd
80 except ImportError:
81 type MultipleCellCoadd = Any # type: ignore[no-redef]
82 try:
83 from lsst.afw.image import PhotoCalib as LegacyPhotoCalib
84 except ImportError:
85 type LegacyPhotoCalib = Any # type: ignore[no-redef]
88# The sky-to-pixel direction of a SkyProjection is a numerical AST inversion
89# of the celestial projection, so it carries a small error floor (~1e-9 px).
90# A relative-only tolerance collapses to exact equality at pixel coordinate
91# zero, so comparisons that can include the origin need a small absolute pixel
92# tolerance to avoid reporting that numerical noise as a difference.
93_AST_INV_ATOL = 1e-7
95# Relative tolerance for comparisons of AST-backed transforms and projections.
96_AST_RTOL = 1e-9
98# Default absolute tolerance for great-circle sky-coordinate comparisons, in
99# the same order of magnitude as _AST_INV_ATOL (which is in pixels). A
100# nonzero default is important because sky coordinates can be near zero
101# (e.g. RA/Dec ~ 0), where a relative tolerance provides no error floor.
102_AST_SKY_ATOL = 1e-7 * u.deg
105class _BaseCmpArg(ABC):
106 """A helper base clas for `assert_values_equal` whose subclasses represent
107 different types of comparisons.
109 Parameters
110 ----------
111 value
112 The coerced operand; type is subclass-dependent..
113 """
115 def __init__(self, value: Any):
116 self.value = value
118 @property
119 @abstractmethod
120 def unit(self) -> u.UnitBase:
121 raise NotImplementedError()
123 def coerce_atol(self, atol: float | u.Quantity | astropy.time.TimeDelta) -> float:
124 try:
125 if isinstance(atol, astropy.time.TimeDelta):
126 return atol.to_value(self.unit)
127 if atol != 0.0:
128 return u.Quantity(atol).to_value(self.unit)
129 return 0.0
130 except u.UnitConversionError as err:
131 err.add_note("Invalid atol for these operands.")
132 raise
134 def get_nan_mask(self) -> np.ndarray:
135 return np.isnan(np.asarray(self.value))
137 @abstractmethod
138 def get_separation(self, other: Self) -> np.ndarray:
139 raise NotImplementedError
141 @abstractmethod
142 def find_mismatches(
143 self,
144 other: Self,
145 *,
146 rtol: float,
147 atol_value: float,
148 equal_nan: bool,
149 ) -> np.ndarray | None:
150 """Search for element-wise mismatches against ``other``.
152 Each comparison family implements its own tolerance semantics.
154 Parameters
155 ----------
156 other
157 The other comparison argument.
158 rtol
159 Relative tolerance.
160 atol_value
161 Absolute tolerance, already coerced to this comparison's unit.
162 equal_nan
163 If `True`, treat NaN as equal to NaN wherever the operand can
164 hold NaN.
166 Returns
167 -------
168 mismatch : `numpy.ndarray` | None
169 Boolean mask marking differing elements, or `None` if all
170 elements agree within tolerance.
171 """
172 raise NotImplementedError()
174 def report_mismatches(
175 self,
176 other: Self,
177 a: Any,
178 b: Any,
179 mismatch: np.ndarray,
180 label: str,
181 ) -> NoReturn:
182 """Raise an `AssertionError` reporting a computed mismatch.
184 Parameters
185 ----------
186 other : _BaseCmpArg
187 The other comparison argument.
188 a
189 First original operand; used for the scalar failure message.
190 b
191 Second original operand; used for the scalar failure message.
192 mismatch
193 Boolean mask marking differing elements, from `find_mismatches`.
194 label
195 Prefix prepended to the failure message.
196 """
197 n_mismatch = int(np.count_nonzero(mismatch))
198 prefix = f"{label}: " if label else ""
199 separation = self.get_separation(other)
200 if mismatch.shape == ():
201 raise AssertionError(f"{prefix}{a!r} != {b!r} (differ by {float(separation)}{self.unit})")
202 nan_mismatch = int(np.count_nonzero(mismatch & (self.get_nan_mask() | other.get_nan_mask())))
203 nan_str = f"; {nan_mismatch} NaN mismatches" if nan_mismatch else ""
204 if np.isnan(separation).all():
205 raise AssertionError(f"{prefix}{n_mismatch}/{mismatch.size} values differ{nan_str}")
206 maxdiff = float(np.nanmax(separation))
207 # Convert to plain ints so the index prints as e.g. (2198,)
208 # rather than (np.int64(2198),).
209 loc = tuple(int(i) for i in np.unravel_index(np.nanargmax(separation), separation.shape))
210 raise AssertionError(
211 f"{prefix}{n_mismatch}/{mismatch.size} values differ; "
212 f"max abs diff {maxdiff}{self.unit} at index {loc} {nan_str}"
213 )
216class _QuantityCmpArg(_BaseCmpArg):
217 """Direct comparison of (possibly dimensionless)
218 `astropy.units.Quantity` instances, used for all non-separation
219 floating-point comparisons.
220 """
222 value: u.Quantity
224 @property
225 def unit(self) -> u.UnitBase:
226 return self.value.unit
228 def get_separation(self, other: Self) -> np.ndarray:
229 a_vals = self.value.to_value(self.unit)
230 b_vals = other.value.to_value(self.unit)
231 return np.abs(a_vals - b_vals)
233 def find_mismatches(
234 self,
235 other: Self,
236 *,
237 rtol: float,
238 atol_value: float,
239 equal_nan: bool,
240 ) -> np.ndarray | None:
241 a_vals = self.value.to_value(self.unit)
242 b_vals = other.value.to_value(self.unit)
243 mismatch = ~np.isclose(a_vals, b_vals, rtol=rtol, atol=atol_value, equal_nan=equal_nan)
244 return mismatch if np.any(mismatch) else None
247class _SeparationCmpArg(_BaseCmpArg):
248 """A base class for operands compared through a pairwise separation.
250 ``rtol`` is meaningless and must be zero for these comparisons.
251 """
253 @property
254 @abstractmethod
255 def unit(self) -> u.UnitBase:
256 raise NotImplementedError
258 def find_mismatches(
259 self,
260 other: Self,
261 *,
262 rtol: float,
263 atol_value: float,
264 equal_nan: bool,
265 ) -> np.ndarray | None:
266 if rtol != 0.0:
267 raise ValueError(f"rtol is meaningless for separations; use atol (in {self.unit}) instead.")
268 # A NaN separation (from a NaN operand element) fails the ``<=``
269 # test and so counts as a mismatch unless equal_nan clears it.
270 separation = np.asarray(self.get_separation(other))
271 mismatch = np.logical_not(separation <= atol_value)
272 if equal_nan:
273 mismatch = mismatch & ~(self.get_nan_mask() & other.get_nan_mask())
274 return mismatch if np.any(mismatch) else None
277class _TimeCmpArg(_SeparationCmpArg):
278 value: astropy.time.Time
280 @property
281 def unit(self) -> u.UnitBase:
282 return u.s
284 def get_separation(self, other: Self) -> np.ndarray:
285 return np.abs((self.value - other.value).sec)
287 def get_nan_mask(self) -> np.ndarray:
288 # astropy.time.Time rejects non-finite values at construction.
289 return np.zeros(self.value.shape, dtype=bool)
292class _SkyCoordCmpArg(_SeparationCmpArg):
293 value: SkyCoord
295 @property
296 def unit(self) -> u.UnitBase:
297 return u.arcsec
299 def get_separation(self, other: Self) -> np.ndarray:
300 return self.value.separation(other.value).to_value(self.unit)
302 def get_nan_mask(self) -> np.ndarray:
303 # A coordinate counts as NaN only when *both* components are NaN;
304 # a half-NaN coordinate is a difference to report, not a NaN match.
305 return np.isnan(self.value.ra) & np.isnan(self.value.dec)
308class _WrappedAngleCmpArg(_SeparationCmpArg):
309 value: u.Quantity
311 @property
312 def unit(self) -> u.UnitBase:
313 return u.rad
315 def get_separation(self, other: Self) -> np.ndarray:
316 a_rad = np.asarray(self.value.to_value(u.rad))
317 b_rad = np.asarray(other.value.to_value(u.rad))
318 return np.abs((a_rad - b_rad + np.pi) % (2.0 * np.pi) - np.pi)
321def _classify_operand(arg: Any, *, wrap_angles: bool = False) -> _BaseCmpArg | None:
322 """Coerce one comparison operand into its comparison-argument wrapper.
324 Parameters
325 ----------
326 arg
327 Operand to classify.
328 wrap_angles
329 If `True`, ``arg`` must be an angular `~astropy.units.Quantity` and
330 is classified for wrapped-angle comparison.
332 Returns
333 -------
334 cmp_arg : _BaseCmpArg | None
335 Wrapper selecting the comparison space for ``arg``, or `None` for
336 operands that cannot hold floating-point values (integers, bools,
337 strings); those are compared exactly by the driver.
338 """
339 if wrap_angles:
340 return _WrappedAngleCmpArg(u.Quantity(arg).to(u.rad))
341 match arg:
342 case astropy.time.Time():
343 return _TimeCmpArg(arg)
344 case astropy.time.TimeDelta():
345 return _QuantityCmpArg(arg.to(u.s))
346 case SkyCoord():
347 return _SkyCoordCmpArg(arg)
348 case u.Quantity():
349 if arg.dtype.kind not in "fc": 349 ↛ 351line 349 didn't jump to line 351 because the condition on line 349 was never true
350 # Coerce integer Quantities to floating point.
351 arg = arg.astype(float)
352 return _QuantityCmpArg(arg)
353 arg = np.asarray(arg)
354 if arg.dtype.kind not in "fc":
355 return None
356 return _QuantityCmpArg(u.Quantity(arg))
359def assert_values_equal(
360 a: (np.ndarray | u.Quantity | float | astropy.time.Time | astropy.time.TimeDelta | SkyCoord),
361 b: (np.ndarray | u.Quantity | float | astropy.time.Time | astropy.time.TimeDelta | SkyCoord),
362 *,
363 rtol: float = 0.0,
364 atol: float | u.Quantity | astropy.time.TimeDelta = 0.0,
365 equal_nan: bool = True,
366 label: str = "",
367 wrap_angles: bool = False,
368) -> None:
369 """Assert that two arrays, quantities, floats, astropy times, or sky
370 coordinates are equal or close.
372 On mismatch this raises an `AssertionError` with a concise report (count,
373 max absolute difference, location, NaN-mismatch count) rather than dumping
374 the arrays. Exact (non-floating-point) mismatches instead report through
375 `numpy.testing.assert_array_equal`.
377 Both operands must be of the same kind: both floating-point numeric
378 (bare or as `~astropy.units.Quantity` or `~astropy.time.TimeDelta`),
379 or both the same astropy type (`~astropy.time.Time` or
380 `astropy.coordinates.SkyCoord`).
381 Operands that cannot hold floating-point values (integers, bools,
382 strings) are compared exactly, with no tolerance or NaN semantics.
383 Tolerances never have their units guessed: a bare-float ``atol`` is
384 only valid for dimensionless operands, and for unitful operands only
385 exactly ``0.0`` (which is unit-independent) is accepted without units.
387 Parameters
388 ----------
389 a
390 First value to compare.
391 b
392 Second value to compare.
393 rtol
394 Relative tolerance, applied to the operand magnitudes. Must be zero
395 for comparisons of separations.
396 atol
397 Absolute tolerance. For dimensionless operands, a bare float. For
398 unitful operands, a `~astropy.units.Quantity` with units compatible
399 with ``a``, or exactly ``0.0``.
400 equal_nan
401 If `True`, treat NaN as equal to NaN. A `astropy.coordinates.SkyCoord`
402 element counts as NaN only when both of its coordinates are NaN.
403 label
404 Prefix prepended to the failure message.
405 wrap_angles
406 If `True`, compare ``a`` and ``b`` (angular quantities) modulo 2 pi.
407 """
408 a_arg = _classify_operand(a, wrap_angles=wrap_angles)
409 b_arg = _classify_operand(b, wrap_angles=wrap_angles)
410 if type(a_arg) is not type(b_arg):
411 raise TypeError(f"cannot compare {type(a).__name__} with {type(b).__name__}")
412 if a_arg is None or b_arg is None:
413 # The type check above means both are None: exact comparison of
414 # operands that cannot hold floating-point values.
415 if rtol != 0.0 or atol != 0.0:
416 raise ValueError("rtol and atol are only supported for floating-point operands.")
417 np.testing.assert_array_equal(
418 np.asarray(a),
419 np.asarray(b),
420 strict=True,
421 err_msg=f"{label}: arrays are not equal." if label else "",
422 )
423 return
424 atol_value = a_arg.coerce_atol(atol)
425 if a_arg.value.shape != b_arg.value.shape:
426 raise AssertionError(f"{label}: shape {a_arg.value.shape} != {b_arg.value.shape}")
427 mismatch = a_arg.find_mismatches(b_arg, rtol=rtol, atol_value=atol_value, equal_nan=equal_nan)
428 if mismatch is not None:
429 a_arg.report_mismatches(b_arg, a, b, mismatch, label)
432@contextmanager
433def annotate_errors(note: str) -> Generator[None]:
434 """Attach a note to any exception raised in the block.
436 Parameters
437 ----------
438 note
439 Note to attach on failure.
440 """
441 try:
442 yield
443 except Exception as err:
444 err.add_note(note)
445 raise
448def assert_images_equal(
449 a: Image,
450 b: Image,
451 *,
452 rtol: float = 0.0,
453 atol: float = 0.0,
454 expect_view: bool | Literal["array"] | None = None,
455) -> None:
456 """Assert that two images are equal or nearly equal.
458 Parameters
459 ----------
460 a
461 First image to compare.
462 b
463 Second image to compare.
464 rtol
465 Relative tolerance for the pixel comparison.
466 atol
467 Absolute tolerance for the pixel comparison.
468 expect_view
469 If not `None`, also assert whether ``b`` shares memory with ``a``
470 (i.e. is a view); ``"array"`` checks only the pixel arrays.
471 """
472 assert a.bbox == b.bbox
473 assert a.unit == b.unit
474 assert_sky_projections_equal(a.sky_projection, b.sky_projection)
475 if expect_view is not None:
476 assert np.may_share_memory(a.array, b.array) == bool(expect_view)
477 if expect_view == "array":
478 assert a.metadata.native == b.metadata.native
479 else:
480 assert (a._metadata is b._metadata) == expect_view
481 if not expect_view:
482 assert_values_equal(a.array, b.array, atol=atol, rtol=rtol)
483 assert a.metadata.native == b.metadata.native
486def _note_mask_difference(
487 err: AssertionError, a: Mask, b: Mask, *, label: str, names: Mapping[str, str] | None = None
488) -> None:
489 """Embed a plane-by-plane `Mask.compare` summary as a note on ``err``.
491 Parameters
492 ----------
493 err
494 Exception to attach the summary to.
495 a
496 First mask compared.
497 b
498 Second mask compared.
499 label
500 Prefix for the failure message.
501 names
502 Optional mapping from mask-plane names used by ``a`` and ``b`` to
503 names to display in the summary (e.g. the legacy names, when
504 comparing against legacy code).
505 """
506 diff = a.compare(b)
507 if not diff: 507 ↛ 508line 507 didn't jump to line 508 because the condition on line 507 was never true
508 return
509 shown = names or {}
510 summary = "\n".join(
511 f"{label}[{shown.get(name, name)}]: +{added} -{removed}" for name, (added, removed) in diff.items()
512 )
513 err.add_note(summary)
516def assert_masks_equal(a: Mask, b: Mask) -> None:
517 """Assert that two masks are equal or nearly equal.
519 Parameters
520 ----------
521 a
522 First mask to compare.
523 b
524 Second mask to compare.
525 """
526 assert a.bbox == b.bbox
527 assert a.schema == b.schema
528 assert a.metadata.native == b.metadata.native
529 assert_sky_projections_equal(a.sky_projection, b.sky_projection)
530 try:
531 assert_values_equal(a.array, b.array, label="mask")
532 except AssertionError as err:
533 _note_mask_difference(err, a, b, label="mask")
534 raise
537def assert_masked_images_equal(
538 a: MaskedImage,
539 b: MaskedImage,
540 *,
541 rtol: float = 0.0,
542 atol: float = 0.0,
543 expect_view: bool | None = None,
544) -> None:
545 """Assert that two masked images are equal or nearly equal.
547 Parameters
548 ----------
549 a
550 First masked image to compare.
551 b
552 Second masked image to compare.
553 rtol
554 Relative tolerance for the pixel comparison.
555 atol
556 Absolute tolerance for the pixel comparison.
557 expect_view
558 If not `None`, also assert whether ``b`` shares memory with ``a``
559 (i.e. is a view).
560 """
561 assert a.metadata.native == b.metadata.native
562 assert_sky_projections_equal(a.sky_projection, b.sky_projection)
563 assert_images_equal(a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view)
564 assert_masks_equal(a.mask, b.mask)
565 assert_images_equal(a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view)
568def assert_psfs_equal(
569 psf1: PointSpreadFunction,
570 psf2: PointSpreadFunction,
571 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
572) -> int:
573 """Compare two PSF objets.
575 Parameters
576 ----------
577 psf1
578 Point-spread function to test.
579 psf2
580 The other point-spread function to test.
581 points
582 Points to evaluate the PSFs at. If not provided, the intersection of
583 the PSF bounding boxes are used to generate points on a grid.
585 Returns
586 -------
587 `int`
588 The number of points actually tested.
589 """
590 if points is None: 590 ↛ 593line 590 didn't jump to line 593 because the condition on line 590 was always true
591 points = psf1.bounds.bbox.intersection(psf2.bounds.bbox).meshgrid(3).map(np.ravel)
593 assert psf1.kernel_bbox == psf2.kernel_bbox
595 n_points_tested: int = 0
596 for x, y in zip(points.x, points.y):
597 # The two PSFs must agree on which points fall inside their input
598 # domain. Querying ``.contains`` directly (rather than relying on
599 # ``compute_kernel_image`` to raise) makes this test tolerant of
600 # implementations that do not raise on out-of-domain points -- in
601 # particular ``CellPointSpreadFunction``, where evaluating in a
602 # missing cell does not always raise ``BoundsError``.
603 contains1 = psf1.bounds.contains(x=x, y=y)
604 contains2 = psf2.bounds.contains(x=x, y=y)
605 assert contains1 == contains2, (
606 f"PSFs disagree on whether ({x}, {y}) is in-bounds: psf1={contains1}, psf2={contains2}"
607 )
608 if not contains1: 608 ↛ 609line 608 didn't jump to line 609 because the condition on line 608 was never true
609 continue
610 assert psf1.compute_kernel_image(x=x, y=y) == psf2.compute_kernel_image(x=x, y=y)
611 assert psf1.compute_stellar_bbox(x=x, y=y) == psf2.compute_stellar_bbox(x=x, y=y)
612 assert psf1.compute_stellar_image(x=x, y=y) == psf2.compute_stellar_image(x=x, y=y)
613 n_points_tested += 1
614 return n_points_tested
617def assert_visit_images_equal(
618 a: VisitImage,
619 b: VisitImage,
620 *,
621 expect_view: bool | None = None,
622) -> None:
623 """Assert that two `.VisitImage` instances carry the same persistent state.
625 Extends `assert_masked_images_equal` with the VisitImage-specific
626 attributes (PSF, filter, observation info, detector, aperture
627 corrections, photometric scaling, backgrounds, polygon bounds,
628 summary stats) so a round-trip check on a `.VisitImage` does not
629 silently miss differences in any of them.
631 Parameters
632 ----------
633 a
634 First visit image to compare.
635 b
636 Second visit image to compare.
637 expect_view
638 If not `None`, also assert whether ``b`` shares memory with ``a``
639 (i.e. is a view).
640 """
641 assert_masked_images_equal(a, b, expect_view=expect_view)
642 assert a.summary_stats == b.summary_stats
643 assert a.physical_filter == b.physical_filter
644 assert a.band == b.band
645 assert a.obs_info == b.obs_info
646 assert a.detector == b.detector
647 assert dict(a.aperture_corrections) == dict(b.aperture_corrections)
648 assert a.photometric_scaling == b.photometric_scaling
649 assert dict(a.backgrounds) == dict(b.backgrounds)
650 assert a.backgrounds.subtracted == b.backgrounds.subtracted
651 assert a.bounds == b.bounds
652 assert_psfs_equal(a.psf, b.psf)
655def assert_cell_coadds_equal(
656 a: CellCoadd,
657 b: CellCoadd,
658 *,
659 expect_view: bool | None = None,
660) -> None:
661 """Assert that two `.CellCoadd` instances carry the same persistent state.
663 Extends the masked-image-style equality check with the
664 CellCoadd-specific attributes (PSF, cell grid, missing cells,
665 backgrounds, patch/tract, band) so a round-trip check does not
666 silently miss differences in any of them.
668 Parameters
669 ----------
670 a
671 First cell coadd to compare.
672 b
673 Second cell coadd to compare.
674 expect_view
675 If not `None`, also assert whether ``b`` shares memory with ``a``
676 (i.e. is a view).
677 """
678 assert_masked_images_equal(a, b, expect_view=expect_view)
679 assert a.band == b.band
680 assert a.patch == b.patch
681 assert a.tract == b.tract
682 assert a.grid == b.grid
683 assert a.bounds.missing == b.bounds.missing
684 assert dict(a.backgrounds) == dict(b.backgrounds)
685 assert a.backgrounds.subtracted == b.backgrounds.subtracted
686 assert_psfs_equal(a.psf, b.psf)
689def compare_image_to_legacy(
690 image: Image, legacy_image: Any, expect_view: bool | None = None, *, label: str = "image"
691) -> None:
692 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object.
694 Parameters
695 ----------
696 image
697 Image to compare.
698 legacy_image
699 Legacy `lsst.afw.image.Image` to compare against.
700 expect_view
701 If not `None`, also assert whether ``image`` shares memory with
702 ``legacy_image`` (i.e. is a view).
703 label
704 Prefix for the failure message.
705 """
706 assert image.bbox == Box.from_legacy(legacy_image.getBBox())
707 if expect_view is not None: 707 ↛ 709line 707 didn't jump to line 709 because the condition on line 707 was always true
708 assert np.may_share_memory(image.array, legacy_image.array) == expect_view
709 if not expect_view: 709 ↛ exitline 709 didn't return from function 'compare_image_to_legacy' because the condition on line 709 was always true
710 assert_values_equal(image.array, legacy_image.array, label=label)
713def compare_mask_to_legacy(
714 mask: Mask,
715 legacy_mask: Any,
716 plane_map: Mapping[str, MaskPlane] | None = None,
717 *,
718 label: str = "mask",
719) -> None:
720 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object.
722 Parameters
723 ----------
724 mask
725 Mask to compare.
726 legacy_mask
727 Legacy `lsst.afw.image.Mask` to compare against.
728 plane_map
729 Mapping from legacy plane name to the new mask plane; defaults to
730 the planes in ``mask.schema``.
731 label
732 Prefix for the failure message; each plane is reported as
733 ``"{label}[<plane>]"``.
734 """
735 assert mask.bbox == Box.from_legacy(legacy_mask.getBBox())
736 if plane_map is None: 736 ↛ 738line 736 didn't jump to line 738 because the condition on line 736 was always true
737 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
738 try:
739 for old_name, new_plane in plane_map.items():
740 assert_values_equal(
741 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
742 mask.get(new_plane.name),
743 label=f"{label}[{old_name}]",
744 )
745 except AssertionError as err:
746 names = {new_plane.name: old_name for old_name, new_plane in plane_map.items()}
747 _note_mask_difference(err, mask, Mask.from_legacy(legacy_mask, plane_map), label=label, names=names)
748 raise
751def compare_masked_image_to_legacy(
752 masked_image: MaskedImage,
753 legacy_masked_image: Any,
754 *,
755 plane_map: Mapping[str, MaskPlane] | None = None,
756 expect_view: bool | None = None,
757 alternates: Mapping[str, Any] | None = None,
758) -> None:
759 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
760 object.
762 Parameters
763 ----------
764 masked_image
765 New image to test.
766 legacy_masked_image
767 Legacy image to test against.
768 plane_map
769 Mapping between new and legacy mask planes.
770 expect_view
771 Whether to test that the image and variance arrays do or do not share
772 memory.
773 alternates
774 A mapping of other versions of one or more (new) components to also
775 check against the legacy versions of those components.
776 """
777 compare_image_to_legacy(
778 masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view, label="image"
779 )
780 compare_mask_to_legacy(
781 masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask"
782 )
783 compare_image_to_legacy(
784 masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view, label="variance"
785 )
786 if alternates: 786 ↛ 787line 786 didn't jump to line 787 because the condition on line 786 was never true
787 if image := alternates.get("image"):
788 compare_image_to_legacy(
789 image, legacy_masked_image.getImage(), expect_view=expect_view, label="image"
790 )
791 if mask := alternates.get("mask"):
792 compare_mask_to_legacy(mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask")
793 if variance := alternates.get("variance"):
794 compare_image_to_legacy(
795 variance, legacy_masked_image.getVariance(), expect_view=expect_view, label="variance"
796 )
799def compare_visit_image_to_legacy(
800 visit_image: VisitImage,
801 legacy_exposure: Any,
802 *,
803 plane_map: Mapping[str, MaskPlane] | None = None,
804 expect_view: bool | None = None,
805 instrument: str,
806 visit: int,
807 detector: int,
808 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
809 alternates: Mapping[str, Any] | None = None,
810 check_photometric_scaling: bool = True,
811) -> None:
812 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
813 object.
815 Parameters
816 ----------
817 visit_image
818 New image to test.
819 legacy_exposure
820 Legacy image to test against.
821 plane_map
822 Mapping between new and legacy mask planes; defaults to the planes
823 used by the legacy (`lsst.afw.image.Exposure`) visit image.
824 expect_view
825 Whether to test that the image and variance arrays do or do not share
826 memory.
827 instrument
828 Expected instrument name.
829 visit
830 Expected visit ID.
831 detector
832 Expected detector ID.
833 applied_legacy_photo_calib
834 Legacy `lsst.afw.image.PhotoCalib` already applied to
835 ``legacy_exposure``, used when comparing photometric scaling.
836 alternates
837 A mapping of other versions of one or more (new) components to also
838 check against the legacy versions of those components.
839 check_photometric_scaling
840 If `False`, skip the photometric-scaling comparison.
841 """
842 if plane_map is None:
843 plane_map = get_legacy_visit_image_mask_planes()
844 compare_masked_image_to_legacy(
845 visit_image,
846 legacy_exposure.getMaskedImage(),
847 plane_map=plane_map,
848 expect_view=expect_view,
849 alternates=alternates,
850 )
851 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
852 with annotate_errors("sky_projection"):
853 compare_sky_projection_to_legacy_wcs(
854 visit_image.sky_projection,
855 legacy_exposure.getWcs(),
856 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
857 visit_image.bbox,
858 )
859 assert visit_image.sky_projection is visit_image.mask.sky_projection
860 assert visit_image.sky_projection is visit_image.variance.sky_projection
861 with annotate_errors("psf"):
862 compare_psf_to_legacy(visit_image.psf, legacy_exposure.getPsf())
863 with annotate_errors("summary_stats"):
864 compare_observation_summary_stats_to_legacy(
865 visit_image.summary_stats, legacy_exposure.info.getSummaryStats()
866 )
867 with annotate_errors("detector"):
868 compare_detector_to_legacy(visit_image.detector, legacy_exposure.getDetector(), is_raw_assembled=True)
869 # Make a tiny box for Field comparisons that need to make arrays; that can
870 # get expensive otherwise.
871 tiny_bbox = detector_bbox.local[2:4, 3:6]
872 with annotate_errors("aperture_corrections"):
873 compare_aperture_corrections_to_legacy(
874 visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
875 )
876 if check_photometric_scaling:
877 with annotate_errors("photometric_scaling"):
878 compare_photo_calib_to_legacy(
879 visit_image.photometric_scaling,
880 legacy_exposure.info.getPhotoCalib(),
881 applied_legacy_photo_calib=applied_legacy_photo_calib,
882 subimage_bbox=tiny_bbox,
883 )
884 if alternates:
885 if (bbox := alternates.get("bbox")) is not None:
886 assert bbox == visit_image.bbox
887 if sky_projection := alternates.get("sky_projection"):
888 with annotate_errors("sky_projection"):
889 compare_sky_projection_to_legacy_wcs(
890 sky_projection,
891 legacy_exposure.getWcs(),
892 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
893 visit_image.bbox,
894 )
895 if psf := alternates.get("psf"):
896 with annotate_errors("psf"):
897 compare_psf_to_legacy(psf, legacy_exposure.getPsf())
898 if summary_stats := alternates.get("summary_stats"):
899 with annotate_errors("summary_stats"):
900 compare_observation_summary_stats_to_legacy(
901 summary_stats, legacy_exposure.info.getSummaryStats()
902 )
903 if detector_obj := alternates.get("detector"):
904 with annotate_errors("detector"):
905 compare_detector_to_legacy(detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True)
906 if obs_info := alternates.get("obs_info"):
907 visitInfo = legacy_exposure.visitInfo
908 assert obs_info.instrument == visitInfo.getInstrumentLabel()
909 if aperture_corrections := alternates.get("aperture_corrections"):
910 with annotate_errors("aperture_corrections"):
911 compare_aperture_corrections_to_legacy(
912 aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
913 )
914 if check_photometric_scaling:
915 if (photometric_scaling := alternates.get("photometric_scaling", ...)) is not ...:
916 with annotate_errors("photometric_scaling"):
917 compare_photo_calib_to_legacy(
918 photometric_scaling,
919 legacy_exposure.info.getPhotoCalib(),
920 applied_legacy_photo_calib=applied_legacy_photo_calib,
921 subimage_bbox=tiny_bbox,
922 )
925def compare_photo_calib_to_legacy(
926 photometric_scaling: BaseField | None,
927 legacy_photo_calib: LegacyPhotoCalib,
928 *,
929 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
930 subimage_bbox: Box,
931) -> None:
932 if legacy_photo_calib._isConstant:
933 if legacy_photo_calib.getCalibrationMean() == 1.0:
934 if applied_legacy_photo_calib is None:
935 assert photometric_scaling is None
936 return
937 else:
938 legacy_photo_calib = applied_legacy_photo_calib
939 if legacy_photo_calib._isConstant:
940 assert isinstance(photometric_scaling, ChebyshevField)
941 assert_values_equal(
942 photometric_scaling.coefficients,
943 np.array([[legacy_photo_calib.getCalibrationMean()]]),
944 # The coefficients come from a numerical fit, so this can be
945 # bigger than just round-off error epislon.
946 rtol=1e-5,
947 )
948 else:
949 assert photometric_scaling is not None
950 compare_field_to_legacy(
951 photometric_scaling / legacy_photo_calib.getCalibrationMean(),
952 legacy_photo_calib.computeScaledCalibration(),
953 subimage_bbox,
954 )
957def compare_cell_coadd_to_legacy(
958 cell_coadd: CellCoadd,
959 legacy_cell_coadd: MultipleCellCoadd,
960 *,
961 tract_bbox: Box,
962 plane_map: Mapping[str, MaskPlane] | None = None,
963 alternates: Mapping[str, Any] | None = None,
964 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
965) -> None:
966 """Compare a `.cells.CellCoadd` object to a legacy
967 `lsst.cell_coadds.MultipleCellCoadd` object.
969 Parameters
970 ----------
971 cell_coadd
972 New coadd to test.
973 legacy_cell_coadd
974 Legacy coadd to test against.
975 tract_bbox
976 Bounding box of the full tract.
977 plane_map
978 Mapping between new and legacy mask planes.
979 alternates
980 A mapping of other versions of one or more (new) components to also
981 check against the legacy versions of those components.
982 psf_points
983 Points to use to compare the PSFs.
984 """
985 legacy_stitched = legacy_cell_coadd.stitch(cell_coadd.bbox.to_legacy())
986 compare_image_to_legacy(cell_coadd.image, legacy_stitched.image, expect_view=False)
987 compare_mask_to_legacy(cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map)
988 compare_image_to_legacy(cell_coadd.variance, legacy_stitched.variance, expect_view=False)
989 if legacy_stitched.mask_fractions is not None:
990 compare_image_to_legacy(
991 cell_coadd.mask_fractions["rejected"], legacy_stitched.mask_fractions, expect_view=False
992 )
993 for n in range(legacy_stitched.n_noise_realizations):
994 compare_image_to_legacy(
995 cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False
996 )
997 assert cell_coadd.skymap == legacy_stitched.identifiers.skymap
998 assert cell_coadd.tract == legacy_stitched.identifiers.tract
999 assert cell_coadd.patch.index.x == legacy_stitched.identifiers.patch.x
1000 assert cell_coadd.patch.index.y == legacy_stitched.identifiers.patch.y
1001 assert cell_coadd.band == legacy_stitched.identifiers.band
1002 assert tract_bbox.contains(cell_coadd.patch.outer_bbox)
1003 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox)
1004 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox)
1005 assert cell_coadd.unit == u.Unit(legacy_cell_coadd.common.units.value)
1006 assert cell_coadd.bounds.bbox.contains(cell_coadd.bbox)
1007 assert cell_coadd.grid.bbox.contains(cell_coadd.bbox)
1008 compare_sky_projection_to_legacy_wcs(
1009 cell_coadd.sky_projection,
1010 legacy_cell_coadd.wcs,
1011 TractFrame(
1012 skymap=legacy_cell_coadd.identifiers.skymap,
1013 tract=legacy_cell_coadd.identifiers.tract,
1014 bbox=tract_bbox,
1015 ),
1016 cell_coadd.bbox,
1017 is_fits=True,
1018 )
1019 assert cell_coadd.sky_projection is cell_coadd.mask.sky_projection
1020 assert cell_coadd.sky_projection is cell_coadd.variance.sky_projection
1021 compare_psf_to_legacy(
1022 cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points
1023 )
1024 compare_aperture_corrections_to_legacy(
1025 cell_coadd.aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
1026 )
1027 compare_cell_coadd_provenance_to_legacy(cell_coadd.provenance, legacy_cell_coadd)
1028 if alternates:
1029 if sky_projection := alternates.get("sky_projection"):
1030 compare_sky_projection_to_legacy_wcs(
1031 sky_projection,
1032 legacy_stitched.wcs,
1033 TractFrame(
1034 skymap=legacy_cell_coadd.identifiers.skymap,
1035 tract=legacy_cell_coadd.identifiers.tract,
1036 bbox=tract_bbox,
1037 ),
1038 cell_coadd.bbox,
1039 is_fits=True,
1040 )
1041 if psf := alternates.get("psf"):
1042 compare_psf_to_legacy(psf, legacy_stitched.psf, points=psf_points)
1043 if aperture_corrections := alternates.get("aperture_corrections"):
1044 compare_aperture_corrections_to_legacy(
1045 aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
1046 )
1047 if provenance := alternates.get("provenance"):
1048 compare_cell_coadd_provenance_to_legacy(provenance, legacy_cell_coadd)
1051def compare_cell_coadd_provenance_to_legacy(
1052 provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd
1053) -> None:
1054 """Compare a `.cells.CoaddProvenance` object to a legacy
1055 `lsst.cell_coadds.MultipleCellCoadd` object.
1057 Parameters
1058 ----------
1059 provenance
1060 New provenance object to test.
1061 legacy_cell_coadd
1062 Legacy coadd to test against.
1063 """
1064 from lsst.cell_coadds import ObservationIdentifiers
1066 for legacy_cell in legacy_cell_coadd.cells.values():
1067 cell_index = CellIJ.from_legacy(legacy_cell.identifiers.cell)
1068 prov = provenance[cell_index]
1069 assert prov is not None, f"Legacy cell {cell_index} has inputs, so its provenance is not empty."
1070 legacy_table = astropy.table.Table(
1071 rows=[
1072 [
1073 ids.instrument,
1074 ids.visit,
1075 ids.detector,
1076 ids.day_obs,
1077 ids.physical_filter,
1078 legacy_input.overlaps_center,
1079 legacy_input.overlap_fraction,
1080 legacy_input.weight,
1081 legacy_input.psf_shape.getIxx(),
1082 legacy_input.psf_shape.getIyy(),
1083 legacy_input.psf_shape.getIxy(),
1084 legacy_input.psf_shape_flag,
1085 ]
1086 for ids, legacy_input in legacy_cell.inputs.items()
1087 ],
1088 dtype=[
1089 np.object_,
1090 np.uint64,
1091 np.uint16,
1092 np.uint32,
1093 np.object_,
1094 np.bool_,
1095 np.float64,
1096 np.float64,
1097 np.float64,
1098 np.float64,
1099 np.float64,
1100 np.bool_,
1101 ],
1102 names=[
1103 "instrument",
1104 "visit",
1105 "detector",
1106 "day_obs",
1107 "physical_filter",
1108 "overlaps_center",
1109 "overlap_fraction",
1110 "weight",
1111 "psf_shape_xx",
1112 "psf_shape_yy",
1113 "psf_shape_xy",
1114 "psf_shape_flag",
1115 ],
1116 )
1117 # For a single cell all 'inputs' are also 'contributions'.
1118 assert len(legacy_cell.inputs) == len(prov.inputs)
1119 assert len(legacy_cell.inputs) == len(prov.contributions)
1120 prov.inputs.sort(["instrument", "visit", "detector"])
1121 prov.contributions.sort(["instrument", "visit", "detector"])
1122 legacy_table.sort(["instrument", "visit", "detector"])
1123 assert_values_equal(prov.inputs["instrument"], prov.contributions["instrument"], label="inputs")
1124 assert_values_equal(prov.inputs["visit"], prov.contributions["visit"], label="inputs")
1125 assert_values_equal(prov.inputs["detector"], prov.contributions["detector"], label="inputs")
1126 assert_values_equal(prov.inputs["instrument"], legacy_table["instrument"], label="inputs")
1127 assert_values_equal(prov.inputs["visit"], legacy_table["visit"], label="inputs")
1128 assert_values_equal(prov.inputs["detector"], legacy_table["detector"], label="inputs")
1129 assert_values_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"], label="inputs")
1130 assert_values_equal(prov.inputs["day_obs"], legacy_table["day_obs"], label="inputs")
1131 assert_values_equal(
1132 prov.contributions["overlaps_center"], legacy_table["overlaps_center"], label="inputs"
1133 )
1134 assert_values_equal(
1135 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"], label="inputs"
1136 )
1137 assert_values_equal(prov.contributions["weight"], legacy_table["weight"], label="inputs")
1138 assert_values_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"], label="inputs")
1139 assert_values_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"], label="inputs")
1140 assert_values_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"], label="inputs")
1141 assert_values_equal(
1142 prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"], label="inputs"
1143 )
1144 for row in prov.inputs:
1145 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"})
1146 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key]
1147 assert legacy_polygon == row["polygon"].to_legacy()
1150def compare_psf_to_legacy(
1151 psf: PointSpreadFunction,
1152 legacy_psf: Any,
1153 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
1154 expect_legacy_raise_on_out_of_bounds: bool = False,
1155) -> int:
1156 """Compare a PSF model object to its legacy interface.
1158 Parameters
1159 ----------
1160 psf
1161 Point-spread function to test.
1162 legacy_psf
1163 Legacy `lsst.afw.detection.Psf` instance to compare with.
1164 points
1165 Points to evaluate the PSFs at. If not provided, the intersection of
1166 the PSF bounding boxes are used to generate points on a grid.
1167 expect_legacy_raise_on_out_of_bounds
1168 If `True`, expect ``legacy_psf`` to raise
1169 `lsst.afw.detection.InvalidPsfError` when evaluated at a position
1170 considered out-of-bounds by ``psf``.
1172 Returns
1173 -------
1174 `int`
1175 The number of points actually tested.
1176 """
1177 from lsst.afw.detection import InvalidPsfError
1179 if points is None:
1180 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel)
1181 legacy_points = arrays_to_legacy_points(points.x, points.y)
1182 n_points_tested: int = 0
1183 for p in legacy_points:
1184 if not psf.bounds.contains(x=p.x, y=p.y):
1185 if expect_legacy_raise_on_out_of_bounds:
1186 with pytest.raises(InvalidPsfError):
1187 legacy_psf.computeKernelImage(p)
1188 continue
1189 assert psf.kernel_bbox == Box.from_legacy(legacy_psf.computeKernelBBox(p))
1190 assert psf.compute_kernel_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeKernelImage(p))
1191 assert psf.compute_stellar_bbox(x=p.x, y=p.y) == Box.from_legacy(legacy_psf.computeImageBBox(p))
1192 assert psf.compute_stellar_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeImage(p))
1193 n_points_tested += 1
1194 return n_points_tested
1197def compare_field_to_legacy(
1198 field: BaseField,
1199 legacy_field: Any,
1200 subimage_bbox: Box,
1201) -> None:
1202 """Test a Field object by comparing it to an equivalent
1203 `lsst.afw.math.BoundedField`.
1205 Parameters
1206 ----------
1207 field
1208 Field to test.
1209 legacy_field : ``lsst.afw.math.BoundedField``
1210 Equivalent legacy bounded field.
1211 subimage_bbox
1212 Bounding box for full-image tests.
1213 """
1214 from lsst.afw.math import BoundedField as LegacyBoundedField
1216 assert field.bounds.bbox == Box.from_legacy(legacy_field.getBBox())
1217 # Pixel coordinates to test the numpy array interface with.
1218 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
1219 if not isinstance(field.bounds, Box): 1219 ↛ 1220line 1219 didn't jump to line 1220 because the condition on line 1219 was never true
1220 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y)
1221 pixel_xy = pixel_xy.map(lambda v: v[mask])
1222 try:
1223 assert_values_equal(
1224 field(x=pixel_xy.x, y=pixel_xy.y),
1225 legacy_field.evaluate(pixel_xy.x, pixel_xy.y),
1226 # The field coefficients come from a numerical fit, so this can
1227 # be bigger than just round-off error epsilon.
1228 rtol=1e-5,
1229 )
1230 except AssertionError as err:
1231 err.add_note(f"evaluated at {pixel_xy}")
1232 raise
1233 if not isinstance(legacy_field, LegacyBoundedField): 1233 ↛ 1236line 1233 didn't jump to line 1236 because the condition on line 1233 was never true
1234 # Legacy StitchedApertureCorrection objects are not true BoundedFields
1235 # and don't have addToImage.
1236 return
1237 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
1238 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
1239 assert_images_equal(
1240 field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13
1241 )
1244def compare_aperture_corrections_to_legacy(
1245 aperture_corrections: Mapping[str, BaseField],
1246 legacy_ap_corr_map: Any,
1247 subimage_bbox: Box,
1248) -> None:
1249 """Test an aperture correction `dict` by comparing it to an equivalent
1250 `lsst.afw.image.ApCorrMap`.
1252 Parameters
1253 ----------
1254 aperture_corrections
1255 Dictionary to test.
1256 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
1257 Equivalent legacy aperture correction map.
1258 subimage_bbox
1259 Bounding box for full-image tests.
1260 """
1261 assert aperture_corrections.keys() == set(legacy_ap_corr_map.keys())
1262 for name, field in aperture_corrections.items():
1263 compare_field_to_legacy(field, legacy_ap_corr_map[name], subimage_bbox)
1266def compare_observation_summary_stats_to_legacy(
1267 summary_stats: ObservationSummaryStats,
1268 legacy_summary_stats: Any,
1269) -> None:
1270 """Test an ObservationSummaryStats object by comparing it to an equivalent
1271 `lsst.afw.image.ExposureSummaryStats`.
1273 Parameters
1274 ----------
1275 summary_stats
1276 Struct to test.
1277 legacy_summary_stats : ``lsst.afw.image.ExposureSummaryStats``
1278 Equivalent legacy struct.
1279 """
1280 for field in dataclasses.fields(legacy_summary_stats):
1281 # Always skip version since that field has no meaning in the new
1282 # type.
1283 if field.name == "version":
1284 continue
1285 a = getattr(legacy_summary_stats, field.name)
1286 b = getattr(summary_stats, field.name)
1287 assert_values_equal(a, b, label=field.name)
1290def compare_sky_projection_to_legacy_wcs[F: Frame](
1291 sky_projection: SkyProjection[F],
1292 legacy_wcs: Any,
1293 pixel_frame: F,
1294 subimage_bbox: Box,
1295 is_fits: bool = False,
1296) -> None:
1297 """Test a Projection object by comparing it to an equivalent
1298 `lsst.afw.geom.SkyWcs`.
1300 Parameters
1301 ----------
1302 sky_projection
1303 Projection to test.
1304 legacy_wcs : ``lsst.afw.geom.SkyWcs``
1305 Equivalent legacy WCS.
1306 pixel_frame
1307 Expected pixel frame for the sky_projection.
1308 subimage_bbox
1309 Bounding box of points to generate for tests.
1310 is_fits
1311 Whether this sky_projection is expected to be exactly representable as
1312 a FITS WCS. If `False` it is assumed to have a FITS approximation
1313 attached instead.
1314 """
1315 # Pixel coordinates to test on over the subimage region of interest:
1316 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
1317 # Array indices of those pixel values (subtract off bbox starts):
1318 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
1319 sky_coords = legacy_coords_to_astropy(
1320 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
1321 )
1322 # Test transforming with the Projection itself, which also tests its
1323 # nested Transform and an Astropy High-Level WCS view with no origin
1324 # change.
1325 check_projection(sky_projection, pixel_xy, sky_coords, pixel_frame, pixel_atol=_AST_INV_ATOL)
1326 # Also test the Astropy High-Level WCS view with an origin change to
1327 # array indices.
1328 check_astropy_wcs_interface(
1329 sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
1330 )
1331 if is_fits: 1331 ↛ 1350line 1331 didn't jump to line 1350 because the condition on line 1331 was always true
1332 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
1333 assert fits_wcs is not None
1334 check_astropy_wcs_interface(fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
1335 # Use that FITS approximation to check that we can make a
1336 # Projection from a FITS WCS, too.
1337 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
1338 check_projection(
1339 fits_projection,
1340 subimage_array_xy,
1341 sky_coords,
1342 pixel_frame,
1343 pixel_atol=1e-5,
1344 )
1345 # We want Projections we create from a FITS WCS to be backed by an
1346 # AST FrameSet so we can convert them into legacy
1347 # `lsst.afw.geom.SkyWcs` objects if desired.
1348 assert "Begin FrameSet" in fits_projection.show()
1349 else:
1350 assert sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False) is None
1351 # The legacy SkyWcs should instead have a FITS approximation
1352 # attached; run the same tests on that.
1353 assert sky_projection.fits_approximation is not None
1354 compare_sky_projection_to_legacy_wcs(
1355 sky_projection.fits_approximation,
1356 legacy_wcs.getFitsApproximation(),
1357 pixel_frame,
1358 subimage_bbox,
1359 is_fits=True,
1360 )
1363def check_transform[I: Frame, O: Frame](
1364 transform: Transform[I, O],
1365 input_xy: XY[np.ndarray],
1366 output_xy: XY[np.ndarray],
1367 in_frame: Frame,
1368 out_frame: Frame,
1369 *,
1370 check_inverted: bool = True,
1371 in_atol: u.Quantity | float | None = None,
1372 out_atol: u.Quantity | float | None = None,
1373 angle_wrapped: bool = False,
1374) -> None:
1375 """Test Transform against known arrays of input and output points.
1377 Parameters
1378 ----------
1379 transform
1380 Transform to test.
1381 input_xy
1382 Arrays of input points.
1383 output_xy
1384 Arrays of output points.
1385 in_frame
1386 Expected input frame.
1387 out_frame
1388 Expected output frame.
1389 check_inverted
1390 If `True`, recurse (once) to test the inverse transform.
1391 in_atol
1392 Expected absolute precision of input points.
1393 out_atol
1394 Expected absolute precision of output points.
1395 angle_wrapped
1396 If `True`, treat comparisons on sky-coordinate axes (as determined
1397 by the frames) as angles modulo 2 pi, so differences that cross
1398 zero or 360 degrees compare correctly.
1399 """
1400 assert transform.in_frame == in_frame
1401 assert transform.out_frame == out_frame
1402 in_atol = 0.0 if in_atol is None else in_atol
1403 out_atol = 0.0 if out_atol is None else out_atol
1404 in_atol_v = in_atol.to_value(in_frame.unit) if isinstance(in_atol, u.Quantity) else float(in_atol)
1405 out_atol_v = out_atol.to_value(out_frame.unit) if isinstance(out_atol, u.Quantity) else float(out_atol)
1406 out_wrapped = angle_wrapped and out_frame == SkyFrame.ICRS
1407 in_wrapped = angle_wrapped and in_frame == SkyFrame.ICRS
1409 def assert_close(
1410 test: np.ndarray | float | u.Quantity,
1411 expected: np.ndarray | float | u.Quantity,
1412 atol: float | u.Quantity,
1413 *,
1414 wrapped: bool,
1415 ) -> None:
1416 """Compare one axis, wrap-aware if that axis is a sky angle.
1418 Parameters
1419 ----------
1420 test
1421 Values to test.
1422 expected
1423 Expected values.
1424 atol
1425 Absolute tolerance, in the units the values are expressed in.
1426 wrapped
1427 If `True`, treat the values as angles modulo 2 pi (radians).
1428 """
1429 if not wrapped:
1430 assert_values_equal(test, expected, atol=atol, rtol=_AST_RTOL)
1431 return
1432 test_values = np.asarray(getattr(test, "value", test))
1433 expected_values = np.asarray(getattr(expected, "value", expected))
1434 atol_rad = atol.to_value(u.rad) if isinstance(atol, u.Quantity) else float(atol)
1435 assert_values_equal(
1436 test_values * u.rad,
1437 expected_values * u.rad,
1438 atol=atol_rad * u.rad,
1439 equal_nan=False,
1440 wrap_angles=True,
1441 )
1443 # Test array interfaces.
1444 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
1445 assert_close(test_output_xy.x, output_xy.x, out_atol_v, wrapped=out_wrapped)
1446 assert_close(test_output_xy.y, output_xy.y, out_atol_v, wrapped=out_wrapped)
1447 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
1448 assert_close(test_input_xy.x, input_xy.x, in_atol_v, wrapped=in_wrapped)
1449 assert_close(test_input_xy.y, input_xy.y, in_atol_v, wrapped=in_wrapped)
1450 # Test scalar interfaces with numpy scalars.
1451 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
1452 forward_xy = transform.apply_forward(x=input_x, y=input_y)
1453 assert_close(forward_xy.x, output_x, out_atol_v, wrapped=out_wrapped)
1454 assert_close(forward_xy.y, output_y, out_atol_v, wrapped=out_wrapped)
1455 inverse_xy = transform.apply_inverse(x=output_x, y=output_y)
1456 assert_close(inverse_xy.x, input_x, in_atol_v, wrapped=in_wrapped)
1457 assert_close(inverse_xy.y, input_y, in_atol_v, wrapped=in_wrapped)
1458 # Test quantity array interfaces.
1459 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
1460 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
1461 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
1462 assert_close(test_output_q_xy.x, output_q_xy.x, out_atol, wrapped=out_wrapped)
1463 assert_close(test_output_q_xy.y, output_q_xy.y, out_atol, wrapped=out_wrapped)
1464 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
1465 assert_close(test_input_q_xy.x, input_q_xy.x, in_atol, wrapped=in_wrapped)
1466 assert_close(test_input_q_xy.y, input_q_xy.y, in_atol, wrapped=in_wrapped)
1467 # Test quantity scalar interfaces.
1468 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
1469 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
1470 ):
1471 forward_q_xy = transform.apply_forward_q(x=input_q_x, y=input_q_y)
1472 assert_close(forward_q_xy.x, output_q_x, out_atol, wrapped=out_wrapped)
1473 assert_close(forward_q_xy.y, output_q_y, out_atol, wrapped=out_wrapped)
1474 inverse_q_xy = transform.apply_inverse_q(x=output_q_x, y=output_q_y)
1475 assert_close(inverse_q_xy.x, input_q_x, in_atol, wrapped=in_wrapped)
1476 assert_close(inverse_q_xy.y, input_q_y, in_atol, wrapped=in_wrapped)
1477 if check_inverted:
1478 # Test the inverse transform.
1479 check_transform(
1480 transform.inverted(),
1481 output_xy,
1482 input_xy,
1483 out_frame,
1484 in_frame,
1485 check_inverted=False,
1486 out_atol=in_atol,
1487 in_atol=out_atol,
1488 angle_wrapped=angle_wrapped,
1489 )
1492def check_projection[P: Frame](
1493 sky_projection: SkyProjection[P],
1494 pixel_xy: XY[np.ndarray],
1495 sky_coords: SkyCoord,
1496 pixel_frame: Frame,
1497 *,
1498 pixel_atol: float | None = None,
1499 sky_atol: u.Quantity | None = None,
1500) -> None:
1501 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1502 coordinates.
1504 Parameters
1505 ----------
1506 sky_projection
1507 Projection to test.
1508 pixel_xy
1509 Arrays of pixel coordinates.
1510 sky_coords
1511 Corresponding sky coordinates.
1512 pixel_frame
1513 Expected pixel frame.
1514 pixel_atol
1515 Expected absolute precision of pixel points.
1516 sky_atol
1517 Expected absolute precision of sky coordinates, as a maximum
1518 great-circle separation.
1519 """
1520 assert sky_projection.pixel_frame == pixel_frame
1521 assert sky_projection.sky_frame == SkyFrame.ICRS
1522 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol)
1523 pixel_atol = 0.0 if pixel_atol is None else pixel_atol
1524 pixel_atol_q = pixel_atol * u.pix
1525 # Test array interfaces.
1526 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1527 assert_values_equal(test_pixel_xy.x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL)
1528 assert_values_equal(test_pixel_xy.y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL)
1529 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1530 assert_values_equal(test_sky_astropy, sky_coords, atol=sky_atol, label="pixel_to_sky")
1531 # Test scalar interfaces.
1532 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1533 assert_values_equal(
1534 sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol, rtol=_AST_RTOL
1535 )
1536 assert_values_equal(
1537 sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol, rtol=_AST_RTOL
1538 )
1539 test_sky_single = sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y)
1540 assert_values_equal(test_sky_single, sky_single, atol=sky_atol, label="pixel_to_sky")
1541 # Test the underlying Transform object.
1542 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1543 check_transform(
1544 sky_projection.pixel_to_sky_transform,
1545 pixel_xy,
1546 sky_xy,
1547 pixel_frame,
1548 SkyFrame.ICRS,
1549 check_inverted=False,
1550 in_atol=pixel_atol_q,
1551 out_atol=sky_atol,
1552 angle_wrapped=True,
1553 )
1554 check_transform(
1555 sky_projection.sky_to_pixel_transform,
1556 sky_xy,
1557 pixel_xy,
1558 SkyFrame.ICRS,
1559 pixel_frame,
1560 check_inverted=False,
1561 in_atol=sky_atol,
1562 out_atol=pixel_atol_q,
1563 angle_wrapped=True,
1564 )
1565 # Test the Astropy interface adapter.
1566 check_astropy_wcs_interface(
1567 sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1568 )
1571def assert_sky_projections_equal(
1572 a: SkyProjection[Any] | None,
1573 b: SkyProjection[Any] | None,
1574 expect_identity: bool | None = None,
1575) -> None:
1576 """Test that two `.SkyProjection` instances are equivalent.
1578 Parameters
1579 ----------
1580 a
1581 First sky projection to compare.
1582 b
1583 Second sky projection to compare.
1584 expect_identity
1585 If not `None`, assert whether ``a`` and ``b`` are the same object.
1586 """
1587 if a is None and b is None:
1588 return
1589 assert a is not None
1590 assert b is not None
1591 match expect_identity:
1592 case True: 1592 ↛ 1593line 1592 didn't jump to line 1593 because the pattern on line 1592 never matched
1593 assert a is b
1594 return
1595 case False:
1596 assert a is not b
1597 case None if a is b:
1598 return
1599 assert a == b
1602def check_astropy_wcs_interface(
1603 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1604 pixel_xy: XY[np.ndarray],
1605 sky_coords: SkyCoord,
1606 *,
1607 pixel_atol: float | None = None,
1608 sky_atol: u.Quantity | None = None,
1609) -> None:
1610 """Test an Astropy WCS instance against known arrays of pixel and
1611 sky coordinates.
1613 Parameters
1614 ----------
1615 wcs
1616 Astropy WCS object to test.
1617 pixel_xy
1618 Arrays of pixel coordinates.
1619 sky_coords
1620 Corresponding sky coordinates.
1621 pixel_atol
1622 Expected absolute precision of pixel points.
1623 sky_atol
1624 Expected absolute precision of sky coordinates, as a maximum
1625 great-circle separation.
1626 """
1627 pixel_atol = 0.0 if pixel_atol is None else pixel_atol
1628 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol)
1629 test_x, test_y = wcs.world_to_pixel(sky_coords)
1630 assert_values_equal(test_x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL)
1631 assert_values_equal(test_y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL)
1632 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
1633 assert_values_equal(test_sky_coords, sky_coords, atol=sky_atol, label="pixel_to_world")
1636def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
1637 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array.
1639 Parameters
1640 ----------
1641 legacy_points
1642 Legacy ``lsst.geom.Point2D`` objects to convert.
1643 """
1644 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
1647def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
1648 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
1649 coordinate object.
1651 Parameters
1652 ----------
1653 legacy_coords
1654 Legacy ``lsst.geom.SpherePoint`` objects to convert.
1655 """
1656 return SkyCoord(
1657 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
1658 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
1659 )
1662def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
1663 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``.
1665 Parameters
1666 ----------
1667 x
1668 X coordinates of the points.
1669 y
1670 Y coordinates of the points.
1671 """
1672 from lsst.geom import Point2D
1674 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
1677def compare_amplifier_to_legacy(
1678 amplifier: Amplifier,
1679 legacy_amplifier: Any,
1680 *,
1681 is_raw_assembled: bool,
1682 expect_nominal_calibrations: bool = True,
1683) -> None:
1684 """Compare an `~.cameras.Amplifier` to a legacy
1685 `lsst.afw.cameraGeom.Amplifier`.
1687 Parameters
1688 ----------
1689 amplifier
1690 Amplifier to compare.
1691 legacy_amplifier
1692 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against.
1693 is_raw_assembled
1694 Whether the raw geometry is expected to be the assembled-raw
1695 geometry (`True`) or the unassembled-raw geometry (`False`).
1696 expect_nominal_calibrations
1697 Whether the amplifier is expected to carry nominal calibrations.
1698 """
1699 assert legacy_amplifier.getName() == amplifier.name
1700 assert Box.from_legacy(legacy_amplifier.getBBox()) == amplifier.bbox
1701 if is_raw_assembled: 1701 ↛ 1704line 1701 didn't jump to line 1704 because the condition on line 1701 was always true
1702 raw_geom = amplifier.assembled_raw_geometry
1703 else:
1704 raw_geom = amplifier.unassembled_raw_geometry
1705 assert raw_geom is not None
1706 assert ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()) == raw_geom.readout_corner
1707 assert Box.from_legacy(legacy_amplifier.getRawBBox()) == raw_geom.bbox
1708 assert Box.from_legacy(legacy_amplifier.getRawDataBBox()) == raw_geom.data_bbox
1709 assert legacy_amplifier.getRawFlipX() == raw_geom.flip_x
1710 assert legacy_amplifier.getRawFlipY() == raw_geom.flip_y
1711 assert legacy_amplifier.getRawXYOffset().getX() == raw_geom.x_offset
1712 assert legacy_amplifier.getRawXYOffset().getY() == raw_geom.y_offset
1713 assert (
1714 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()) == raw_geom.horizontal_overscan_bbox
1715 )
1716 assert Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()) == raw_geom.vertical_overscan_bbox
1717 if raw_geom.horizontal_prescan_bbox is None:
1718 # `Box` requires a positive size, so an amplifier with no prescan
1719 # region carries `None` here.
1720 assert legacy_amplifier.getRawPrescanBBox().isEmpty()
1721 else:
1722 assert Box.from_legacy(legacy_amplifier.getRawPrescanBBox()) == raw_geom.horizontal_prescan_bbox
1723 if expect_nominal_calibrations: 1723 ↛ exitline 1723 didn't return from function 'compare_amplifier_to_legacy' because the condition on line 1723 was always true
1724 assert amplifier.nominal_calibrations is not None
1725 assert_values_equal(
1726 legacy_amplifier.getGain(),
1727 amplifier.nominal_calibrations.gain,
1728 equal_nan=True,
1729 label="gain",
1730 )
1731 assert_values_equal(
1732 legacy_amplifier.getReadNoise(),
1733 amplifier.nominal_calibrations.read_noise,
1734 equal_nan=True,
1735 label="read_noise",
1736 )
1737 assert_values_equal(
1738 legacy_amplifier.getSaturation(),
1739 amplifier.nominal_calibrations.saturation,
1740 equal_nan=True,
1741 label="saturation",
1742 )
1743 assert_values_equal(
1744 legacy_amplifier.getSuspectLevel(),
1745 amplifier.nominal_calibrations.suspect_level,
1746 equal_nan=True,
1747 label="suspect_level",
1748 )
1749 assert_values_equal(
1750 legacy_amplifier.getLinearityCoeffs(),
1751 amplifier.nominal_calibrations.linearity_coefficients,
1752 label="linearity_coefficients",
1753 )
1754 assert legacy_amplifier.getLinearityType() == amplifier.nominal_calibrations.linearity_type
1757def compare_detector_to_legacy(
1758 detector: Detector,
1759 legacy_detector: Any,
1760 *,
1761 is_raw_assembled: bool,
1762 expect_nominal_calibrations: bool = True,
1763) -> None:
1764 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`.
1766 Parameters
1767 ----------
1768 detector
1769 Detector to compare.
1770 legacy_detector
1771 Legacy `lsst.afw.cameraGeom.Detector` to compare against.
1772 is_raw_assembled
1773 Whether the raw geometry is expected to be the assembled-raw
1774 geometry (`True`) or the unassembled-raw geometry (`False`).
1775 expect_nominal_calibrations
1776 Whether the detector's amplifiers are expected to carry nominal
1777 calibrations.
1778 """
1779 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
1781 assert legacy_detector.getName() == detector.name
1782 assert legacy_detector.getId() == detector.id
1783 assert DetectorType.from_legacy(legacy_detector.getType()) == detector.type
1784 assert Box.from_legacy(legacy_detector.getBBox()) == detector.bbox
1785 assert legacy_detector.getSerial() == detector.serial
1786 legacy_orientation = legacy_detector.getOrientation()
1787 assert legacy_orientation.getFpPosition3().getX() == detector.orientation.focal_plane_x
1788 assert legacy_orientation.getFpPosition3().getY() == detector.orientation.focal_plane_y
1789 assert legacy_orientation.getFpPosition3().getZ() == detector.orientation.focal_plane_z
1790 assert legacy_orientation.getReferencePoint().getX() == detector.orientation.pixel_reference_x
1791 assert legacy_orientation.getReferencePoint().getY() == detector.orientation.pixel_reference_y
1792 assert legacy_orientation.getYaw().asRadians() == detector.orientation.yaw.to_value(u.rad)
1793 assert legacy_orientation.getPitch().asRadians() == detector.orientation.pitch.to_value(u.rad)
1794 assert legacy_orientation.getRoll().asRadians() == detector.orientation.roll.to_value(u.rad)
1795 assert legacy_detector.getPixelSize().getX() == detector.pixel_size
1796 assert legacy_detector.getPhysicalType() == detector.physical_type
1797 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
1798 compare_amplifier_to_legacy(
1799 amplifier,
1800 legacy_amplifier,
1801 is_raw_assembled=is_raw_assembled,
1802 expect_nominal_calibrations=expect_nominal_calibrations,
1803 )
1804 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
1805 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
1806 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
1807 check_transform(
1808 detector.to_focal_plane,
1809 pixel_xy,
1810 legacy_points_to_xy_array(fp_legacy_points),
1811 detector.frame,
1812 detector.to_focal_plane.out_frame,
1813 in_atol=1e-9 * u.pix,
1814 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
1815 )
1816 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
1817 check_transform(
1818 detector.to_field_angle,
1819 pixel_xy,
1820 legacy_points_to_xy_array(fa_legacy_points),
1821 detector.frame,
1822 detector.to_field_angle.out_frame,
1823 in_atol=1e-9 * u.pix,
1824 out_atol=1e-7 * u.arcsec,
1825 )
1828def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
1829 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
1831 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
1832 classes. Importing this module already imports every ``lsst.images``
1833 module that defines a subclass, so the tree is fully populated by the time
1834 this is called.
1836 This discovery is deliberately separate from
1837 `check_archive_tree_class_invariants` so that the per-class check stays
1838 usable on a single class even if this metaprogramming is removed later.
1839 """
1840 seen: set[type] = set()
1841 stack: list[type] = [ArchiveTree]
1842 while stack:
1843 kls = stack.pop()
1844 sub: type
1845 for sub in kls.__subclasses__():
1846 if sub in seen: 1846 ↛ 1847line 1846 didn't jump to line 1847 because the condition on line 1846 was never true
1847 continue
1848 seen.add(sub)
1849 stack.append(sub)
1850 if not getattr(sub, "__abstractmethods__", None): 1850 ↛ 1845line 1850 didn't jump to line 1845 because the condition on line 1850 was always true
1851 yield sub
1854def check_bounds_contains_broadcasting(bounds: Bounds) -> None:
1855 """Verify that `~lsst.images.Bounds.contains` accepts array-like inputs.
1857 Uses the scalar overload as the reference and checks that 1-D arrays,
1858 list inputs, mixed scalar-plus-array inputs, 2-D broadcast inputs, and
1859 `.XY` / `.YX` positional-argument forms all produce results consistent
1860 with calling the scalar overload on each ``(x, y)`` pair individually.
1862 Parameters
1863 ----------
1864 bounds
1865 The `~lsst.images.Bounds` implementation to exercise.
1866 """
1867 bbox = bounds.bbox
1868 # One point outside each boundary, one on each boundary, one in the
1869 # interior. The boundary points (start and stop) straddle inside/outside
1870 # (start is inside, stop is outside) so we are guaranteed a mix of True
1871 # and False results without hard-coding expected values.
1872 ys = np.array(
1873 [
1874 bbox.y.start - 1,
1875 bbox.y.start,
1876 (bbox.y.start + bbox.y.stop) // 2,
1877 bbox.y.stop,
1878 ]
1879 )
1880 xs = np.array(
1881 [
1882 bbox.x.start - 1,
1883 bbox.x.start,
1884 (bbox.x.start + bbox.x.stop) // 2,
1885 bbox.x.stop,
1886 ]
1887 )
1888 # 2-D reference: expected[i, j] == contains(x=xs[j], y=ys[i]).
1889 expected = np.array(
1890 [[bounds.contains(x=int(xi), y=int(yi)) for xi in xs] for yi in ys],
1891 dtype=bool,
1892 )
1893 # 1-D ndarray (diagonal pairs).
1894 assert_values_equal(
1895 bounds.contains(x=xs, y=ys),
1896 np.diagonal(expected),
1897 label="bounds.contains",
1898 )
1899 # list inputs (array-like).
1900 assert_values_equal(
1901 bounds.contains(x=xs.tolist(), y=ys.tolist()),
1902 np.diagonal(expected),
1903 label="bounds.contains",
1904 )
1905 # Mixed: scalar y, 1-D array x — produces a 1-D result.
1906 fixed_yi = 1 # index into ys; ys[1] == bbox.y.start (on the boundary)
1907 assert_values_equal(
1908 bounds.contains(x=xs, y=int(ys[fixed_yi])),
1909 expected[fixed_yi, :],
1910 label="bounds.contains",
1911 )
1912 # Mixed: 1-D array y, scalar x.
1913 fixed_xi = 1 # index into xs; xs[1] == bbox.x.start (on the boundary)
1914 assert_values_equal(
1915 bounds.contains(x=int(xs[fixed_xi]), y=ys),
1916 expected[:, fixed_xi],
1917 label="bounds.contains",
1918 )
1919 # Float scalars: results must match the int-scalar reference.
1920 assert bounds.contains(x=float(xs[fixed_xi]), y=float(ys[fixed_yi])) == expected[fixed_yi, fixed_xi]
1921 # XY / YX scalar: results must match the keyword-scalar reference.
1922 assert bounds.contains(XY(x=int(xs[fixed_xi]), y=int(ys[fixed_yi]))) == expected[fixed_yi, fixed_xi]
1923 assert bounds.contains(YX(y=int(ys[fixed_yi]), x=int(xs[fixed_xi]))) == expected[fixed_yi, fixed_xi]
1924 # XY / YX array: results must match the 1-D ndarray reference.
1925 assert_values_equal(bounds.contains(XY(x=xs, y=ys)), np.diagonal(expected), label="bounds.contains")
1926 assert_values_equal(bounds.contains(YX(y=ys, x=xs)), np.diagonal(expected), label="bounds.contains")
1927 # 2-D broadcast: y shape (N, 1) × x shape (1, M) → (N, M).
1928 assert_values_equal(
1929 bounds.contains(x=xs.reshape(1, -1), y=ys.reshape(-1, 1)),
1930 expected,
1931 label="bounds.contains",
1932 )
1935def check_archive_tree_class_invariants(tree_type: type[ArchiveTree]) -> None:
1936 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
1937 well-formed schema-version constants and an in-memory type.
1939 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
1940 ``PUBLIC_TYPE`` are present and well-typed, that the version is
1941 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
1942 schema major.
1944 Parameters
1945 ----------
1946 tree_type
1947 The concrete `.serialization.ArchiveTree` subclass to check.
1948 """
1949 assert hasattr(tree_type, "SCHEMA_NAME"), f"{tree_type.__name__} lacks SCHEMA_NAME"
1950 assert hasattr(tree_type, "SCHEMA_VERSION"), f"{tree_type.__name__} lacks SCHEMA_VERSION"
1951 assert hasattr(tree_type, "MIN_READ_VERSION"), f"{tree_type.__name__} lacks MIN_READ_VERSION"
1952 assert hasattr(tree_type, "PUBLIC_TYPE"), f"{tree_type.__name__} lacks PUBLIC_TYPE"
1953 assert isinstance(tree_type.SCHEMA_NAME, str)
1954 assert len(tree_type.SCHEMA_NAME) > 0
1955 # Allow an optional PEP 440 development-release suffix (e.g. 1.0.0.dev0)
1956 # for schemas still in development.
1957 assert re.fullmatch(r"^\d+\.\d+\.\d+(\.dev\d+)?$", tree_type.SCHEMA_VERSION)
1958 assert isinstance(tree_type.MIN_READ_VERSION, int)
1959 assert tree_type.MIN_READ_VERSION >= 1
1960 assert isinstance(tree_type.PUBLIC_TYPE, type)
1961 major = int(tree_type.SCHEMA_VERSION.split(".")[0])
1962 assert tree_type.MIN_READ_VERSION <= major
1965def _angle_to_astropy(angle: Any) -> u.Quantity:
1966 """Convert an `lsst.geom.Angle` to an Astropy angular quantity."""
1967 return u.Quantity(angle.asRadians(), u.rad)
1970def _observatory_to_location(observatory: Any) -> coords.EarthLocation:
1971 """Convert an `lsst.afw.coord.Observatory` to an Astropy
1972 `~astropy.coordinates.EarthLocation`.
1973 """
1974 return coords.EarthLocation.from_geodetic(
1975 lon=observatory.getLongitude().asRadians() * u.rad,
1976 lat=observatory.getLatitude().asRadians() * u.rad,
1977 height=observatory.getElevation() * u.m,
1978 )
1981def assert_obs_metadata_fields_equal(actual: Any, expected: Any, *, label: str = "") -> None:
1982 """Assert that two ObservationInfo or VisitInfo field values are equal.
1984 Dispatches on the type of ``expected`` so that values affected by afw's
1985 float64 date storage, unit conversions, and angle wrapping are compared
1986 with appropriate tolerances; anything else is compared with ``==``.
1988 Parameters
1989 ----------
1990 actual : `object`
1991 The field value to check.
1992 expected : `object`
1993 The field value to check it against.
1994 label : `str`, optional
1995 Prefix for the failure message.
1996 """
1997 from lsst.afw.coord import Observatory, Weather
1998 from lsst.daf.base import DateTime
1999 from lsst.geom import Angle
2001 if isinstance(expected, DateTime):
2002 assert_obs_metadata_fields_equal(actual.toAstropy(), expected.toAstropy(), label=label)
2003 elif isinstance(expected, astropy.time.Time):
2004 # afw stores VisitInfo dates as float64 MJDs, which limits precision.
2005 assert_values_equal(actual, expected, atol=1e-9 * u.s, label=label)
2006 elif isinstance(expected, coords.AltAz):
2007 assert_values_equal(actual.az, expected.az, wrap_angles=True, label=f"{label}.az")
2008 assert_values_equal(actual.alt, expected.alt, atol=1e-3 * u.arcsec, label=f"{label}.alt")
2009 elif isinstance(expected, coords.SkyCoord):
2010 assert_values_equal(actual, expected, atol=1e-3 * u.arcsec, label=label)
2011 elif isinstance(expected, coords.EarthLocation):
2012 assert_values_equal(
2013 u.Quantity([actual.x, actual.y, actual.z]),
2014 u.Quantity([expected.x, expected.y, expected.z]),
2015 rtol=1e-9,
2016 label=label,
2017 )
2018 elif isinstance(expected, Observatory):
2019 assert_obs_metadata_fields_equal(
2020 _observatory_to_location(actual), _observatory_to_location(expected), label=label
2021 )
2022 elif isinstance(expected, Weather):
2023 assert_obs_metadata_fields_equal(
2024 u.Quantity(actual.getAirTemperature(), u.deg_C),
2025 u.Quantity(expected.getAirTemperature(), u.deg_C),
2026 label=f"{label}.temperature",
2027 )
2028 assert_obs_metadata_fields_equal(
2029 actual.getAirPressure(), expected.getAirPressure(), label=f"{label}.pressure"
2030 )
2031 assert_obs_metadata_fields_equal(
2032 actual.getHumidity(), expected.getHumidity(), label=f"{label}.humidity"
2033 )
2034 elif isinstance(expected, Angle):
2035 assert_values_equal(
2036 _angle_to_astropy(actual),
2037 _angle_to_astropy(expected),
2038 atol=1e-12 * u.rad,
2039 wrap_angles=True,
2040 label=label,
2041 )
2042 elif (
2043 isinstance(expected, u.Quantity)
2044 and isinstance(actual, u.Quantity)
2045 and expected.unit.physical_type == "temperature"
2046 and actual.unit.physical_type == "temperature"
2047 ):
2048 assert_values_equal(
2049 actual.to(u.K, equivalencies=u.temperature()),
2050 expected.to(u.K, equivalencies=u.temperature()),
2051 atol=1e-6 * u.K,
2052 label=label,
2053 )
2054 elif isinstance(expected, u.Quantity | float):
2055 assert_values_equal(actual, expected, rtol=1e-9, label=label)
2056 elif actual != expected:
2057 raise AssertionError(f"{label}: {actual!r} != {expected!r}")