Coverage for python/lsst/images/tests/_checks.py: 69%
674 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 19:06 +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``. Schema planes the map does not name
731 are compared under their own names. A plane the legacy mask does
732 not define must have no pixels set in ``mask``.
733 label
734 Prefix for the failure message; each plane is reported as
735 ``"{label}[<plane>]"``.
736 """
737 assert mask.bbox == Box.from_legacy(legacy_mask.getBBox())
738 if plane_map is None:
739 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
740 legacy_planes = legacy_mask.getMaskPlaneDict()
741 # Compare schema planes the map does not name (e.g. the planes of
742 # `get_legacy_optional_mask_planes`) under their own names.
743 plane_map = dict(plane_map)
744 mapped = {plane.name for plane in plane_map.values()}
745 plane_map.update(
746 {plane.name: plane for plane in mask.schema if plane is not None and plane.name not in mapped}
747 )
748 try:
749 for old_name, new_plane in plane_map.items():
750 pixels = mask.get(new_plane.name)
751 if old_name not in legacy_planes: 751 ↛ 753line 751 didn't jump to line 753 because the condition on line 751 was never true
752 # A plane the legacy mask does not define must be empty.
753 assert_values_equal(np.zeros_like(pixels), pixels, label=f"{label}[{old_name}]")
754 continue
755 assert_values_equal(
756 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
757 pixels,
758 label=f"{label}[{old_name}]",
759 )
760 except AssertionError as err:
761 names = {new_plane.name: old_name for old_name, new_plane in plane_map.items()}
762 _note_mask_difference(err, mask, Mask.from_legacy(legacy_mask, plane_map), label=label, names=names)
763 raise
766def compare_masked_image_to_legacy(
767 masked_image: MaskedImage,
768 legacy_masked_image: Any,
769 *,
770 plane_map: Mapping[str, MaskPlane] | None = None,
771 expect_view: bool | None = None,
772 alternates: Mapping[str, Any] | None = None,
773) -> None:
774 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
775 object.
777 Parameters
778 ----------
779 masked_image
780 New image to test.
781 legacy_masked_image
782 Legacy image to test against.
783 plane_map
784 Mapping between new and legacy mask planes.
785 expect_view
786 Whether to test that the image and variance arrays do or do not share
787 memory.
788 alternates
789 A mapping of other versions of one or more (new) components to also
790 check against the legacy versions of those components.
791 """
792 compare_image_to_legacy(
793 masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view, label="image"
794 )
795 compare_mask_to_legacy(
796 masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask"
797 )
798 compare_image_to_legacy(
799 masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view, label="variance"
800 )
801 if alternates: 801 ↛ 802line 801 didn't jump to line 802 because the condition on line 801 was never true
802 if image := alternates.get("image"):
803 compare_image_to_legacy(
804 image, legacy_masked_image.getImage(), expect_view=expect_view, label="image"
805 )
806 if mask := alternates.get("mask"):
807 compare_mask_to_legacy(mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask")
808 if variance := alternates.get("variance"):
809 compare_image_to_legacy(
810 variance, legacy_masked_image.getVariance(), expect_view=expect_view, label="variance"
811 )
814def compare_visit_image_to_legacy(
815 visit_image: VisitImage,
816 legacy_exposure: Any,
817 *,
818 plane_map: Mapping[str, MaskPlane] | None = None,
819 expect_view: bool | None = None,
820 instrument: str,
821 visit: int,
822 detector: int,
823 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
824 alternates: Mapping[str, Any] | None = None,
825 check_photometric_scaling: bool = True,
826) -> None:
827 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
828 object.
830 Parameters
831 ----------
832 visit_image
833 New image to test.
834 legacy_exposure
835 Legacy image to test against.
836 plane_map
837 Mapping between new and legacy mask planes; defaults to the planes
838 used by the legacy (`lsst.afw.image.Exposure`) visit image.
839 expect_view
840 Whether to test that the image and variance arrays do or do not share
841 memory.
842 instrument
843 Expected instrument name.
844 visit
845 Expected visit ID.
846 detector
847 Expected detector ID.
848 applied_legacy_photo_calib
849 Legacy `lsst.afw.image.PhotoCalib` already applied to
850 ``legacy_exposure``, used when comparing photometric scaling.
851 alternates
852 A mapping of other versions of one or more (new) components to also
853 check against the legacy versions of those components.
854 check_photometric_scaling
855 If `False`, skip the photometric-scaling comparison.
856 """
857 if plane_map is None:
858 plane_map = get_legacy_visit_image_mask_planes()
859 compare_masked_image_to_legacy(
860 visit_image,
861 legacy_exposure.getMaskedImage(),
862 plane_map=plane_map,
863 expect_view=expect_view,
864 alternates=alternates,
865 )
866 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
867 with annotate_errors("sky_projection"):
868 compare_sky_projection_to_legacy_wcs(
869 visit_image.sky_projection,
870 legacy_exposure.getWcs(),
871 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
872 visit_image.bbox,
873 )
874 assert visit_image.sky_projection is visit_image.mask.sky_projection
875 assert visit_image.sky_projection is visit_image.variance.sky_projection
876 with annotate_errors("psf"):
877 compare_psf_to_legacy(visit_image.psf, legacy_exposure.getPsf())
878 with annotate_errors("summary_stats"):
879 compare_observation_summary_stats_to_legacy(
880 visit_image.summary_stats, legacy_exposure.info.getSummaryStats()
881 )
882 with annotate_errors("detector"):
883 compare_detector_to_legacy(visit_image.detector, legacy_exposure.getDetector(), is_raw_assembled=True)
884 # Make a tiny box for Field comparisons that need to make arrays; that can
885 # get expensive otherwise.
886 tiny_bbox = detector_bbox.local[2:4, 3:6]
887 with annotate_errors("aperture_corrections"):
888 compare_aperture_corrections_to_legacy(
889 visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
890 )
891 if check_photometric_scaling:
892 with annotate_errors("photometric_scaling"):
893 compare_photo_calib_to_legacy(
894 visit_image.photometric_scaling,
895 legacy_exposure.info.getPhotoCalib(),
896 applied_legacy_photo_calib=applied_legacy_photo_calib,
897 subimage_bbox=tiny_bbox,
898 )
899 if alternates:
900 if (bbox := alternates.get("bbox")) is not None:
901 assert bbox == visit_image.bbox
902 if sky_projection := alternates.get("sky_projection"):
903 with annotate_errors("sky_projection"):
904 compare_sky_projection_to_legacy_wcs(
905 sky_projection,
906 legacy_exposure.getWcs(),
907 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
908 visit_image.bbox,
909 )
910 if psf := alternates.get("psf"):
911 with annotate_errors("psf"):
912 compare_psf_to_legacy(psf, legacy_exposure.getPsf())
913 if summary_stats := alternates.get("summary_stats"):
914 with annotate_errors("summary_stats"):
915 compare_observation_summary_stats_to_legacy(
916 summary_stats, legacy_exposure.info.getSummaryStats()
917 )
918 if detector_obj := alternates.get("detector"):
919 with annotate_errors("detector"):
920 compare_detector_to_legacy(detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True)
921 if obs_info := alternates.get("obs_info"):
922 visitInfo = legacy_exposure.visitInfo
923 assert obs_info.instrument == visitInfo.getInstrumentLabel()
924 if aperture_corrections := alternates.get("aperture_corrections"):
925 with annotate_errors("aperture_corrections"):
926 compare_aperture_corrections_to_legacy(
927 aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
928 )
929 if check_photometric_scaling:
930 if (photometric_scaling := alternates.get("photometric_scaling", ...)) is not ...:
931 with annotate_errors("photometric_scaling"):
932 compare_photo_calib_to_legacy(
933 photometric_scaling,
934 legacy_exposure.info.getPhotoCalib(),
935 applied_legacy_photo_calib=applied_legacy_photo_calib,
936 subimage_bbox=tiny_bbox,
937 )
940def compare_photo_calib_to_legacy(
941 photometric_scaling: BaseField | None,
942 legacy_photo_calib: LegacyPhotoCalib,
943 *,
944 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
945 subimage_bbox: Box,
946) -> None:
947 if legacy_photo_calib._isConstant:
948 if legacy_photo_calib.getCalibrationMean() == 1.0:
949 if applied_legacy_photo_calib is None:
950 assert photometric_scaling is None
951 return
952 else:
953 legacy_photo_calib = applied_legacy_photo_calib
954 if legacy_photo_calib._isConstant:
955 assert isinstance(photometric_scaling, ChebyshevField)
956 assert_values_equal(
957 photometric_scaling.coefficients,
958 np.array([[legacy_photo_calib.getCalibrationMean()]]),
959 # The coefficients come from a numerical fit, so this can be
960 # bigger than just round-off error epislon.
961 rtol=1e-5,
962 )
963 else:
964 assert photometric_scaling is not None
965 compare_field_to_legacy(
966 photometric_scaling / legacy_photo_calib.getCalibrationMean(),
967 legacy_photo_calib.computeScaledCalibration(),
968 subimage_bbox,
969 )
972def compare_cell_coadd_to_legacy(
973 cell_coadd: CellCoadd,
974 legacy_cell_coadd: MultipleCellCoadd,
975 *,
976 tract_bbox: Box,
977 plane_map: Mapping[str, MaskPlane] | None = None,
978 alternates: Mapping[str, Any] | None = None,
979 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
980) -> None:
981 """Compare a `.cells.CellCoadd` object to a legacy
982 `lsst.cell_coadds.MultipleCellCoadd` object.
984 Parameters
985 ----------
986 cell_coadd
987 New coadd to test.
988 legacy_cell_coadd
989 Legacy coadd to test against.
990 tract_bbox
991 Bounding box of the full tract.
992 plane_map
993 Mapping between new and legacy mask planes.
994 alternates
995 A mapping of other versions of one or more (new) components to also
996 check against the legacy versions of those components.
997 psf_points
998 Points to use to compare the PSFs.
999 """
1000 legacy_stitched = legacy_cell_coadd.stitch(cell_coadd.bbox.to_legacy())
1001 compare_image_to_legacy(cell_coadd.image, legacy_stitched.image, expect_view=False)
1002 compare_mask_to_legacy(cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map)
1003 compare_image_to_legacy(cell_coadd.variance, legacy_stitched.variance, expect_view=False)
1004 if legacy_stitched.mask_fractions is not None:
1005 compare_image_to_legacy(
1006 cell_coadd.mask_fractions["rejected"], legacy_stitched.mask_fractions, expect_view=False
1007 )
1008 for n in range(legacy_stitched.n_noise_realizations):
1009 compare_image_to_legacy(
1010 cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False
1011 )
1012 assert cell_coadd.skymap == legacy_stitched.identifiers.skymap
1013 assert cell_coadd.tract == legacy_stitched.identifiers.tract
1014 assert cell_coadd.patch.index.x == legacy_stitched.identifiers.patch.x
1015 assert cell_coadd.patch.index.y == legacy_stitched.identifiers.patch.y
1016 assert cell_coadd.band == legacy_stitched.identifiers.band
1017 assert tract_bbox.contains(cell_coadd.patch.outer_bbox)
1018 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox)
1019 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox)
1020 assert cell_coadd.unit == u.Unit(legacy_cell_coadd.common.units.value)
1021 assert cell_coadd.bounds.bbox.contains(cell_coadd.bbox)
1022 assert cell_coadd.grid.bbox.contains(cell_coadd.bbox)
1023 compare_sky_projection_to_legacy_wcs(
1024 cell_coadd.sky_projection,
1025 legacy_cell_coadd.wcs,
1026 TractFrame(
1027 skymap=legacy_cell_coadd.identifiers.skymap,
1028 tract=legacy_cell_coadd.identifiers.tract,
1029 bbox=tract_bbox,
1030 ),
1031 cell_coadd.bbox,
1032 is_fits=True,
1033 )
1034 assert cell_coadd.sky_projection is cell_coadd.mask.sky_projection
1035 assert cell_coadd.sky_projection is cell_coadd.variance.sky_projection
1036 compare_psf_to_legacy(
1037 cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points
1038 )
1039 compare_aperture_corrections_to_legacy(
1040 cell_coadd.aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
1041 )
1042 compare_cell_coadd_provenance_to_legacy(cell_coadd.provenance, legacy_cell_coadd)
1043 if alternates:
1044 if sky_projection := alternates.get("sky_projection"):
1045 compare_sky_projection_to_legacy_wcs(
1046 sky_projection,
1047 legacy_stitched.wcs,
1048 TractFrame(
1049 skymap=legacy_cell_coadd.identifiers.skymap,
1050 tract=legacy_cell_coadd.identifiers.tract,
1051 bbox=tract_bbox,
1052 ),
1053 cell_coadd.bbox,
1054 is_fits=True,
1055 )
1056 if psf := alternates.get("psf"):
1057 compare_psf_to_legacy(psf, legacy_stitched.psf, points=psf_points)
1058 if aperture_corrections := alternates.get("aperture_corrections"):
1059 compare_aperture_corrections_to_legacy(
1060 aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
1061 )
1062 if provenance := alternates.get("provenance"):
1063 compare_cell_coadd_provenance_to_legacy(provenance, legacy_cell_coadd)
1066def compare_cell_coadd_provenance_to_legacy(
1067 provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd
1068) -> None:
1069 """Compare a `.cells.CoaddProvenance` object to a legacy
1070 `lsst.cell_coadds.MultipleCellCoadd` object.
1072 Parameters
1073 ----------
1074 provenance
1075 New provenance object to test.
1076 legacy_cell_coadd
1077 Legacy coadd to test against.
1078 """
1079 from lsst.cell_coadds import ObservationIdentifiers
1081 for legacy_cell in legacy_cell_coadd.cells.values():
1082 cell_index = CellIJ.from_legacy(legacy_cell.identifiers.cell)
1083 prov = provenance[cell_index]
1084 assert prov is not None, f"Legacy cell {cell_index} has inputs, so its provenance is not empty."
1085 legacy_table = astropy.table.Table(
1086 rows=[
1087 [
1088 ids.instrument,
1089 ids.visit,
1090 ids.detector,
1091 ids.day_obs,
1092 ids.physical_filter,
1093 legacy_input.overlaps_center,
1094 legacy_input.overlap_fraction,
1095 legacy_input.weight,
1096 legacy_input.psf_shape.getIxx(),
1097 legacy_input.psf_shape.getIyy(),
1098 legacy_input.psf_shape.getIxy(),
1099 legacy_input.psf_shape_flag,
1100 ]
1101 for ids, legacy_input in legacy_cell.inputs.items()
1102 ],
1103 dtype=[
1104 np.object_,
1105 np.uint64,
1106 np.uint16,
1107 np.uint32,
1108 np.object_,
1109 np.bool_,
1110 np.float64,
1111 np.float64,
1112 np.float64,
1113 np.float64,
1114 np.float64,
1115 np.bool_,
1116 ],
1117 names=[
1118 "instrument",
1119 "visit",
1120 "detector",
1121 "day_obs",
1122 "physical_filter",
1123 "overlaps_center",
1124 "overlap_fraction",
1125 "weight",
1126 "psf_shape_xx",
1127 "psf_shape_yy",
1128 "psf_shape_xy",
1129 "psf_shape_flag",
1130 ],
1131 )
1132 # For a single cell all 'inputs' are also 'contributions'.
1133 assert len(legacy_cell.inputs) == len(prov.inputs)
1134 assert len(legacy_cell.inputs) == len(prov.contributions)
1135 prov.inputs.sort(["instrument", "visit", "detector"])
1136 prov.contributions.sort(["instrument", "visit", "detector"])
1137 legacy_table.sort(["instrument", "visit", "detector"])
1138 assert_values_equal(prov.inputs["instrument"], prov.contributions["instrument"], label="inputs")
1139 assert_values_equal(prov.inputs["visit"], prov.contributions["visit"], label="inputs")
1140 assert_values_equal(prov.inputs["detector"], prov.contributions["detector"], label="inputs")
1141 assert_values_equal(prov.inputs["instrument"], legacy_table["instrument"], label="inputs")
1142 assert_values_equal(prov.inputs["visit"], legacy_table["visit"], label="inputs")
1143 assert_values_equal(prov.inputs["detector"], legacy_table["detector"], label="inputs")
1144 assert_values_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"], label="inputs")
1145 assert_values_equal(prov.inputs["day_obs"], legacy_table["day_obs"], label="inputs")
1146 assert_values_equal(
1147 prov.contributions["overlaps_center"], legacy_table["overlaps_center"], label="inputs"
1148 )
1149 assert_values_equal(
1150 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"], label="inputs"
1151 )
1152 assert_values_equal(prov.contributions["weight"], legacy_table["weight"], label="inputs")
1153 assert_values_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"], label="inputs")
1154 assert_values_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"], label="inputs")
1155 assert_values_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"], label="inputs")
1156 assert_values_equal(
1157 prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"], label="inputs"
1158 )
1159 for row in prov.inputs:
1160 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"})
1161 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key]
1162 assert legacy_polygon == row["polygon"].to_legacy()
1165def compare_psf_to_legacy(
1166 psf: PointSpreadFunction,
1167 legacy_psf: Any,
1168 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
1169 expect_legacy_raise_on_out_of_bounds: bool = False,
1170) -> int:
1171 """Compare a PSF model object to its legacy interface.
1173 Parameters
1174 ----------
1175 psf
1176 Point-spread function to test.
1177 legacy_psf
1178 Legacy `lsst.afw.detection.Psf` instance to compare with.
1179 points
1180 Points to evaluate the PSFs at. If not provided, the intersection of
1181 the PSF bounding boxes are used to generate points on a grid.
1182 expect_legacy_raise_on_out_of_bounds
1183 If `True`, expect ``legacy_psf`` to raise
1184 `lsst.afw.detection.InvalidPsfError` when evaluated at a position
1185 considered out-of-bounds by ``psf``.
1187 Returns
1188 -------
1189 `int`
1190 The number of points actually tested.
1191 """
1192 from lsst.afw.detection import InvalidPsfError
1194 if points is None:
1195 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel)
1196 legacy_points = arrays_to_legacy_points(points.x, points.y)
1197 n_points_tested: int = 0
1198 for p in legacy_points:
1199 if not psf.bounds.contains(x=p.x, y=p.y):
1200 if expect_legacy_raise_on_out_of_bounds:
1201 with pytest.raises(InvalidPsfError):
1202 legacy_psf.computeKernelImage(p)
1203 continue
1204 assert psf.kernel_bbox == Box.from_legacy(legacy_psf.computeKernelBBox(p))
1205 assert psf.compute_kernel_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeKernelImage(p))
1206 assert psf.compute_stellar_bbox(x=p.x, y=p.y) == Box.from_legacy(legacy_psf.computeImageBBox(p))
1207 assert psf.compute_stellar_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeImage(p))
1208 n_points_tested += 1
1209 return n_points_tested
1212def compare_field_to_legacy(
1213 field: BaseField,
1214 legacy_field: Any,
1215 subimage_bbox: Box,
1216) -> None:
1217 """Test a Field object by comparing it to an equivalent
1218 `lsst.afw.math.BoundedField`.
1220 Parameters
1221 ----------
1222 field
1223 Field to test.
1224 legacy_field : ``lsst.afw.math.BoundedField``
1225 Equivalent legacy bounded field.
1226 subimage_bbox
1227 Bounding box for full-image tests.
1228 """
1229 from lsst.afw.math import BoundedField as LegacyBoundedField
1231 assert field.bounds.bbox == Box.from_legacy(legacy_field.getBBox())
1232 # Pixel coordinates to test the numpy array interface with.
1233 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
1234 if not isinstance(field.bounds, Box): 1234 ↛ 1235line 1234 didn't jump to line 1235 because the condition on line 1234 was never true
1235 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y)
1236 pixel_xy = pixel_xy.map(lambda v: v[mask])
1237 try:
1238 assert_values_equal(
1239 field(x=pixel_xy.x, y=pixel_xy.y),
1240 legacy_field.evaluate(pixel_xy.x, pixel_xy.y),
1241 # The field coefficients come from a numerical fit, so this can
1242 # be bigger than just round-off error epsilon.
1243 rtol=1e-5,
1244 )
1245 except AssertionError as err:
1246 err.add_note(f"evaluated at {pixel_xy}")
1247 raise
1248 if not isinstance(legacy_field, LegacyBoundedField): 1248 ↛ 1251line 1248 didn't jump to line 1251 because the condition on line 1248 was never true
1249 # Legacy StitchedApertureCorrection objects are not true BoundedFields
1250 # and don't have addToImage.
1251 return
1252 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
1253 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
1254 assert_images_equal(
1255 field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13
1256 )
1259def compare_aperture_corrections_to_legacy(
1260 aperture_corrections: Mapping[str, BaseField],
1261 legacy_ap_corr_map: Any,
1262 subimage_bbox: Box,
1263) -> None:
1264 """Test an aperture correction `dict` by comparing it to an equivalent
1265 `lsst.afw.image.ApCorrMap`.
1267 Parameters
1268 ----------
1269 aperture_corrections
1270 Dictionary to test.
1271 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
1272 Equivalent legacy aperture correction map.
1273 subimage_bbox
1274 Bounding box for full-image tests.
1275 """
1276 assert aperture_corrections.keys() == set(legacy_ap_corr_map.keys())
1277 for name, field in aperture_corrections.items():
1278 compare_field_to_legacy(field, legacy_ap_corr_map[name], subimage_bbox)
1281def compare_observation_summary_stats_to_legacy(
1282 summary_stats: ObservationSummaryStats,
1283 legacy_summary_stats: Any,
1284) -> None:
1285 """Test an ObservationSummaryStats object by comparing it to an equivalent
1286 `lsst.afw.image.ExposureSummaryStats`.
1288 Parameters
1289 ----------
1290 summary_stats
1291 Struct to test.
1292 legacy_summary_stats : ``lsst.afw.image.ExposureSummaryStats``
1293 Equivalent legacy struct.
1294 """
1295 for field in dataclasses.fields(legacy_summary_stats):
1296 # Always skip version since that field has no meaning in the new
1297 # type.
1298 if field.name == "version":
1299 continue
1300 a = getattr(legacy_summary_stats, field.name)
1301 b = getattr(summary_stats, field.name)
1302 assert_values_equal(a, b, label=field.name)
1305def compare_sky_projection_to_legacy_wcs[F: Frame](
1306 sky_projection: SkyProjection[F],
1307 legacy_wcs: Any,
1308 pixel_frame: F,
1309 subimage_bbox: Box,
1310 is_fits: bool = False,
1311) -> None:
1312 """Test a Projection object by comparing it to an equivalent
1313 `lsst.afw.geom.SkyWcs`.
1315 Parameters
1316 ----------
1317 sky_projection
1318 Projection to test.
1319 legacy_wcs : ``lsst.afw.geom.SkyWcs``
1320 Equivalent legacy WCS.
1321 pixel_frame
1322 Expected pixel frame for the sky_projection.
1323 subimage_bbox
1324 Bounding box of points to generate for tests.
1325 is_fits
1326 Whether this sky_projection is expected to be exactly representable as
1327 a FITS WCS. If `False` it is assumed to have a FITS approximation
1328 attached instead.
1329 """
1330 # Pixel coordinates to test on over the subimage region of interest:
1331 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
1332 # Array indices of those pixel values (subtract off bbox starts):
1333 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
1334 sky_coords = legacy_coords_to_astropy(
1335 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
1336 )
1337 # Test transforming with the Projection itself, which also tests its
1338 # nested Transform and an Astropy High-Level WCS view with no origin
1339 # change.
1340 check_projection(sky_projection, pixel_xy, sky_coords, pixel_frame, pixel_atol=_AST_INV_ATOL)
1341 # Also test the Astropy High-Level WCS view with an origin change to
1342 # array indices.
1343 check_astropy_wcs_interface(
1344 sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
1345 )
1346 if is_fits: 1346 ↛ 1365line 1346 didn't jump to line 1365 because the condition on line 1346 was always true
1347 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
1348 assert fits_wcs is not None
1349 check_astropy_wcs_interface(fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
1350 # Use that FITS approximation to check that we can make a
1351 # Projection from a FITS WCS, too.
1352 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
1353 check_projection(
1354 fits_projection,
1355 subimage_array_xy,
1356 sky_coords,
1357 pixel_frame,
1358 pixel_atol=1e-5,
1359 )
1360 # We want Projections we create from a FITS WCS to be backed by an
1361 # AST FrameSet so we can convert them into legacy
1362 # `lsst.afw.geom.SkyWcs` objects if desired.
1363 assert "Begin FrameSet" in fits_projection.show()
1364 else:
1365 assert sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False) is None
1366 # The legacy SkyWcs should instead have a FITS approximation
1367 # attached; run the same tests on that.
1368 assert sky_projection.fits_approximation is not None
1369 compare_sky_projection_to_legacy_wcs(
1370 sky_projection.fits_approximation,
1371 legacy_wcs.getFitsApproximation(),
1372 pixel_frame,
1373 subimage_bbox,
1374 is_fits=True,
1375 )
1378def check_transform[I: Frame, O: Frame](
1379 transform: Transform[I, O],
1380 input_xy: XY[np.ndarray],
1381 output_xy: XY[np.ndarray],
1382 in_frame: Frame,
1383 out_frame: Frame,
1384 *,
1385 check_inverted: bool = True,
1386 in_atol: u.Quantity | float | None = None,
1387 out_atol: u.Quantity | float | None = None,
1388 angle_wrapped: bool = False,
1389) -> None:
1390 """Test Transform against known arrays of input and output points.
1392 Parameters
1393 ----------
1394 transform
1395 Transform to test.
1396 input_xy
1397 Arrays of input points.
1398 output_xy
1399 Arrays of output points.
1400 in_frame
1401 Expected input frame.
1402 out_frame
1403 Expected output frame.
1404 check_inverted
1405 If `True`, recurse (once) to test the inverse transform.
1406 in_atol
1407 Expected absolute precision of input points.
1408 out_atol
1409 Expected absolute precision of output points.
1410 angle_wrapped
1411 If `True`, treat comparisons on sky-coordinate axes (as determined
1412 by the frames) as angles modulo 2 pi, so differences that cross
1413 zero or 360 degrees compare correctly.
1414 """
1415 assert transform.in_frame == in_frame
1416 assert transform.out_frame == out_frame
1417 in_atol = 0.0 if in_atol is None else in_atol
1418 out_atol = 0.0 if out_atol is None else out_atol
1419 in_atol_v = in_atol.to_value(in_frame.unit) if isinstance(in_atol, u.Quantity) else float(in_atol)
1420 out_atol_v = out_atol.to_value(out_frame.unit) if isinstance(out_atol, u.Quantity) else float(out_atol)
1421 out_wrapped = angle_wrapped and out_frame == SkyFrame.ICRS
1422 in_wrapped = angle_wrapped and in_frame == SkyFrame.ICRS
1424 def assert_close(
1425 test: np.ndarray | float | u.Quantity,
1426 expected: np.ndarray | float | u.Quantity,
1427 atol: float | u.Quantity,
1428 *,
1429 wrapped: bool,
1430 ) -> None:
1431 """Compare one axis, wrap-aware if that axis is a sky angle.
1433 Parameters
1434 ----------
1435 test
1436 Values to test.
1437 expected
1438 Expected values.
1439 atol
1440 Absolute tolerance, in the units the values are expressed in.
1441 wrapped
1442 If `True`, treat the values as angles modulo 2 pi (radians).
1443 """
1444 if not wrapped:
1445 assert_values_equal(test, expected, atol=atol, rtol=_AST_RTOL)
1446 return
1447 test_values = np.asarray(getattr(test, "value", test))
1448 expected_values = np.asarray(getattr(expected, "value", expected))
1449 atol_rad = atol.to_value(u.rad) if isinstance(atol, u.Quantity) else float(atol)
1450 assert_values_equal(
1451 test_values * u.rad,
1452 expected_values * u.rad,
1453 atol=atol_rad * u.rad,
1454 equal_nan=False,
1455 wrap_angles=True,
1456 )
1458 # Test array interfaces.
1459 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
1460 assert_close(test_output_xy.x, output_xy.x, out_atol_v, wrapped=out_wrapped)
1461 assert_close(test_output_xy.y, output_xy.y, out_atol_v, wrapped=out_wrapped)
1462 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
1463 assert_close(test_input_xy.x, input_xy.x, in_atol_v, wrapped=in_wrapped)
1464 assert_close(test_input_xy.y, input_xy.y, in_atol_v, wrapped=in_wrapped)
1465 # Test scalar interfaces with numpy scalars.
1466 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
1467 forward_xy = transform.apply_forward(x=input_x, y=input_y)
1468 assert_close(forward_xy.x, output_x, out_atol_v, wrapped=out_wrapped)
1469 assert_close(forward_xy.y, output_y, out_atol_v, wrapped=out_wrapped)
1470 inverse_xy = transform.apply_inverse(x=output_x, y=output_y)
1471 assert_close(inverse_xy.x, input_x, in_atol_v, wrapped=in_wrapped)
1472 assert_close(inverse_xy.y, input_y, in_atol_v, wrapped=in_wrapped)
1473 # Test quantity array interfaces.
1474 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
1475 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
1476 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
1477 assert_close(test_output_q_xy.x, output_q_xy.x, out_atol, wrapped=out_wrapped)
1478 assert_close(test_output_q_xy.y, output_q_xy.y, out_atol, wrapped=out_wrapped)
1479 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
1480 assert_close(test_input_q_xy.x, input_q_xy.x, in_atol, wrapped=in_wrapped)
1481 assert_close(test_input_q_xy.y, input_q_xy.y, in_atol, wrapped=in_wrapped)
1482 # Test quantity scalar interfaces.
1483 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
1484 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
1485 ):
1486 forward_q_xy = transform.apply_forward_q(x=input_q_x, y=input_q_y)
1487 assert_close(forward_q_xy.x, output_q_x, out_atol, wrapped=out_wrapped)
1488 assert_close(forward_q_xy.y, output_q_y, out_atol, wrapped=out_wrapped)
1489 inverse_q_xy = transform.apply_inverse_q(x=output_q_x, y=output_q_y)
1490 assert_close(inverse_q_xy.x, input_q_x, in_atol, wrapped=in_wrapped)
1491 assert_close(inverse_q_xy.y, input_q_y, in_atol, wrapped=in_wrapped)
1492 if check_inverted:
1493 # Test the inverse transform.
1494 check_transform(
1495 transform.inverted(),
1496 output_xy,
1497 input_xy,
1498 out_frame,
1499 in_frame,
1500 check_inverted=False,
1501 out_atol=in_atol,
1502 in_atol=out_atol,
1503 angle_wrapped=angle_wrapped,
1504 )
1507def check_projection[P: Frame](
1508 sky_projection: SkyProjection[P],
1509 pixel_xy: XY[np.ndarray],
1510 sky_coords: SkyCoord,
1511 pixel_frame: Frame,
1512 *,
1513 pixel_atol: float | None = None,
1514 sky_atol: u.Quantity | None = None,
1515) -> None:
1516 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1517 coordinates.
1519 Parameters
1520 ----------
1521 sky_projection
1522 Projection to test.
1523 pixel_xy
1524 Arrays of pixel coordinates.
1525 sky_coords
1526 Corresponding sky coordinates.
1527 pixel_frame
1528 Expected pixel frame.
1529 pixel_atol
1530 Expected absolute precision of pixel points.
1531 sky_atol
1532 Expected absolute precision of sky coordinates, as a maximum
1533 great-circle separation.
1534 """
1535 assert sky_projection.pixel_frame == pixel_frame
1536 assert sky_projection.sky_frame == SkyFrame.ICRS
1537 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol)
1538 pixel_atol = 0.0 if pixel_atol is None else pixel_atol
1539 pixel_atol_q = pixel_atol * u.pix
1540 # Test array interfaces.
1541 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1542 assert_values_equal(test_pixel_xy.x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL)
1543 assert_values_equal(test_pixel_xy.y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL)
1544 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1545 assert_values_equal(test_sky_astropy, sky_coords, atol=sky_atol, label="pixel_to_sky")
1546 # Test scalar interfaces.
1547 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1548 assert_values_equal(
1549 sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol, rtol=_AST_RTOL
1550 )
1551 assert_values_equal(
1552 sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol, rtol=_AST_RTOL
1553 )
1554 test_sky_single = sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y)
1555 assert_values_equal(test_sky_single, sky_single, atol=sky_atol, label="pixel_to_sky")
1556 # Test the underlying Transform object.
1557 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1558 check_transform(
1559 sky_projection.pixel_to_sky_transform,
1560 pixel_xy,
1561 sky_xy,
1562 pixel_frame,
1563 SkyFrame.ICRS,
1564 check_inverted=False,
1565 in_atol=pixel_atol_q,
1566 out_atol=sky_atol,
1567 angle_wrapped=True,
1568 )
1569 check_transform(
1570 sky_projection.sky_to_pixel_transform,
1571 sky_xy,
1572 pixel_xy,
1573 SkyFrame.ICRS,
1574 pixel_frame,
1575 check_inverted=False,
1576 in_atol=sky_atol,
1577 out_atol=pixel_atol_q,
1578 angle_wrapped=True,
1579 )
1580 # Test the Astropy interface adapter.
1581 check_astropy_wcs_interface(
1582 sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1583 )
1586def assert_sky_projections_equal(
1587 a: SkyProjection[Any] | None,
1588 b: SkyProjection[Any] | None,
1589 expect_identity: bool | None = None,
1590) -> None:
1591 """Test that two `.SkyProjection` instances are equivalent.
1593 Parameters
1594 ----------
1595 a
1596 First sky projection to compare.
1597 b
1598 Second sky projection to compare.
1599 expect_identity
1600 If not `None`, assert whether ``a`` and ``b`` are the same object.
1601 """
1602 if a is None and b is None:
1603 return
1604 assert a is not None
1605 assert b is not None
1606 match expect_identity:
1607 case True: 1607 ↛ 1608line 1607 didn't jump to line 1608 because the pattern on line 1607 never matched
1608 assert a is b
1609 return
1610 case False:
1611 assert a is not b
1612 case None if a is b:
1613 return
1614 assert a == b
1617def check_astropy_wcs_interface(
1618 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1619 pixel_xy: XY[np.ndarray],
1620 sky_coords: SkyCoord,
1621 *,
1622 pixel_atol: float | None = None,
1623 sky_atol: u.Quantity | None = None,
1624) -> None:
1625 """Test an Astropy WCS instance against known arrays of pixel and
1626 sky coordinates.
1628 Parameters
1629 ----------
1630 wcs
1631 Astropy WCS object to test.
1632 pixel_xy
1633 Arrays of pixel coordinates.
1634 sky_coords
1635 Corresponding sky coordinates.
1636 pixel_atol
1637 Expected absolute precision of pixel points.
1638 sky_atol
1639 Expected absolute precision of sky coordinates, as a maximum
1640 great-circle separation.
1641 """
1642 pixel_atol = 0.0 if pixel_atol is None else pixel_atol
1643 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol)
1644 test_x, test_y = wcs.world_to_pixel(sky_coords)
1645 assert_values_equal(test_x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL)
1646 assert_values_equal(test_y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL)
1647 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
1648 assert_values_equal(test_sky_coords, sky_coords, atol=sky_atol, label="pixel_to_world")
1651def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
1652 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array.
1654 Parameters
1655 ----------
1656 legacy_points
1657 Legacy ``lsst.geom.Point2D`` objects to convert.
1658 """
1659 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
1662def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
1663 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
1664 coordinate object.
1666 Parameters
1667 ----------
1668 legacy_coords
1669 Legacy ``lsst.geom.SpherePoint`` objects to convert.
1670 """
1671 return SkyCoord(
1672 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
1673 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
1674 )
1677def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
1678 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``.
1680 Parameters
1681 ----------
1682 x
1683 X coordinates of the points.
1684 y
1685 Y coordinates of the points.
1686 """
1687 from lsst.geom import Point2D
1689 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
1692def compare_amplifier_to_legacy(
1693 amplifier: Amplifier,
1694 legacy_amplifier: Any,
1695 *,
1696 is_raw_assembled: bool,
1697 expect_nominal_calibrations: bool = True,
1698) -> None:
1699 """Compare an `~.cameras.Amplifier` to a legacy
1700 `lsst.afw.cameraGeom.Amplifier`.
1702 Parameters
1703 ----------
1704 amplifier
1705 Amplifier to compare.
1706 legacy_amplifier
1707 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against.
1708 is_raw_assembled
1709 Whether the raw geometry is expected to be the assembled-raw
1710 geometry (`True`) or the unassembled-raw geometry (`False`).
1711 expect_nominal_calibrations
1712 Whether the amplifier is expected to carry nominal calibrations.
1713 """
1714 assert legacy_amplifier.getName() == amplifier.name
1715 assert Box.from_legacy(legacy_amplifier.getBBox()) == amplifier.bbox
1716 if is_raw_assembled: 1716 ↛ 1719line 1716 didn't jump to line 1719 because the condition on line 1716 was always true
1717 raw_geom = amplifier.assembled_raw_geometry
1718 else:
1719 raw_geom = amplifier.unassembled_raw_geometry
1720 assert raw_geom is not None
1721 assert ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()) == raw_geom.readout_corner
1722 assert Box.from_legacy(legacy_amplifier.getRawBBox()) == raw_geom.bbox
1723 assert Box.from_legacy(legacy_amplifier.getRawDataBBox()) == raw_geom.data_bbox
1724 assert legacy_amplifier.getRawFlipX() == raw_geom.flip_x
1725 assert legacy_amplifier.getRawFlipY() == raw_geom.flip_y
1726 assert legacy_amplifier.getRawXYOffset().getX() == raw_geom.x_offset
1727 assert legacy_amplifier.getRawXYOffset().getY() == raw_geom.y_offset
1728 assert (
1729 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()) == raw_geom.horizontal_overscan_bbox
1730 )
1731 assert Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()) == raw_geom.vertical_overscan_bbox
1732 if raw_geom.horizontal_prescan_bbox is None:
1733 # `Box` requires a positive size, so an amplifier with no prescan
1734 # region carries `None` here.
1735 assert legacy_amplifier.getRawPrescanBBox().isEmpty()
1736 else:
1737 assert Box.from_legacy(legacy_amplifier.getRawPrescanBBox()) == raw_geom.horizontal_prescan_bbox
1738 if expect_nominal_calibrations: 1738 ↛ exitline 1738 didn't return from function 'compare_amplifier_to_legacy' because the condition on line 1738 was always true
1739 assert amplifier.nominal_calibrations is not None
1740 assert_values_equal(
1741 legacy_amplifier.getGain(),
1742 amplifier.nominal_calibrations.gain,
1743 equal_nan=True,
1744 label="gain",
1745 )
1746 assert_values_equal(
1747 legacy_amplifier.getReadNoise(),
1748 amplifier.nominal_calibrations.read_noise,
1749 equal_nan=True,
1750 label="read_noise",
1751 )
1752 assert_values_equal(
1753 legacy_amplifier.getSaturation(),
1754 amplifier.nominal_calibrations.saturation,
1755 equal_nan=True,
1756 label="saturation",
1757 )
1758 assert_values_equal(
1759 legacy_amplifier.getSuspectLevel(),
1760 amplifier.nominal_calibrations.suspect_level,
1761 equal_nan=True,
1762 label="suspect_level",
1763 )
1764 assert_values_equal(
1765 legacy_amplifier.getLinearityCoeffs(),
1766 amplifier.nominal_calibrations.linearity_coefficients,
1767 label="linearity_coefficients",
1768 )
1769 assert legacy_amplifier.getLinearityType() == amplifier.nominal_calibrations.linearity_type
1772def compare_detector_to_legacy(
1773 detector: Detector,
1774 legacy_detector: Any,
1775 *,
1776 is_raw_assembled: bool,
1777 expect_nominal_calibrations: bool = True,
1778) -> None:
1779 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`.
1781 Parameters
1782 ----------
1783 detector
1784 Detector to compare.
1785 legacy_detector
1786 Legacy `lsst.afw.cameraGeom.Detector` to compare against.
1787 is_raw_assembled
1788 Whether the raw geometry is expected to be the assembled-raw
1789 geometry (`True`) or the unassembled-raw geometry (`False`).
1790 expect_nominal_calibrations
1791 Whether the detector's amplifiers are expected to carry nominal
1792 calibrations.
1793 """
1794 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
1796 assert legacy_detector.getName() == detector.name
1797 assert legacy_detector.getId() == detector.id
1798 assert DetectorType.from_legacy(legacy_detector.getType()) == detector.type
1799 assert Box.from_legacy(legacy_detector.getBBox()) == detector.bbox
1800 assert legacy_detector.getSerial() == detector.serial
1801 legacy_orientation = legacy_detector.getOrientation()
1802 assert legacy_orientation.getFpPosition3().getX() == detector.orientation.focal_plane_x
1803 assert legacy_orientation.getFpPosition3().getY() == detector.orientation.focal_plane_y
1804 assert legacy_orientation.getFpPosition3().getZ() == detector.orientation.focal_plane_z
1805 assert legacy_orientation.getReferencePoint().getX() == detector.orientation.pixel_reference_x
1806 assert legacy_orientation.getReferencePoint().getY() == detector.orientation.pixel_reference_y
1807 assert legacy_orientation.getYaw().asRadians() == detector.orientation.yaw.to_value(u.rad)
1808 assert legacy_orientation.getPitch().asRadians() == detector.orientation.pitch.to_value(u.rad)
1809 assert legacy_orientation.getRoll().asRadians() == detector.orientation.roll.to_value(u.rad)
1810 assert legacy_detector.getPixelSize().getX() == detector.pixel_size
1811 assert legacy_detector.getPhysicalType() == detector.physical_type
1812 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
1813 compare_amplifier_to_legacy(
1814 amplifier,
1815 legacy_amplifier,
1816 is_raw_assembled=is_raw_assembled,
1817 expect_nominal_calibrations=expect_nominal_calibrations,
1818 )
1819 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
1820 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
1821 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
1822 check_transform(
1823 detector.to_focal_plane,
1824 pixel_xy,
1825 legacy_points_to_xy_array(fp_legacy_points),
1826 detector.frame,
1827 detector.to_focal_plane.out_frame,
1828 in_atol=1e-9 * u.pix,
1829 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
1830 )
1831 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
1832 check_transform(
1833 detector.to_field_angle,
1834 pixel_xy,
1835 legacy_points_to_xy_array(fa_legacy_points),
1836 detector.frame,
1837 detector.to_field_angle.out_frame,
1838 in_atol=1e-9 * u.pix,
1839 out_atol=1e-7 * u.arcsec,
1840 )
1843def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
1844 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
1846 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
1847 classes. Importing this module already imports every ``lsst.images``
1848 module that defines a subclass, so the tree is fully populated by the time
1849 this is called.
1851 This discovery is deliberately separate from
1852 `check_archive_tree_class_invariants` so that the per-class check stays
1853 usable on a single class even if this metaprogramming is removed later.
1854 """
1855 seen: set[type] = set()
1856 stack: list[type] = [ArchiveTree]
1857 while stack:
1858 kls = stack.pop()
1859 sub: type
1860 for sub in kls.__subclasses__():
1861 if sub in seen: 1861 ↛ 1862line 1861 didn't jump to line 1862 because the condition on line 1861 was never true
1862 continue
1863 seen.add(sub)
1864 stack.append(sub)
1865 if not getattr(sub, "__abstractmethods__", None): 1865 ↛ 1860line 1865 didn't jump to line 1860 because the condition on line 1865 was always true
1866 yield sub
1869def check_bounds_contains_broadcasting(bounds: Bounds) -> None:
1870 """Verify that `~lsst.images.Bounds.contains` accepts array-like inputs.
1872 Uses the scalar overload as the reference and checks that 1-D arrays,
1873 list inputs, mixed scalar-plus-array inputs, 2-D broadcast inputs, and
1874 `.XY` / `.YX` positional-argument forms all produce results consistent
1875 with calling the scalar overload on each ``(x, y)`` pair individually.
1877 Parameters
1878 ----------
1879 bounds
1880 The `~lsst.images.Bounds` implementation to exercise.
1881 """
1882 bbox = bounds.bbox
1883 # One point outside each boundary, one on each boundary, one in the
1884 # interior. The boundary points (start and stop) straddle inside/outside
1885 # (start is inside, stop is outside) so we are guaranteed a mix of True
1886 # and False results without hard-coding expected values.
1887 ys = np.array(
1888 [
1889 bbox.y.start - 1,
1890 bbox.y.start,
1891 (bbox.y.start + bbox.y.stop) // 2,
1892 bbox.y.stop,
1893 ]
1894 )
1895 xs = np.array(
1896 [
1897 bbox.x.start - 1,
1898 bbox.x.start,
1899 (bbox.x.start + bbox.x.stop) // 2,
1900 bbox.x.stop,
1901 ]
1902 )
1903 # 2-D reference: expected[i, j] == contains(x=xs[j], y=ys[i]).
1904 expected = np.array(
1905 [[bounds.contains(x=int(xi), y=int(yi)) for xi in xs] for yi in ys],
1906 dtype=bool,
1907 )
1908 # 1-D ndarray (diagonal pairs).
1909 assert_values_equal(
1910 bounds.contains(x=xs, y=ys),
1911 np.diagonal(expected),
1912 label="bounds.contains",
1913 )
1914 # list inputs (array-like).
1915 assert_values_equal(
1916 bounds.contains(x=xs.tolist(), y=ys.tolist()),
1917 np.diagonal(expected),
1918 label="bounds.contains",
1919 )
1920 # Mixed: scalar y, 1-D array x — produces a 1-D result.
1921 fixed_yi = 1 # index into ys; ys[1] == bbox.y.start (on the boundary)
1922 assert_values_equal(
1923 bounds.contains(x=xs, y=int(ys[fixed_yi])),
1924 expected[fixed_yi, :],
1925 label="bounds.contains",
1926 )
1927 # Mixed: 1-D array y, scalar x.
1928 fixed_xi = 1 # index into xs; xs[1] == bbox.x.start (on the boundary)
1929 assert_values_equal(
1930 bounds.contains(x=int(xs[fixed_xi]), y=ys),
1931 expected[:, fixed_xi],
1932 label="bounds.contains",
1933 )
1934 # Float scalars: results must match the int-scalar reference.
1935 assert bounds.contains(x=float(xs[fixed_xi]), y=float(ys[fixed_yi])) == expected[fixed_yi, fixed_xi]
1936 # XY / YX scalar: results must match the keyword-scalar reference.
1937 assert bounds.contains(XY(x=int(xs[fixed_xi]), y=int(ys[fixed_yi]))) == expected[fixed_yi, fixed_xi]
1938 assert bounds.contains(YX(y=int(ys[fixed_yi]), x=int(xs[fixed_xi]))) == expected[fixed_yi, fixed_xi]
1939 # XY / YX array: results must match the 1-D ndarray reference.
1940 assert_values_equal(bounds.contains(XY(x=xs, y=ys)), np.diagonal(expected), label="bounds.contains")
1941 assert_values_equal(bounds.contains(YX(y=ys, x=xs)), np.diagonal(expected), label="bounds.contains")
1942 # 2-D broadcast: y shape (N, 1) × x shape (1, M) → (N, M).
1943 assert_values_equal(
1944 bounds.contains(x=xs.reshape(1, -1), y=ys.reshape(-1, 1)),
1945 expected,
1946 label="bounds.contains",
1947 )
1950def check_archive_tree_class_invariants(tree_type: type[ArchiveTree]) -> None:
1951 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
1952 well-formed schema-version constants and an in-memory type.
1954 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
1955 ``PUBLIC_TYPE`` are present and well-typed, that the version is
1956 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
1957 schema major.
1959 Parameters
1960 ----------
1961 tree_type
1962 The concrete `.serialization.ArchiveTree` subclass to check.
1963 """
1964 assert hasattr(tree_type, "SCHEMA_NAME"), f"{tree_type.__name__} lacks SCHEMA_NAME"
1965 assert hasattr(tree_type, "SCHEMA_VERSION"), f"{tree_type.__name__} lacks SCHEMA_VERSION"
1966 assert hasattr(tree_type, "MIN_READ_VERSION"), f"{tree_type.__name__} lacks MIN_READ_VERSION"
1967 assert hasattr(tree_type, "PUBLIC_TYPE"), f"{tree_type.__name__} lacks PUBLIC_TYPE"
1968 assert isinstance(tree_type.SCHEMA_NAME, str)
1969 assert len(tree_type.SCHEMA_NAME) > 0
1970 # Allow an optional PEP 440 development-release suffix (e.g. 1.0.0.dev0)
1971 # for schemas still in development.
1972 assert re.fullmatch(r"^\d+\.\d+\.\d+(\.dev\d+)?$", tree_type.SCHEMA_VERSION)
1973 assert isinstance(tree_type.MIN_READ_VERSION, int)
1974 assert tree_type.MIN_READ_VERSION >= 1
1975 assert isinstance(tree_type.PUBLIC_TYPE, type)
1976 major = int(tree_type.SCHEMA_VERSION.split(".")[0])
1977 assert tree_type.MIN_READ_VERSION <= major
1980def _angle_to_astropy(angle: Any) -> u.Quantity:
1981 """Convert an `lsst.geom.Angle` to an Astropy angular quantity."""
1982 return u.Quantity(angle.asRadians(), u.rad)
1985def _observatory_to_location(observatory: Any) -> coords.EarthLocation:
1986 """Convert an `lsst.afw.coord.Observatory` to an Astropy
1987 `~astropy.coordinates.EarthLocation`.
1988 """
1989 return coords.EarthLocation.from_geodetic(
1990 lon=observatory.getLongitude().asRadians() * u.rad,
1991 lat=observatory.getLatitude().asRadians() * u.rad,
1992 height=observatory.getElevation() * u.m,
1993 )
1996def assert_obs_metadata_fields_equal(actual: Any, expected: Any, *, label: str = "") -> None:
1997 """Assert that two ObservationInfo or VisitInfo field values are equal.
1999 Dispatches on the type of ``expected`` so that values affected by afw's
2000 float64 date storage, unit conversions, and angle wrapping are compared
2001 with appropriate tolerances; anything else is compared with ``==``.
2003 Parameters
2004 ----------
2005 actual : `object`
2006 The field value to check.
2007 expected : `object`
2008 The field value to check it against.
2009 label : `str`, optional
2010 Prefix for the failure message.
2011 """
2012 from lsst.afw.coord import Observatory, Weather
2013 from lsst.daf.base import DateTime
2014 from lsst.geom import Angle
2016 if isinstance(expected, DateTime):
2017 assert_obs_metadata_fields_equal(actual.toAstropy(), expected.toAstropy(), label=label)
2018 elif isinstance(expected, astropy.time.Time):
2019 # afw stores VisitInfo dates as float64 MJDs, which limits precision.
2020 assert_values_equal(actual, expected, atol=1e-9 * u.s, label=label)
2021 elif isinstance(expected, coords.AltAz):
2022 assert_values_equal(actual.az, expected.az, wrap_angles=True, label=f"{label}.az")
2023 assert_values_equal(actual.alt, expected.alt, atol=1e-3 * u.arcsec, label=f"{label}.alt")
2024 elif isinstance(expected, coords.SkyCoord):
2025 assert_values_equal(actual, expected, atol=1e-3 * u.arcsec, label=label)
2026 elif isinstance(expected, coords.EarthLocation):
2027 assert_values_equal(
2028 u.Quantity([actual.x, actual.y, actual.z]),
2029 u.Quantity([expected.x, expected.y, expected.z]),
2030 rtol=1e-9,
2031 label=label,
2032 )
2033 elif isinstance(expected, Observatory):
2034 assert_obs_metadata_fields_equal(
2035 _observatory_to_location(actual), _observatory_to_location(expected), label=label
2036 )
2037 elif isinstance(expected, Weather):
2038 assert_obs_metadata_fields_equal(
2039 u.Quantity(actual.getAirTemperature(), u.deg_C),
2040 u.Quantity(expected.getAirTemperature(), u.deg_C),
2041 label=f"{label}.temperature",
2042 )
2043 assert_obs_metadata_fields_equal(
2044 actual.getAirPressure(), expected.getAirPressure(), label=f"{label}.pressure"
2045 )
2046 assert_obs_metadata_fields_equal(
2047 actual.getHumidity(), expected.getHumidity(), label=f"{label}.humidity"
2048 )
2049 elif isinstance(expected, Angle):
2050 assert_values_equal(
2051 _angle_to_astropy(actual),
2052 _angle_to_astropy(expected),
2053 atol=1e-12 * u.rad,
2054 wrap_angles=True,
2055 label=label,
2056 )
2057 elif (
2058 isinstance(expected, u.Quantity)
2059 and isinstance(actual, u.Quantity)
2060 and expected.unit.physical_type == "temperature"
2061 and actual.unit.physical_type == "temperature"
2062 ):
2063 assert_values_equal(
2064 actual.to(u.K, equivalencies=u.temperature()),
2065 expected.to(u.K, equivalencies=u.temperature()),
2066 atol=1e-6 * u.K,
2067 label=label,
2068 )
2069 elif isinstance(expected, u.Quantity | float):
2070 assert_values_equal(actual, expected, rtol=1e-9, label=label)
2071 elif actual != expected:
2072 raise AssertionError(f"{label}: {actual!r} != {expected!r}")