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