Coverage for python/lsst/images/tests/_checks.py: 49%
455 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-16 15:24 -0700
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 a = getattr(legacy_summary_stats, field.name)
855 b = getattr(summary_stats, field.name)
856 if isinstance(b, tuple):
857 for ai, bi in zip(a, b):
858 assert ai == bi or (math.isnan(ai) and math.isnan(bi)), f"{field.name}: {a} != {b}"
859 else:
860 assert a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}"
863def compare_sky_projection_to_legacy_wcs[F: Frame](
864 sky_projection: SkyProjection[F],
865 legacy_wcs: Any,
866 pixel_frame: F,
867 subimage_bbox: Box,
868 is_fits: bool = False,
869) -> None:
870 """Test a Projection object by comparing it to an equivalent
871 `lsst.afw.geom.SkyWcs`.
873 Parameters
874 ----------
875 sky_projection
876 Projection to test.
877 legacy_wcs : ``lsst.afw.geom.SkyWcs``
878 Equivalent legacy WCS.
879 pixel_frame
880 Expected pixel frame for the sky_projection.
881 subimage_bbox
882 Bounding box of points to generate for tests.
883 is_fits
884 Whether this sky_projection is expected to be exactly representable as
885 a FITS WCS. If `False` it is assumed to have a FITS approximation
886 attached instead.
887 """
888 # Pixel coordinates to test on over the subimage region of interest:
889 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
890 # Array indices of those pixel values (subtract off bbox starts):
891 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
892 sky_coords = legacy_coords_to_astropy(
893 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
894 )
895 # Test transforming with the Projection itself, which also tests its
896 # nested Transform and an Astropy High-Level WCS view with no origin
897 # change.
898 check_projection(sky_projection, pixel_xy, sky_coords, pixel_frame)
899 # Also test the Astropy High-Level WCS view with an origin change to
900 # array indices.
901 check_astropy_wcs_interface(
902 sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
903 )
904 if is_fits: 904 ↛ 923line 904 didn't jump to line 923 because the condition on line 904 was always true
905 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
906 assert fits_wcs is not None
907 check_astropy_wcs_interface(fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
908 # Use that FITS approximation to check that we can make a
909 # Projection from a FITS WCS, too.
910 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
911 check_projection(
912 fits_projection,
913 subimage_array_xy,
914 sky_coords,
915 pixel_frame,
916 pixel_atol=1e-5,
917 )
918 # We want Projections we create from a FITS WCS to be backed by an
919 # AST FrameSet so we can convert them into legacy
920 # `lsst.afw.geom.SkyWcs` objects if desired.
921 assert "Begin FrameSet" in fits_projection.show()
922 else:
923 assert sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False) is None
924 # The legacy SkyWcs should instead have a FITS approximation
925 # attached; run the same tests on that.
926 assert sky_projection.fits_approximation is not None
927 compare_sky_projection_to_legacy_wcs(
928 sky_projection.fits_approximation,
929 legacy_wcs.getFitsApproximation(),
930 pixel_frame,
931 subimage_bbox,
932 is_fits=True,
933 )
936def check_transform[I: Frame, O: Frame](
937 transform: Transform[I, O],
938 input_xy: XY[np.ndarray],
939 output_xy: XY[np.ndarray],
940 in_frame: Frame,
941 out_frame: Frame,
942 *,
943 check_inverted: bool = True,
944 in_atol: u.Quantity | None = None,
945 out_atol: u.Quantity | None = None,
946) -> None:
947 """Test Transform against known arrays of input and output points.
949 Parameters
950 ----------
951 transform
952 Transform to test.
953 input_xy
954 Arrays of input points.
955 output_xy
956 Arrays of output points.
957 in_frame
958 Expected input frame.
959 out_frame
960 Expected output frame.
961 check_inverted
962 If `True`, recurse (once) to test the inverse transform.
963 in_atol
964 Expected absolute precision of input points.
965 out_atol
966 Expected absolute precision of output points.
967 """
968 assert transform.in_frame == in_frame
969 assert transform.out_frame == out_frame
970 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None
971 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None
972 # Test array interfaces.
973 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
974 assert_close(test_output_xy.x, output_xy.x, atol=out_atol_v)
975 assert_close(test_output_xy.y, output_xy.y, atol=out_atol_v)
976 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
977 assert_close(test_input_xy.x, input_xy.x, atol=in_atol_v)
978 assert_close(test_input_xy.y, input_xy.y, atol=in_atol_v)
979 # Test scalar interfaces with numpy scalars.
980 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
981 assert_close(transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v)
982 assert_close(transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v)
983 assert_close(transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v)
984 assert_close(transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v)
985 # Test quantity array interfaces.
986 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
987 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
988 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
989 assert_close(test_output_q_xy.x, output_q_xy.x, atol=out_atol)
990 assert_close(test_output_q_xy.y, output_q_xy.y, atol=out_atol)
991 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
992 assert_close(test_input_q_xy.x, input_q_xy.x, atol=in_atol)
993 assert_close(test_input_q_xy.y, input_q_xy.y, atol=in_atol)
994 # Test quantity scalar interfaces.
995 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
996 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
997 ):
998 assert_close(transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol)
999 assert_close(transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol)
1000 assert_close(transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol)
1001 assert_close(transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol)
1002 if check_inverted:
1003 # Test the inverse transform.
1004 check_transform(
1005 transform.inverted(),
1006 output_xy,
1007 input_xy,
1008 out_frame,
1009 in_frame,
1010 check_inverted=False,
1011 out_atol=in_atol,
1012 in_atol=out_atol,
1013 )
1016def check_projection[P: Frame](
1017 sky_projection: SkyProjection[P],
1018 pixel_xy: XY[np.ndarray],
1019 sky_coords: SkyCoord,
1020 pixel_frame: Frame,
1021 *,
1022 pixel_atol: float | None = None,
1023 sky_atol: u.Quantity | None = None,
1024) -> None:
1025 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1026 coordinates.
1028 Parameters
1029 ----------
1030 sky_projection
1031 Projection to test.
1032 pixel_xy
1033 Arrays of pixel coordinates.
1034 sky_coords
1035 Corresponding sky coordinates.
1036 pixel_frame
1037 Expected pixel frame.
1038 pixel_atol
1039 Expected absolute precision of pixel points.
1040 sky_atol
1041 Expected absolute precision of sky coordinates.
1042 """
1043 assert sky_projection.pixel_frame == pixel_frame
1044 assert sky_projection.sky_frame == SkyFrame.ICRS
1045 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None
1046 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None
1047 # Test array interfaces.
1048 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1049 assert_close(test_pixel_xy.x, pixel_xy.x, atol=pixel_atol)
1050 assert_close(test_pixel_xy.y, pixel_xy.y, atol=pixel_atol)
1051 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1052 assert_close(test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v)
1053 assert_close(test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v)
1054 # Test scalar interfaces.
1055 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1056 assert_close(sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol)
1057 assert_close(sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol)
1058 assert_close(sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v)
1059 assert_close(sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v)
1060 # Test the underlying Transform object.
1061 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1062 check_transform(
1063 sky_projection.pixel_to_sky_transform,
1064 pixel_xy,
1065 sky_xy,
1066 pixel_frame,
1067 SkyFrame.ICRS,
1068 check_inverted=False,
1069 in_atol=pixel_atol_q,
1070 out_atol=sky_atol,
1071 )
1072 check_transform(
1073 sky_projection.sky_to_pixel_transform,
1074 sky_xy,
1075 pixel_xy,
1076 SkyFrame.ICRS,
1077 pixel_frame,
1078 check_inverted=False,
1079 in_atol=sky_atol,
1080 out_atol=pixel_atol_q,
1081 )
1082 # Test the Astropy interface adapter.
1083 check_astropy_wcs_interface(
1084 sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1085 )
1088def assert_sky_projections_equal(
1089 a: SkyProjection[Any] | None,
1090 b: SkyProjection[Any] | None,
1091 expect_identity: bool | None = None,
1092) -> None:
1093 """Test that two `.SkyProjection` instances are equivalent.
1095 Parameters
1096 ----------
1097 a
1098 First sky projection to compare.
1099 b
1100 Second sky projection to compare.
1101 expect_identity
1102 If not `None`, assert whether ``a`` and ``b`` are the same object.
1103 """
1104 if a is None and b is None:
1105 return
1106 assert a is not None and b is not None
1107 match expect_identity:
1108 case True: 1108 ↛ 1109line 1108 didn't jump to line 1109 because the pattern on line 1108 never matched
1109 assert a is b
1110 return
1111 case False:
1112 assert a is not b
1113 case None if a is b:
1114 return
1115 assert a == b
1118def check_astropy_wcs_interface(
1119 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1120 pixel_xy: XY[np.ndarray],
1121 sky_coords: SkyCoord,
1122 *,
1123 pixel_atol: float | None = None,
1124 sky_atol: u.Quantity | None = None,
1125) -> None:
1126 """Test an Astropy WCS instance against known arrays of pixel and
1127 sky coordinates.
1129 Parameters
1130 ----------
1131 wcs
1132 Astropy WCS object to test.
1133 pixel_xy
1134 Arrays of pixel coordinates.
1135 sky_coords
1136 Corresponding sky coordinates.
1137 pixel_atol
1138 Expected absolute precision of pixel points.
1139 sky_atol
1140 Expected absolute precision of sky coordinates.
1141 """
1142 test_x, test_y = wcs.world_to_pixel(sky_coords)
1143 assert_close(test_x, pixel_xy.x, atol=pixel_atol)
1144 assert_close(test_y, pixel_xy.y, atol=pixel_atol)
1145 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
1146 assert_close(test_sky_coords.ra, sky_coords.ra, atol=sky_atol)
1147 assert_close(test_sky_coords.dec, sky_coords.dec, atol=sky_atol)
1150def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
1151 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array.
1153 Parameters
1154 ----------
1155 legacy_points
1156 Legacy ``lsst.geom.Point2D`` objects to convert.
1157 """
1158 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
1161def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
1162 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
1163 coordinate object.
1165 Parameters
1166 ----------
1167 legacy_coords
1168 Legacy ``lsst.geom.SpherePoint`` objects to convert.
1169 """
1170 return SkyCoord(
1171 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
1172 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
1173 )
1176def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
1177 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``.
1179 Parameters
1180 ----------
1181 x
1182 X coordinates of the points.
1183 y
1184 Y coordinates of the points.
1185 """
1186 from lsst.geom import Point2D
1188 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
1191def compare_amplifier_to_legacy(
1192 amplifier: Amplifier,
1193 legacy_amplifier: Any,
1194 *,
1195 is_raw_assembled: bool,
1196 expect_nominal_calibrations: bool = True,
1197) -> None:
1198 """Compare an `~.cameras.Amplifier` to a legacy
1199 `lsst.afw.cameraGeom.Amplifier`.
1201 Parameters
1202 ----------
1203 amplifier
1204 Amplifier to compare.
1205 legacy_amplifier
1206 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against.
1207 is_raw_assembled
1208 Whether the raw geometry is expected to be the assembled-raw
1209 geometry (`True`) or the unassembled-raw geometry (`False`).
1210 expect_nominal_calibrations
1211 Whether the amplifier is expected to carry nominal calibrations.
1212 """
1213 assert legacy_amplifier.getName() == amplifier.name
1214 assert Box.from_legacy(legacy_amplifier.getBBox()) == amplifier.bbox
1215 if is_raw_assembled:
1216 raw_geom = amplifier.assembled_raw_geometry
1217 else:
1218 raw_geom = amplifier.unassembled_raw_geometry
1219 assert raw_geom is not None
1220 assert ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()) == raw_geom.readout_corner
1221 assert Box.from_legacy(legacy_amplifier.getRawBBox()) == raw_geom.bbox
1222 assert Box.from_legacy(legacy_amplifier.getRawDataBBox()) == raw_geom.data_bbox
1223 assert legacy_amplifier.getRawFlipX() == raw_geom.flip_x
1224 assert legacy_amplifier.getRawFlipY() == raw_geom.flip_y
1225 assert legacy_amplifier.getRawXYOffset().getX() == raw_geom.x_offset
1226 assert legacy_amplifier.getRawXYOffset().getY() == raw_geom.y_offset
1227 assert (
1228 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()) == raw_geom.horizontal_overscan_bbox
1229 )
1230 assert Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()) == raw_geom.vertical_overscan_bbox
1231 assert Box.from_legacy(legacy_amplifier.getRawPrescanBBox()) == raw_geom.horizontal_prescan_bbox
1232 if expect_nominal_calibrations:
1233 assert amplifier.nominal_calibrations is not None
1234 assert_equal_allow_nan(legacy_amplifier.getGain(), amplifier.nominal_calibrations.gain)
1235 assert_equal_allow_nan(legacy_amplifier.getReadNoise(), amplifier.nominal_calibrations.read_noise)
1236 assert_equal_allow_nan(legacy_amplifier.getSaturation(), amplifier.nominal_calibrations.saturation)
1237 assert_equal_allow_nan(
1238 legacy_amplifier.getSuspectLevel(), amplifier.nominal_calibrations.suspect_level
1239 )
1240 np.testing.assert_array_equal(
1241 legacy_amplifier.getLinearityCoeffs(), amplifier.nominal_calibrations.linearity_coefficients
1242 )
1243 assert legacy_amplifier.getLinearityType() == amplifier.nominal_calibrations.linearity_type
1246def compare_detector_to_legacy(
1247 detector: Detector,
1248 legacy_detector: Any,
1249 *,
1250 is_raw_assembled: bool,
1251 expect_nominal_calibrations: bool = True,
1252) -> None:
1253 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`.
1255 Parameters
1256 ----------
1257 detector
1258 Detector to compare.
1259 legacy_detector
1260 Legacy `lsst.afw.cameraGeom.Detector` to compare against.
1261 is_raw_assembled
1262 Whether the raw geometry is expected to be the assembled-raw
1263 geometry (`True`) or the unassembled-raw geometry (`False`).
1264 expect_nominal_calibrations
1265 Whether the detector's amplifiers are expected to carry nominal
1266 calibrations.
1267 """
1268 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
1270 assert legacy_detector.getName() == detector.name
1271 assert legacy_detector.getId() == detector.id
1272 assert DetectorType.from_legacy(legacy_detector.getType()) == detector.type
1273 assert Box.from_legacy(legacy_detector.getBBox()) == detector.bbox
1274 assert legacy_detector.getSerial() == detector.serial
1275 legacy_orientation = legacy_detector.getOrientation()
1276 assert legacy_orientation.getFpPosition3().getX() == detector.orientation.focal_plane_x
1277 assert legacy_orientation.getFpPosition3().getY() == detector.orientation.focal_plane_y
1278 assert legacy_orientation.getFpPosition3().getZ() == detector.orientation.focal_plane_z
1279 assert legacy_orientation.getReferencePoint().getX() == detector.orientation.pixel_reference_x
1280 assert legacy_orientation.getReferencePoint().getY() == detector.orientation.pixel_reference_y
1281 assert legacy_orientation.getYaw().asRadians() == detector.orientation.yaw.to_value(u.rad)
1282 assert legacy_orientation.getPitch().asRadians() == detector.orientation.pitch.to_value(u.rad)
1283 assert legacy_orientation.getRoll().asRadians() == detector.orientation.roll.to_value(u.rad)
1284 assert legacy_detector.getPixelSize().getX() == detector.pixel_size
1285 assert legacy_detector.getPhysicalType() == detector.physical_type
1286 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
1287 compare_amplifier_to_legacy(
1288 amplifier,
1289 legacy_amplifier,
1290 is_raw_assembled=is_raw_assembled,
1291 expect_nominal_calibrations=expect_nominal_calibrations,
1292 )
1293 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
1294 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
1295 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
1296 check_transform(
1297 detector.to_focal_plane,
1298 pixel_xy,
1299 legacy_points_to_xy_array(fp_legacy_points),
1300 detector.frame,
1301 detector.to_focal_plane.out_frame,
1302 in_atol=1e-9 * u.pix,
1303 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
1304 )
1305 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
1306 check_transform(
1307 detector.to_field_angle,
1308 pixel_xy,
1309 legacy_points_to_xy_array(fa_legacy_points),
1310 detector.frame,
1311 detector.to_field_angle.out_frame,
1312 in_atol=1e-9 * u.pix,
1313 out_atol=1e-7 * u.arcsec,
1314 )
1317def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
1318 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
1320 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
1321 classes. Importing this module already imports every ``lsst.images``
1322 module that defines a subclass, so the tree is fully populated by the time
1323 this is called.
1325 This discovery is deliberately separate from
1326 `check_archive_tree_class_invariants` so that the per-class check stays
1327 usable on a single class even if this metaprogramming is removed later.
1328 """
1329 seen: set[type] = set()
1330 stack: list[type] = [ArchiveTree]
1331 while stack:
1332 kls = stack.pop()
1333 for sub in kls.__subclasses__():
1334 if sub in seen: 1334 ↛ 1335line 1334 didn't jump to line 1335 because the condition on line 1334 was never true
1335 continue
1336 seen.add(sub)
1337 stack.append(sub)
1338 if not getattr(sub, "__abstractmethods__", None): 1338 ↛ 1333line 1338 didn't jump to line 1333 because the condition on line 1338 was always true
1339 yield sub
1342def check_bounds_contains_broadcasting(bounds: Bounds) -> None:
1343 """Verify that `~lsst.images.Bounds.contains` accepts array-like inputs.
1345 Uses the scalar overload as the reference and checks that 1-D arrays,
1346 list inputs, mixed scalar-plus-array inputs, 2-D broadcast inputs, and
1347 `.XY` / `.YX` positional-argument forms all produce results consistent
1348 with calling the scalar overload on each ``(x, y)`` pair individually.
1350 Parameters
1351 ----------
1352 bounds
1353 The `~lsst.images.Bounds` implementation to exercise.
1354 """
1355 bbox = bounds.bbox
1356 # One point outside each boundary, one on each boundary, one in the
1357 # interior. The boundary points (start and stop) straddle inside/outside
1358 # (start is inside, stop is outside) so we are guaranteed a mix of True
1359 # and False results without hard-coding expected values.
1360 ys = np.array(
1361 [
1362 bbox.y.start - 1,
1363 bbox.y.start,
1364 (bbox.y.start + bbox.y.stop) // 2,
1365 bbox.y.stop,
1366 ]
1367 )
1368 xs = np.array(
1369 [
1370 bbox.x.start - 1,
1371 bbox.x.start,
1372 (bbox.x.start + bbox.x.stop) // 2,
1373 bbox.x.stop,
1374 ]
1375 )
1376 # 2-D reference: expected[i, j] == contains(x=xs[j], y=ys[i]).
1377 expected = np.array(
1378 [[bounds.contains(x=int(xi), y=int(yi)) for xi in xs] for yi in ys],
1379 dtype=bool,
1380 )
1381 # 1-D ndarray (diagonal pairs).
1382 np.testing.assert_array_equal(
1383 bounds.contains(x=xs, y=ys),
1384 np.diagonal(expected),
1385 )
1386 # list inputs (array-like).
1387 np.testing.assert_array_equal(
1388 bounds.contains(x=xs.tolist(), y=ys.tolist()),
1389 np.diagonal(expected),
1390 )
1391 # Mixed: scalar y, 1-D array x — produces a 1-D result.
1392 fixed_yi = 1 # index into ys; ys[1] == bbox.y.start (on the boundary)
1393 np.testing.assert_array_equal(
1394 bounds.contains(x=xs, y=int(ys[fixed_yi])),
1395 expected[fixed_yi, :],
1396 )
1397 # Mixed: 1-D array y, scalar x.
1398 fixed_xi = 1 # index into xs; xs[1] == bbox.x.start (on the boundary)
1399 np.testing.assert_array_equal(
1400 bounds.contains(x=int(xs[fixed_xi]), y=ys),
1401 expected[:, fixed_xi],
1402 )
1403 # Float scalars: results must match the int-scalar reference.
1404 assert bounds.contains(x=float(xs[fixed_xi]), y=float(ys[fixed_yi])) == expected[fixed_yi, fixed_xi]
1405 # XY / YX scalar: results must match the keyword-scalar reference.
1406 assert bounds.contains(XY(x=int(xs[fixed_xi]), y=int(ys[fixed_yi]))) == expected[fixed_yi, fixed_xi]
1407 assert bounds.contains(YX(y=int(ys[fixed_yi]), x=int(xs[fixed_xi]))) == expected[fixed_yi, fixed_xi]
1408 # XY / YX array: results must match the 1-D ndarray reference.
1409 np.testing.assert_array_equal(bounds.contains(XY(x=xs, y=ys)), np.diagonal(expected))
1410 np.testing.assert_array_equal(bounds.contains(YX(y=ys, x=xs)), np.diagonal(expected))
1411 # 2-D broadcast: y shape (N, 1) × x shape (1, M) → (N, M).
1412 np.testing.assert_array_equal(
1413 bounds.contains(x=xs.reshape(1, -1), y=ys.reshape(-1, 1)),
1414 expected,
1415 )
1418def check_archive_tree_class_invariants(tree_type: type[ArchiveTree]) -> None:
1419 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
1420 well-formed schema-version constants and an in-memory type.
1422 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
1423 ``PUBLIC_TYPE`` are present and well-typed, that the version is
1424 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
1425 schema major.
1427 Parameters
1428 ----------
1429 tree_type
1430 The concrete `.serialization.ArchiveTree` subclass to check.
1431 """
1432 assert hasattr(tree_type, "SCHEMA_NAME"), f"{tree_type.__name__} lacks SCHEMA_NAME"
1433 assert hasattr(tree_type, "SCHEMA_VERSION"), f"{tree_type.__name__} lacks SCHEMA_VERSION"
1434 assert hasattr(tree_type, "MIN_READ_VERSION"), f"{tree_type.__name__} lacks MIN_READ_VERSION"
1435 assert hasattr(tree_type, "PUBLIC_TYPE"), f"{tree_type.__name__} lacks PUBLIC_TYPE"
1436 assert isinstance(tree_type.SCHEMA_NAME, str)
1437 assert len(tree_type.SCHEMA_NAME) > 0
1438 assert re.fullmatch(r"^\d+\.\d+\.\d+$", tree_type.SCHEMA_VERSION)
1439 assert isinstance(tree_type.MIN_READ_VERSION, int)
1440 assert tree_type.MIN_READ_VERSION >= 1
1441 assert isinstance(tree_type.PUBLIC_TYPE, type)
1442 major = int(tree_type.SCHEMA_VERSION.split(".")[0])
1443 assert tree_type.MIN_READ_VERSION <= major