Coverage for python/lsst/images/tests/_checks.py: 49%
459 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-08-30 09:30 +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 assert prov is not None, f"Legacy cell {cell_index} has inputs, so its provenance is not empty."
650 legacy_table = astropy.table.Table(
651 rows=[
652 [
653 ids.instrument,
654 ids.visit,
655 ids.detector,
656 ids.day_obs,
657 ids.physical_filter,
658 legacy_input.overlaps_center,
659 legacy_input.overlap_fraction,
660 legacy_input.weight,
661 legacy_input.psf_shape.getIxx(),
662 legacy_input.psf_shape.getIyy(),
663 legacy_input.psf_shape.getIxy(),
664 legacy_input.psf_shape_flag,
665 ]
666 for ids, legacy_input in legacy_cell.inputs.items()
667 ],
668 dtype=[
669 np.object_,
670 np.uint64,
671 np.uint16,
672 np.uint32,
673 np.object_,
674 np.bool_,
675 np.float64,
676 np.float64,
677 np.float64,
678 np.float64,
679 np.float64,
680 np.bool_,
681 ],
682 names=[
683 "instrument",
684 "visit",
685 "detector",
686 "day_obs",
687 "physical_filter",
688 "overlaps_center",
689 "overlap_fraction",
690 "weight",
691 "psf_shape_xx",
692 "psf_shape_yy",
693 "psf_shape_xy",
694 "psf_shape_flag",
695 ],
696 )
697 # For a single cell all 'inputs' are also 'contributions'.
698 assert len(legacy_cell.inputs) == len(prov.inputs)
699 assert len(legacy_cell.inputs) == len(prov.contributions)
700 prov.inputs.sort(["instrument", "visit", "detector"])
701 prov.contributions.sort(["instrument", "visit", "detector"])
702 legacy_table.sort(["instrument", "visit", "detector"])
703 np.testing.assert_array_equal(prov.inputs["instrument"], prov.contributions["instrument"])
704 np.testing.assert_array_equal(prov.inputs["visit"], prov.contributions["visit"])
705 np.testing.assert_array_equal(prov.inputs["detector"], prov.contributions["detector"])
706 np.testing.assert_array_equal(prov.inputs["instrument"], legacy_table["instrument"])
707 np.testing.assert_array_equal(prov.inputs["visit"], legacy_table["visit"])
708 np.testing.assert_array_equal(prov.inputs["detector"], legacy_table["detector"])
709 np.testing.assert_array_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"])
710 np.testing.assert_array_equal(prov.inputs["day_obs"], legacy_table["day_obs"])
711 np.testing.assert_array_equal(prov.contributions["overlaps_center"], legacy_table["overlaps_center"])
712 np.testing.assert_array_equal(
713 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"]
714 )
715 np.testing.assert_array_equal(prov.contributions["weight"], legacy_table["weight"])
716 np.testing.assert_array_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"])
717 np.testing.assert_array_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"])
718 np.testing.assert_array_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"])
719 np.testing.assert_array_equal(prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"])
720 for row in prov.inputs:
721 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"})
722 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key]
723 assert legacy_polygon == row["polygon"].to_legacy()
726def compare_psf_to_legacy(
727 psf: PointSpreadFunction,
728 legacy_psf: Any,
729 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
730 expect_legacy_raise_on_out_of_bounds: bool = False,
731) -> int:
732 """Compare a PSF model object to its legacy interface.
734 Parameters
735 ----------
736 psf
737 Point-spread function to test.
738 legacy_psf
739 Legacy `lsst.afw.detection.Psf` instance to compare with.
740 points
741 Points to evaluate the PSFs at. If not provided, the intersection of
742 the PSF bounding boxes are used to generate points on a grid.
743 expect_legacy_raise_on_out_of_bounds
744 If `True`, expect ``legacy_psf`` to raise
745 `lsst.afw.detection.InvalidPsfError` when evaluated at a position
746 considered out-of-bounds by ``psf``.
748 Returns
749 -------
750 `int`
751 The number of points actually tested.
752 """
753 from lsst.afw.detection import InvalidPsfError
755 if points is None:
756 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel)
757 legacy_points = arrays_to_legacy_points(points.x, points.y)
758 n_points_tested: int = 0
759 for p in legacy_points:
760 if not psf.bounds.contains(x=p.x, y=p.y):
761 if expect_legacy_raise_on_out_of_bounds:
762 with pytest.raises(InvalidPsfError):
763 legacy_psf.computeKernelImage(p)
764 continue
765 assert psf.kernel_bbox == Box.from_legacy(legacy_psf.computeKernelBBox(p))
766 assert psf.compute_kernel_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeKernelImage(p))
767 assert psf.compute_stellar_bbox(x=p.x, y=p.y) == Box.from_legacy(legacy_psf.computeImageBBox(p))
768 assert psf.compute_stellar_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeImage(p))
769 n_points_tested += 1
770 return n_points_tested
773def compare_field_to_legacy(
774 field: BaseField,
775 legacy_field: Any,
776 subimage_bbox: Box,
777) -> None:
778 """Test a Field object by comparing it to an equivalent
779 `lsst.afw.math.BoundedField`.
781 Parameters
782 ----------
783 field
784 Field to test.
785 legacy_field : ``lsst.afw.math.BoundedField``
786 Equivalent legacy bounded field.
787 subimage_bbox
788 Bounding box for full-image tests.
789 """
790 from lsst.afw.math import BoundedField as LegacyBoundedField
792 assert field.bounds.bbox == Box.from_legacy(legacy_field.getBBox())
793 # Pixel coordinates to test the numpy array interface with.
794 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
795 if not isinstance(field.bounds, Box): 795 ↛ 796line 795 didn't jump to line 796 because the condition on line 795 was never true
796 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y)
797 pixel_xy = pixel_xy.map(lambda v: v[mask])
798 try:
799 assert_close(
800 field(x=pixel_xy.x, y=pixel_xy.y),
801 legacy_field.evaluate(pixel_xy.x, pixel_xy.y),
802 equal_nan=True,
803 )
804 except AssertionError as err:
805 err.add_note(f"evaluated at {pixel_xy}")
806 raise
807 if not isinstance(legacy_field, LegacyBoundedField): 807 ↛ 810line 807 didn't jump to line 810 because the condition on line 807 was never true
808 # Legacy StitchedApertureCorrection objects are not true BoundedFields
809 # and don't have addToImage.
810 return
811 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
812 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
813 assert_images_equal(
814 field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13
815 )
818def compare_aperture_corrections_to_legacy(
819 aperture_corrections: Mapping[str, BaseField],
820 legacy_ap_corr_map: Any,
821 subimage_bbox: Box,
822) -> None:
823 """Test an aperture correction `dict` by comparing it to an equivalent
824 `lsst.afw.image.ApCorrMap`.
826 Parameters
827 ----------
828 aperture_corrections
829 Dictionary to test.
830 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
831 Equivalent legacy aperture correction map.
832 subimage_bbox
833 Bounding box for full-image tests.
834 """
835 assert aperture_corrections.keys() == set(legacy_ap_corr_map.keys())
836 for name, field in aperture_corrections.items():
837 compare_field_to_legacy(field, legacy_ap_corr_map[name], subimage_bbox)
840def compare_observation_summary_stats_to_legacy(
841 summary_stats: ObservationSummaryStats,
842 legacy_summary_stats: Any,
843) -> None:
844 """Test an ObservationSummaryStats object by comparing it to an equivalent
845 `lsst.afw.image.ExposureSummaryStats`.
847 Parameters
848 ----------
849 summary_stats
850 Struct to test.
851 legacy_summary_stats : ``lsst.afw.image.ExposureSummaryStats``
852 Equivalent legacy struct.
853 """
854 for field in dataclasses.fields(legacy_summary_stats):
855 # Always skip version since that field has no meaning in the new
856 # type.
857 if field.name == "version":
858 continue
859 a = getattr(legacy_summary_stats, field.name)
860 b = getattr(summary_stats, field.name)
861 if isinstance(b, tuple):
862 for ai, bi in zip(a, b):
863 assert ai == bi or (math.isnan(ai) and math.isnan(bi)), f"{field.name}: {a} != {b}"
864 else:
865 assert a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}"
868def compare_sky_projection_to_legacy_wcs[F: Frame](
869 sky_projection: SkyProjection[F],
870 legacy_wcs: Any,
871 pixel_frame: F,
872 subimage_bbox: Box,
873 is_fits: bool = False,
874) -> None:
875 """Test a Projection object by comparing it to an equivalent
876 `lsst.afw.geom.SkyWcs`.
878 Parameters
879 ----------
880 sky_projection
881 Projection to test.
882 legacy_wcs : ``lsst.afw.geom.SkyWcs``
883 Equivalent legacy WCS.
884 pixel_frame
885 Expected pixel frame for the sky_projection.
886 subimage_bbox
887 Bounding box of points to generate for tests.
888 is_fits
889 Whether this sky_projection is expected to be exactly representable as
890 a FITS WCS. If `False` it is assumed to have a FITS approximation
891 attached instead.
892 """
893 # Pixel coordinates to test on over the subimage region of interest:
894 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
895 # Array indices of those pixel values (subtract off bbox starts):
896 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
897 sky_coords = legacy_coords_to_astropy(
898 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
899 )
900 # Test transforming with the Projection itself, which also tests its
901 # nested Transform and an Astropy High-Level WCS view with no origin
902 # change.
903 check_projection(sky_projection, pixel_xy, sky_coords, pixel_frame)
904 # Also test the Astropy High-Level WCS view with an origin change to
905 # array indices.
906 check_astropy_wcs_interface(
907 sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
908 )
909 if is_fits: 909 ↛ 928line 909 didn't jump to line 928 because the condition on line 909 was always true
910 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
911 assert fits_wcs is not None
912 check_astropy_wcs_interface(fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
913 # Use that FITS approximation to check that we can make a
914 # Projection from a FITS WCS, too.
915 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
916 check_projection(
917 fits_projection,
918 subimage_array_xy,
919 sky_coords,
920 pixel_frame,
921 pixel_atol=1e-5,
922 )
923 # We want Projections we create from a FITS WCS to be backed by an
924 # AST FrameSet so we can convert them into legacy
925 # `lsst.afw.geom.SkyWcs` objects if desired.
926 assert "Begin FrameSet" in fits_projection.show()
927 else:
928 assert sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False) is None
929 # The legacy SkyWcs should instead have a FITS approximation
930 # attached; run the same tests on that.
931 assert sky_projection.fits_approximation is not None
932 compare_sky_projection_to_legacy_wcs(
933 sky_projection.fits_approximation,
934 legacy_wcs.getFitsApproximation(),
935 pixel_frame,
936 subimage_bbox,
937 is_fits=True,
938 )
941def check_transform[I: Frame, O: Frame](
942 transform: Transform[I, O],
943 input_xy: XY[np.ndarray],
944 output_xy: XY[np.ndarray],
945 in_frame: Frame,
946 out_frame: Frame,
947 *,
948 check_inverted: bool = True,
949 in_atol: u.Quantity | None = None,
950 out_atol: u.Quantity | None = None,
951) -> None:
952 """Test Transform against known arrays of input and output points.
954 Parameters
955 ----------
956 transform
957 Transform to test.
958 input_xy
959 Arrays of input points.
960 output_xy
961 Arrays of output points.
962 in_frame
963 Expected input frame.
964 out_frame
965 Expected output frame.
966 check_inverted
967 If `True`, recurse (once) to test the inverse transform.
968 in_atol
969 Expected absolute precision of input points.
970 out_atol
971 Expected absolute precision of output points.
972 """
973 assert transform.in_frame == in_frame
974 assert transform.out_frame == out_frame
975 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None
976 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None
977 # Test array interfaces.
978 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
979 assert_close(test_output_xy.x, output_xy.x, atol=out_atol_v)
980 assert_close(test_output_xy.y, output_xy.y, atol=out_atol_v)
981 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
982 assert_close(test_input_xy.x, input_xy.x, atol=in_atol_v)
983 assert_close(test_input_xy.y, input_xy.y, atol=in_atol_v)
984 # Test scalar interfaces with numpy scalars.
985 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
986 assert_close(transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v)
987 assert_close(transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v)
988 assert_close(transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v)
989 assert_close(transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v)
990 # Test quantity array interfaces.
991 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
992 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
993 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
994 assert_close(test_output_q_xy.x, output_q_xy.x, atol=out_atol)
995 assert_close(test_output_q_xy.y, output_q_xy.y, atol=out_atol)
996 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
997 assert_close(test_input_q_xy.x, input_q_xy.x, atol=in_atol)
998 assert_close(test_input_q_xy.y, input_q_xy.y, atol=in_atol)
999 # Test quantity scalar interfaces.
1000 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
1001 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
1002 ):
1003 assert_close(transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol)
1004 assert_close(transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol)
1005 assert_close(transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol)
1006 assert_close(transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol)
1007 if check_inverted:
1008 # Test the inverse transform.
1009 check_transform(
1010 transform.inverted(),
1011 output_xy,
1012 input_xy,
1013 out_frame,
1014 in_frame,
1015 check_inverted=False,
1016 out_atol=in_atol,
1017 in_atol=out_atol,
1018 )
1021def check_projection[P: Frame](
1022 sky_projection: SkyProjection[P],
1023 pixel_xy: XY[np.ndarray],
1024 sky_coords: SkyCoord,
1025 pixel_frame: Frame,
1026 *,
1027 pixel_atol: float | None = None,
1028 sky_atol: u.Quantity | None = None,
1029) -> None:
1030 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1031 coordinates.
1033 Parameters
1034 ----------
1035 sky_projection
1036 Projection to test.
1037 pixel_xy
1038 Arrays of pixel coordinates.
1039 sky_coords
1040 Corresponding sky coordinates.
1041 pixel_frame
1042 Expected pixel frame.
1043 pixel_atol
1044 Expected absolute precision of pixel points.
1045 sky_atol
1046 Expected absolute precision of sky coordinates.
1047 """
1048 assert sky_projection.pixel_frame == pixel_frame
1049 assert sky_projection.sky_frame == SkyFrame.ICRS
1050 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None
1051 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None
1052 # Test array interfaces.
1053 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1054 assert_close(test_pixel_xy.x, pixel_xy.x, atol=pixel_atol)
1055 assert_close(test_pixel_xy.y, pixel_xy.y, atol=pixel_atol)
1056 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1057 assert_close(test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v)
1058 assert_close(test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v)
1059 # Test scalar interfaces.
1060 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1061 assert_close(sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol)
1062 assert_close(sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol)
1063 assert_close(sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v)
1064 assert_close(sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v)
1065 # Test the underlying Transform object.
1066 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1067 check_transform(
1068 sky_projection.pixel_to_sky_transform,
1069 pixel_xy,
1070 sky_xy,
1071 pixel_frame,
1072 SkyFrame.ICRS,
1073 check_inverted=False,
1074 in_atol=pixel_atol_q,
1075 out_atol=sky_atol,
1076 )
1077 check_transform(
1078 sky_projection.sky_to_pixel_transform,
1079 sky_xy,
1080 pixel_xy,
1081 SkyFrame.ICRS,
1082 pixel_frame,
1083 check_inverted=False,
1084 in_atol=sky_atol,
1085 out_atol=pixel_atol_q,
1086 )
1087 # Test the Astropy interface adapter.
1088 check_astropy_wcs_interface(
1089 sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1090 )
1093def assert_sky_projections_equal(
1094 a: SkyProjection[Any] | None,
1095 b: SkyProjection[Any] | None,
1096 expect_identity: bool | None = None,
1097) -> None:
1098 """Test that two `.SkyProjection` instances are equivalent.
1100 Parameters
1101 ----------
1102 a
1103 First sky projection to compare.
1104 b
1105 Second sky projection to compare.
1106 expect_identity
1107 If not `None`, assert whether ``a`` and ``b`` are the same object.
1108 """
1109 if a is None and b is None:
1110 return
1111 assert a is not None
1112 assert b is not None
1113 match expect_identity:
1114 case True: 1114 ↛ 1115line 1114 didn't jump to line 1115 because the pattern on line 1114 never matched
1115 assert a is b
1116 return
1117 case False:
1118 assert a is not b
1119 case None if a is b:
1120 return
1121 assert a == b
1124def check_astropy_wcs_interface(
1125 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1126 pixel_xy: XY[np.ndarray],
1127 sky_coords: SkyCoord,
1128 *,
1129 pixel_atol: float | None = None,
1130 sky_atol: u.Quantity | None = None,
1131) -> None:
1132 """Test an Astropy WCS instance against known arrays of pixel and
1133 sky coordinates.
1135 Parameters
1136 ----------
1137 wcs
1138 Astropy WCS object to test.
1139 pixel_xy
1140 Arrays of pixel coordinates.
1141 sky_coords
1142 Corresponding sky coordinates.
1143 pixel_atol
1144 Expected absolute precision of pixel points.
1145 sky_atol
1146 Expected absolute precision of sky coordinates.
1147 """
1148 test_x, test_y = wcs.world_to_pixel(sky_coords)
1149 assert_close(test_x, pixel_xy.x, atol=pixel_atol)
1150 assert_close(test_y, pixel_xy.y, atol=pixel_atol)
1151 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
1152 assert_close(test_sky_coords.ra, sky_coords.ra, atol=sky_atol)
1153 assert_close(test_sky_coords.dec, sky_coords.dec, atol=sky_atol)
1156def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
1157 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array.
1159 Parameters
1160 ----------
1161 legacy_points
1162 Legacy ``lsst.geom.Point2D`` objects to convert.
1163 """
1164 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
1167def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
1168 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
1169 coordinate object.
1171 Parameters
1172 ----------
1173 legacy_coords
1174 Legacy ``lsst.geom.SpherePoint`` objects to convert.
1175 """
1176 return SkyCoord(
1177 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
1178 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
1179 )
1182def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
1183 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``.
1185 Parameters
1186 ----------
1187 x
1188 X coordinates of the points.
1189 y
1190 Y coordinates of the points.
1191 """
1192 from lsst.geom import Point2D
1194 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
1197def compare_amplifier_to_legacy(
1198 amplifier: Amplifier,
1199 legacy_amplifier: Any,
1200 *,
1201 is_raw_assembled: bool,
1202 expect_nominal_calibrations: bool = True,
1203) -> None:
1204 """Compare an `~.cameras.Amplifier` to a legacy
1205 `lsst.afw.cameraGeom.Amplifier`.
1207 Parameters
1208 ----------
1209 amplifier
1210 Amplifier to compare.
1211 legacy_amplifier
1212 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against.
1213 is_raw_assembled
1214 Whether the raw geometry is expected to be the assembled-raw
1215 geometry (`True`) or the unassembled-raw geometry (`False`).
1216 expect_nominal_calibrations
1217 Whether the amplifier is expected to carry nominal calibrations.
1218 """
1219 assert legacy_amplifier.getName() == amplifier.name
1220 assert Box.from_legacy(legacy_amplifier.getBBox()) == amplifier.bbox
1221 if is_raw_assembled:
1222 raw_geom = amplifier.assembled_raw_geometry
1223 else:
1224 raw_geom = amplifier.unassembled_raw_geometry
1225 assert raw_geom is not None
1226 assert ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()) == raw_geom.readout_corner
1227 assert Box.from_legacy(legacy_amplifier.getRawBBox()) == raw_geom.bbox
1228 assert Box.from_legacy(legacy_amplifier.getRawDataBBox()) == raw_geom.data_bbox
1229 assert legacy_amplifier.getRawFlipX() == raw_geom.flip_x
1230 assert legacy_amplifier.getRawFlipY() == raw_geom.flip_y
1231 assert legacy_amplifier.getRawXYOffset().getX() == raw_geom.x_offset
1232 assert legacy_amplifier.getRawXYOffset().getY() == raw_geom.y_offset
1233 assert (
1234 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()) == raw_geom.horizontal_overscan_bbox
1235 )
1236 assert Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()) == raw_geom.vertical_overscan_bbox
1237 assert Box.from_legacy(legacy_amplifier.getRawPrescanBBox()) == raw_geom.horizontal_prescan_bbox
1238 if expect_nominal_calibrations:
1239 assert amplifier.nominal_calibrations is not None
1240 assert_equal_allow_nan(legacy_amplifier.getGain(), amplifier.nominal_calibrations.gain)
1241 assert_equal_allow_nan(legacy_amplifier.getReadNoise(), amplifier.nominal_calibrations.read_noise)
1242 assert_equal_allow_nan(legacy_amplifier.getSaturation(), amplifier.nominal_calibrations.saturation)
1243 assert_equal_allow_nan(
1244 legacy_amplifier.getSuspectLevel(), amplifier.nominal_calibrations.suspect_level
1245 )
1246 np.testing.assert_array_equal(
1247 legacy_amplifier.getLinearityCoeffs(), amplifier.nominal_calibrations.linearity_coefficients
1248 )
1249 assert legacy_amplifier.getLinearityType() == amplifier.nominal_calibrations.linearity_type
1252def compare_detector_to_legacy(
1253 detector: Detector,
1254 legacy_detector: Any,
1255 *,
1256 is_raw_assembled: bool,
1257 expect_nominal_calibrations: bool = True,
1258) -> None:
1259 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`.
1261 Parameters
1262 ----------
1263 detector
1264 Detector to compare.
1265 legacy_detector
1266 Legacy `lsst.afw.cameraGeom.Detector` to compare against.
1267 is_raw_assembled
1268 Whether the raw geometry is expected to be the assembled-raw
1269 geometry (`True`) or the unassembled-raw geometry (`False`).
1270 expect_nominal_calibrations
1271 Whether the detector's amplifiers are expected to carry nominal
1272 calibrations.
1273 """
1274 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
1276 assert legacy_detector.getName() == detector.name
1277 assert legacy_detector.getId() == detector.id
1278 assert DetectorType.from_legacy(legacy_detector.getType()) == detector.type
1279 assert Box.from_legacy(legacy_detector.getBBox()) == detector.bbox
1280 assert legacy_detector.getSerial() == detector.serial
1281 legacy_orientation = legacy_detector.getOrientation()
1282 assert legacy_orientation.getFpPosition3().getX() == detector.orientation.focal_plane_x
1283 assert legacy_orientation.getFpPosition3().getY() == detector.orientation.focal_plane_y
1284 assert legacy_orientation.getFpPosition3().getZ() == detector.orientation.focal_plane_z
1285 assert legacy_orientation.getReferencePoint().getX() == detector.orientation.pixel_reference_x
1286 assert legacy_orientation.getReferencePoint().getY() == detector.orientation.pixel_reference_y
1287 assert legacy_orientation.getYaw().asRadians() == detector.orientation.yaw.to_value(u.rad)
1288 assert legacy_orientation.getPitch().asRadians() == detector.orientation.pitch.to_value(u.rad)
1289 assert legacy_orientation.getRoll().asRadians() == detector.orientation.roll.to_value(u.rad)
1290 assert legacy_detector.getPixelSize().getX() == detector.pixel_size
1291 assert legacy_detector.getPhysicalType() == detector.physical_type
1292 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
1293 compare_amplifier_to_legacy(
1294 amplifier,
1295 legacy_amplifier,
1296 is_raw_assembled=is_raw_assembled,
1297 expect_nominal_calibrations=expect_nominal_calibrations,
1298 )
1299 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
1300 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
1301 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
1302 check_transform(
1303 detector.to_focal_plane,
1304 pixel_xy,
1305 legacy_points_to_xy_array(fp_legacy_points),
1306 detector.frame,
1307 detector.to_focal_plane.out_frame,
1308 in_atol=1e-9 * u.pix,
1309 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
1310 )
1311 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
1312 check_transform(
1313 detector.to_field_angle,
1314 pixel_xy,
1315 legacy_points_to_xy_array(fa_legacy_points),
1316 detector.frame,
1317 detector.to_field_angle.out_frame,
1318 in_atol=1e-9 * u.pix,
1319 out_atol=1e-7 * u.arcsec,
1320 )
1323def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
1324 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
1326 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
1327 classes. Importing this module already imports every ``lsst.images``
1328 module that defines a subclass, so the tree is fully populated by the time
1329 this is called.
1331 This discovery is deliberately separate from
1332 `check_archive_tree_class_invariants` so that the per-class check stays
1333 usable on a single class even if this metaprogramming is removed later.
1334 """
1335 seen: set[type] = set()
1336 stack: list[type] = [ArchiveTree]
1337 while stack:
1338 kls = stack.pop()
1339 for sub in kls.__subclasses__():
1340 if sub in seen: 1340 ↛ 1341line 1340 didn't jump to line 1341 because the condition on line 1340 was never true
1341 continue
1342 seen.add(sub)
1343 stack.append(sub)
1344 if not getattr(sub, "__abstractmethods__", None): 1344 ↛ 1339line 1344 didn't jump to line 1339 because the condition on line 1344 was always true
1345 yield sub
1348def check_bounds_contains_broadcasting(bounds: Bounds) -> None:
1349 """Verify that `~lsst.images.Bounds.contains` accepts array-like inputs.
1351 Uses the scalar overload as the reference and checks that 1-D arrays,
1352 list inputs, mixed scalar-plus-array inputs, 2-D broadcast inputs, and
1353 `.XY` / `.YX` positional-argument forms all produce results consistent
1354 with calling the scalar overload on each ``(x, y)`` pair individually.
1356 Parameters
1357 ----------
1358 bounds
1359 The `~lsst.images.Bounds` implementation to exercise.
1360 """
1361 bbox = bounds.bbox
1362 # One point outside each boundary, one on each boundary, one in the
1363 # interior. The boundary points (start and stop) straddle inside/outside
1364 # (start is inside, stop is outside) so we are guaranteed a mix of True
1365 # and False results without hard-coding expected values.
1366 ys = np.array(
1367 [
1368 bbox.y.start - 1,
1369 bbox.y.start,
1370 (bbox.y.start + bbox.y.stop) // 2,
1371 bbox.y.stop,
1372 ]
1373 )
1374 xs = np.array(
1375 [
1376 bbox.x.start - 1,
1377 bbox.x.start,
1378 (bbox.x.start + bbox.x.stop) // 2,
1379 bbox.x.stop,
1380 ]
1381 )
1382 # 2-D reference: expected[i, j] == contains(x=xs[j], y=ys[i]).
1383 expected = np.array(
1384 [[bounds.contains(x=int(xi), y=int(yi)) for xi in xs] for yi in ys],
1385 dtype=bool,
1386 )
1387 # 1-D ndarray (diagonal pairs).
1388 np.testing.assert_array_equal(
1389 bounds.contains(x=xs, y=ys),
1390 np.diagonal(expected),
1391 )
1392 # list inputs (array-like).
1393 np.testing.assert_array_equal(
1394 bounds.contains(x=xs.tolist(), y=ys.tolist()),
1395 np.diagonal(expected),
1396 )
1397 # Mixed: scalar y, 1-D array x — produces a 1-D result.
1398 fixed_yi = 1 # index into ys; ys[1] == bbox.y.start (on the boundary)
1399 np.testing.assert_array_equal(
1400 bounds.contains(x=xs, y=int(ys[fixed_yi])),
1401 expected[fixed_yi, :],
1402 )
1403 # Mixed: 1-D array y, scalar x.
1404 fixed_xi = 1 # index into xs; xs[1] == bbox.x.start (on the boundary)
1405 np.testing.assert_array_equal(
1406 bounds.contains(x=int(xs[fixed_xi]), y=ys),
1407 expected[:, fixed_xi],
1408 )
1409 # Float scalars: results must match the int-scalar reference.
1410 assert bounds.contains(x=float(xs[fixed_xi]), y=float(ys[fixed_yi])) == expected[fixed_yi, fixed_xi]
1411 # XY / YX scalar: results must match the keyword-scalar reference.
1412 assert bounds.contains(XY(x=int(xs[fixed_xi]), y=int(ys[fixed_yi]))) == expected[fixed_yi, fixed_xi]
1413 assert bounds.contains(YX(y=int(ys[fixed_yi]), x=int(xs[fixed_xi]))) == expected[fixed_yi, fixed_xi]
1414 # XY / YX array: results must match the 1-D ndarray reference.
1415 np.testing.assert_array_equal(bounds.contains(XY(x=xs, y=ys)), np.diagonal(expected))
1416 np.testing.assert_array_equal(bounds.contains(YX(y=ys, x=xs)), np.diagonal(expected))
1417 # 2-D broadcast: y shape (N, 1) × x shape (1, M) → (N, M).
1418 np.testing.assert_array_equal(
1419 bounds.contains(x=xs.reshape(1, -1), y=ys.reshape(-1, 1)),
1420 expected,
1421 )
1424def check_archive_tree_class_invariants(tree_type: type[ArchiveTree]) -> None:
1425 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
1426 well-formed schema-version constants and an in-memory type.
1428 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
1429 ``PUBLIC_TYPE`` are present and well-typed, that the version is
1430 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
1431 schema major.
1433 Parameters
1434 ----------
1435 tree_type
1436 The concrete `.serialization.ArchiveTree` subclass to check.
1437 """
1438 assert hasattr(tree_type, "SCHEMA_NAME"), f"{tree_type.__name__} lacks SCHEMA_NAME"
1439 assert hasattr(tree_type, "SCHEMA_VERSION"), f"{tree_type.__name__} lacks SCHEMA_VERSION"
1440 assert hasattr(tree_type, "MIN_READ_VERSION"), f"{tree_type.__name__} lacks MIN_READ_VERSION"
1441 assert hasattr(tree_type, "PUBLIC_TYPE"), f"{tree_type.__name__} lacks PUBLIC_TYPE"
1442 assert isinstance(tree_type.SCHEMA_NAME, str)
1443 assert len(tree_type.SCHEMA_NAME) > 0
1444 # Allow an optional PEP 440 development-release suffix (e.g. 1.0.0.dev0)
1445 # for schemas still in development.
1446 assert re.fullmatch(r"^\d+\.\d+\.\d+(\.dev\d+)?$", tree_type.SCHEMA_VERSION)
1447 assert isinstance(tree_type.MIN_READ_VERSION, int)
1448 assert tree_type.MIN_READ_VERSION >= 1
1449 assert isinstance(tree_type.PUBLIC_TYPE, type)
1450 major = int(tree_type.SCHEMA_VERSION.split(".")[0])
1451 assert tree_type.MIN_READ_VERSION <= major