Coverage for python/lsst/images/tests/_checks.py: 61%
778 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-10 09:42 +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 "cell_center_points",
27 "check_archive_tree_class_invariants",
28 "check_astropy_wcs_interface",
29 "check_bounds_contains_broadcasting",
30 "check_projection",
31 "check_transform",
32 "compare_amplifier_to_legacy",
33 "compare_aperture_corrections_to_legacy",
34 "compare_cell_coadd_provenance_to_legacy_lazy",
35 "compare_cell_coadd_to_legacy",
36 "compare_cell_coadd_to_legacy_lazy",
37 "compare_detector_to_legacy",
38 "compare_field_to_legacy",
39 "compare_image_to_legacy",
40 "compare_mask_to_legacy",
41 "compare_masked_image_to_legacy",
42 "compare_observation_summary_stats_to_legacy",
43 "compare_photo_calib_to_legacy",
44 "compare_psf_to_legacy",
45 "compare_sky_projection_to_legacy_wcs",
46 "compare_visit_image_to_legacy",
47 "iter_concrete_archive_tree_subclasses",
48 "legacy_coords_to_astropy",
49 "legacy_points_to_xy_array",
50)
52import dataclasses
53import re
54from abc import ABC, abstractmethod
55from collections.abc import Generator, Iterator, Mapping
56from contextlib import contextmanager
57from typing import TYPE_CHECKING, Any, Literal, NoReturn, Self, cast
59import astropy.coordinates as coords
60import astropy.time
61import astropy.units as u
62import astropy.wcs.wcsapi
63import numpy as np
64import pytest
65from astropy.coordinates import SkyCoord
67from .._geom import XY, YX, Bounds, Box
68from .._image import Image
69from .._mask import (
70 Mask,
71 MaskPlane,
72 get_legacy_non_cell_coadd_mask_planes,
73 get_legacy_visit_image_mask_planes,
74)
75from .._masked_image import MaskedImage
76from .._observation_summary_stats import ObservationSummaryStats
77from .._transforms import DetectorFrame, Frame, SkyFrame, SkyProjection, TractFrame, Transform
78from .._visit_image import VisitImage
79from ..cameras import Amplifier, Detector, DetectorType, ReadoutCorner
80from ..cells import CellCoadd, CellGridBounds, CellIJ, CoaddProvenance
81from ..fields import BaseField, ChebyshevField
82from ..psfs import PointSpreadFunction
83from ..serialization import ArchiveTree
85if TYPE_CHECKING:
86 try:
87 from lsst.cell_coadds import MultipleCellCoadd
88 except ImportError:
89 type MultipleCellCoadd = Any # type: ignore[no-redef]
90 try:
91 from lsst.afw.image import PhotoCalib as LegacyPhotoCalib
92 except ImportError:
93 type LegacyPhotoCalib = Any # type: ignore[no-redef]
96# The sky-to-pixel direction of a SkyProjection is a numerical AST inversion
97# of the celestial projection, so it carries a small error floor (~1e-9 px).
98# A relative-only tolerance collapses to exact equality at pixel coordinate
99# zero, so comparisons that can include the origin need a small absolute pixel
100# tolerance to avoid reporting that numerical noise as a difference.
101_AST_INV_ATOL = 1e-7
103# Relative tolerance for comparisons of AST-backed transforms and projections.
104_AST_RTOL = 1e-9
106# Default absolute tolerance for great-circle sky-coordinate comparisons, in
107# the same order of magnitude as _AST_INV_ATOL (which is in pixels). A
108# nonzero default is important because sky coordinates can be near zero
109# (e.g. RA/Dec ~ 0), where a relative tolerance provides no error floor.
110_AST_SKY_ATOL = 1e-7 * u.deg
113class _BaseCmpArg(ABC):
114 """A helper base clas for `assert_values_equal` whose subclasses represent
115 different types of comparisons.
117 Parameters
118 ----------
119 value
120 The coerced operand; type is subclass-dependent..
121 """
123 def __init__(self, value: Any):
124 self.value = value
126 @property
127 @abstractmethod
128 def unit(self) -> u.UnitBase:
129 raise NotImplementedError()
131 def coerce_atol(self, atol: float | u.Quantity | astropy.time.TimeDelta) -> float:
132 try:
133 if isinstance(atol, astropy.time.TimeDelta):
134 return atol.to_value(self.unit)
135 if atol != 0.0:
136 return u.Quantity(atol).to_value(self.unit)
137 return 0.0
138 except u.UnitConversionError as err:
139 err.add_note("Invalid atol for these operands.")
140 raise
142 def get_nan_mask(self) -> np.ndarray:
143 return np.isnan(np.asarray(self.value))
145 @abstractmethod
146 def get_separation(self, other: Self) -> np.ndarray:
147 raise NotImplementedError
149 @abstractmethod
150 def find_mismatches(
151 self,
152 other: Self,
153 *,
154 rtol: float,
155 atol_value: float,
156 equal_nan: bool,
157 ) -> np.ndarray | None:
158 """Search for element-wise mismatches against ``other``.
160 Each comparison family implements its own tolerance semantics.
162 Parameters
163 ----------
164 other
165 The other comparison argument.
166 rtol
167 Relative tolerance.
168 atol_value
169 Absolute tolerance, already coerced to this comparison's unit.
170 equal_nan
171 If `True`, treat NaN as equal to NaN wherever the operand can
172 hold NaN.
174 Returns
175 -------
176 mismatch : `numpy.ndarray` | None
177 Boolean mask marking differing elements, or `None` if all
178 elements agree within tolerance.
179 """
180 raise NotImplementedError()
182 def report_mismatches(
183 self,
184 other: Self,
185 a: Any,
186 b: Any,
187 mismatch: np.ndarray,
188 label: str,
189 ) -> NoReturn:
190 """Raise an `AssertionError` reporting a computed mismatch.
192 Parameters
193 ----------
194 other : _BaseCmpArg
195 The other comparison argument.
196 a
197 First original operand; used for the scalar failure message.
198 b
199 Second original operand; used for the scalar failure message.
200 mismatch
201 Boolean mask marking differing elements, from `find_mismatches`.
202 label
203 Prefix prepended to the failure message.
204 """
205 n_mismatch = int(np.count_nonzero(mismatch))
206 prefix = f"{label}: " if label else ""
207 separation = self.get_separation(other)
208 if mismatch.shape == ():
209 raise AssertionError(f"{prefix}{a!r} != {b!r} (differ by {float(separation)}{self.unit})")
210 nan_mismatch = int(np.count_nonzero(mismatch & (self.get_nan_mask() | other.get_nan_mask())))
211 nan_str = f"; {nan_mismatch} NaN mismatches" if nan_mismatch else ""
212 if np.isnan(separation).all():
213 raise AssertionError(f"{prefix}{n_mismatch}/{mismatch.size} values differ{nan_str}")
214 maxdiff = float(np.nanmax(separation))
215 # Convert to plain ints so the index prints as e.g. (2198,)
216 # rather than (np.int64(2198),).
217 loc = tuple(int(i) for i in np.unravel_index(np.nanargmax(separation), separation.shape))
218 raise AssertionError(
219 f"{prefix}{n_mismatch}/{mismatch.size} values differ; "
220 f"max abs diff {maxdiff}{self.unit} at index {loc} {nan_str}"
221 )
224class _QuantityCmpArg(_BaseCmpArg):
225 """Direct comparison of (possibly dimensionless)
226 `astropy.units.Quantity` instances, used for all non-separation
227 floating-point comparisons.
228 """
230 value: u.Quantity
232 @property
233 def unit(self) -> u.UnitBase:
234 return self.value.unit
236 def get_separation(self, other: Self) -> np.ndarray:
237 a_vals = self.value.to_value(self.unit)
238 b_vals = other.value.to_value(self.unit)
239 return np.abs(a_vals - b_vals)
241 def find_mismatches(
242 self,
243 other: Self,
244 *,
245 rtol: float,
246 atol_value: float,
247 equal_nan: bool,
248 ) -> np.ndarray | None:
249 a_vals = self.value.to_value(self.unit)
250 b_vals = other.value.to_value(self.unit)
251 mismatch = ~np.isclose(a_vals, b_vals, rtol=rtol, atol=atol_value, equal_nan=equal_nan)
252 return mismatch if np.any(mismatch) else None
255class _SeparationCmpArg(_BaseCmpArg):
256 """A base class for operands compared through a pairwise separation.
258 ``rtol`` is meaningless and must be zero for these comparisons.
259 """
261 @property
262 @abstractmethod
263 def unit(self) -> u.UnitBase:
264 raise NotImplementedError
266 def find_mismatches(
267 self,
268 other: Self,
269 *,
270 rtol: float,
271 atol_value: float,
272 equal_nan: bool,
273 ) -> np.ndarray | None:
274 if rtol != 0.0:
275 raise ValueError(f"rtol is meaningless for separations; use atol (in {self.unit}) instead.")
276 # A NaN separation (from a NaN operand element) fails the ``<=``
277 # test and so counts as a mismatch unless equal_nan clears it.
278 separation = np.asarray(self.get_separation(other))
279 mismatch = np.logical_not(separation <= atol_value)
280 if equal_nan:
281 mismatch = mismatch & ~(self.get_nan_mask() & other.get_nan_mask())
282 return mismatch if np.any(mismatch) else None
285class _TimeCmpArg(_SeparationCmpArg):
286 value: astropy.time.Time
288 @property
289 def unit(self) -> u.UnitBase:
290 return u.s
292 def get_separation(self, other: Self) -> np.ndarray:
293 return np.abs((self.value - other.value).sec)
295 def get_nan_mask(self) -> np.ndarray:
296 # astropy.time.Time rejects non-finite values at construction.
297 return np.zeros(self.value.shape, dtype=bool)
300class _SkyCoordCmpArg(_SeparationCmpArg):
301 value: SkyCoord
303 @property
304 def unit(self) -> u.UnitBase:
305 return u.arcsec
307 def get_separation(self, other: Self) -> np.ndarray:
308 return self.value.separation(other.value).to_value(self.unit)
310 def get_nan_mask(self) -> np.ndarray:
311 # A coordinate counts as NaN only when *both* components are NaN;
312 # a half-NaN coordinate is a difference to report, not a NaN match.
313 return np.isnan(self.value.ra) & np.isnan(self.value.dec)
316class _WrappedAngleCmpArg(_SeparationCmpArg):
317 value: u.Quantity
319 @property
320 def unit(self) -> u.UnitBase:
321 return u.rad
323 def get_separation(self, other: Self) -> np.ndarray:
324 a_rad = np.asarray(self.value.to_value(u.rad))
325 b_rad = np.asarray(other.value.to_value(u.rad))
326 return np.abs((a_rad - b_rad + np.pi) % (2.0 * np.pi) - np.pi)
329def _classify_operand(arg: Any, *, wrap_angles: bool = False) -> _BaseCmpArg | None:
330 """Coerce one comparison operand into its comparison-argument wrapper.
332 Parameters
333 ----------
334 arg
335 Operand to classify.
336 wrap_angles
337 If `True`, ``arg`` must be an angular `~astropy.units.Quantity` and
338 is classified for wrapped-angle comparison.
340 Returns
341 -------
342 cmp_arg : _BaseCmpArg | None
343 Wrapper selecting the comparison space for ``arg``, or `None` for
344 operands that cannot hold floating-point values (integers, bools,
345 strings); those are compared exactly by the driver.
346 """
347 if wrap_angles:
348 return _WrappedAngleCmpArg(u.Quantity(arg).to(u.rad))
349 match arg:
350 case astropy.time.Time():
351 return _TimeCmpArg(arg)
352 case astropy.time.TimeDelta():
353 return _QuantityCmpArg(arg.to(u.s))
354 case SkyCoord():
355 return _SkyCoordCmpArg(arg)
356 case u.Quantity():
357 if arg.dtype.kind not in "fc": 357 ↛ 359line 357 didn't jump to line 359 because the condition on line 357 was never true
358 # Coerce integer Quantities to floating point.
359 arg = arg.astype(float)
360 return _QuantityCmpArg(arg)
361 arg = np.asarray(arg)
362 if arg.dtype.kind not in "fc":
363 return None
364 return _QuantityCmpArg(u.Quantity(arg))
367def assert_values_equal(
368 a: (np.ndarray | u.Quantity | float | astropy.time.Time | astropy.time.TimeDelta | SkyCoord),
369 b: (np.ndarray | u.Quantity | float | astropy.time.Time | astropy.time.TimeDelta | SkyCoord),
370 *,
371 rtol: float = 0.0,
372 atol: float | u.Quantity | astropy.time.TimeDelta = 0.0,
373 equal_nan: bool = True,
374 label: str = "",
375 wrap_angles: bool = False,
376) -> None:
377 """Assert that two arrays, quantities, floats, astropy times, or sky
378 coordinates are equal or close.
380 On mismatch this raises an `AssertionError` with a concise report (count,
381 max absolute difference, location, NaN-mismatch count) rather than dumping
382 the arrays. Exact (non-floating-point) mismatches instead report through
383 `numpy.testing.assert_array_equal`.
385 Both operands must be of the same kind: both floating-point numeric
386 (bare or as `~astropy.units.Quantity` or `~astropy.time.TimeDelta`),
387 or both the same astropy type (`~astropy.time.Time` or
388 `astropy.coordinates.SkyCoord`).
389 Operands that cannot hold floating-point values (integers, bools,
390 strings) are compared exactly, with no tolerance or NaN semantics.
391 Tolerances never have their units guessed: a bare-float ``atol`` is
392 only valid for dimensionless operands, and for unitful operands only
393 exactly ``0.0`` (which is unit-independent) is accepted without units.
395 Parameters
396 ----------
397 a
398 First value to compare.
399 b
400 Second value to compare.
401 rtol
402 Relative tolerance, applied to the operand magnitudes. Must be zero
403 for comparisons of separations.
404 atol
405 Absolute tolerance. For dimensionless operands, a bare float. For
406 unitful operands, a `~astropy.units.Quantity` with units compatible
407 with ``a``, or exactly ``0.0``.
408 equal_nan
409 If `True`, treat NaN as equal to NaN. A `astropy.coordinates.SkyCoord`
410 element counts as NaN only when both of its coordinates are NaN.
411 label
412 Prefix prepended to the failure message.
413 wrap_angles
414 If `True`, compare ``a`` and ``b`` (angular quantities) modulo 2 pi.
415 """
416 a_arg = _classify_operand(a, wrap_angles=wrap_angles)
417 b_arg = _classify_operand(b, wrap_angles=wrap_angles)
418 if type(a_arg) is not type(b_arg):
419 raise TypeError(f"cannot compare {type(a).__name__} with {type(b).__name__}")
420 if a_arg is None or b_arg is None:
421 # The type check above means both are None: exact comparison of
422 # operands that cannot hold floating-point values.
423 if rtol != 0.0 or atol != 0.0:
424 raise ValueError("rtol and atol are only supported for floating-point operands.")
425 np.testing.assert_array_equal(
426 np.asarray(a),
427 np.asarray(b),
428 strict=True,
429 err_msg=f"{label}: arrays are not equal." if label else "",
430 )
431 return
432 atol_value = a_arg.coerce_atol(atol)
433 if a_arg.value.shape != b_arg.value.shape:
434 raise AssertionError(f"{label}: shape {a_arg.value.shape} != {b_arg.value.shape}")
435 mismatch = a_arg.find_mismatches(b_arg, rtol=rtol, atol_value=atol_value, equal_nan=equal_nan)
436 if mismatch is not None:
437 a_arg.report_mismatches(b_arg, a, b, mismatch, label)
440@contextmanager
441def annotate_errors(note: str) -> Generator[None]:
442 """Attach a note to any exception raised in the block.
444 Parameters
445 ----------
446 note
447 Note to attach on failure.
448 """
449 try:
450 yield
451 except Exception as err:
452 err.add_note(note)
453 raise
456def assert_images_equal(
457 a: Image,
458 b: Image,
459 *,
460 rtol: float = 0.0,
461 atol: float = 0.0,
462 expect_view: bool | Literal["array"] | None = None,
463) -> None:
464 """Assert that two images are equal or nearly equal.
466 Parameters
467 ----------
468 a
469 First image to compare.
470 b
471 Second image to compare.
472 rtol
473 Relative tolerance for the pixel comparison.
474 atol
475 Absolute tolerance for the pixel comparison.
476 expect_view
477 If not `None`, also assert whether ``b`` shares memory with ``a``
478 (i.e. is a view); ``"array"`` checks only the pixel arrays.
479 """
480 assert a.bbox == b.bbox
481 assert a.unit == b.unit
482 assert_sky_projections_equal(a.sky_projection, b.sky_projection)
483 if expect_view is not None:
484 assert np.may_share_memory(a.array, b.array) == bool(expect_view)
485 if expect_view == "array":
486 assert a.metadata.native == b.metadata.native
487 else:
488 assert (a._metadata is b._metadata) == expect_view
489 if not expect_view:
490 assert_values_equal(a.array, b.array, atol=atol, rtol=rtol)
491 assert a.metadata.native == b.metadata.native
494def _note_mask_difference(
495 err: AssertionError, a: Mask, b: Mask, *, label: str, names: Mapping[str, str] | None = None
496) -> None:
497 """Embed a plane-by-plane `Mask.compare` summary as a note on ``err``.
499 Parameters
500 ----------
501 err
502 Exception to attach the summary to.
503 a
504 First mask compared.
505 b
506 Second mask compared.
507 label
508 Prefix for the failure message.
509 names
510 Optional mapping from mask-plane names used by ``a`` and ``b`` to
511 names to display in the summary (e.g. the legacy names, when
512 comparing against legacy code).
513 """
514 diff = a.compare(b)
515 if not diff: 515 ↛ 516line 515 didn't jump to line 516 because the condition on line 515 was never true
516 return
517 shown = names or {}
518 summary = "\n".join(
519 f"{label}[{shown.get(name, name)}]: +{added} -{removed}" for name, (added, removed) in diff.items()
520 )
521 err.add_note(summary)
524def assert_masks_equal(a: Mask, b: Mask) -> None:
525 """Assert that two masks are equal or nearly equal.
527 Parameters
528 ----------
529 a
530 First mask to compare.
531 b
532 Second mask to compare.
533 """
534 assert a.bbox == b.bbox
535 assert a.schema == b.schema
536 assert a.metadata.native == b.metadata.native
537 assert_sky_projections_equal(a.sky_projection, b.sky_projection)
538 try:
539 assert_values_equal(a.array, b.array, label="mask")
540 except AssertionError as err:
541 _note_mask_difference(err, a, b, label="mask")
542 raise
545def assert_masked_images_equal(
546 a: MaskedImage,
547 b: MaskedImage,
548 *,
549 rtol: float = 0.0,
550 atol: float = 0.0,
551 expect_view: bool | None = None,
552) -> None:
553 """Assert that two masked images are equal or nearly equal.
555 Parameters
556 ----------
557 a
558 First masked image to compare.
559 b
560 Second masked image to compare.
561 rtol
562 Relative tolerance for the pixel comparison.
563 atol
564 Absolute tolerance for the pixel comparison.
565 expect_view
566 If not `None`, also assert whether ``b`` shares memory with ``a``
567 (i.e. is a view).
568 """
569 assert a.metadata.native == b.metadata.native
570 assert_sky_projections_equal(a.sky_projection, b.sky_projection)
571 assert_images_equal(a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view)
572 assert_masks_equal(a.mask, b.mask)
573 assert_images_equal(a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view)
576def assert_psfs_equal(
577 psf1: PointSpreadFunction,
578 psf2: PointSpreadFunction,
579 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
580) -> int:
581 """Compare two PSF objets.
583 Parameters
584 ----------
585 psf1
586 Point-spread function to test.
587 psf2
588 The other point-spread function to test.
589 points
590 Points to evaluate the PSFs at. If not provided, the intersection of
591 the PSF bounding boxes are used to generate points on a grid.
593 Returns
594 -------
595 `int`
596 The number of points actually tested.
597 """
598 if points is None: 598 ↛ 601line 598 didn't jump to line 601 because the condition on line 598 was always true
599 points = psf1.bounds.bbox.intersection(psf2.bounds.bbox).meshgrid(3).map(np.ravel)
601 assert psf1.kernel_bbox == psf2.kernel_bbox
603 n_points_tested: int = 0
604 for x, y in zip(points.x, points.y):
605 # The two PSFs must agree on which points fall inside their input
606 # domain. Querying ``.contains`` directly (rather than relying on
607 # ``compute_kernel_image`` to raise) makes this test tolerant of
608 # implementations that do not raise on out-of-domain points -- in
609 # particular ``CellPointSpreadFunction``, where evaluating in a
610 # missing cell does not always raise ``BoundsError``.
611 contains1 = psf1.bounds.contains(x=x, y=y)
612 contains2 = psf2.bounds.contains(x=x, y=y)
613 assert contains1 == contains2, (
614 f"PSFs disagree on whether ({x}, {y}) is in-bounds: psf1={contains1}, psf2={contains2}"
615 )
616 if not contains1: 616 ↛ 617line 616 didn't jump to line 617 because the condition on line 616 was never true
617 continue
618 assert psf1.compute_kernel_image(x=x, y=y) == psf2.compute_kernel_image(x=x, y=y)
619 assert psf1.compute_stellar_bbox(x=x, y=y) == psf2.compute_stellar_bbox(x=x, y=y)
620 assert psf1.compute_stellar_image(x=x, y=y) == psf2.compute_stellar_image(x=x, y=y)
621 n_points_tested += 1
622 return n_points_tested
625def assert_visit_images_equal(
626 a: VisitImage,
627 b: VisitImage,
628 *,
629 expect_view: bool | None = None,
630) -> None:
631 """Assert that two `.VisitImage` instances carry the same persistent state.
633 Extends `assert_masked_images_equal` with the VisitImage-specific
634 attributes (PSF, filter, observation info, detector, aperture
635 corrections, photometric scaling, backgrounds, polygon bounds,
636 summary stats) so a round-trip check on a `.VisitImage` does not
637 silently miss differences in any of them.
639 Parameters
640 ----------
641 a
642 First visit image to compare.
643 b
644 Second visit image to compare.
645 expect_view
646 If not `None`, also assert whether ``b`` shares memory with ``a``
647 (i.e. is a view).
648 """
649 assert_masked_images_equal(a, b, expect_view=expect_view)
650 assert a.summary_stats == b.summary_stats
651 assert a.physical_filter == b.physical_filter
652 assert a.band == b.band
653 assert a.obs_info == b.obs_info
654 assert a.detector == b.detector
655 assert dict(a.aperture_corrections) == dict(b.aperture_corrections)
656 assert a.photometric_scaling == b.photometric_scaling
657 assert dict(a.backgrounds) == dict(b.backgrounds)
658 assert a.backgrounds.subtracted == b.backgrounds.subtracted
659 assert a.bounds == b.bounds
660 assert_psfs_equal(a.psf, b.psf)
663def assert_cell_coadds_equal(
664 a: CellCoadd,
665 b: CellCoadd,
666 *,
667 expect_view: bool | None = None,
668) -> None:
669 """Assert that two `.CellCoadd` instances carry the same persistent state.
671 Extends the masked-image-style equality check with the
672 CellCoadd-specific attributes (PSF, cell grid, missing cells,
673 backgrounds, patch/tract, band) so a round-trip check does not
674 silently miss differences in any of them.
676 Parameters
677 ----------
678 a
679 First cell coadd to compare.
680 b
681 Second cell coadd to compare.
682 expect_view
683 If not `None`, also assert whether ``b`` shares memory with ``a``
684 (i.e. is a view).
685 """
686 assert_masked_images_equal(a, b, expect_view=expect_view)
687 assert a.band == b.band
688 assert a.patch == b.patch
689 assert a.tract == b.tract
690 assert a.grid == b.grid
691 assert a.bounds.missing == b.bounds.missing
692 assert dict(a.backgrounds) == dict(b.backgrounds)
693 assert a.backgrounds.subtracted == b.backgrounds.subtracted
694 assert_psfs_equal(a.psf, b.psf)
697def compare_image_to_legacy(
698 image: Image,
699 legacy_image: Any,
700 expect_view: bool | None = None,
701 *,
702 label: str = "image",
703 atol: float = 0.0,
704) -> None:
705 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object.
707 Parameters
708 ----------
709 image
710 Image to compare.
711 legacy_image
712 Legacy `lsst.afw.image.Image` to compare against.
713 expect_view
714 If not `None`, also assert whether ``image`` shares memory with
715 ``legacy_image`` (i.e. is a view).
716 label
717 Prefix for the failure message.
718 atol
719 Absolute tolerance for the pixel-value comparison.
720 """
721 assert image.bbox == Box.from_legacy(legacy_image.getBBox())
722 if expect_view is not None: 722 ↛ 724line 722 didn't jump to line 724 because the condition on line 722 was always true
723 assert np.may_share_memory(image.array, legacy_image.array) == expect_view
724 if not expect_view: 724 ↛ exitline 724 didn't return from function 'compare_image_to_legacy' because the condition on line 724 was always true
725 assert_values_equal(image.array, legacy_image.array, atol=atol, label=label)
728def compare_mask_to_legacy(
729 mask: Mask,
730 legacy_mask: Any,
731 plane_map: Mapping[str, MaskPlane] | None = None,
732 *,
733 label: str = "mask",
734) -> None:
735 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object.
737 Parameters
738 ----------
739 mask
740 Mask to compare.
741 legacy_mask
742 Legacy `lsst.afw.image.Mask` to compare against.
743 plane_map
744 Mapping from legacy plane name to the new mask plane; defaults to
745 the planes in ``mask.schema``. Schema planes the map does not name
746 are compared under their own names. A plane the legacy mask does
747 not define must have no pixels set in ``mask``.
748 label
749 Prefix for the failure message; each plane is reported as
750 ``"{label}[<plane>]"``.
751 """
752 assert mask.bbox == Box.from_legacy(legacy_mask.getBBox())
753 if plane_map is None:
754 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
755 legacy_planes = legacy_mask.getMaskPlaneDict()
756 # Compare schema planes the map does not name (e.g. the planes of
757 # `get_legacy_optional_mask_planes`) under their own names.
758 plane_map = dict(plane_map)
759 mapped = {plane.name for plane in plane_map.values()}
760 plane_map.update(
761 {plane.name: plane for plane in mask.schema if plane is not None and plane.name not in mapped}
762 )
763 try:
764 for old_name, new_plane in plane_map.items():
765 pixels = mask.get(new_plane.name)
766 if old_name not in legacy_planes: 766 ↛ 768line 766 didn't jump to line 768 because the condition on line 766 was never true
767 # A plane the legacy mask does not define must be empty.
768 assert_values_equal(np.zeros_like(pixels), pixels, label=f"{label}[{old_name}]")
769 continue
770 assert_values_equal(
771 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
772 pixels,
773 label=f"{label}[{old_name}]",
774 )
775 except AssertionError as err:
776 names = {new_plane.name: old_name for old_name, new_plane in plane_map.items()}
777 _note_mask_difference(err, mask, Mask.from_legacy(legacy_mask, plane_map), label=label, names=names)
778 raise
781def compare_masked_image_to_legacy(
782 masked_image: MaskedImage,
783 legacy_masked_image: Any,
784 *,
785 plane_map: Mapping[str, MaskPlane] | None = None,
786 expect_view: bool | None = None,
787 alternates: Mapping[str, Any] | None = None,
788 image_atol: float = 0.0,
789 variance_atol: float = 0.0,
790) -> None:
791 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
792 object.
794 Parameters
795 ----------
796 masked_image
797 New image to test.
798 legacy_masked_image
799 Legacy image to test against.
800 plane_map
801 Mapping between new and legacy mask planes.
802 expect_view
803 Whether to test that the image and variance arrays do or do not share
804 memory.
805 alternates
806 A mapping of other versions of one or more (new) components to also
807 check against the legacy versions of those components.
808 image_atol
809 Absolute tolerance for the image pixel comparison.
810 variance_atol
811 Absolute tolerance for the variance pixel comparison. Separate from
812 ``image_atol`` because the two have different units. The mask
813 comparison is always exact.
814 """
815 compare_image_to_legacy(
816 masked_image.image,
817 legacy_masked_image.getImage(),
818 expect_view=expect_view,
819 label="image",
820 atol=image_atol,
821 )
822 compare_mask_to_legacy(
823 masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask"
824 )
825 compare_image_to_legacy(
826 masked_image.variance,
827 legacy_masked_image.getVariance(),
828 expect_view=expect_view,
829 label="variance",
830 atol=variance_atol,
831 )
832 if alternates: 832 ↛ 833line 832 didn't jump to line 833 because the condition on line 832 was never true
833 if image := alternates.get("image"):
834 compare_image_to_legacy(
835 image,
836 legacy_masked_image.getImage(),
837 expect_view=expect_view,
838 label="image",
839 atol=image_atol,
840 )
841 if mask := alternates.get("mask"):
842 compare_mask_to_legacy(mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask")
843 if variance := alternates.get("variance"):
844 compare_image_to_legacy(
845 variance,
846 legacy_masked_image.getVariance(),
847 expect_view=expect_view,
848 label="variance",
849 atol=variance_atol,
850 )
853def compare_visit_image_to_legacy(
854 visit_image: VisitImage,
855 legacy_exposure: Any,
856 *,
857 plane_map: Mapping[str, MaskPlane] | None = None,
858 expect_view: bool | None = None,
859 instrument: str,
860 visit: int,
861 detector: int,
862 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
863 alternates: Mapping[str, Any] | None = None,
864 check_photometric_scaling: bool = True,
865) -> None:
866 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
867 object.
869 Parameters
870 ----------
871 visit_image
872 New image to test.
873 legacy_exposure
874 Legacy image to test against.
875 plane_map
876 Mapping between new and legacy mask planes; defaults to the planes
877 used by the legacy (`lsst.afw.image.Exposure`) visit image.
878 expect_view
879 Whether to test that the image and variance arrays do or do not share
880 memory.
881 instrument
882 Expected instrument name.
883 visit
884 Expected visit ID.
885 detector
886 Expected detector ID.
887 applied_legacy_photo_calib
888 Legacy `lsst.afw.image.PhotoCalib` already applied to
889 ``legacy_exposure``, used when comparing photometric scaling.
890 alternates
891 A mapping of other versions of one or more (new) components to also
892 check against the legacy versions of those components.
893 check_photometric_scaling
894 If `False`, skip the photometric-scaling comparison.
895 """
896 if plane_map is None:
897 plane_map = get_legacy_visit_image_mask_planes()
898 compare_masked_image_to_legacy(
899 visit_image,
900 legacy_exposure.getMaskedImage(),
901 plane_map=plane_map,
902 expect_view=expect_view,
903 alternates=alternates,
904 )
905 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
906 with annotate_errors("sky_projection"):
907 compare_sky_projection_to_legacy_wcs(
908 visit_image.sky_projection,
909 legacy_exposure.getWcs(),
910 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
911 visit_image.bbox,
912 )
913 assert visit_image.sky_projection is visit_image.mask.sky_projection
914 assert visit_image.sky_projection is visit_image.variance.sky_projection
915 with annotate_errors("psf"):
916 compare_psf_to_legacy(visit_image.psf, legacy_exposure.getPsf())
917 with annotate_errors("summary_stats"):
918 compare_observation_summary_stats_to_legacy(
919 visit_image.summary_stats, legacy_exposure.info.getSummaryStats()
920 )
921 with annotate_errors("detector"):
922 compare_detector_to_legacy(visit_image.detector, legacy_exposure.getDetector(), is_raw_assembled=True)
923 # Make a tiny box for Field comparisons that need to make arrays; that can
924 # get expensive otherwise.
925 tiny_bbox = detector_bbox.local[2:4, 3:6]
926 with annotate_errors("aperture_corrections"):
927 compare_aperture_corrections_to_legacy(
928 visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
929 )
930 if check_photometric_scaling:
931 with annotate_errors("photometric_scaling"):
932 compare_photo_calib_to_legacy(
933 visit_image.photometric_scaling,
934 legacy_exposure.info.getPhotoCalib(),
935 applied_legacy_photo_calib=applied_legacy_photo_calib,
936 subimage_bbox=tiny_bbox,
937 )
938 if alternates:
939 if (bbox := alternates.get("bbox")) is not None:
940 assert bbox == visit_image.bbox
941 if sky_projection := alternates.get("sky_projection"):
942 with annotate_errors("sky_projection"):
943 compare_sky_projection_to_legacy_wcs(
944 sky_projection,
945 legacy_exposure.getWcs(),
946 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
947 visit_image.bbox,
948 )
949 if psf := alternates.get("psf"):
950 with annotate_errors("psf"):
951 compare_psf_to_legacy(psf, legacy_exposure.getPsf())
952 if summary_stats := alternates.get("summary_stats"):
953 with annotate_errors("summary_stats"):
954 compare_observation_summary_stats_to_legacy(
955 summary_stats, legacy_exposure.info.getSummaryStats()
956 )
957 if detector_obj := alternates.get("detector"):
958 with annotate_errors("detector"):
959 compare_detector_to_legacy(detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True)
960 if obs_info := alternates.get("obs_info"):
961 visitInfo = legacy_exposure.visitInfo
962 assert obs_info.instrument == visitInfo.getInstrumentLabel()
963 if aperture_corrections := alternates.get("aperture_corrections"):
964 with annotate_errors("aperture_corrections"):
965 compare_aperture_corrections_to_legacy(
966 aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
967 )
968 if check_photometric_scaling:
969 if (photometric_scaling := alternates.get("photometric_scaling", ...)) is not ...:
970 with annotate_errors("photometric_scaling"):
971 compare_photo_calib_to_legacy(
972 photometric_scaling,
973 legacy_exposure.info.getPhotoCalib(),
974 applied_legacy_photo_calib=applied_legacy_photo_calib,
975 subimage_bbox=tiny_bbox,
976 )
979def compare_photo_calib_to_legacy(
980 photometric_scaling: BaseField | None,
981 legacy_photo_calib: LegacyPhotoCalib,
982 *,
983 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
984 subimage_bbox: Box,
985) -> None:
986 if legacy_photo_calib._isConstant:
987 if legacy_photo_calib.getCalibrationMean() == 1.0:
988 if applied_legacy_photo_calib is None:
989 assert photometric_scaling is None
990 return
991 else:
992 legacy_photo_calib = applied_legacy_photo_calib
993 if legacy_photo_calib._isConstant:
994 assert isinstance(photometric_scaling, ChebyshevField)
995 assert_values_equal(
996 photometric_scaling.coefficients,
997 np.array([[legacy_photo_calib.getCalibrationMean()]]),
998 # The coefficients come from a numerical fit, so this can be
999 # bigger than just round-off error epislon.
1000 rtol=1e-5,
1001 )
1002 else:
1003 assert photometric_scaling is not None
1004 compare_field_to_legacy(
1005 photometric_scaling / legacy_photo_calib.getCalibrationMean(),
1006 legacy_photo_calib.computeScaledCalibration(),
1007 subimage_bbox,
1008 )
1011def _check_cell_coadd_geometry(cell_coadd: CellCoadd, tract_bbox: Box) -> None:
1012 """Check the geometric self-consistency of a cell coadd's patch,
1013 bounds, grid, and bounding box.
1015 Parameters
1016 ----------
1017 cell_coadd
1018 Coadd to test.
1019 tract_bbox
1020 Bounding box of the full tract.
1021 """
1022 assert tract_bbox.contains(cell_coadd.patch.outer_bbox)
1023 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox)
1024 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox)
1025 # bounds.bbox is snapped to cell boundaries, cell_coadd.bbox is not.
1026 assert cell_coadd.bounds.bbox.contains(cell_coadd.bbox)
1027 assert cell_coadd.grid.bbox.contains(cell_coadd.bbox)
1030def cell_center_points(bounds: CellGridBounds) -> YX[np.ndarray]:
1031 """Return the centers of the non-missing cells of cell grid bounds.
1033 These are the points at which approximations made by conversions from
1034 legacy (non-cell) coadds are expected to be accurate.
1036 Parameters
1037 ----------
1038 bounds
1039 Cell grid bounds whose cells' centers are returned.
1041 Returns
1042 -------
1043 `.YX` [ `numpy.ndarray` ]
1044 Centers of the non-missing cells within ``bounds``.
1045 """
1046 xs: list[float] = []
1047 ys: list[float] = []
1048 for cell in bounds.cell_indices():
1049 bbox = bounds.grid.bbox_of(cell)
1050 xs.append(bbox.x.center)
1051 ys.append(bbox.y.center)
1052 return YX(y=np.array(ys), x=np.array(xs))
1055def compare_cell_coadd_to_legacy(
1056 cell_coadd: CellCoadd,
1057 legacy_cell_coadd: MultipleCellCoadd,
1058 *,
1059 tract_bbox: Box,
1060 plane_map: Mapping[str, MaskPlane] | None = None,
1061 alternates: Mapping[str, Any] | None = None,
1062 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
1063) -> None:
1064 """Compare a `.cells.CellCoadd` object to a legacy
1065 `lsst.cell_coadds.MultipleCellCoadd` object.
1067 Parameters
1068 ----------
1069 cell_coadd
1070 New coadd to test.
1071 legacy_cell_coadd
1072 Legacy coadd to test against.
1073 tract_bbox
1074 Bounding box of the full tract.
1075 plane_map
1076 Mapping between new and legacy mask planes.
1077 alternates
1078 A mapping of other versions of one or more (new) components to also
1079 check against the legacy versions of those components.
1080 psf_points
1081 Points to use to compare the PSFs.
1082 """
1083 legacy_stitched = legacy_cell_coadd.stitch(cell_coadd.bbox.to_legacy())
1084 compare_image_to_legacy(cell_coadd.image, legacy_stitched.image, expect_view=False)
1085 compare_mask_to_legacy(cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map)
1086 compare_image_to_legacy(cell_coadd.variance, legacy_stitched.variance, expect_view=False)
1087 if legacy_stitched.mask_fractions is not None:
1088 compare_image_to_legacy(
1089 cell_coadd.mask_fractions["rejected"], legacy_stitched.mask_fractions, expect_view=False
1090 )
1091 for n in range(legacy_stitched.n_noise_realizations):
1092 compare_image_to_legacy(
1093 cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False
1094 )
1095 assert cell_coadd.skymap == legacy_stitched.identifiers.skymap
1096 assert cell_coadd.tract == legacy_stitched.identifiers.tract
1097 assert cell_coadd.patch.index.x == legacy_stitched.identifiers.patch.x
1098 assert cell_coadd.patch.index.y == legacy_stitched.identifiers.patch.y
1099 assert cell_coadd.band == legacy_stitched.identifiers.band
1100 _check_cell_coadd_geometry(cell_coadd, tract_bbox)
1101 assert cell_coadd.unit == u.Unit(legacy_cell_coadd.common.units.value)
1102 compare_sky_projection_to_legacy_wcs(
1103 cell_coadd.sky_projection,
1104 legacy_cell_coadd.wcs,
1105 TractFrame(
1106 skymap=legacy_cell_coadd.identifiers.skymap,
1107 tract=legacy_cell_coadd.identifiers.tract,
1108 bbox=tract_bbox,
1109 ),
1110 cell_coadd.bbox,
1111 is_fits=True,
1112 )
1113 assert cell_coadd.sky_projection is cell_coadd.mask.sky_projection
1114 assert cell_coadd.sky_projection is cell_coadd.variance.sky_projection
1115 compare_psf_to_legacy(
1116 cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points
1117 )
1118 compare_aperture_corrections_to_legacy(
1119 cell_coadd.aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
1120 )
1121 compare_cell_coadd_provenance_to_legacy(cell_coadd.provenance, legacy_cell_coadd)
1122 if alternates:
1123 if sky_projection := alternates.get("sky_projection"):
1124 compare_sky_projection_to_legacy_wcs(
1125 sky_projection,
1126 legacy_stitched.wcs,
1127 TractFrame(
1128 skymap=legacy_cell_coadd.identifiers.skymap,
1129 tract=legacy_cell_coadd.identifiers.tract,
1130 bbox=tract_bbox,
1131 ),
1132 cell_coadd.bbox,
1133 is_fits=True,
1134 )
1135 if psf := alternates.get("psf"):
1136 compare_psf_to_legacy(psf, legacy_stitched.psf, points=psf_points)
1137 if aperture_corrections := alternates.get("aperture_corrections"):
1138 compare_aperture_corrections_to_legacy(
1139 aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
1140 )
1141 if provenance := alternates.get("provenance"):
1142 compare_cell_coadd_provenance_to_legacy(provenance, legacy_cell_coadd)
1145def compare_cell_coadd_provenance_to_legacy(
1146 provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd
1147) -> None:
1148 """Compare a `.cells.CoaddProvenance` object to a legacy
1149 `lsst.cell_coadds.MultipleCellCoadd` object.
1151 Parameters
1152 ----------
1153 provenance
1154 New provenance object to test.
1155 legacy_cell_coadd
1156 Legacy coadd to test against.
1157 """
1158 from lsst.cell_coadds import ObservationIdentifiers
1160 for legacy_cell in legacy_cell_coadd.cells.values():
1161 cell_index = CellIJ.from_legacy(legacy_cell.identifiers.cell)
1162 prov = provenance[cell_index]
1163 assert prov is not None, f"Legacy cell {cell_index} has inputs, so its provenance is not empty."
1164 legacy_table = astropy.table.Table(
1165 rows=[
1166 [
1167 ids.instrument,
1168 ids.visit,
1169 ids.detector,
1170 ids.day_obs,
1171 ids.physical_filter,
1172 legacy_input.overlaps_center,
1173 legacy_input.overlap_fraction,
1174 legacy_input.weight,
1175 legacy_input.psf_shape.getIxx(),
1176 legacy_input.psf_shape.getIyy(),
1177 legacy_input.psf_shape.getIxy(),
1178 legacy_input.psf_shape_flag,
1179 ]
1180 for ids, legacy_input in legacy_cell.inputs.items()
1181 ],
1182 dtype=[
1183 np.object_,
1184 np.uint64,
1185 np.uint16,
1186 np.uint32,
1187 np.object_,
1188 np.bool_,
1189 np.float64,
1190 np.float64,
1191 np.float64,
1192 np.float64,
1193 np.float64,
1194 np.bool_,
1195 ],
1196 names=[
1197 "instrument",
1198 "visit",
1199 "detector",
1200 "day_obs",
1201 "physical_filter",
1202 "overlaps_center",
1203 "overlap_fraction",
1204 "weight",
1205 "psf_shape_xx",
1206 "psf_shape_yy",
1207 "psf_shape_xy",
1208 "psf_shape_flag",
1209 ],
1210 )
1211 # For a single cell all 'inputs' are also 'contributions'.
1212 assert len(legacy_cell.inputs) == len(prov.inputs)
1213 assert len(legacy_cell.inputs) == len(prov.contributions)
1214 prov.inputs.sort(["instrument", "visit", "detector"])
1215 prov.contributions.sort(["instrument", "visit", "detector"])
1216 legacy_table.sort(["instrument", "visit", "detector"])
1217 assert_values_equal(prov.inputs["instrument"], prov.contributions["instrument"], label="inputs")
1218 assert_values_equal(prov.inputs["visit"], prov.contributions["visit"], label="inputs")
1219 assert_values_equal(prov.inputs["detector"], prov.contributions["detector"], label="inputs")
1220 assert_values_equal(prov.inputs["instrument"], legacy_table["instrument"], label="inputs")
1221 assert_values_equal(prov.inputs["visit"], legacy_table["visit"], label="inputs")
1222 assert_values_equal(prov.inputs["detector"], legacy_table["detector"], label="inputs")
1223 assert_values_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"], label="inputs")
1224 assert_values_equal(prov.inputs["day_obs"], legacy_table["day_obs"], label="inputs")
1225 assert_values_equal(
1226 prov.contributions["overlaps_center"], legacy_table["overlaps_center"], label="inputs"
1227 )
1228 assert_values_equal(
1229 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"], label="inputs"
1230 )
1231 assert_values_equal(prov.contributions["weight"], legacy_table["weight"], label="inputs")
1232 assert_values_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"], label="inputs")
1233 assert_values_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"], label="inputs")
1234 assert_values_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"], label="inputs")
1235 assert_values_equal(
1236 prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"], label="inputs"
1237 )
1238 for row in prov.inputs:
1239 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"})
1240 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key]
1241 assert legacy_polygon == row["polygon"].to_legacy()
1244def compare_cell_coadd_to_legacy_lazy(
1245 cell_coadd: CellCoadd,
1246 legacy_exposure: Any,
1247 *,
1248 skymap: str,
1249 tract: int,
1250 patch: int,
1251 band: str,
1252 tract_bbox: Box,
1253 plane_map: Mapping[str, MaskPlane] | None = None,
1254 unit: u.UnitBase = u.nJy,
1255 expect_view: bool | None = None,
1256 alternates: Mapping[str, Any] | None = None,
1257 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
1258 image_atol: float = 0.0,
1259 variance_atol: float = 0.0,
1260) -> None:
1261 """Compare a `.cells.CellCoadd` object to the legacy
1262 `lsst.afw.image.Exposure` instance it was converted from
1263 by ``CellCoadd.from_legacy_lazy_coadd``.
1265 Parameters
1266 ----------
1267 cell_coadd
1268 New coadd to test.
1269 legacy_exposure
1270 Legacy `lsst.afw.image.Exposure` to test against.
1271 skymap
1272 Expected skymap name, which the legacy exposure does not carry.
1273 tract
1274 Expected tract ID.
1275 patch
1276 Expected patch sequential ID.
1277 band
1278 Expected band name.
1279 tract_bbox
1280 Bounding box of the full tract.
1281 plane_map
1282 Mapping between new and legacy mask planes; defaults to the planes
1283 used by legacy (non-cell) coadds like DP2's ``template_coadd``.
1284 unit
1285 Expected pixel unit of the coadd.
1286 expect_view
1287 Whether to test that the image and variance arrays do or do not
1288 share memory with ``legacy_exposure``.
1289 alternates
1290 A mapping of other versions of one or more (new) components to also
1291 check against the legacy versions of those components.
1292 psf_points
1293 Points to use to compare the PSFs. Defaults to the centers of the
1294 coadd's non-missing cells.
1295 image_atol
1296 Absolute tolerance for the image pixel comparison.
1297 variance_atol
1298 Absolute tolerance for the variance pixel comparison.
1299 """
1300 if plane_map is None:
1301 plane_map = get_legacy_non_cell_coadd_mask_planes()
1302 compare_masked_image_to_legacy(
1303 cell_coadd,
1304 legacy_exposure.getMaskedImage(),
1305 plane_map=plane_map,
1306 expect_view=expect_view,
1307 image_atol=image_atol,
1308 variance_atol=variance_atol,
1309 )
1310 # A legacy (non-cell) coadd has no place for any of these components,
1311 # so a faithful conversion must not invent them.
1312 assert cell_coadd.mask_fractions == {}
1313 assert cell_coadd.noise_realizations == []
1314 assert cell_coadd.backgrounds.keys() == set()
1315 assert cell_coadd.backgrounds.subtracted is None
1316 assert cell_coadd.skymap == skymap
1317 assert cell_coadd.tract == tract
1318 assert cell_coadd.patch.id == patch
1319 assert cell_coadd.band == band
1320 assert cell_coadd.unit == unit
1321 _check_cell_coadd_geometry(cell_coadd, tract_bbox)
1322 assert cell_coadd.sky_projection is cell_coadd.mask.sky_projection
1323 assert cell_coadd.sky_projection is cell_coadd.variance.sky_projection
1324 compare_sky_projection_to_legacy_wcs(
1325 cell_coadd.sky_projection,
1326 legacy_exposure.getWcs(),
1327 TractFrame(skymap=skymap, tract=tract, bbox=tract_bbox),
1328 cell_coadd.bbox,
1329 is_fits=True,
1330 )
1331 if psf_points is None:
1332 psf_points = cell_center_points(cell_coadd.bounds)
1333 n_psf_points = compare_psf_to_legacy(
1334 cell_coadd.psf, legacy_exposure.getPsf(), points=psf_points, rtol=1e-5
1335 )
1336 assert n_psf_points > 0, "The legacy PSF could not be evaluated at any of the PSF comparison points."
1337 compare_aperture_corrections_to_legacy(
1338 cell_coadd.aperture_corrections,
1339 legacy_exposure.info.getApCorrMap(),
1340 cell_coadd.bbox,
1341 eval_points=psf_points,
1342 expect_subset_bounds=True,
1343 )
1344 compare_cell_coadd_provenance_to_legacy_lazy(cell_coadd.provenance, legacy_exposure, cell_coadd.bounds)
1345 if alternates:
1346 if sky_projection := alternates.get("sky_projection"):
1347 compare_sky_projection_to_legacy_wcs(
1348 sky_projection,
1349 legacy_exposure.getWcs(),
1350 TractFrame(skymap=skymap, tract=tract, bbox=tract_bbox),
1351 cell_coadd.bbox,
1352 is_fits=True,
1353 )
1354 if psf := alternates.get("psf"):
1355 assert compare_psf_to_legacy(psf, legacy_exposure.getPsf(), points=psf_points, rtol=1e-5) > 0
1356 if aperture_corrections := alternates.get("aperture_corrections"):
1357 compare_aperture_corrections_to_legacy(
1358 aperture_corrections,
1359 legacy_exposure.info.getApCorrMap(),
1360 cell_coadd.bbox,
1361 eval_points=psf_points,
1362 expect_subset_bounds=True,
1363 )
1364 if provenance := alternates.get("provenance"):
1365 compare_cell_coadd_provenance_to_legacy_lazy(provenance, legacy_exposure, cell_coadd.bounds)
1368def compare_cell_coadd_provenance_to_legacy_lazy(
1369 provenance: CoaddProvenance,
1370 legacy_exposure: Any,
1371 bounds: CellGridBounds,
1372) -> None:
1373 """Compare a `.cells.CoaddProvenance` object to the
1374 `lsst.afw.image.CoaddInputs` of the legacy (non-cell) coadd it was
1375 converted from.
1377 Parameters
1378 ----------
1379 provenance
1380 New provenance object to test.
1381 legacy_exposure
1382 Legacy `lsst.afw.image.Exposure` whose ``info.getCoaddInputs()``
1383 the provenance should be consistent with.
1384 bounds
1385 The cell bounds of the coadd this provenance belongs to.
1386 """
1387 coadd_inputs = legacy_exposure.info.getCoaddInputs()
1388 assert coadd_inputs is not None
1389 wcs = legacy_exposure.getWcs()
1390 legacy_by_key: dict[tuple[int, int], Any] = {(row["visit"], row["ccd"]): row for row in coadd_inputs.ccds}
1391 inputs = provenance.inputs
1392 prov_keys = {(row["visit"], row["detector"]) for row in inputs}
1393 # A coadd converted from a larger CoaddInputs (e.g. a postage stamp cut
1394 # from a whole-tract template) legitimately drops inputs that intersect
1395 # its bounds nowhere, so this is subset, not equality.
1396 assert prov_keys <= set(legacy_by_key), "Provenance inputs must be a subset of CoaddInputs ccds."
1397 instruments = {row["instrument"] for row in inputs}
1398 assert len(instruments) == 1, "A coadd's inputs must all share one instrument."
1399 for row in inputs:
1400 legacy_row = legacy_by_key[(row["visit"], row["detector"])]
1401 assert row["physical_filter"] == legacy_row["filter"]
1402 assert row["polygon"].area > 0
1403 # Group contribution rows by input, checking weights are carried over
1404 # exactly, and then group them by cell for the geometric checks below.
1405 contributions_by_input: dict[tuple[int, int], list[Any]] = {}
1406 contributions_by_cell: dict[CellIJ, list[Any]] = {}
1407 for row in provenance.contributions:
1408 key = (row["visit"], row["detector"])
1409 assert key in prov_keys
1410 contributions_by_input.setdefault(key, []).append(row)
1411 contributions_by_cell.setdefault(CellIJ(i=row["cell_i"], j=row["cell_j"]), []).append(row)
1412 for key, rows in contributions_by_input.items():
1413 assert rows, f"Input (visit={key[0]}, detector={key[1]}) must contribute somewhere."
1414 assert_values_equal(
1415 np.array([row["weight"] for row in rows]),
1416 np.full(len(rows), legacy_by_key[key]["weight"]),
1417 rtol=1e-12,
1418 label="contributions[weight]",
1419 )
1420 overlap = np.asarray(provenance.contributions["overlap_fraction"], dtype=np.float64)
1421 assert np.all((overlap > 0.0) & (overlap <= 1.0)), "Overlap fractions must be valid cell fractions."
1422 assert set(contributions_by_cell) <= set(bounds.cell_indices()), (
1423 "Contributions must be on non-missing cells."
1424 )
1425 for cell, rows in contributions_by_cell.items():
1426 bbox = bounds.grid.bbox_of(cell)
1427 overlapping = {
1428 (row["visit"], row["ccd"])
1429 for row in coadd_inputs.ccds.subsetContaining(
1430 bbox.center.to_legacy_float_point(), wcs, includeValidPolygon=True
1431 )
1432 }
1433 keys = {(row["visit"], row["detector"]) for row in rows}
1434 assert overlapping <= keys, "Inputs overlapping a cell center must contribute to that cell."
1435 centered = {(row["visit"], row["detector"]) for row in rows if row["overlaps_center"]}
1436 assert centered == overlapping, "overlaps_center must match which inputs cover the cell center."
1439def compare_psf_to_legacy(
1440 psf: PointSpreadFunction,
1441 legacy_psf: Any,
1442 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
1443 expect_legacy_raise_on_out_of_bounds: bool = False,
1444 *,
1445 rtol: float = 0.0,
1446 atol: float = 0.0,
1447) -> int:
1448 """Compare a PSF model object to its legacy interface.
1450 Parameters
1451 ----------
1452 psf
1453 Point-spread function to test.
1454 legacy_psf
1455 Legacy `lsst.afw.detection.Psf` instance to compare with.
1456 points
1457 Points to evaluate the PSFs at. If not provided, the intersection of
1458 the PSF bounding boxes are used to generate points on a grid.
1459 expect_legacy_raise_on_out_of_bounds
1460 If `True`, expect ``legacy_psf`` to raise
1461 `lsst.afw.detection.InvalidPsfError` when evaluated at a position
1462 considered out-of-bounds by ``psf``.
1463 rtol
1464 Relative tolerance for kernel and stellar image pixel values.
1465 The default of zero (with ``atol`` also zero) requires exact
1466 equality, appropriate when both PSFs were converted from the same
1467 representation. Nonzero values are for comparing approximations,
1468 such as cell-based evaluations of a spatially varying legacy PSF.
1469 atol
1470 Absolute tolerance for kernel and stellar image pixel values.
1472 Returns
1473 -------
1474 `int`
1475 The number of points actually tested.
1476 """
1477 from lsst.afw.detection import InvalidPsfError
1479 approximate = rtol != 0.0 or atol != 0.0
1480 if points is None:
1481 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel)
1482 legacy_points = arrays_to_legacy_points(points.x, points.y)
1483 n_points_tested: int = 0
1484 for p in legacy_points:
1485 if not psf.bounds.contains(x=p.x, y=p.y):
1486 if expect_legacy_raise_on_out_of_bounds:
1487 with pytest.raises(InvalidPsfError):
1488 legacy_psf.computeKernelImage(p)
1489 continue
1490 try:
1491 legacy_kernel = Image.from_legacy(legacy_psf.computeKernelImage(p))
1492 except InvalidPsfError:
1493 if not approximate:
1494 raise
1495 # An approximate conversion samples a cell's corners when the
1496 # legacy PSF cannot be evaluated at its center, so there is
1497 # nothing to compare at this point.
1498 continue
1499 legacy_kernel_bbox = Box.from_legacy(legacy_psf.computeKernelBBox(p))
1500 legacy_stellar_bbox = Box.from_legacy(legacy_psf.computeImageBBox(p))
1501 legacy_stellar = Image.from_legacy(legacy_psf.computeImage(p))
1502 kernel = psf.compute_kernel_image(x=p.x, y=p.y)
1503 stellar_bbox = psf.compute_stellar_bbox(x=p.x, y=p.y)
1504 stellar = psf.compute_stellar_image(x=p.x, y=p.y)
1505 if approximate:
1506 # An approximate conversion stores every cell's stamp in one
1507 # kernel bbox large enough for the legacy kernel at every sampled
1508 # point, so the legacy kernel must match the stamp's sub-box and
1509 # the rest of the stamp must be empty. Resampling the stamp into
1510 # the stellar image spreads flux into that padding, so only the
1511 # sub-box is comparable there.
1512 assert psf.kernel_bbox.contains(legacy_kernel_bbox)
1513 expected_kernel = Image(0.0, bbox=kernel.bbox, dtype=np.float64)
1514 expected_kernel[legacy_kernel.bbox] = legacy_kernel
1515 with annotate_errors("psf kernel image"):
1516 assert_images_equal(kernel, expected_kernel, rtol=rtol, atol=atol)
1517 assert stellar_bbox.contains(legacy_stellar_bbox)
1518 with annotate_errors("psf stellar image"):
1519 assert_images_equal(stellar[legacy_stellar.bbox], legacy_stellar, rtol=rtol, atol=atol)
1520 else:
1521 assert psf.kernel_bbox == legacy_kernel_bbox
1522 assert kernel == legacy_kernel
1523 assert stellar_bbox == legacy_stellar_bbox
1524 assert stellar == legacy_stellar
1525 n_points_tested += 1
1526 return n_points_tested
1529def compare_field_to_legacy(
1530 field: BaseField,
1531 legacy_field: Any,
1532 subimage_bbox: Box,
1533 eval_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
1534 expect_subset_bounds: bool = False,
1535 *,
1536 rtol: float = 1e-5,
1537) -> None:
1538 """Test a Field object by comparing it to an equivalent
1539 `lsst.afw.math.BoundedField`.
1541 Parameters
1542 ----------
1543 field
1544 Field to test.
1545 legacy_field : ``lsst.afw.math.BoundedField``
1546 Equivalent legacy bounded field.
1547 subimage_bbox
1548 Bounding box for full-image tests.
1549 eval_points
1550 Points at which to compare field values, instead of testing over a
1551 grid covering ``field.bounds``. This also disables image-rendering
1552 comparison checks.
1553 expect_subset_bounds
1554 If `True`, expect ``field.bounds`` to be contained by
1555 ``legacy_field.getBBox()``, not equal to it.
1556 rtol
1557 Relative tolerance for the evaluation-point field-value comparison.
1558 """
1559 from lsst.afw.math import BoundedField as LegacyBoundedField
1561 if expect_subset_bounds: 1561 ↛ 1562line 1561 didn't jump to line 1562 because the condition on line 1561 was never true
1562 assert Box.from_legacy(legacy_field.getBBox()).contains(field.bounds.bbox)
1563 else:
1564 assert field.bounds.bbox == Box.from_legacy(legacy_field.getBBox())
1565 if eval_points is not None: 1565 ↛ 1566line 1565 didn't jump to line 1566 because the condition on line 1565 was never true
1566 pixel_xy = eval_points.map(np.ravel)
1567 else:
1568 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
1569 if not isinstance(field.bounds, Box): 1569 ↛ 1570line 1569 didn't jump to line 1570 because the condition on line 1569 was never true
1570 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y)
1571 pixel_xy = pixel_xy.map(lambda v: v[mask])
1572 try:
1573 assert_values_equal(
1574 field(x=pixel_xy.x, y=pixel_xy.y),
1575 legacy_field.evaluate(pixel_xy.x, pixel_xy.y),
1576 # The field coefficients come from a numerical fit, so this can
1577 # be bigger than just round-off error epsilon.
1578 rtol=rtol,
1579 )
1580 except AssertionError as err:
1581 err.add_note(f"evaluated at {pixel_xy}")
1582 raise
1583 if not isinstance(legacy_field, LegacyBoundedField): 1583 ↛ 1586line 1583 didn't jump to line 1586 because the condition on line 1583 was never true
1584 # Legacy StitchedApertureCorrection objects are not true BoundedFields
1585 # and don't have addToImage.
1586 return
1587 if eval_points is not None: 1587 ↛ 1588line 1587 didn't jump to line 1588 because the condition on line 1587 was never true
1588 return
1589 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
1590 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
1591 assert_images_equal(
1592 field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13
1593 )
1596def compare_aperture_corrections_to_legacy(
1597 aperture_corrections: Mapping[str, BaseField],
1598 legacy_ap_corr_map: Any,
1599 subimage_bbox: Box,
1600 *,
1601 eval_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
1602 expect_subset_bounds: bool = False,
1603 rtol: float = 1e-5,
1604) -> None:
1605 """Test an aperture correction `dict` by comparing it to an equivalent
1606 `lsst.afw.image.ApCorrMap`.
1608 Parameters
1609 ----------
1610 aperture_corrections
1611 Dictionary to test.
1612 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
1613 Equivalent legacy aperture correction map.
1614 subimage_bbox
1615 Bounding box for full-image tests.
1616 eval_points
1617 Points at which to compare field values, forwarded to
1618 `compare_field_to_legacy`. This also disables image-rendering
1619 comparison checks.
1620 expect_subset_bounds
1621 If `True`, expect ``field.bounds`` to be contained by
1622 ``legacy_field.getBBox()``, not equal to it.
1623 rtol
1624 Relative tolerance for the evaluation-point field-value comparison,
1625 forwarded to `compare_field_to_legacy`.
1626 """
1627 assert aperture_corrections.keys() == set(legacy_ap_corr_map.keys())
1628 for name, field in aperture_corrections.items():
1629 compare_field_to_legacy(
1630 field,
1631 legacy_ap_corr_map[name],
1632 subimage_bbox,
1633 eval_points=eval_points,
1634 expect_subset_bounds=expect_subset_bounds,
1635 rtol=rtol,
1636 )
1639def compare_observation_summary_stats_to_legacy(
1640 summary_stats: ObservationSummaryStats,
1641 legacy_summary_stats: Any,
1642) -> None:
1643 """Test an ObservationSummaryStats object by comparing it to an equivalent
1644 `lsst.afw.image.ExposureSummaryStats`.
1646 Parameters
1647 ----------
1648 summary_stats
1649 Struct to test.
1650 legacy_summary_stats : ``lsst.afw.image.ExposureSummaryStats``
1651 Equivalent legacy struct.
1652 """
1653 for field in dataclasses.fields(legacy_summary_stats):
1654 # Always skip version since that field has no meaning in the new
1655 # type.
1656 if field.name == "version":
1657 continue
1658 a = getattr(legacy_summary_stats, field.name)
1659 b = getattr(summary_stats, field.name)
1660 assert_values_equal(a, b, label=field.name)
1663def compare_sky_projection_to_legacy_wcs[F: Frame](
1664 sky_projection: SkyProjection[F],
1665 legacy_wcs: Any,
1666 pixel_frame: F,
1667 subimage_bbox: Box,
1668 is_fits: bool = False,
1669) -> None:
1670 """Test a Projection object by comparing it to an equivalent
1671 `lsst.afw.geom.SkyWcs`.
1673 Parameters
1674 ----------
1675 sky_projection
1676 Projection to test.
1677 legacy_wcs : ``lsst.afw.geom.SkyWcs``
1678 Equivalent legacy WCS.
1679 pixel_frame
1680 Expected pixel frame for the sky_projection.
1681 subimage_bbox
1682 Bounding box of points to generate for tests.
1683 is_fits
1684 Whether this sky_projection is expected to be exactly representable as
1685 a FITS WCS. If `False` it is assumed to have a FITS approximation
1686 attached instead.
1687 """
1688 # Pixel coordinates to test on over the subimage region of interest:
1689 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
1690 # Array indices of those pixel values (subtract off bbox starts):
1691 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
1692 sky_coords = legacy_coords_to_astropy(
1693 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
1694 )
1695 # Test transforming with the Projection itself, which also tests its
1696 # nested Transform and an Astropy High-Level WCS view with no origin
1697 # change.
1698 check_projection(sky_projection, pixel_xy, sky_coords, pixel_frame, pixel_atol=_AST_INV_ATOL)
1699 # Also test the Astropy High-Level WCS view with an origin change to
1700 # array indices.
1701 check_astropy_wcs_interface(
1702 sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
1703 )
1704 if is_fits: 1704 ↛ 1723line 1704 didn't jump to line 1723 because the condition on line 1704 was always true
1705 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
1706 assert fits_wcs is not None
1707 check_astropy_wcs_interface(fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
1708 # Use that FITS approximation to check that we can make a
1709 # Projection from a FITS WCS, too.
1710 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
1711 check_projection(
1712 fits_projection,
1713 subimage_array_xy,
1714 sky_coords,
1715 pixel_frame,
1716 pixel_atol=1e-5,
1717 )
1718 # We want Projections we create from a FITS WCS to be backed by an
1719 # AST FrameSet so we can convert them into legacy
1720 # `lsst.afw.geom.SkyWcs` objects if desired.
1721 assert "Begin FrameSet" in fits_projection.show()
1722 else:
1723 assert sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False) is None
1724 # The legacy SkyWcs should instead have a FITS approximation
1725 # attached; run the same tests on that.
1726 assert sky_projection.fits_approximation is not None
1727 compare_sky_projection_to_legacy_wcs(
1728 sky_projection.fits_approximation,
1729 legacy_wcs.getFitsApproximation(),
1730 pixel_frame,
1731 subimage_bbox,
1732 is_fits=True,
1733 )
1736def check_transform[I: Frame, O: Frame](
1737 transform: Transform[I, O],
1738 input_xy: XY[np.ndarray],
1739 output_xy: XY[np.ndarray],
1740 in_frame: Frame,
1741 out_frame: Frame,
1742 *,
1743 check_inverted: bool = True,
1744 in_atol: u.Quantity | float | None = None,
1745 out_atol: u.Quantity | float | None = None,
1746 angle_wrapped: bool = False,
1747) -> None:
1748 """Test Transform against known arrays of input and output points.
1750 Parameters
1751 ----------
1752 transform
1753 Transform to test.
1754 input_xy
1755 Arrays of input points.
1756 output_xy
1757 Arrays of output points.
1758 in_frame
1759 Expected input frame.
1760 out_frame
1761 Expected output frame.
1762 check_inverted
1763 If `True`, recurse (once) to test the inverse transform.
1764 in_atol
1765 Expected absolute precision of input points.
1766 out_atol
1767 Expected absolute precision of output points.
1768 angle_wrapped
1769 If `True`, treat comparisons on sky-coordinate axes (as determined
1770 by the frames) as angles modulo 2 pi, so differences that cross
1771 zero or 360 degrees compare correctly.
1772 """
1773 assert transform.in_frame == in_frame
1774 assert transform.out_frame == out_frame
1775 in_atol = 0.0 if in_atol is None else in_atol
1776 out_atol = 0.0 if out_atol is None else out_atol
1777 in_atol_v = in_atol.to_value(in_frame.unit) if isinstance(in_atol, u.Quantity) else float(in_atol)
1778 out_atol_v = out_atol.to_value(out_frame.unit) if isinstance(out_atol, u.Quantity) else float(out_atol)
1779 out_wrapped = angle_wrapped and out_frame == SkyFrame.ICRS
1780 in_wrapped = angle_wrapped and in_frame == SkyFrame.ICRS
1782 def assert_close(
1783 test: np.ndarray | float | u.Quantity,
1784 expected: np.ndarray | float | u.Quantity,
1785 atol: float | u.Quantity,
1786 *,
1787 wrapped: bool,
1788 ) -> None:
1789 """Compare one axis, wrap-aware if that axis is a sky angle.
1791 Parameters
1792 ----------
1793 test
1794 Values to test.
1795 expected
1796 Expected values.
1797 atol
1798 Absolute tolerance, in the units the values are expressed in.
1799 wrapped
1800 If `True`, treat the values as angles modulo 2 pi (radians).
1801 """
1802 if not wrapped:
1803 assert_values_equal(test, expected, atol=atol, rtol=_AST_RTOL)
1804 return
1805 test_values = np.asarray(getattr(test, "value", test))
1806 expected_values = np.asarray(getattr(expected, "value", expected))
1807 atol_rad = atol.to_value(u.rad) if isinstance(atol, u.Quantity) else float(atol)
1808 assert_values_equal(
1809 test_values * u.rad,
1810 expected_values * u.rad,
1811 atol=atol_rad * u.rad,
1812 equal_nan=False,
1813 wrap_angles=True,
1814 )
1816 # Test array interfaces.
1817 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
1818 assert_close(test_output_xy.x, output_xy.x, out_atol_v, wrapped=out_wrapped)
1819 assert_close(test_output_xy.y, output_xy.y, out_atol_v, wrapped=out_wrapped)
1820 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
1821 assert_close(test_input_xy.x, input_xy.x, in_atol_v, wrapped=in_wrapped)
1822 assert_close(test_input_xy.y, input_xy.y, in_atol_v, wrapped=in_wrapped)
1823 # Test scalar interfaces with numpy scalars.
1824 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
1825 forward_xy = transform.apply_forward(x=input_x, y=input_y)
1826 assert_close(forward_xy.x, output_x, out_atol_v, wrapped=out_wrapped)
1827 assert_close(forward_xy.y, output_y, out_atol_v, wrapped=out_wrapped)
1828 inverse_xy = transform.apply_inverse(x=output_x, y=output_y)
1829 assert_close(inverse_xy.x, input_x, in_atol_v, wrapped=in_wrapped)
1830 assert_close(inverse_xy.y, input_y, in_atol_v, wrapped=in_wrapped)
1831 # Test quantity array interfaces.
1832 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
1833 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
1834 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
1835 assert_close(test_output_q_xy.x, output_q_xy.x, out_atol, wrapped=out_wrapped)
1836 assert_close(test_output_q_xy.y, output_q_xy.y, out_atol, wrapped=out_wrapped)
1837 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
1838 assert_close(test_input_q_xy.x, input_q_xy.x, in_atol, wrapped=in_wrapped)
1839 assert_close(test_input_q_xy.y, input_q_xy.y, in_atol, wrapped=in_wrapped)
1840 # Test quantity scalar interfaces.
1841 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
1842 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
1843 ):
1844 forward_q_xy = transform.apply_forward_q(x=input_q_x, y=input_q_y)
1845 assert_close(forward_q_xy.x, output_q_x, out_atol, wrapped=out_wrapped)
1846 assert_close(forward_q_xy.y, output_q_y, out_atol, wrapped=out_wrapped)
1847 inverse_q_xy = transform.apply_inverse_q(x=output_q_x, y=output_q_y)
1848 assert_close(inverse_q_xy.x, input_q_x, in_atol, wrapped=in_wrapped)
1849 assert_close(inverse_q_xy.y, input_q_y, in_atol, wrapped=in_wrapped)
1850 if check_inverted:
1851 # Test the inverse transform.
1852 check_transform(
1853 transform.inverted(),
1854 output_xy,
1855 input_xy,
1856 out_frame,
1857 in_frame,
1858 check_inverted=False,
1859 out_atol=in_atol,
1860 in_atol=out_atol,
1861 angle_wrapped=angle_wrapped,
1862 )
1865def check_projection[P: Frame](
1866 sky_projection: SkyProjection[P],
1867 pixel_xy: XY[np.ndarray],
1868 sky_coords: SkyCoord,
1869 pixel_frame: Frame,
1870 *,
1871 pixel_atol: float | None = None,
1872 sky_atol: u.Quantity | None = None,
1873) -> None:
1874 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1875 coordinates.
1877 Parameters
1878 ----------
1879 sky_projection
1880 Projection to test.
1881 pixel_xy
1882 Arrays of pixel coordinates.
1883 sky_coords
1884 Corresponding sky coordinates.
1885 pixel_frame
1886 Expected pixel frame.
1887 pixel_atol
1888 Expected absolute precision of pixel points.
1889 sky_atol
1890 Expected absolute precision of sky coordinates, as a maximum
1891 great-circle separation.
1892 """
1893 assert sky_projection.pixel_frame == pixel_frame
1894 assert sky_projection.sky_frame == SkyFrame.ICRS
1895 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol)
1896 pixel_atol = 0.0 if pixel_atol is None else pixel_atol
1897 pixel_atol_q = pixel_atol * u.pix
1898 # Test array interfaces.
1899 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1900 assert_values_equal(test_pixel_xy.x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL)
1901 assert_values_equal(test_pixel_xy.y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL)
1902 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1903 assert_values_equal(test_sky_astropy, sky_coords, atol=sky_atol, label="pixel_to_sky")
1904 # Test scalar interfaces.
1905 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1906 assert_values_equal(
1907 sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol, rtol=_AST_RTOL
1908 )
1909 assert_values_equal(
1910 sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol, rtol=_AST_RTOL
1911 )
1912 test_sky_single = sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y)
1913 assert_values_equal(test_sky_single, sky_single, atol=sky_atol, label="pixel_to_sky")
1914 # Test the underlying Transform object.
1915 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1916 check_transform(
1917 sky_projection.pixel_to_sky_transform,
1918 pixel_xy,
1919 sky_xy,
1920 pixel_frame,
1921 SkyFrame.ICRS,
1922 check_inverted=False,
1923 in_atol=pixel_atol_q,
1924 out_atol=sky_atol,
1925 angle_wrapped=True,
1926 )
1927 check_transform(
1928 sky_projection.sky_to_pixel_transform,
1929 sky_xy,
1930 pixel_xy,
1931 SkyFrame.ICRS,
1932 pixel_frame,
1933 check_inverted=False,
1934 in_atol=sky_atol,
1935 out_atol=pixel_atol_q,
1936 angle_wrapped=True,
1937 )
1938 # Test the Astropy interface adapter.
1939 check_astropy_wcs_interface(
1940 sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1941 )
1944def assert_sky_projections_equal(
1945 a: SkyProjection[Any] | None,
1946 b: SkyProjection[Any] | None,
1947 expect_identity: bool | None = None,
1948) -> None:
1949 """Test that two `.SkyProjection` instances are equivalent.
1951 Parameters
1952 ----------
1953 a
1954 First sky projection to compare.
1955 b
1956 Second sky projection to compare.
1957 expect_identity
1958 If not `None`, assert whether ``a`` and ``b`` are the same object.
1959 """
1960 if a is None and b is None:
1961 return
1962 assert a is not None
1963 assert b is not None
1964 match expect_identity:
1965 case True: 1965 ↛ 1966line 1965 didn't jump to line 1966 because the pattern on line 1965 never matched
1966 assert a is b
1967 return
1968 case False:
1969 assert a is not b
1970 case None if a is b:
1971 return
1972 assert a == b
1975def check_astropy_wcs_interface(
1976 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1977 pixel_xy: XY[np.ndarray],
1978 sky_coords: SkyCoord,
1979 *,
1980 pixel_atol: float | None = None,
1981 sky_atol: u.Quantity | None = None,
1982) -> None:
1983 """Test an Astropy WCS instance against known arrays of pixel and
1984 sky coordinates.
1986 Parameters
1987 ----------
1988 wcs
1989 Astropy WCS object to test.
1990 pixel_xy
1991 Arrays of pixel coordinates.
1992 sky_coords
1993 Corresponding sky coordinates.
1994 pixel_atol
1995 Expected absolute precision of pixel points.
1996 sky_atol
1997 Expected absolute precision of sky coordinates, as a maximum
1998 great-circle separation.
1999 """
2000 pixel_atol = 0.0 if pixel_atol is None else pixel_atol
2001 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol)
2002 test_x, test_y = wcs.world_to_pixel(sky_coords)
2003 assert_values_equal(test_x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL)
2004 assert_values_equal(test_y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL)
2005 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
2006 assert_values_equal(test_sky_coords, sky_coords, atol=sky_atol, label="pixel_to_world")
2009def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
2010 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array.
2012 Parameters
2013 ----------
2014 legacy_points
2015 Legacy ``lsst.geom.Point2D`` objects to convert.
2016 """
2017 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
2020def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
2021 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
2022 coordinate object.
2024 Parameters
2025 ----------
2026 legacy_coords
2027 Legacy ``lsst.geom.SpherePoint`` objects to convert.
2028 """
2029 return SkyCoord(
2030 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
2031 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
2032 )
2035def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
2036 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``.
2038 Parameters
2039 ----------
2040 x
2041 X coordinates of the points.
2042 y
2043 Y coordinates of the points.
2044 """
2045 from lsst.geom import Point2D
2047 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
2050def compare_amplifier_to_legacy(
2051 amplifier: Amplifier,
2052 legacy_amplifier: Any,
2053 *,
2054 is_raw_assembled: bool,
2055 expect_nominal_calibrations: bool = True,
2056) -> None:
2057 """Compare an `~.cameras.Amplifier` to a legacy
2058 `lsst.afw.cameraGeom.Amplifier`.
2060 Parameters
2061 ----------
2062 amplifier
2063 Amplifier to compare.
2064 legacy_amplifier
2065 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against.
2066 is_raw_assembled
2067 Whether the raw geometry is expected to be the assembled-raw
2068 geometry (`True`) or the unassembled-raw geometry (`False`).
2069 expect_nominal_calibrations
2070 Whether the amplifier is expected to carry nominal calibrations.
2071 """
2072 assert legacy_amplifier.getName() == amplifier.name
2073 assert Box.from_legacy(legacy_amplifier.getBBox()) == amplifier.bbox
2074 if is_raw_assembled: 2074 ↛ 2077line 2074 didn't jump to line 2077 because the condition on line 2074 was always true
2075 raw_geom = amplifier.assembled_raw_geometry
2076 else:
2077 raw_geom = amplifier.unassembled_raw_geometry
2078 assert raw_geom is not None
2079 assert ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()) == raw_geom.readout_corner
2080 assert Box.from_legacy(legacy_amplifier.getRawBBox()) == raw_geom.bbox
2081 assert Box.from_legacy(legacy_amplifier.getRawDataBBox()) == raw_geom.data_bbox
2082 assert legacy_amplifier.getRawFlipX() == raw_geom.flip_x
2083 assert legacy_amplifier.getRawFlipY() == raw_geom.flip_y
2084 assert legacy_amplifier.getRawXYOffset().getX() == raw_geom.x_offset
2085 assert legacy_amplifier.getRawXYOffset().getY() == raw_geom.y_offset
2086 assert (
2087 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()) == raw_geom.horizontal_overscan_bbox
2088 )
2089 assert Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()) == raw_geom.vertical_overscan_bbox
2090 if raw_geom.horizontal_prescan_bbox is None:
2091 # `Box` requires a positive size, so an amplifier with no prescan
2092 # region carries `None` here.
2093 assert legacy_amplifier.getRawPrescanBBox().isEmpty()
2094 else:
2095 assert Box.from_legacy(legacy_amplifier.getRawPrescanBBox()) == raw_geom.horizontal_prescan_bbox
2096 if expect_nominal_calibrations: 2096 ↛ exitline 2096 didn't return from function 'compare_amplifier_to_legacy' because the condition on line 2096 was always true
2097 assert amplifier.nominal_calibrations is not None
2098 assert_values_equal(
2099 legacy_amplifier.getGain(),
2100 amplifier.nominal_calibrations.gain,
2101 equal_nan=True,
2102 label="gain",
2103 )
2104 assert_values_equal(
2105 legacy_amplifier.getReadNoise(),
2106 amplifier.nominal_calibrations.read_noise,
2107 equal_nan=True,
2108 label="read_noise",
2109 )
2110 assert_values_equal(
2111 legacy_amplifier.getSaturation(),
2112 amplifier.nominal_calibrations.saturation,
2113 equal_nan=True,
2114 label="saturation",
2115 )
2116 assert_values_equal(
2117 legacy_amplifier.getSuspectLevel(),
2118 amplifier.nominal_calibrations.suspect_level,
2119 equal_nan=True,
2120 label="suspect_level",
2121 )
2122 assert_values_equal(
2123 legacy_amplifier.getLinearityCoeffs(),
2124 amplifier.nominal_calibrations.linearity_coefficients,
2125 label="linearity_coefficients",
2126 )
2127 assert legacy_amplifier.getLinearityType() == amplifier.nominal_calibrations.linearity_type
2130def compare_detector_to_legacy(
2131 detector: Detector,
2132 legacy_detector: Any,
2133 *,
2134 is_raw_assembled: bool,
2135 expect_nominal_calibrations: bool = True,
2136) -> None:
2137 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`.
2139 Parameters
2140 ----------
2141 detector
2142 Detector to compare.
2143 legacy_detector
2144 Legacy `lsst.afw.cameraGeom.Detector` to compare against.
2145 is_raw_assembled
2146 Whether the raw geometry is expected to be the assembled-raw
2147 geometry (`True`) or the unassembled-raw geometry (`False`).
2148 expect_nominal_calibrations
2149 Whether the detector's amplifiers are expected to carry nominal
2150 calibrations.
2151 """
2152 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
2154 assert legacy_detector.getName() == detector.name
2155 assert legacy_detector.getId() == detector.id
2156 assert DetectorType.from_legacy(legacy_detector.getType()) == detector.type
2157 assert Box.from_legacy(legacy_detector.getBBox()) == detector.bbox
2158 assert legacy_detector.getSerial() == detector.serial
2159 legacy_orientation = legacy_detector.getOrientation()
2160 assert legacy_orientation.getFpPosition3().getX() == detector.orientation.focal_plane_x
2161 assert legacy_orientation.getFpPosition3().getY() == detector.orientation.focal_plane_y
2162 assert legacy_orientation.getFpPosition3().getZ() == detector.orientation.focal_plane_z
2163 assert legacy_orientation.getReferencePoint().getX() == detector.orientation.pixel_reference_x
2164 assert legacy_orientation.getReferencePoint().getY() == detector.orientation.pixel_reference_y
2165 assert legacy_orientation.getYaw().asRadians() == detector.orientation.yaw.to_value(u.rad)
2166 assert legacy_orientation.getPitch().asRadians() == detector.orientation.pitch.to_value(u.rad)
2167 assert legacy_orientation.getRoll().asRadians() == detector.orientation.roll.to_value(u.rad)
2168 assert legacy_detector.getPixelSize().getX() == detector.pixel_size
2169 assert legacy_detector.getPhysicalType() == detector.physical_type
2170 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
2171 compare_amplifier_to_legacy(
2172 amplifier,
2173 legacy_amplifier,
2174 is_raw_assembled=is_raw_assembled,
2175 expect_nominal_calibrations=expect_nominal_calibrations,
2176 )
2177 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
2178 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
2179 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
2180 check_transform(
2181 detector.to_focal_plane,
2182 pixel_xy,
2183 legacy_points_to_xy_array(fp_legacy_points),
2184 detector.frame,
2185 detector.to_focal_plane.out_frame,
2186 in_atol=1e-9 * u.pix,
2187 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
2188 )
2189 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
2190 check_transform(
2191 detector.to_field_angle,
2192 pixel_xy,
2193 legacy_points_to_xy_array(fa_legacy_points),
2194 detector.frame,
2195 detector.to_field_angle.out_frame,
2196 in_atol=1e-9 * u.pix,
2197 out_atol=1e-7 * u.arcsec,
2198 )
2201def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
2202 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
2204 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
2205 classes. Importing this module already imports every ``lsst.images``
2206 module that defines a subclass, so the tree is fully populated by the time
2207 this is called.
2209 This discovery is deliberately separate from
2210 `check_archive_tree_class_invariants` so that the per-class check stays
2211 usable on a single class even if this metaprogramming is removed later.
2212 """
2213 seen: set[type] = set()
2214 stack: list[type] = [ArchiveTree]
2215 while stack:
2216 kls = stack.pop()
2217 sub: type
2218 for sub in kls.__subclasses__():
2219 if sub in seen: 2219 ↛ 2220line 2219 didn't jump to line 2220 because the condition on line 2219 was never true
2220 continue
2221 seen.add(sub)
2222 stack.append(sub)
2223 if not getattr(sub, "__abstractmethods__", None): 2223 ↛ 2218line 2223 didn't jump to line 2218 because the condition on line 2223 was always true
2224 yield sub
2227def check_bounds_contains_broadcasting(bounds: Bounds) -> None:
2228 """Verify that `~lsst.images.Bounds.contains` accepts array-like inputs.
2230 Uses the scalar overload as the reference and checks that 1-D arrays,
2231 list inputs, mixed scalar-plus-array inputs, 2-D broadcast inputs, and
2232 `.XY` / `.YX` positional-argument forms all produce results consistent
2233 with calling the scalar overload on each ``(x, y)`` pair individually.
2235 Parameters
2236 ----------
2237 bounds
2238 The `~lsst.images.Bounds` implementation to exercise.
2239 """
2240 bbox = bounds.bbox
2241 # One point outside each boundary, one on each boundary, one in the
2242 # interior. The boundary points (start and stop) straddle inside/outside
2243 # (start is inside, stop is outside) so we are guaranteed a mix of True
2244 # and False results without hard-coding expected values.
2245 ys = np.array(
2246 [
2247 bbox.y.start - 1,
2248 bbox.y.start,
2249 (bbox.y.start + bbox.y.stop) // 2,
2250 bbox.y.stop,
2251 ]
2252 )
2253 xs = np.array(
2254 [
2255 bbox.x.start - 1,
2256 bbox.x.start,
2257 (bbox.x.start + bbox.x.stop) // 2,
2258 bbox.x.stop,
2259 ]
2260 )
2261 # 2-D reference: expected[i, j] == contains(x=xs[j], y=ys[i]).
2262 expected = np.array(
2263 [[bounds.contains(x=int(xi), y=int(yi)) for xi in xs] for yi in ys],
2264 dtype=bool,
2265 )
2266 # 1-D ndarray (diagonal pairs).
2267 assert_values_equal(
2268 bounds.contains(x=xs, y=ys),
2269 np.diagonal(expected),
2270 label="bounds.contains",
2271 )
2272 # list inputs (array-like).
2273 assert_values_equal(
2274 bounds.contains(x=xs.tolist(), y=ys.tolist()),
2275 np.diagonal(expected),
2276 label="bounds.contains",
2277 )
2278 # Mixed: scalar y, 1-D array x — produces a 1-D result.
2279 fixed_yi = 1 # index into ys; ys[1] == bbox.y.start (on the boundary)
2280 assert_values_equal(
2281 bounds.contains(x=xs, y=int(ys[fixed_yi])),
2282 expected[fixed_yi, :],
2283 label="bounds.contains",
2284 )
2285 # Mixed: 1-D array y, scalar x.
2286 fixed_xi = 1 # index into xs; xs[1] == bbox.x.start (on the boundary)
2287 assert_values_equal(
2288 bounds.contains(x=int(xs[fixed_xi]), y=ys),
2289 expected[:, fixed_xi],
2290 label="bounds.contains",
2291 )
2292 # Float scalars: results must match the int-scalar reference.
2293 assert bounds.contains(x=float(xs[fixed_xi]), y=float(ys[fixed_yi])) == expected[fixed_yi, fixed_xi]
2294 # XY / YX scalar: results must match the keyword-scalar reference.
2295 assert bounds.contains(XY(x=int(xs[fixed_xi]), y=int(ys[fixed_yi]))) == expected[fixed_yi, fixed_xi]
2296 assert bounds.contains(YX(y=int(ys[fixed_yi]), x=int(xs[fixed_xi]))) == expected[fixed_yi, fixed_xi]
2297 # XY / YX array: results must match the 1-D ndarray reference.
2298 assert_values_equal(bounds.contains(XY(x=xs, y=ys)), np.diagonal(expected), label="bounds.contains")
2299 assert_values_equal(bounds.contains(YX(y=ys, x=xs)), np.diagonal(expected), label="bounds.contains")
2300 # 2-D broadcast: y shape (N, 1) × x shape (1, M) → (N, M).
2301 assert_values_equal(
2302 bounds.contains(x=xs.reshape(1, -1), y=ys.reshape(-1, 1)),
2303 expected,
2304 label="bounds.contains",
2305 )
2308def check_archive_tree_class_invariants(tree_type: type[ArchiveTree]) -> None:
2309 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
2310 well-formed schema-version constants and an in-memory type.
2312 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
2313 ``PUBLIC_TYPE`` are present and well-typed, that the version is
2314 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
2315 schema major.
2317 Parameters
2318 ----------
2319 tree_type
2320 The concrete `.serialization.ArchiveTree` subclass to check.
2321 """
2322 assert hasattr(tree_type, "SCHEMA_NAME"), f"{tree_type.__name__} lacks SCHEMA_NAME"
2323 assert hasattr(tree_type, "SCHEMA_VERSION"), f"{tree_type.__name__} lacks SCHEMA_VERSION"
2324 assert hasattr(tree_type, "MIN_READ_VERSION"), f"{tree_type.__name__} lacks MIN_READ_VERSION"
2325 assert hasattr(tree_type, "PUBLIC_TYPE"), f"{tree_type.__name__} lacks PUBLIC_TYPE"
2326 assert isinstance(tree_type.SCHEMA_NAME, str)
2327 assert len(tree_type.SCHEMA_NAME) > 0
2328 # Allow an optional PEP 440 development-release suffix (e.g. 1.0.0.dev0)
2329 # for schemas still in development.
2330 assert re.fullmatch(r"^\d+\.\d+\.\d+(\.dev\d+)?$", tree_type.SCHEMA_VERSION)
2331 assert isinstance(tree_type.MIN_READ_VERSION, int)
2332 assert tree_type.MIN_READ_VERSION >= 1
2333 assert isinstance(tree_type.PUBLIC_TYPE, type)
2334 major = int(tree_type.SCHEMA_VERSION.split(".")[0])
2335 assert tree_type.MIN_READ_VERSION <= major
2338def _angle_to_astropy(angle: Any) -> u.Quantity:
2339 """Convert an `lsst.geom.Angle` to an Astropy angular quantity."""
2340 return u.Quantity(angle.asRadians(), u.rad)
2343def _observatory_to_location(observatory: Any) -> coords.EarthLocation:
2344 """Convert an `lsst.afw.coord.Observatory` to an Astropy
2345 `~astropy.coordinates.EarthLocation`.
2346 """
2347 return coords.EarthLocation.from_geodetic(
2348 lon=observatory.getLongitude().asRadians() * u.rad,
2349 lat=observatory.getLatitude().asRadians() * u.rad,
2350 height=observatory.getElevation() * u.m,
2351 )
2354def assert_obs_metadata_fields_equal(actual: Any, expected: Any, *, label: str = "") -> None:
2355 """Assert that two ObservationInfo or VisitInfo field values are equal.
2357 Dispatches on the type of ``expected`` so that values affected by afw's
2358 float64 date storage, unit conversions, and angle wrapping are compared
2359 with appropriate tolerances; anything else is compared with ``==``.
2361 Parameters
2362 ----------
2363 actual : `object`
2364 The field value to check.
2365 expected : `object`
2366 The field value to check it against.
2367 label : `str`, optional
2368 Prefix for the failure message.
2369 """
2370 from lsst.afw.coord import Observatory, Weather
2371 from lsst.daf.base import DateTime
2372 from lsst.geom import Angle
2374 if isinstance(expected, DateTime):
2375 assert_obs_metadata_fields_equal(actual.toAstropy(), expected.toAstropy(), label=label)
2376 elif isinstance(expected, astropy.time.Time):
2377 # afw stores VisitInfo dates as float64 MJDs, which limits precision.
2378 assert_values_equal(actual, expected, atol=1e-9 * u.s, label=label)
2379 elif isinstance(expected, coords.AltAz):
2380 assert_values_equal(actual.az, expected.az, wrap_angles=True, label=f"{label}.az")
2381 assert_values_equal(actual.alt, expected.alt, atol=1e-3 * u.arcsec, label=f"{label}.alt")
2382 elif isinstance(expected, coords.SkyCoord):
2383 assert_values_equal(actual, expected, atol=1e-3 * u.arcsec, label=label)
2384 elif isinstance(expected, coords.EarthLocation):
2385 assert_values_equal(
2386 u.Quantity([actual.x, actual.y, actual.z]),
2387 u.Quantity([expected.x, expected.y, expected.z]),
2388 rtol=1e-9,
2389 label=label,
2390 )
2391 elif isinstance(expected, Observatory):
2392 assert_obs_metadata_fields_equal(
2393 _observatory_to_location(actual), _observatory_to_location(expected), label=label
2394 )
2395 elif isinstance(expected, Weather):
2396 assert_obs_metadata_fields_equal(
2397 u.Quantity(actual.getAirTemperature(), u.deg_C),
2398 u.Quantity(expected.getAirTemperature(), u.deg_C),
2399 label=f"{label}.temperature",
2400 )
2401 assert_obs_metadata_fields_equal(
2402 actual.getAirPressure(), expected.getAirPressure(), label=f"{label}.pressure"
2403 )
2404 assert_obs_metadata_fields_equal(
2405 actual.getHumidity(), expected.getHumidity(), label=f"{label}.humidity"
2406 )
2407 elif isinstance(expected, Angle):
2408 assert_values_equal(
2409 _angle_to_astropy(actual),
2410 _angle_to_astropy(expected),
2411 atol=1e-12 * u.rad,
2412 wrap_angles=True,
2413 label=label,
2414 )
2415 elif (
2416 isinstance(expected, u.Quantity)
2417 and isinstance(actual, u.Quantity)
2418 and expected.unit.physical_type == "temperature"
2419 and actual.unit.physical_type == "temperature"
2420 ):
2421 assert_values_equal(
2422 actual.to(u.K, equivalencies=u.temperature()),
2423 expected.to(u.K, equivalencies=u.temperature()),
2424 atol=1e-6 * u.K,
2425 label=label,
2426 )
2427 elif isinstance(expected, u.Quantity | float):
2428 assert_values_equal(actual, expected, rtol=1e-9, label=label)
2429 elif actual != expected:
2430 raise AssertionError(f"{label}: {actual!r} != {expected!r}")