Coverage for python/lsst/images/tests/_checks.py: 40%
441 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-25 08:24 +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_projection",
28 "check_transform",
29 "compare_amplifier_to_legacy",
30 "compare_aperture_corrections_to_legacy",
31 "compare_cell_coadd_to_legacy",
32 "compare_detector_to_legacy",
33 "compare_field_to_legacy",
34 "compare_image_to_legacy",
35 "compare_mask_to_legacy",
36 "compare_masked_image_to_legacy",
37 "compare_observation_summary_stats_to_legacy",
38 "compare_photo_calib_to_legacy",
39 "compare_psf_to_legacy",
40 "compare_sky_projection_to_legacy_wcs",
41 "compare_visit_image_to_legacy",
42 "iter_concrete_archive_tree_subclasses",
43 "legacy_coords_to_astropy",
44 "legacy_points_to_xy_array",
45)
47import dataclasses
48import math
49import unittest
50from collections.abc import Iterator, Mapping
51from typing import TYPE_CHECKING, Any, Literal, cast
53import astropy.units as u
54import astropy.wcs.wcsapi
55import numpy as np
56from astropy.coordinates import SkyCoord
58from .._geom import XY, YX, Box
59from .._image import Image
60from .._mask import Mask, MaskPlane
61from .._masked_image import MaskedImage
62from .._observation_summary_stats import ObservationSummaryStats
63from .._transforms import DetectorFrame, Frame, SkyFrame, SkyProjection, TractFrame, Transform
64from .._visit_image import VisitImage
65from ..cameras import Amplifier, Detector, DetectorType, ReadoutCorner
66from ..cells import CellCoadd, CellIJ, CoaddProvenance
67from ..fields import BaseField, ChebyshevField
68from ..psfs import PointSpreadFunction
69from ..serialization import ArchiveTree
71if TYPE_CHECKING:
72 try:
73 from lsst.cell_coadds import MultipleCellCoadd
74 except ImportError:
75 type MultipleCellCoadd = Any # type: ignore[no-redef]
76 try:
77 from lsst.afw.image import PhotoCalib as LegacyPhotoCalib
78 except ImportError:
79 type LegacyPhotoCalib = Any # type: ignore[no-redef]
82def assert_close(
83 tc: unittest.TestCase,
84 a: np.ndarray | u.Quantity | float,
85 b: np.ndarray | u.Quantity | float,
86 **kwargs: Any,
87) -> None:
88 """Test that two arrays, floats, or quantities are almost equal.
90 Parameters
91 ----------
92 tc
93 Test case object with assert methods to use.
94 a
95 Array, float, or quantity to compare.
96 b
97 Array, float, or quantity to compare.
98 **kwargs
99 Forwarded to `astropy.units.allclose`.
100 """
101 tc.assertTrue(u.allclose(a, b, **kwargs), msg=f"{a} != {b}")
104def assert_equal_allow_nan(tc: unittest.TestCase, a: float, b: float) -> None:
105 """Test that two floating point values are equal, with nan == nan."""
106 try:
107 tc.assertEqual(a, b)
108 except AssertionError:
109 if not (math.isnan(a) and math.isnan(b)):
110 raise
113def assert_images_equal(
114 tc: unittest.TestCase,
115 a: Image,
116 b: Image,
117 *,
118 rtol: float = 0.0,
119 atol: float = 0.0,
120 expect_view: bool | Literal["array"] | None = None,
121) -> None:
122 """Assert that two images are equal or nearly equal."""
123 tc.assertEqual(a.bbox, b.bbox)
124 tc.assertEqual(a.unit, b.unit)
125 assert_sky_projections_equal(tc, a.sky_projection, b.sky_projection)
126 if expect_view is not None:
127 tc.assertEqual(np.may_share_memory(a.array, b.array), bool(expect_view))
128 if expect_view == "array":
129 tc.assertEqual(a.metadata, b.metadata)
130 else:
131 tc.assertEqual(a.metadata is b.metadata, expect_view)
132 if not expect_view:
133 assert_close(tc, a.array, b.array, atol=atol, rtol=rtol)
134 tc.assertEqual(a.metadata, b.metadata)
137def assert_masks_equal(tc: unittest.TestCase, a: Mask, b: Mask) -> None:
138 """Assert that two masks are equal or nearly equal."""
139 tc.assertEqual(a.bbox, b.bbox)
140 tc.assertEqual(a.schema, b.schema)
141 tc.assertEqual(a.metadata, b.metadata)
142 assert_sky_projections_equal(tc, a.sky_projection, b.sky_projection)
143 np.testing.assert_array_equal(a.array, b.array)
146def assert_masked_images_equal(
147 tc: unittest.TestCase,
148 a: MaskedImage,
149 b: MaskedImage,
150 *,
151 rtol: float = 0.0,
152 atol: float = 0.0,
153 expect_view: bool | None = None,
154) -> None:
155 """Assert that two masked images are equal or nearly equal."""
156 tc.assertEqual(a.metadata, b.metadata)
157 assert_sky_projections_equal(tc, a.sky_projection, b.sky_projection)
158 assert_images_equal(tc, a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view)
159 assert_masks_equal(tc, a.mask, b.mask)
160 assert_images_equal(tc, a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view)
163def assert_psfs_equal(
164 tc: unittest.TestCase,
165 psf1: PointSpreadFunction,
166 psf2: PointSpreadFunction,
167 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
168) -> int:
169 """Compare two PSF objets.
171 Parameters
172 ----------
173 tc
174 Test case object with assert methods to use.
175 psf1
176 Point-spread function to test.
177 psf2
178 The other point-spread function to test.
179 points
180 Points to evaluate the PSFs at. If not provided, the intersection of
181 the PSF bounding boxes are used to generate points on a grid.
183 Returns
184 -------
185 `int`
186 The number of points actually tested.
187 """
188 if points is None: 188 ↛ 191line 188 didn't jump to line 191 because the condition on line 188 was always true
189 points = psf1.bounds.bbox.intersection(psf2.bounds.bbox).meshgrid(3).map(np.ravel)
191 tc.assertEqual(psf1.kernel_bbox, psf2.kernel_bbox)
193 n_points_tested: int = 0
194 for x, y in zip(points.x, points.y):
195 # The two PSFs must agree on which points fall inside their input
196 # domain. Querying ``.contains`` directly (rather than relying on
197 # ``compute_kernel_image`` to raise) makes this test tolerant of
198 # implementations that do not raise on out-of-domain points -- in
199 # particular ``CellPointSpreadFunction``, where evaluating in a
200 # missing cell does not always raise ``BoundsError``.
201 contains1 = psf1.bounds.contains(x=x, y=y)
202 contains2 = psf2.bounds.contains(x=x, y=y)
203 tc.assertEqual(
204 contains1,
205 contains2,
206 f"PSFs disagree on whether ({x}, {y}) is in-bounds: psf1={contains1}, psf2={contains2}",
207 )
208 if not contains1: 208 ↛ 209line 208 didn't jump to line 209 because the condition on line 208 was never true
209 continue
210 tc.assertEqual(psf1.compute_kernel_image(x=x, y=y), psf2.compute_kernel_image(x=x, y=y))
211 tc.assertEqual(psf1.compute_stellar_bbox(x=x, y=y), psf2.compute_stellar_bbox(x=x, y=y))
212 tc.assertEqual(psf1.compute_stellar_image(x=x, y=y), psf2.compute_stellar_image(x=x, y=y))
213 n_points_tested += 1
214 return n_points_tested
217def assert_visit_images_equal(
218 tc: unittest.TestCase,
219 a: VisitImage,
220 b: VisitImage,
221 *,
222 expect_view: bool | None = None,
223) -> None:
224 """Assert that two `.VisitImage` instances carry the same persistent state.
226 Extends `assert_masked_images_equal` with the VisitImage-specific
227 attributes (PSF, filter, observation info, detector, aperture
228 corrections, photometric scaling, backgrounds, polygon bounds,
229 summary stats) so a round-trip check on a `.VisitImage` does not
230 silently miss differences in any of them.
231 """
232 assert_masked_images_equal(tc, a, b, expect_view=expect_view)
233 tc.assertEqual(a.summary_stats, b.summary_stats)
234 tc.assertEqual(a.physical_filter, b.physical_filter)
235 tc.assertEqual(a.band, b.band)
236 tc.assertEqual(a.obs_info, b.obs_info)
237 tc.assertEqual(a.detector, b.detector)
238 tc.assertEqual(dict(a.aperture_corrections), dict(b.aperture_corrections))
239 tc.assertEqual(a.photometric_scaling, b.photometric_scaling)
240 tc.assertEqual(dict(a.backgrounds), dict(b.backgrounds))
241 tc.assertEqual(a.backgrounds.subtracted, b.backgrounds.subtracted)
242 tc.assertEqual(a.bounds, b.bounds)
243 assert_psfs_equal(tc, a.psf, b.psf)
246def assert_cell_coadds_equal(
247 tc: unittest.TestCase,
248 a: CellCoadd,
249 b: CellCoadd,
250 *,
251 expect_view: bool | None = None,
252) -> None:
253 """Assert that two `.CellCoadd` instances carry the same persistent state.
255 Extends the masked-image-style equality check with the
256 CellCoadd-specific attributes (PSF, cell grid, missing cells,
257 backgrounds, patch/tract, band) so a round-trip check does not
258 silently miss differences in any of them.
259 """
260 assert_masked_images_equal(tc, a, b, expect_view=expect_view)
261 tc.assertEqual(a.band, b.band)
262 tc.assertEqual(a.patch, b.patch)
263 tc.assertEqual(a.tract, b.tract)
264 tc.assertEqual(a.grid, b.grid)
265 tc.assertEqual(a.bounds.missing, b.bounds.missing)
266 tc.assertEqual(dict(a.backgrounds), dict(b.backgrounds))
267 tc.assertEqual(a.backgrounds.subtracted, b.backgrounds.subtracted)
268 assert_psfs_equal(tc, a.psf, b.psf)
271def compare_image_to_legacy(
272 tc: unittest.TestCase, image: Image, legacy_image: Any, expect_view: bool | None = None
273) -> None:
274 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object."""
275 tc.assertEqual(image.bbox, Box.from_legacy(legacy_image.getBBox()))
276 if expect_view is not None: 276 ↛ 278line 276 didn't jump to line 278 because the condition on line 276 was always true
277 tc.assertEqual(np.may_share_memory(image.array, legacy_image.array), expect_view)
278 if not expect_view: 278 ↛ exitline 278 didn't return from function 'compare_image_to_legacy' because the condition on line 278 was always true
279 np.testing.assert_array_equal(image.array, legacy_image.array)
282def compare_mask_to_legacy(
283 tc: unittest.TestCase, mask: Mask, legacy_mask: Any, plane_map: Mapping[str, MaskPlane] | None = None
284) -> None:
285 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object."""
286 tc.assertEqual(mask.bbox, Box.from_legacy(legacy_mask.getBBox()))
287 if plane_map is None: 287 ↛ 289line 287 didn't jump to line 289 because the condition on line 287 was always true
288 plane_map = {plane.name: plane for plane in mask.schema if plane is not None}
289 for old_name, new_plane in plane_map.items():
290 np.testing.assert_array_equal(
291 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool),
292 mask.get(new_plane.name),
293 )
296def compare_masked_image_to_legacy(
297 tc: unittest.TestCase,
298 masked_image: MaskedImage,
299 legacy_masked_image: Any,
300 *,
301 plane_map: Mapping[str, MaskPlane] | None = None,
302 expect_view: bool | None = None,
303 alternates: Mapping[str, Any] | None = None,
304) -> None:
305 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage`
306 object.
308 Parameters
309 ----------
310 tc
311 Test case to use for asserts.
312 masked_image
313 New image to test.
314 legacy_masked_image
315 Legacy image to test against.
316 plane_map
317 Mapping between new and legacy mask planes.
318 expect_view
319 Whether to test that the image and variance arrays do or do not share
320 memory.
321 alternates
322 A mapping of other versions of one or more (new) components to also
323 check against the legacy versions of those components.
324 """
325 compare_image_to_legacy(tc, masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view)
326 compare_mask_to_legacy(tc, masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map)
327 compare_image_to_legacy(
328 tc, masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view
329 )
330 if alternates: 330 ↛ 331line 330 didn't jump to line 331 because the condition on line 330 was never true
331 if image := alternates.get("image"):
332 compare_image_to_legacy(tc, image, legacy_masked_image.getImage(), expect_view=expect_view)
333 if mask := alternates.get("mask"):
334 compare_mask_to_legacy(tc, mask, legacy_masked_image.getMask(), plane_map=plane_map)
335 if variance := alternates.get("variance"):
336 compare_image_to_legacy(tc, variance, legacy_masked_image.getVariance(), expect_view=expect_view)
339def compare_visit_image_to_legacy(
340 tc: unittest.TestCase,
341 visit_image: VisitImage,
342 legacy_exposure: Any,
343 *,
344 plane_map: Mapping[str, MaskPlane] | None = None,
345 expect_view: bool | None = None,
346 instrument: str,
347 visit: int,
348 detector: int,
349 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
350 alternates: Mapping[str, Any] | None = None,
351) -> None:
352 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure`
353 object.
355 Parameters
356 ----------
357 tc
358 Test case to use for asserts.
359 visit_image
360 New image to test.
361 legacy_exposure
362 Legacy image to test against.
363 plane_map
364 Mapping between new and legacy mask planes.
365 expect_view
366 Whether to test that the image and variance arrays do or do not share
367 memory.
368 instrument
369 Expected instrument name.
370 visit
371 Expected visit ID.
372 detector
373 Expected detector ID.
374 alternates
375 A mapping of other versions of one or more (new) components to also
376 check against the legacy versions of those components.
377 """
378 compare_masked_image_to_legacy(
379 tc,
380 visit_image,
381 legacy_exposure.getMaskedImage(),
382 plane_map=plane_map,
383 expect_view=expect_view,
384 alternates=alternates,
385 )
386 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox())
387 compare_sky_projection_to_legacy_wcs(
388 tc,
389 visit_image.sky_projection,
390 legacy_exposure.getWcs(),
391 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
392 visit_image.bbox,
393 )
394 tc.assertIs(visit_image.sky_projection, visit_image.mask.sky_projection)
395 tc.assertIs(visit_image.sky_projection, visit_image.variance.sky_projection)
396 compare_psf_to_legacy(tc, visit_image.psf, legacy_exposure.getPsf())
397 compare_observation_summary_stats_to_legacy(
398 tc, visit_image.summary_stats, legacy_exposure.info.getSummaryStats()
399 )
400 compare_detector_to_legacy(tc, visit_image.detector, legacy_exposure.getDetector(), is_raw_assembled=True)
401 # Make a tiny box for Field comparisons that need to make arrays; that can
402 # get expensive otherwisre.
403 tiny_bbox = detector_bbox.local[2:4, 3:6]
404 compare_aperture_corrections_to_legacy(
405 tc, visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
406 )
407 compare_photo_calib_to_legacy(
408 tc,
409 visit_image.photometric_scaling,
410 legacy_exposure.info.getPhotoCalib(),
411 applied_legacy_photo_calib=applied_legacy_photo_calib,
412 subimage_bbox=tiny_bbox,
413 )
414 if alternates:
415 if (bbox := alternates.get("bbox")) is not None:
416 tc.assertEqual(bbox, visit_image.bbox)
417 if sky_projection := alternates.get("sky_projection"):
418 compare_sky_projection_to_legacy_wcs(
419 tc,
420 sky_projection,
421 legacy_exposure.getWcs(),
422 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox),
423 visit_image.bbox,
424 )
425 if psf := alternates.get("psf"):
426 compare_psf_to_legacy(tc, psf, legacy_exposure.getPsf())
427 if summary_stats := alternates.get("summary_stats"):
428 compare_observation_summary_stats_to_legacy(
429 tc, summary_stats, legacy_exposure.info.getSummaryStats()
430 )
431 if detector_obj := alternates.get("detector"):
432 compare_detector_to_legacy(tc, detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True)
433 if obs_info := alternates.get("obs_info"):
434 visitInfo = legacy_exposure.visitInfo
435 tc.assertEqual(obs_info.instrument, visitInfo.getInstrumentLabel())
436 if aperture_corrections := alternates.get("aperture_corrections"):
437 compare_aperture_corrections_to_legacy(
438 tc, aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox
439 )
440 if (photometric_scaling := alternates.get("photometic_scaling", ...)) is not ...:
441 compare_photo_calib_to_legacy(
442 tc,
443 photometric_scaling,
444 legacy_exposure.info.getPhotoCalib(),
445 applied_legacy_photo_calib=applied_legacy_photo_calib,
446 subimage_bbox=tiny_bbox,
447 )
450def compare_photo_calib_to_legacy(
451 tc: unittest.TestCase,
452 photometric_scaling: BaseField | None,
453 legacy_photo_calib: LegacyPhotoCalib,
454 *,
455 applied_legacy_photo_calib: LegacyPhotoCalib | None = None,
456 subimage_bbox: Box,
457) -> None:
458 if legacy_photo_calib._isConstant:
459 if legacy_photo_calib.getCalibrationMean() == 1.0:
460 if applied_legacy_photo_calib is None:
461 tc.assertIsNone(photometric_scaling)
462 return
463 else:
464 legacy_photo_calib = applied_legacy_photo_calib
465 if legacy_photo_calib._isConstant:
466 assert isinstance(photometric_scaling, ChebyshevField)
467 assert_close(
468 tc, photometric_scaling.coefficients, np.array([[legacy_photo_calib.getCalibrationMean()]])
469 )
470 else:
471 assert photometric_scaling is not None
472 compare_field_to_legacy(
473 tc,
474 photometric_scaling / legacy_photo_calib.getCalibrationMean(),
475 legacy_photo_calib.computeScaledCalibration(),
476 subimage_bbox,
477 )
480def compare_cell_coadd_to_legacy(
481 tc: unittest.TestCase,
482 cell_coadd: CellCoadd,
483 legacy_cell_coadd: MultipleCellCoadd,
484 *,
485 tract_bbox: Box,
486 plane_map: Mapping[str, MaskPlane] | None = None,
487 alternates: Mapping[str, Any] | None = None,
488 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None,
489) -> None:
490 """Compare a `.cells.CellCoadd` object to a legacy
491 `lsst.cell_coadds.MultipleCellCoadd` object.
493 Parameters
494 ----------
495 tc
496 Test case to use for asserts.
497 cell_coadd
498 New coadd to test.
499 legacy_cell_coadd
500 Legacy coadd to test against.
501 tract_bbox
502 Bounding box of the full tract.
503 psf_points
504 Points to use to compare the PSFs.
505 plane_map
506 Mapping between new and legacy mask planes.
507 alternates
508 A mapping of other versions of one or more (new) components to also
509 check against the legacy versions of those components.
510 """
511 legacy_stitched = legacy_cell_coadd.stitch(cell_coadd.bbox.to_legacy())
512 compare_image_to_legacy(tc, cell_coadd.image, legacy_stitched.image, expect_view=False)
513 compare_mask_to_legacy(tc, cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map)
514 compare_image_to_legacy(tc, cell_coadd.variance, legacy_stitched.variance, expect_view=False)
515 if legacy_stitched.mask_fractions is not None:
516 compare_image_to_legacy(
517 tc, cell_coadd.mask_fractions["rejected"], legacy_stitched.mask_fractions, expect_view=False
518 )
519 for n in range(legacy_stitched.n_noise_realizations):
520 compare_image_to_legacy(
521 tc, cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False
522 )
523 tc.assertEqual(cell_coadd.skymap, legacy_stitched.identifiers.skymap)
524 tc.assertEqual(cell_coadd.tract, legacy_stitched.identifiers.tract)
525 tc.assertEqual(cell_coadd.patch.index.x, legacy_stitched.identifiers.patch.x)
526 tc.assertEqual(cell_coadd.patch.index.y, legacy_stitched.identifiers.patch.y)
527 tc.assertEqual(cell_coadd.band, legacy_stitched.identifiers.band)
528 tc.assertTrue(tract_bbox.contains(cell_coadd.patch.outer_bbox))
529 tc.assertTrue(cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox))
530 tc.assertTrue(cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox))
531 tc.assertEqual(cell_coadd.unit, u.Unit(legacy_cell_coadd.common.units.value))
532 tc.assertTrue(cell_coadd.bounds.bbox.contains(cell_coadd.bbox))
533 tc.assertTrue(cell_coadd.grid.bbox.contains(cell_coadd.bbox))
534 compare_sky_projection_to_legacy_wcs(
535 tc,
536 cell_coadd.sky_projection,
537 legacy_cell_coadd.wcs,
538 TractFrame(
539 skymap=legacy_cell_coadd.identifiers.skymap,
540 tract=legacy_cell_coadd.identifiers.tract,
541 bbox=tract_bbox,
542 ),
543 cell_coadd.bbox,
544 is_fits=True,
545 )
546 tc.assertIs(cell_coadd.sky_projection, cell_coadd.mask.sky_projection)
547 tc.assertIs(cell_coadd.sky_projection, cell_coadd.variance.sky_projection)
548 compare_psf_to_legacy(
549 tc, cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points
550 )
551 compare_aperture_corrections_to_legacy(
552 tc, cell_coadd.aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
553 )
554 compare_cell_coadd_provenance_to_legacy(tc, cell_coadd.provenance, legacy_cell_coadd)
555 if alternates:
556 if sky_projection := alternates.get("sky_projection"):
557 compare_sky_projection_to_legacy_wcs(
558 tc,
559 sky_projection,
560 legacy_stitched.wcs,
561 TractFrame(
562 skymap=legacy_cell_coadd.identifiers.skymap,
563 tract=legacy_cell_coadd.identifiers.tract,
564 bbox=tract_bbox,
565 ),
566 cell_coadd.bbox,
567 is_fits=True,
568 )
569 if psf := alternates.get("psf"):
570 compare_psf_to_legacy(tc, psf, legacy_stitched.psf, points=psf_points)
571 if aperture_corrections := alternates.get("aperture_corrections"):
572 compare_aperture_corrections_to_legacy(
573 tc, aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox
574 )
575 if provenance := alternates.get("provenance"):
576 compare_cell_coadd_provenance_to_legacy(tc, provenance, legacy_cell_coadd)
579def compare_cell_coadd_provenance_to_legacy(
580 tc: unittest.TestCase, provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd
581) -> None:
582 """Compare a `.cells.CoaddProvenance` object to a legacy
583 `lsst.cell_coadds.MultipleCellCoadd` object.
585 Parameters
586 ----------
587 tc
588 Test case to use for asserts.
589 provenance
590 New provenance object to test.
591 legacy_cell_coadd
592 Legacy coadd to test against.
593 """
594 from lsst.cell_coadds import ObservationIdentifiers
596 for legacy_cell in legacy_cell_coadd.cells.values():
597 cell_index = CellIJ.from_legacy(legacy_cell.identifiers.cell)
598 prov = provenance[cell_index]
599 legacy_table = astropy.table.Table(
600 rows=[
601 [
602 ids.instrument,
603 ids.visit,
604 ids.detector,
605 ids.day_obs,
606 ids.physical_filter,
607 legacy_input.overlaps_center,
608 legacy_input.overlap_fraction,
609 legacy_input.weight,
610 legacy_input.psf_shape.getIxx(),
611 legacy_input.psf_shape.getIyy(),
612 legacy_input.psf_shape.getIxy(),
613 legacy_input.psf_shape_flag,
614 ]
615 for ids, legacy_input in legacy_cell.inputs.items()
616 ],
617 dtype=[
618 np.object_,
619 np.uint64,
620 np.uint16,
621 np.uint32,
622 np.object_,
623 np.bool_,
624 np.float64,
625 np.float64,
626 np.float64,
627 np.float64,
628 np.float64,
629 np.bool_,
630 ],
631 names=[
632 "instrument",
633 "visit",
634 "detector",
635 "day_obs",
636 "physical_filter",
637 "overlaps_center",
638 "overlap_fraction",
639 "weight",
640 "psf_shape_xx",
641 "psf_shape_yy",
642 "psf_shape_xy",
643 "psf_shape_flag",
644 ],
645 )
646 # For a single cell all 'inputs' are also 'contributions'.
647 tc.assertEqual(len(legacy_cell.inputs), len(prov.inputs))
648 tc.assertEqual(len(legacy_cell.inputs), len(prov.contributions))
649 prov.inputs.sort(["instrument", "visit", "detector"])
650 prov.contributions.sort(["instrument", "visit", "detector"])
651 legacy_table.sort(["instrument", "visit", "detector"])
652 np.testing.assert_array_equal(prov.inputs["instrument"], prov.contributions["instrument"])
653 np.testing.assert_array_equal(prov.inputs["visit"], prov.contributions["visit"])
654 np.testing.assert_array_equal(prov.inputs["detector"], prov.contributions["detector"])
655 np.testing.assert_array_equal(prov.inputs["instrument"], legacy_table["instrument"])
656 np.testing.assert_array_equal(prov.inputs["visit"], legacy_table["visit"])
657 np.testing.assert_array_equal(prov.inputs["detector"], legacy_table["detector"])
658 np.testing.assert_array_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"])
659 np.testing.assert_array_equal(prov.inputs["day_obs"], legacy_table["day_obs"])
660 np.testing.assert_array_equal(prov.contributions["overlaps_center"], legacy_table["overlaps_center"])
661 np.testing.assert_array_equal(
662 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"]
663 )
664 np.testing.assert_array_equal(prov.contributions["weight"], legacy_table["weight"])
665 np.testing.assert_array_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"])
666 np.testing.assert_array_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"])
667 np.testing.assert_array_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"])
668 np.testing.assert_array_equal(prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"])
669 for row in prov.inputs:
670 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"})
671 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key]
672 tc.assertEqual(legacy_polygon, row["polygon"].to_legacy())
675def compare_psf_to_legacy(
676 tc: unittest.TestCase,
677 psf: PointSpreadFunction,
678 legacy_psf: Any,
679 points: YX[np.ndarray] | XY[np.ndarray] | None = None,
680 expect_legacy_raise_on_out_of_bounds: bool = False,
681) -> int:
682 """Compare a PSF model object to its legacy interface.
684 Parameters
685 ----------
686 tc
687 Test case object with assert methods to use.
688 psf
689 Point-spread function to test.
690 legacy_psf
691 Legacy `lsst.afw.detection.Psf` instance to compare with.
692 points
693 Points to evaluate the PSFs at. If not provided, the intersection of
694 the PSF bounding boxes are used to generate points on a grid.
695 expect_legacy_raise_on_out_of_bounds
696 If `True`, expect ``legacy_psf`` to raise
697 `lsst.afw.detection.InvalidPsfError` when evaluated at a position
698 considered out-of-bounds by ``psf``.
700 Returns
701 -------
702 `int`
703 The number of points actually tested.
704 """
705 from lsst.afw.detection import InvalidPsfError
707 if points is None:
708 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel)
709 legacy_points = arrays_to_legacy_points(points.x, points.y)
710 n_points_tested: int = 0
711 for p in legacy_points:
712 if not psf.bounds.contains(x=p.x, y=p.y):
713 if expect_legacy_raise_on_out_of_bounds:
714 with tc.assertRaises(InvalidPsfError):
715 legacy_psf.computeKernelImage(p)
716 continue
717 tc.assertEqual(psf.kernel_bbox, Box.from_legacy(legacy_psf.computeKernelBBox(p)))
718 tc.assertEqual(
719 psf.compute_kernel_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeKernelImage(p))
720 )
721 tc.assertEqual(
722 psf.compute_stellar_bbox(x=p.x, y=p.y), Box.from_legacy(legacy_psf.computeImageBBox(p))
723 )
724 tc.assertEqual(psf.compute_stellar_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeImage(p)))
725 n_points_tested += 1
726 return n_points_tested
729def compare_field_to_legacy(
730 tc: unittest.TestCase,
731 field: BaseField,
732 legacy_field: Any,
733 subimage_bbox: Box,
734) -> None:
735 """Test a Field object by comparing it to an equivalent
736 `lsst.afw.math.BoundedField`.
738 Parameters
739 ----------
740 tc
741 Test case object with assert methods to use.
742 field
743 Field to test.
744 legacy_field : ``lsst.afw.math.BoundedField``
745 Equivalent legacy bounded field.
746 subimage_bbox
747 Bounding box for full-image tests.
748 """
749 from lsst.afw.math import BoundedField as LegacyBoundedField
751 tc.assertEqual(field.bounds.bbox, Box.from_legacy(legacy_field.getBBox()))
752 # Pixel coordinates to test the numpy array interface with.
753 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel)
754 if not isinstance(field.bounds, Box): 754 ↛ 755line 754 didn't jump to line 755 because the condition on line 754 was never true
755 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y)
756 pixel_xy = pixel_xy.map(lambda v: v[mask])
757 try:
758 assert_close(
759 tc,
760 field(x=pixel_xy.x, y=pixel_xy.y),
761 legacy_field.evaluate(pixel_xy.x, pixel_xy.y),
762 equal_nan=True,
763 )
764 except AssertionError as err:
765 err.add_note(f"evaluated at {pixel_xy}")
766 raise
767 if not isinstance(legacy_field, LegacyBoundedField): 767 ↛ 770line 767 didn't jump to line 770 because the condition on line 767 was never true
768 # Legacy StitchedApertureCorrection objects are not true BoundedFields
769 # and don't have addToImage.
770 return
771 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy()
772 legacy_field.addToImage(legacy_image_1, overlapOnly=True)
773 assert_images_equal(
774 tc, field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13
775 )
778def compare_aperture_corrections_to_legacy(
779 tc: unittest.TestCase,
780 aperture_corrections: Mapping[str, BaseField],
781 legacy_ap_corr_map: Any,
782 subimage_bbox: Box,
783) -> None:
784 """Test an aperture correction `dict` by comparing it to an equivalent
785 `lsst.afw.image.ApCorrMap`.
787 Parameters
788 ----------
789 tc
790 Test case object with assert methods to use.
791 aperture_corrections
792 Dictionary to test.
793 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap``
794 Equivalent legacy aperture correction map.
795 subimage_bbox
796 Bounding box for full-image tests.
797 """
798 tc.assertEqual(aperture_corrections.keys(), set(legacy_ap_corr_map.keys()))
799 for name, field in aperture_corrections.items():
800 compare_field_to_legacy(tc, field, legacy_ap_corr_map[name], subimage_bbox)
803def compare_observation_summary_stats_to_legacy(
804 tc: unittest.TestCase,
805 summary_stats: ObservationSummaryStats,
806 legacy_summary_stats: Any,
807) -> None:
808 """Test an ObservationSummaryStats object by comparing it to an equivalent
809 `lsst.afw.image.ExposureSummaryStats`.
811 Parameters
812 ----------
813 tc
814 Test case object with assert methods to use.
815 summary_stats
816 Struct to test.
817 legacy : ``lsst.afw.image.ExposureSummaryStats``
818 Equivalent legacy struct.
819 """
820 for field in dataclasses.fields(legacy_summary_stats):
821 a = getattr(legacy_summary_stats, field.name)
822 b = getattr(summary_stats, field.name)
823 if isinstance(b, tuple):
824 for ai, bi in zip(a, b):
825 tc.assertTrue(ai == bi or (math.isnan(ai) and math.isnan(bi)), f"{field.name}: {a} != {b}")
826 else:
827 tc.assertTrue(a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}")
830def compare_sky_projection_to_legacy_wcs[F: Frame](
831 tc: unittest.TestCase,
832 sky_projection: SkyProjection[F],
833 legacy_wcs: Any,
834 pixel_frame: F,
835 subimage_bbox: Box,
836 is_fits: bool = False,
837) -> None:
838 """Test a Projection object by comparing it to an equivalent
839 `lsst.afw.geom.SkyWcs`.
841 Parameters
842 ----------
843 tc
844 Test case object with assert methods to use.
845 sky_projection
846 Projection to test.
847 legacy_wcs : ``lsst.afw.geom.SkyWcs``
848 Equivalent legacy WCS.
849 pixel_frame
850 Expected pixel frame for the sky_projection.
851 subimage_bbox
852 Bounding box of points to generate for tests.
853 is_fits
854 Whether this sky_projection is expected to be exactly representable as
855 a FITS WCS. If `False` it is assumed to have a FITS approximation
856 attached instead.
857 """
858 # Pixel coordinates to test on over the subimage region of interest:
859 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel)
860 # Array indices of those pixel values (subtract off bbox starts):
861 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start)
862 sky_coords = legacy_coords_to_astropy(
863 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y))
864 )
865 # Test transforming with the Projection itself, which also tests its
866 # nested Transform and an Astropy High-Level WCS view with no origin
867 # change.
868 check_projection(tc, sky_projection, pixel_xy, sky_coords, pixel_frame)
869 # Also test the Astropy High-Level WCS view with an origin change to
870 # array indices.
871 check_astropy_wcs_interface(
872 tc, sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5
873 )
874 if is_fits:
875 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True)
876 assert fits_wcs is not None
877 check_astropy_wcs_interface(tc, fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5)
878 # Use that FITS approximation to check that we can make a
879 # Projection from a FITS WCS, too.
880 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame)
881 check_projection(
882 tc,
883 fits_projection,
884 subimage_array_xy,
885 sky_coords,
886 pixel_frame,
887 pixel_atol=1e-5,
888 )
889 # We want Projections we create from a FITS WCS to be backed by an
890 # AST FrameSet so we can convert them into legacy
891 # `lsst.afw.geom.SkyWcs` objects if desired.
892 tc.assertIn("Begin FrameSet", fits_projection.show())
893 else:
894 tc.assertIsNone(sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False))
895 # The legacy SkyWcs should instead have a FITS approximation
896 # attached; run the same tests on that.
897 assert sky_projection.fits_approximation is not None
898 compare_sky_projection_to_legacy_wcs(
899 tc,
900 sky_projection.fits_approximation,
901 legacy_wcs.getFitsApproximation(),
902 pixel_frame,
903 subimage_bbox,
904 is_fits=True,
905 )
908def check_transform[I: Frame, O: Frame](
909 tc: unittest.TestCase,
910 transform: Transform[I, O],
911 input_xy: XY[np.ndarray],
912 output_xy: XY[np.ndarray],
913 in_frame: Frame,
914 out_frame: Frame,
915 *,
916 check_inverted: bool = True,
917 in_atol: u.Quantity | None = None,
918 out_atol: u.Quantity | None = None,
919) -> None:
920 """Test Transform against known arrays of input and output points.
922 Parameters
923 ----------
924 tc
925 Test case object with assert methods to use.
926 transform
927 Transform to test.
928 input_xy
929 Arrays of input points.
930 output_xy
931 Arrays of output points.
932 in_frame
933 Expected input frame.
934 out_frame
935 Expected output frame.
936 check_inverted
937 If `True`, recurse (once) to test the inverse transform.
938 in_atol
939 Expected absolute precision of input points.
940 out_atol
941 Expected absolute precision of output points.
942 """
943 tc.assertEqual(transform.in_frame, in_frame)
944 tc.assertEqual(transform.out_frame, out_frame)
945 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None
946 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None
947 # Test array interfaces.
948 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y)
949 assert_close(tc, test_output_xy.x, output_xy.x, atol=out_atol_v)
950 assert_close(tc, test_output_xy.y, output_xy.y, atol=out_atol_v)
951 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y)
952 assert_close(tc, test_input_xy.x, input_xy.x, atol=in_atol_v)
953 assert_close(tc, test_input_xy.y, input_xy.y, atol=in_atol_v)
954 # Test scalar interfaces with numpy scalars.
955 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y):
956 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v)
957 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v)
958 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v)
959 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v)
960 # Test quantity array interfaces.
961 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit)
962 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit)
963 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y)
964 assert_close(tc, test_output_q_xy.x, output_q_xy.x, atol=out_atol)
965 assert_close(tc, test_output_q_xy.y, output_q_xy.y, atol=out_atol)
966 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y)
967 assert_close(tc, test_input_q_xy.x, input_q_xy.x, atol=in_atol)
968 assert_close(tc, test_input_q_xy.y, input_q_xy.y, atol=in_atol)
969 # Test quantity scalar interfaces.
970 for input_q_x, input_q_y, output_q_x, output_q_y in zip(
971 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y
972 ):
973 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol)
974 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol)
975 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol)
976 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol)
977 if check_inverted:
978 # Test the inverse transform.
979 check_transform(
980 tc,
981 transform.inverted(),
982 output_xy,
983 input_xy,
984 out_frame,
985 in_frame,
986 check_inverted=False,
987 out_atol=in_atol,
988 in_atol=out_atol,
989 )
992def check_projection[P: Frame](
993 tc: unittest.TestCase,
994 sky_projection: SkyProjection[P],
995 pixel_xy: XY[np.ndarray],
996 sky_coords: SkyCoord,
997 pixel_frame: Frame,
998 *,
999 pixel_atol: float | None = None,
1000 sky_atol: u.Quantity | None = None,
1001) -> None:
1002 """Test a `.SkyProjection` instance against known arrays of pixel and sky
1003 coordinates.
1005 Parameters
1006 ----------
1007 tc
1008 Test case object with assert methods to use.
1009 sky_projection
1010 Projection to test.
1011 pixel_xy
1012 Arrays of pixel coordinates.
1013 sky_coords
1014 Corresponding sky coordinates.
1015 pixel_frame
1016 Expected pixel frame.
1017 pixel_atol
1018 Expected absolute precision of pixel points.
1019 sky_atol
1020 Expected absolute precision of sky coordinates.
1021 """
1022 tc.assertEqual(sky_projection.pixel_frame, pixel_frame)
1023 tc.assertEqual(sky_projection.sky_frame, SkyFrame.ICRS)
1024 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None
1025 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None
1026 # Test array interfaces.
1027 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords))
1028 assert_close(tc, test_pixel_xy.x, pixel_xy.x, atol=pixel_atol)
1029 assert_close(tc, test_pixel_xy.y, pixel_xy.y, atol=pixel_atol)
1030 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y)
1031 assert_close(tc, test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v)
1032 assert_close(tc, test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v)
1033 # Test scalar interfaces.
1034 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords):
1035 assert_close(tc, sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol)
1036 assert_close(tc, sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol)
1037 assert_close(tc, sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v)
1038 assert_close(
1039 tc, sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v
1040 )
1041 # Test the underlying Transform object.
1042 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad))
1043 check_transform(
1044 tc,
1045 sky_projection.pixel_to_sky_transform,
1046 pixel_xy,
1047 sky_xy,
1048 pixel_frame,
1049 SkyFrame.ICRS,
1050 check_inverted=False,
1051 in_atol=pixel_atol_q,
1052 out_atol=sky_atol,
1053 )
1054 check_transform(
1055 tc,
1056 sky_projection.sky_to_pixel_transform,
1057 sky_xy,
1058 pixel_xy,
1059 SkyFrame.ICRS,
1060 pixel_frame,
1061 check_inverted=False,
1062 in_atol=sky_atol,
1063 out_atol=pixel_atol_q,
1064 )
1065 # Test the Astropy interface adapter.
1066 check_astropy_wcs_interface(
1067 tc, sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol
1068 )
1071def assert_sky_projections_equal(
1072 tc: unittest.TestCase,
1073 a: SkyProjection[Any] | None,
1074 b: SkyProjection[Any] | None,
1075 expect_identity: bool | None = None,
1076) -> None:
1077 """Test that two `.SkyProjection` instances are equivalent."""
1078 if a is None and b is None:
1079 return
1080 assert a is not None and b is not None
1081 match expect_identity:
1082 case True:
1083 tc.assertIs(a, b)
1084 return
1085 case False:
1086 tc.assertIsNot(a, b)
1087 case None if a is b:
1088 return
1089 tc.assertEqual(a, b)
1092def check_astropy_wcs_interface(
1093 tc: unittest.TestCase,
1094 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS,
1095 pixel_xy: XY[np.ndarray],
1096 sky_coords: SkyCoord,
1097 *,
1098 pixel_atol: float | None = None,
1099 sky_atol: u.Quantity | None = None,
1100) -> None:
1101 """Test an Astropy WCS instance against known arrays of pixel and
1102 sky coordinates.
1104 Parameters
1105 ----------
1106 tc
1107 Test case object with assert methods to use.
1108 wcs
1109 Astropy WCS object to test.
1110 pixel_xy
1111 Arrays of pixel coordinates.
1112 sky_coords
1113 Corresponding sky coordinates.
1114 pixel_atol
1115 Expected absolute precision of pixel points.
1116 sky_atol
1117 Expected absolute precision of sky coordinates.
1118 """
1119 test_x, test_y = wcs.world_to_pixel(sky_coords)
1120 assert_close(tc, test_x, pixel_xy.x, atol=pixel_atol)
1121 assert_close(tc, test_y, pixel_xy.y, atol=pixel_atol)
1122 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y)
1123 assert_close(tc, test_sky_coords.ra, sky_coords.ra, atol=sky_atol)
1124 assert_close(tc, test_sky_coords.dec, sky_coords.dec, atol=sky_atol)
1127def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]:
1128 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array."""
1129 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points]))
1132def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord:
1133 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy
1134 coordinate object.
1135 """
1136 return SkyCoord(
1137 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad,
1138 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad,
1139 )
1142def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]:
1143 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``."""
1144 from lsst.geom import Point2D
1146 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)]
1149def compare_amplifier_to_legacy(
1150 tc: unittest.TestCase,
1151 amplifier: Amplifier,
1152 legacy_amplifier: Any,
1153 *,
1154 is_raw_assembled: bool,
1155 expect_nominal_calibrations: bool = True,
1156) -> None:
1157 """Compare an `~.cameras.Amplifier` to a legacy
1158 `lsst.afw.cameraGeom.Amplifier`.
1159 """
1160 tc.assertEqual(legacy_amplifier.getName(), amplifier.name)
1161 tc.assertEqual(Box.from_legacy(legacy_amplifier.getBBox()), amplifier.bbox)
1162 if is_raw_assembled:
1163 raw_geom = amplifier.assembled_raw_geometry
1164 else:
1165 raw_geom = amplifier.unassembled_raw_geometry
1166 assert raw_geom is not None
1167 tc.assertEqual(ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()), raw_geom.readout_corner)
1168 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawBBox()), raw_geom.bbox)
1169 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawDataBBox()), raw_geom.data_bbox)
1170 tc.assertEqual(legacy_amplifier.getRawFlipX(), raw_geom.flip_x)
1171 tc.assertEqual(legacy_amplifier.getRawFlipY(), raw_geom.flip_y)
1172 tc.assertEqual(legacy_amplifier.getRawXYOffset().getX(), raw_geom.x_offset)
1173 tc.assertEqual(legacy_amplifier.getRawXYOffset().getY(), raw_geom.y_offset)
1174 tc.assertEqual(
1175 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()), raw_geom.horizontal_overscan_bbox
1176 )
1177 tc.assertEqual(
1178 Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()), raw_geom.vertical_overscan_bbox
1179 )
1180 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawPrescanBBox()), raw_geom.horizontal_prescan_bbox)
1181 if expect_nominal_calibrations:
1182 assert amplifier.nominal_calibrations is not None
1183 assert_equal_allow_nan(tc, legacy_amplifier.getGain(), amplifier.nominal_calibrations.gain)
1184 assert_equal_allow_nan(tc, legacy_amplifier.getReadNoise(), amplifier.nominal_calibrations.read_noise)
1185 assert_equal_allow_nan(
1186 tc, legacy_amplifier.getSaturation(), amplifier.nominal_calibrations.saturation
1187 )
1188 assert_equal_allow_nan(
1189 tc, legacy_amplifier.getSuspectLevel(), amplifier.nominal_calibrations.suspect_level
1190 )
1191 np.testing.assert_array_equal(
1192 legacy_amplifier.getLinearityCoeffs(), amplifier.nominal_calibrations.linearity_coefficients
1193 )
1194 tc.assertEqual(legacy_amplifier.getLinearityType(), amplifier.nominal_calibrations.linearity_type)
1197def compare_detector_to_legacy(
1198 tc: unittest.TestCase,
1199 detector: Detector,
1200 legacy_detector: Any,
1201 *,
1202 is_raw_assembled: bool,
1203 expect_nominal_calibrations: bool = True,
1204) -> None:
1205 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`."""
1206 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS
1208 tc.assertEqual(legacy_detector.getName(), detector.name)
1209 tc.assertEqual(legacy_detector.getId(), detector.id)
1210 tc.assertEqual(DetectorType.from_legacy(legacy_detector.getType()), detector.type)
1211 tc.assertEqual(Box.from_legacy(legacy_detector.getBBox()), detector.bbox)
1212 tc.assertEqual(legacy_detector.getSerial(), detector.serial)
1213 legacy_orientation = legacy_detector.getOrientation()
1214 tc.assertEqual(legacy_orientation.getFpPosition3().getX(), detector.orientation.focal_plane_x)
1215 tc.assertEqual(legacy_orientation.getFpPosition3().getY(), detector.orientation.focal_plane_y)
1216 tc.assertEqual(legacy_orientation.getFpPosition3().getZ(), detector.orientation.focal_plane_z)
1217 tc.assertEqual(legacy_orientation.getReferencePoint().getX(), detector.orientation.pixel_reference_x)
1218 tc.assertEqual(legacy_orientation.getReferencePoint().getY(), detector.orientation.pixel_reference_y)
1219 tc.assertEqual(legacy_orientation.getYaw().asRadians(), detector.orientation.yaw.to_value(u.rad))
1220 tc.assertEqual(legacy_orientation.getPitch().asRadians(), detector.orientation.pitch.to_value(u.rad))
1221 tc.assertEqual(legacy_orientation.getRoll().asRadians(), detector.orientation.roll.to_value(u.rad))
1222 tc.assertEqual(legacy_detector.getPixelSize().getX(), detector.pixel_size)
1223 tc.assertEqual(legacy_detector.getPhysicalType(), detector.physical_type)
1224 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True):
1225 compare_amplifier_to_legacy(
1226 tc,
1227 amplifier,
1228 legacy_amplifier,
1229 is_raw_assembled=is_raw_assembled,
1230 expect_nominal_calibrations=expect_nominal_calibrations,
1231 )
1232 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64))
1233 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x)
1234 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE)
1235 check_transform(
1236 tc,
1237 detector.to_focal_plane,
1238 pixel_xy,
1239 legacy_points_to_xy_array(fp_legacy_points),
1240 detector.frame,
1241 detector.to_focal_plane.out_frame,
1242 in_atol=1e-9 * u.pix,
1243 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit,
1244 )
1245 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE)
1246 check_transform(
1247 tc,
1248 detector.to_field_angle,
1249 pixel_xy,
1250 legacy_points_to_xy_array(fa_legacy_points),
1251 detector.frame,
1252 detector.to_field_angle.out_frame,
1253 in_atol=1e-9 * u.pix,
1254 out_atol=1e-7 * u.arcsec,
1255 )
1258def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]:
1259 """Yield every importable concrete `.serialization.ArchiveTree` subclass.
1261 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract
1262 classes. Importing this module already imports every ``lsst.images``
1263 module that defines a subclass, so the tree is fully populated by the time
1264 this is called.
1266 This discovery is deliberately separate from
1267 `check_archive_tree_class_invariants` so that the per-class check stays
1268 usable on a single class even if this metaprogramming is removed later.
1269 """
1270 seen: set[type] = set()
1271 stack: list[type] = [ArchiveTree]
1272 while stack:
1273 kls = stack.pop()
1274 for sub in kls.__subclasses__():
1275 if sub in seen: 1275 ↛ 1276line 1275 didn't jump to line 1276 because the condition on line 1275 was never true
1276 continue
1277 seen.add(sub)
1278 stack.append(sub)
1279 if not getattr(sub, "__abstractmethods__", None): 1279 ↛ 1274line 1279 didn't jump to line 1274 because the condition on line 1279 was always true
1280 yield sub
1283def check_archive_tree_class_invariants(tc: unittest.TestCase, cls: type[ArchiveTree]) -> None:
1284 """Assert that one concrete `.serialization.ArchiveTree` subclass declares
1285 well-formed schema-version constants and an in-memory type.
1287 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and
1288 ``PUBLIC_TYPE`` are present and well-typed, that the version is
1289 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the
1290 schema major.
1292 Parameters
1293 ----------
1294 tc
1295 Test case object with assert methods to use.
1296 cls
1297 The concrete `.serialization.ArchiveTree` subclass to check.
1298 """
1299 tc.assertTrue(hasattr(cls, "SCHEMA_NAME"), f"{cls.__name__} lacks SCHEMA_NAME")
1300 tc.assertTrue(hasattr(cls, "SCHEMA_VERSION"), f"{cls.__name__} lacks SCHEMA_VERSION")
1301 tc.assertTrue(hasattr(cls, "MIN_READ_VERSION"), f"{cls.__name__} lacks MIN_READ_VERSION")
1302 tc.assertTrue(hasattr(cls, "PUBLIC_TYPE"), f"{cls.__name__} lacks PUBLIC_TYPE")
1303 tc.assertIsInstance(cls.SCHEMA_NAME, str)
1304 tc.assertGreater(len(cls.SCHEMA_NAME), 0)
1305 tc.assertRegex(cls.SCHEMA_VERSION, r"^\d+\.\d+\.\d+$")
1306 tc.assertIsInstance(cls.MIN_READ_VERSION, int)
1307 tc.assertGreaterEqual(cls.MIN_READ_VERSION, 1)
1308 tc.assertIsInstance(cls.PUBLIC_TYPE, type)
1309 major = int(cls.SCHEMA_VERSION.split(".")[0])
1310 tc.assertLessEqual(cls.MIN_READ_VERSION, major)