Coverage for python/lsst/images/tests/_checks.py: 53%

565 statements  

« prev     ^ index     » next       coverage.py v7.16.0, created at 2026-09-22 03:11 -0700

1# This file is part of lsst-images. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# Use of this source code is governed by a 3-clause BSD-style 

10# license that can be found in the LICENSE file. 

11 

12from __future__ import annotations 

13 

14__all__ = ( 

15 "annotate_errors", 

16 "arrays_to_legacy_points", 

17 "assert_cell_coadds_equal", 

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_coords_close", 

24 "assert_sky_projections_equal", 

25 "assert_values_equal", 

26 "assert_visit_images_equal", 

27 "check_archive_tree_class_invariants", 

28 "check_astropy_wcs_interface", 

29 "check_bounds_contains_broadcasting", 

30 "check_projection", 

31 "check_transform", 

32 "compare_amplifier_to_legacy", 

33 "compare_aperture_corrections_to_legacy", 

34 "compare_cell_coadd_to_legacy", 

35 "compare_detector_to_legacy", 

36 "compare_field_to_legacy", 

37 "compare_image_to_legacy", 

38 "compare_mask_to_legacy", 

39 "compare_masked_image_to_legacy", 

40 "compare_observation_summary_stats_to_legacy", 

41 "compare_photo_calib_to_legacy", 

42 "compare_psf_to_legacy", 

43 "compare_sky_projection_to_legacy_wcs", 

44 "compare_visit_image_to_legacy", 

45 "iter_concrete_archive_tree_subclasses", 

46 "legacy_coords_to_astropy", 

47 "legacy_points_to_xy_array", 

48) 

49 

50import dataclasses 

51import math 

52import re 

53from collections.abc import Generator, Iterator, Mapping 

54from contextlib import contextmanager 

55from typing import TYPE_CHECKING, Any, Literal, cast 

56 

57import astropy.units as u 

58import astropy.wcs.wcsapi 

59import numpy as np 

60import pytest 

61from astropy.coordinates import SkyCoord 

62 

63from .._geom import XY, YX, Bounds, Box 

64from .._image import Image 

65from .._mask import Mask, MaskPlane, get_legacy_visit_image_mask_planes 

66from .._masked_image import MaskedImage 

67from .._observation_summary_stats import ObservationSummaryStats 

68from .._transforms import DetectorFrame, Frame, SkyFrame, SkyProjection, TractFrame, Transform 

69from .._visit_image import VisitImage 

70from ..cameras import Amplifier, Detector, DetectorType, ReadoutCorner 

71from ..cells import CellCoadd, CellIJ, CoaddProvenance 

72from ..fields import BaseField, ChebyshevField 

73from ..psfs import PointSpreadFunction 

74from ..serialization import ArchiveTree 

75 

76if TYPE_CHECKING: 

77 try: 

78 from lsst.cell_coadds import MultipleCellCoadd 

79 except ImportError: 

80 type MultipleCellCoadd = Any # type: ignore[no-redef] 

81 try: 

82 from lsst.afw.image import PhotoCalib as LegacyPhotoCalib 

83 except ImportError: 

84 type LegacyPhotoCalib = Any # type: ignore[no-redef] 

85 

86 

87# The sky-to-pixel direction of a SkyProjection is a numerical AST inversion 

88# of the celestial projection, so it carries a small error floor (~1e-9 px). 

89# A relative-only tolerance collapses to exact equality at pixel coordinate 

90# zero, so comparisons that can include the origin need a small absolute pixel 

91# tolerance to avoid reporting that numerical noise as a difference. 

92_AST_INV_ATOL = 1e-7 

93 

94# Relative tolerance for comparisons of AST-backed transforms and projections. 

95_AST_RTOL = 1e-9 

96 

97# Default absolute tolerance for great-circle sky-coordinate comparisons, in 

98# the same order of magnitude as _AST_INV_ATOL (which is in pixels). A 

99# nonzero default is important because sky coordinates can be near zero 

100# (e.g. RA/Dec ~ 0), where a relative tolerance provides no error floor. 

101_AST_SKY_ATOL = 1e-7 * u.deg 

102 

103 

104def assert_values_equal( 

105 a: np.ndarray | u.Quantity | float, 

106 b: np.ndarray | u.Quantity | float, 

107 *, 

108 rtol: float = 0.0, 

109 atol: float | u.Quantity = 0.0, 

110 equal_nan: bool = True, 

111 label: str = "", 

112) -> None: 

113 """Assert that two arrays, quantities, or floats are equal or close. 

114 

115 On mismatch this raises an `AssertionError` with a concise report (count, 

116 max absolute difference, location, NaN-mismatch count) rather than dumping 

117 the arrays. 

118 

119 Parameters 

120 ---------- 

121 a 

122 First value to compare. 

123 b 

124 Second value to compare. 

125 rtol 

126 Relative tolerance. 

127 atol 

128 Absolute tolerance; a `~astropy.units.Quantity` is converted to the 

129 unit of ``a``, or of the quantity-valued operand when only one 

130 operand has a unit. 

131 equal_nan 

132 If `True`, treat NaN as equal to NaN. 

133 label 

134 Prefix prepended to the failure message. 

135 """ 

136 if isinstance(a, u.Quantity) and isinstance(b, u.Quantity): 

137 unit = a.unit 

138 a_vals = np.asarray(a.value) 

139 b_vals = np.asarray(b.to_value(unit)) 

140 if isinstance(atol, u.Quantity): 

141 atol = atol.to_value(unit) 

142 elif isinstance(a, u.Quantity) or isinstance(b, u.Quantity): 142 ↛ 143line 142 didn't jump to line 143 because the condition on line 142 was never true

143 unit = getattr(a, "unit", None) or getattr(b, "unit", None) 

144 a_vals = np.asarray(getattr(a, "value", a)) 

145 b_vals = np.asarray(getattr(b, "value", b)) 

146 if isinstance(atol, u.Quantity): 

147 atol = atol.to_value(unit) 

148 else: 

149 unit = None 

150 a_vals = np.asarray(a) 

151 b_vals = np.asarray(b) 

152 if a_vals.shape != b_vals.shape: 

153 raise AssertionError(f"{label}: shape {a_vals.shape} != {b_vals.shape}") 

154 if rtol == 0.0 and atol == 0.0: 

155 mismatch = a_vals != b_vals 

156 if equal_nan and a_vals.dtype.kind in "fc": 

157 mismatch = np.logical_and( 

158 mismatch, 

159 np.logical_not(np.logical_and(np.isnan(a_vals), np.isnan(b_vals))), 

160 ) 

161 else: 

162 mismatch = ~np.isclose(a_vals, b_vals, rtol=rtol, atol=atol, equal_nan=equal_nan) 

163 n_mismatch = int(np.count_nonzero(mismatch)) 

164 if n_mismatch == 0: 

165 return 

166 prefix = f"{label}: " if label else "" 

167 if a_vals.shape == (): 

168 raise AssertionError(f"{prefix}{a!r} != {b!r}") 

169 if a_vals.dtype.kind in "fc": 

170 diff = np.abs(a_vals - b_vals) 

171 nan_mismatch = int( 

172 np.count_nonzero(np.logical_and(mismatch, np.logical_or(np.isnan(a_vals), np.isnan(b_vals)))) 

173 ) 

174 unit_str = f" {unit}" if unit is not None else "" 

175 nan_str = f"; {nan_mismatch} NaN mismatches" if nan_mismatch else "" 

176 if np.isnan(diff).all(): 

177 raise AssertionError(f"{prefix}{n_mismatch}/{a_vals.size} values differ{nan_str}{unit_str}") 

178 maxdiff = float(np.nanmax(diff)) 

179 # Convert to plain ints so the index prints as e.g. (2198,) rather than 

180 # (np.int64(2198),). 

181 loc = tuple(int(i) for i in np.unravel_index(np.nanargmax(diff), diff.shape)) 

182 raise AssertionError( 

183 f"{prefix}{n_mismatch}/{a_vals.size} values differ; " 

184 f"max abs diff {maxdiff} at index {loc}{nan_str}{unit_str}" 

185 ) 

186 else: 

187 first = np.argwhere(mismatch)[:5].tolist() 

188 raise AssertionError( 

189 f"{prefix}{n_mismatch}/{a_vals.size} values differ; first differing indices: {first}" 

190 ) 

191 

192 

193@contextmanager 

194def annotate_errors(note: str) -> Generator[None]: 

195 """Attach a note to any exception raised in the block. 

196 

197 Parameters 

198 ---------- 

199 note 

200 Note to attach on failure. 

201 """ 

202 try: 

203 yield 

204 except Exception as err: 

205 err.add_note(note) 

206 raise 

207 

208 

209def assert_equal_allow_nan(a: float, b: float) -> None: 

210 """Test that two floating point values are equal, with nan == nan. 

211 

212 Parameters 

213 ---------- 

214 a 

215 First value to compare. 

216 b 

217 Second value to compare. 

218 """ 

219 if not (a == b or (math.isnan(a) and math.isnan(b))): 

220 raise AssertionError(f"{a!r} != {b!r}") 

221 

222 

223def assert_sky_coords_close( 

224 test_sky: SkyCoord, expected_sky: SkyCoord, atol: u.Quantity, *, label: str = "sky_coords" 

225) -> None: 

226 """Assert that two astropy sky-coordinate sets are close. 

227 

228 Great-circle separations are used instead of per-axis RA/Dec 

229 comparisons, so the comparison is correct across the RA = 0 meridian 

230 and near the poles. 

231 

232 Parameters 

233 ---------- 

234 test_sky 

235 Sky coordinates to test, of the same shape as ``expected_sky``. 

236 expected_sky 

237 Expected sky coordinates. 

238 atol 

239 Maximum allowed great-circle separation. 

240 label 

241 Prefix prepended to the failure message. 

242 """ 

243 separation = expected_sky.separation(test_sky) 

244 if np.all(separation <= atol): 244 ↛ 246line 244 didn't jump to line 246 because the condition on line 244 was always true

245 return 

246 if separation.shape: 

247 worst = int(np.argmax(separation.value)) 

248 detail = f"separation {separation[worst]} > atol {atol} at index {worst}" 

249 else: 

250 detail = f"separation {separation} > atol {atol}" 

251 raise AssertionError(f"{label}: {detail}") 

252 

253 

254def _assert_wrapped_angles_close( 

255 a: np.ndarray | float, b: np.ndarray | float, *, atol: float, label: str = "" 

256) -> None: 

257 """Assert that angle arrays (in radians) are equal modulo 2 pi. 

258 

259 Parameters 

260 ---------- 

261 a 

262 First angle or angle array. 

263 b 

264 Second angle or angle array. 

265 atol 

266 Absolute tolerance, in radians. 

267 label 

268 Prefix prepended to the failure message. 

269 """ 

270 delta = (np.asarray(a) - np.asarray(b) + np.pi) % (2.0 * np.pi) - np.pi 

271 assert_values_equal(delta, np.zeros_like(delta), atol=atol, equal_nan=False, label=label) 

272 

273 

274def assert_images_equal( 

275 a: Image, 

276 b: Image, 

277 *, 

278 rtol: float = 0.0, 

279 atol: float = 0.0, 

280 expect_view: bool | Literal["array"] | None = None, 

281) -> None: 

282 """Assert that two images are equal or nearly equal. 

283 

284 Parameters 

285 ---------- 

286 a 

287 First image to compare. 

288 b 

289 Second image to compare. 

290 rtol 

291 Relative tolerance for the pixel comparison. 

292 atol 

293 Absolute tolerance for the pixel comparison. 

294 expect_view 

295 If not `None`, also assert whether ``b`` shares memory with ``a`` 

296 (i.e. is a view); ``"array"`` checks only the pixel arrays. 

297 """ 

298 assert a.bbox == b.bbox 

299 assert a.unit == b.unit 

300 assert_sky_projections_equal(a.sky_projection, b.sky_projection) 

301 if expect_view is not None: 

302 assert np.may_share_memory(a.array, b.array) == bool(expect_view) 

303 if expect_view == "array": 

304 assert a.metadata == b.metadata 

305 else: 

306 assert (a.metadata is b.metadata) == expect_view 

307 if not expect_view: 

308 assert_values_equal(a.array, b.array, atol=atol, rtol=rtol) 

309 assert a.metadata == b.metadata 

310 

311 

312def _note_mask_difference( 

313 err: AssertionError, a: Mask, b: Mask, *, label: str, names: Mapping[str, str] | None = None 

314) -> None: 

315 """Embed a plane-by-plane `Mask.compare` summary as a note on ``err``. 

316 

317 Parameters 

318 ---------- 

319 err 

320 Exception to attach the summary to. 

321 a 

322 First mask compared. 

323 b 

324 Second mask compared. 

325 label 

326 Prefix for the failure message. 

327 names 

328 Optional mapping from mask-plane names used by ``a`` and ``b`` to 

329 names to display in the summary (e.g. the legacy names, when 

330 comparing against legacy code). 

331 """ 

332 diff = a.compare(b) 

333 if not diff: 333 ↛ 334line 333 didn't jump to line 334 because the condition on line 333 was never true

334 return 

335 shown = names or {} 

336 summary = "\n".join( 

337 f"{label}[{shown.get(name, name)}]: +{added} -{removed}" for name, (added, removed) in diff.items() 

338 ) 

339 err.add_note(summary) 

340 

341 

342def assert_masks_equal(a: Mask, b: Mask) -> None: 

343 """Assert that two masks are equal or nearly equal. 

344 

345 Parameters 

346 ---------- 

347 a 

348 First mask to compare. 

349 b 

350 Second mask to compare. 

351 """ 

352 assert a.bbox == b.bbox 

353 assert a.schema == b.schema 

354 assert a.metadata == b.metadata 

355 assert_sky_projections_equal(a.sky_projection, b.sky_projection) 

356 try: 

357 assert_values_equal(a.array, b.array, label="mask") 

358 except AssertionError as err: 

359 _note_mask_difference(err, a, b, label="mask") 

360 raise 

361 

362 

363def assert_masked_images_equal( 

364 a: MaskedImage, 

365 b: MaskedImage, 

366 *, 

367 rtol: float = 0.0, 

368 atol: float = 0.0, 

369 expect_view: bool | None = None, 

370) -> None: 

371 """Assert that two masked images are equal or nearly equal. 

372 

373 Parameters 

374 ---------- 

375 a 

376 First masked image to compare. 

377 b 

378 Second masked image to compare. 

379 rtol 

380 Relative tolerance for the pixel comparison. 

381 atol 

382 Absolute tolerance for the pixel comparison. 

383 expect_view 

384 If not `None`, also assert whether ``b`` shares memory with ``a`` 

385 (i.e. is a view). 

386 """ 

387 assert a.metadata == b.metadata 

388 assert_sky_projections_equal(a.sky_projection, b.sky_projection) 

389 assert_images_equal(a.image, b.image, rtol=rtol, atol=atol, expect_view=expect_view) 

390 assert_masks_equal(a.mask, b.mask) 

391 assert_images_equal(a.variance, b.variance, rtol=rtol, atol=atol, expect_view=expect_view) 

392 

393 

394def assert_psfs_equal( 

395 psf1: PointSpreadFunction, 

396 psf2: PointSpreadFunction, 

397 points: YX[np.ndarray] | XY[np.ndarray] | None = None, 

398) -> int: 

399 """Compare two PSF objets. 

400 

401 Parameters 

402 ---------- 

403 psf1 

404 Point-spread function to test. 

405 psf2 

406 The other point-spread function to test. 

407 points 

408 Points to evaluate the PSFs at. If not provided, the intersection of 

409 the PSF bounding boxes are used to generate points on a grid. 

410 

411 Returns 

412 ------- 

413 `int` 

414 The number of points actually tested. 

415 """ 

416 if points is None: 416 ↛ 419line 416 didn't jump to line 419 because the condition on line 416 was always true

417 points = psf1.bounds.bbox.intersection(psf2.bounds.bbox).meshgrid(3).map(np.ravel) 

418 

419 assert psf1.kernel_bbox == psf2.kernel_bbox 

420 

421 n_points_tested: int = 0 

422 for x, y in zip(points.x, points.y): 

423 # The two PSFs must agree on which points fall inside their input 

424 # domain. Querying ``.contains`` directly (rather than relying on 

425 # ``compute_kernel_image`` to raise) makes this test tolerant of 

426 # implementations that do not raise on out-of-domain points -- in 

427 # particular ``CellPointSpreadFunction``, where evaluating in a 

428 # missing cell does not always raise ``BoundsError``. 

429 contains1 = psf1.bounds.contains(x=x, y=y) 

430 contains2 = psf2.bounds.contains(x=x, y=y) 

431 assert contains1 == contains2, ( 

432 f"PSFs disagree on whether ({x}, {y}) is in-bounds: psf1={contains1}, psf2={contains2}" 

433 ) 

434 if not contains1: 434 ↛ 435line 434 didn't jump to line 435 because the condition on line 434 was never true

435 continue 

436 assert psf1.compute_kernel_image(x=x, y=y) == psf2.compute_kernel_image(x=x, y=y) 

437 assert psf1.compute_stellar_bbox(x=x, y=y) == psf2.compute_stellar_bbox(x=x, y=y) 

438 assert psf1.compute_stellar_image(x=x, y=y) == psf2.compute_stellar_image(x=x, y=y) 

439 n_points_tested += 1 

440 return n_points_tested 

441 

442 

443def assert_visit_images_equal( 

444 a: VisitImage, 

445 b: VisitImage, 

446 *, 

447 expect_view: bool | None = None, 

448) -> None: 

449 """Assert that two `.VisitImage` instances carry the same persistent state. 

450 

451 Extends `assert_masked_images_equal` with the VisitImage-specific 

452 attributes (PSF, filter, observation info, detector, aperture 

453 corrections, photometric scaling, backgrounds, polygon bounds, 

454 summary stats) so a round-trip check on a `.VisitImage` does not 

455 silently miss differences in any of them. 

456 

457 Parameters 

458 ---------- 

459 a 

460 First visit image to compare. 

461 b 

462 Second visit image to compare. 

463 expect_view 

464 If not `None`, also assert whether ``b`` shares memory with ``a`` 

465 (i.e. is a view). 

466 """ 

467 assert_masked_images_equal(a, b, expect_view=expect_view) 

468 assert a.summary_stats == b.summary_stats 

469 assert a.physical_filter == b.physical_filter 

470 assert a.band == b.band 

471 assert a.obs_info == b.obs_info 

472 assert a.detector == b.detector 

473 assert dict(a.aperture_corrections) == dict(b.aperture_corrections) 

474 assert a.photometric_scaling == b.photometric_scaling 

475 assert dict(a.backgrounds) == dict(b.backgrounds) 

476 assert a.backgrounds.subtracted == b.backgrounds.subtracted 

477 assert a.bounds == b.bounds 

478 assert_psfs_equal(a.psf, b.psf) 

479 

480 

481def assert_cell_coadds_equal( 

482 a: CellCoadd, 

483 b: CellCoadd, 

484 *, 

485 expect_view: bool | None = None, 

486) -> None: 

487 """Assert that two `.CellCoadd` instances carry the same persistent state. 

488 

489 Extends the masked-image-style equality check with the 

490 CellCoadd-specific attributes (PSF, cell grid, missing cells, 

491 backgrounds, patch/tract, band) so a round-trip check does not 

492 silently miss differences in any of them. 

493 

494 Parameters 

495 ---------- 

496 a 

497 First cell coadd to compare. 

498 b 

499 Second cell coadd to compare. 

500 expect_view 

501 If not `None`, also assert whether ``b`` shares memory with ``a`` 

502 (i.e. is a view). 

503 """ 

504 assert_masked_images_equal(a, b, expect_view=expect_view) 

505 assert a.band == b.band 

506 assert a.patch == b.patch 

507 assert a.tract == b.tract 

508 assert a.grid == b.grid 

509 assert a.bounds.missing == b.bounds.missing 

510 assert dict(a.backgrounds) == dict(b.backgrounds) 

511 assert a.backgrounds.subtracted == b.backgrounds.subtracted 

512 assert_psfs_equal(a.psf, b.psf) 

513 

514 

515def compare_image_to_legacy( 

516 image: Image, legacy_image: Any, expect_view: bool | None = None, *, label: str = "image" 

517) -> None: 

518 """Compare an `.Image` object to a legacy `lsst.afw.image.Image` object. 

519 

520 Parameters 

521 ---------- 

522 image 

523 Image to compare. 

524 legacy_image 

525 Legacy `lsst.afw.image.Image` to compare against. 

526 expect_view 

527 If not `None`, also assert whether ``image`` shares memory with 

528 ``legacy_image`` (i.e. is a view). 

529 label 

530 Prefix for the failure message. 

531 """ 

532 assert image.bbox == Box.from_legacy(legacy_image.getBBox()) 

533 if expect_view is not None: 533 ↛ 535line 533 didn't jump to line 535 because the condition on line 533 was always true

534 assert np.may_share_memory(image.array, legacy_image.array) == expect_view 

535 if not expect_view: 535 ↛ exitline 535 didn't return from function 'compare_image_to_legacy' because the condition on line 535 was always true

536 assert_values_equal(image.array, legacy_image.array, label=label) 

537 

538 

539def compare_mask_to_legacy( 

540 mask: Mask, 

541 legacy_mask: Any, 

542 plane_map: Mapping[str, MaskPlane] | None = None, 

543 *, 

544 label: str = "mask", 

545) -> None: 

546 """Compare a `.Mask` object to a legacy `lsst.afw.image.Mask` object. 

547 

548 Parameters 

549 ---------- 

550 mask 

551 Mask to compare. 

552 legacy_mask 

553 Legacy `lsst.afw.image.Mask` to compare against. 

554 plane_map 

555 Mapping from legacy plane name to the new mask plane; defaults to 

556 the planes in ``mask.schema``. 

557 label 

558 Prefix for the failure message; each plane is reported as 

559 ``"{label}[<plane>]"``. 

560 """ 

561 assert mask.bbox == Box.from_legacy(legacy_mask.getBBox()) 

562 if plane_map is None: 562 ↛ 564line 562 didn't jump to line 564 because the condition on line 562 was always true

563 plane_map = {plane.name: plane for plane in mask.schema if plane is not None} 

564 try: 

565 for old_name, new_plane in plane_map.items(): 

566 assert_values_equal( 

567 (legacy_mask.array & legacy_mask.getPlaneBitMask(old_name)).astype(bool), 

568 mask.get(new_plane.name), 

569 label=f"{label}[{old_name}]", 

570 ) 

571 except AssertionError as err: 

572 names = {new_plane.name: old_name for old_name, new_plane in plane_map.items()} 

573 _note_mask_difference(err, mask, Mask.from_legacy(legacy_mask, plane_map), label=label, names=names) 

574 raise 

575 

576 

577def compare_masked_image_to_legacy( 

578 masked_image: MaskedImage, 

579 legacy_masked_image: Any, 

580 *, 

581 plane_map: Mapping[str, MaskPlane] | None = None, 

582 expect_view: bool | None = None, 

583 alternates: Mapping[str, Any] | None = None, 

584) -> None: 

585 """Compare a `.MaskedImage` object to a legacy `lsst.afw.image.MaskedImage` 

586 object. 

587 

588 Parameters 

589 ---------- 

590 masked_image 

591 New image to test. 

592 legacy_masked_image 

593 Legacy image to test against. 

594 plane_map 

595 Mapping between new and legacy mask planes. 

596 expect_view 

597 Whether to test that the image and variance arrays do or do not share 

598 memory. 

599 alternates 

600 A mapping of other versions of one or more (new) components to also 

601 check against the legacy versions of those components. 

602 """ 

603 compare_image_to_legacy( 

604 masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view, label="image" 

605 ) 

606 compare_mask_to_legacy( 

607 masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask" 

608 ) 

609 compare_image_to_legacy( 

610 masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view, label="variance" 

611 ) 

612 if alternates: 612 ↛ 613line 612 didn't jump to line 613 because the condition on line 612 was never true

613 if image := alternates.get("image"): 

614 compare_image_to_legacy( 

615 image, legacy_masked_image.getImage(), expect_view=expect_view, label="image" 

616 ) 

617 if mask := alternates.get("mask"): 

618 compare_mask_to_legacy(mask, legacy_masked_image.getMask(), plane_map=plane_map, label="mask") 

619 if variance := alternates.get("variance"): 

620 compare_image_to_legacy( 

621 variance, legacy_masked_image.getVariance(), expect_view=expect_view, label="variance" 

622 ) 

623 

624 

625def compare_visit_image_to_legacy( 

626 visit_image: VisitImage, 

627 legacy_exposure: Any, 

628 *, 

629 plane_map: Mapping[str, MaskPlane] | None = None, 

630 expect_view: bool | None = None, 

631 instrument: str, 

632 visit: int, 

633 detector: int, 

634 applied_legacy_photo_calib: LegacyPhotoCalib | None = None, 

635 alternates: Mapping[str, Any] | None = None, 

636 check_photometric_scaling: bool = True, 

637) -> None: 

638 """Compare a `.VisitImage` object to a legacy `lsst.afw.image.Exposure` 

639 object. 

640 

641 Parameters 

642 ---------- 

643 visit_image 

644 New image to test. 

645 legacy_exposure 

646 Legacy image to test against. 

647 plane_map 

648 Mapping between new and legacy mask planes; defaults to the planes 

649 used by the legacy (`lsst.afw.image.Exposure`) visit image. 

650 expect_view 

651 Whether to test that the image and variance arrays do or do not share 

652 memory. 

653 instrument 

654 Expected instrument name. 

655 visit 

656 Expected visit ID. 

657 detector 

658 Expected detector ID. 

659 applied_legacy_photo_calib 

660 Legacy `lsst.afw.image.PhotoCalib` already applied to 

661 ``legacy_exposure``, used when comparing photometric scaling. 

662 alternates 

663 A mapping of other versions of one or more (new) components to also 

664 check against the legacy versions of those components. 

665 check_photometric_scaling 

666 If `False`, skip the photometric-scaling comparison. 

667 """ 

668 if plane_map is None: 

669 plane_map = get_legacy_visit_image_mask_planes() 

670 compare_masked_image_to_legacy( 

671 visit_image, 

672 legacy_exposure.getMaskedImage(), 

673 plane_map=plane_map, 

674 expect_view=expect_view, 

675 alternates=alternates, 

676 ) 

677 detector_bbox = Box.from_legacy(legacy_exposure.getDetector().getBBox()) 

678 with annotate_errors("sky_projection"): 

679 compare_sky_projection_to_legacy_wcs( 

680 visit_image.sky_projection, 

681 legacy_exposure.getWcs(), 

682 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox), 

683 visit_image.bbox, 

684 ) 

685 assert visit_image.sky_projection is visit_image.mask.sky_projection 

686 assert visit_image.sky_projection is visit_image.variance.sky_projection 

687 with annotate_errors("psf"): 

688 compare_psf_to_legacy(visit_image.psf, legacy_exposure.getPsf()) 

689 with annotate_errors("summary_stats"): 

690 compare_observation_summary_stats_to_legacy( 

691 visit_image.summary_stats, legacy_exposure.info.getSummaryStats() 

692 ) 

693 with annotate_errors("detector"): 

694 compare_detector_to_legacy(visit_image.detector, legacy_exposure.getDetector(), is_raw_assembled=True) 

695 # Make a tiny box for Field comparisons that need to make arrays; that can 

696 # get expensive otherwise. 

697 tiny_bbox = detector_bbox.local[2:4, 3:6] 

698 with annotate_errors("aperture_corrections"): 

699 compare_aperture_corrections_to_legacy( 

700 visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox 

701 ) 

702 if check_photometric_scaling: 

703 with annotate_errors("photometric_scaling"): 

704 compare_photo_calib_to_legacy( 

705 visit_image.photometric_scaling, 

706 legacy_exposure.info.getPhotoCalib(), 

707 applied_legacy_photo_calib=applied_legacy_photo_calib, 

708 subimage_bbox=tiny_bbox, 

709 ) 

710 if alternates: 

711 if (bbox := alternates.get("bbox")) is not None: 

712 assert bbox == visit_image.bbox 

713 if sky_projection := alternates.get("sky_projection"): 

714 with annotate_errors("sky_projection"): 

715 compare_sky_projection_to_legacy_wcs( 

716 sky_projection, 

717 legacy_exposure.getWcs(), 

718 DetectorFrame(instrument=instrument, visit=visit, detector=detector, bbox=detector_bbox), 

719 visit_image.bbox, 

720 ) 

721 if psf := alternates.get("psf"): 

722 with annotate_errors("psf"): 

723 compare_psf_to_legacy(psf, legacy_exposure.getPsf()) 

724 if summary_stats := alternates.get("summary_stats"): 

725 with annotate_errors("summary_stats"): 

726 compare_observation_summary_stats_to_legacy( 

727 summary_stats, legacy_exposure.info.getSummaryStats() 

728 ) 

729 if detector_obj := alternates.get("detector"): 

730 with annotate_errors("detector"): 

731 compare_detector_to_legacy(detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True) 

732 if obs_info := alternates.get("obs_info"): 

733 visitInfo = legacy_exposure.visitInfo 

734 assert obs_info.instrument == visitInfo.getInstrumentLabel() 

735 if aperture_corrections := alternates.get("aperture_corrections"): 

736 with annotate_errors("aperture_corrections"): 

737 compare_aperture_corrections_to_legacy( 

738 aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox 

739 ) 

740 if check_photometric_scaling: 

741 if (photometric_scaling := alternates.get("photometric_scaling", ...)) is not ...: 

742 with annotate_errors("photometric_scaling"): 

743 compare_photo_calib_to_legacy( 

744 photometric_scaling, 

745 legacy_exposure.info.getPhotoCalib(), 

746 applied_legacy_photo_calib=applied_legacy_photo_calib, 

747 subimage_bbox=tiny_bbox, 

748 ) 

749 

750 

751def compare_photo_calib_to_legacy( 

752 photometric_scaling: BaseField | None, 

753 legacy_photo_calib: LegacyPhotoCalib, 

754 *, 

755 applied_legacy_photo_calib: LegacyPhotoCalib | None = None, 

756 subimage_bbox: Box, 

757) -> None: 

758 if legacy_photo_calib._isConstant: 

759 if legacy_photo_calib.getCalibrationMean() == 1.0: 

760 if applied_legacy_photo_calib is None: 

761 assert photometric_scaling is None 

762 return 

763 else: 

764 legacy_photo_calib = applied_legacy_photo_calib 

765 if legacy_photo_calib._isConstant: 

766 assert isinstance(photometric_scaling, ChebyshevField) 

767 assert_values_equal( 

768 photometric_scaling.coefficients, 

769 np.array([[legacy_photo_calib.getCalibrationMean()]]), 

770 # The coefficients come from a numerical fit, so this can be 

771 # bigger than just round-off error epislon. 

772 rtol=1e-5, 

773 ) 

774 else: 

775 assert photometric_scaling is not None 

776 compare_field_to_legacy( 

777 photometric_scaling / legacy_photo_calib.getCalibrationMean(), 

778 legacy_photo_calib.computeScaledCalibration(), 

779 subimage_bbox, 

780 ) 

781 

782 

783def compare_cell_coadd_to_legacy( 

784 cell_coadd: CellCoadd, 

785 legacy_cell_coadd: MultipleCellCoadd, 

786 *, 

787 tract_bbox: Box, 

788 plane_map: Mapping[str, MaskPlane] | None = None, 

789 alternates: Mapping[str, Any] | None = None, 

790 psf_points: XY[np.ndarray] | YX[np.ndarray] | None = None, 

791) -> None: 

792 """Compare a `.cells.CellCoadd` object to a legacy 

793 `lsst.cell_coadds.MultipleCellCoadd` object. 

794 

795 Parameters 

796 ---------- 

797 cell_coadd 

798 New coadd to test. 

799 legacy_cell_coadd 

800 Legacy coadd to test against. 

801 tract_bbox 

802 Bounding box of the full tract. 

803 plane_map 

804 Mapping between new and legacy mask planes. 

805 alternates 

806 A mapping of other versions of one or more (new) components to also 

807 check against the legacy versions of those components. 

808 psf_points 

809 Points to use to compare the PSFs. 

810 """ 

811 legacy_stitched = legacy_cell_coadd.stitch(cell_coadd.bbox.to_legacy()) 

812 compare_image_to_legacy(cell_coadd.image, legacy_stitched.image, expect_view=False) 

813 compare_mask_to_legacy(cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map) 

814 compare_image_to_legacy(cell_coadd.variance, legacy_stitched.variance, expect_view=False) 

815 if legacy_stitched.mask_fractions is not None: 

816 compare_image_to_legacy( 

817 cell_coadd.mask_fractions["rejected"], legacy_stitched.mask_fractions, expect_view=False 

818 ) 

819 for n in range(legacy_stitched.n_noise_realizations): 

820 compare_image_to_legacy( 

821 cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False 

822 ) 

823 assert cell_coadd.skymap == legacy_stitched.identifiers.skymap 

824 assert cell_coadd.tract == legacy_stitched.identifiers.tract 

825 assert cell_coadd.patch.index.x == legacy_stitched.identifiers.patch.x 

826 assert cell_coadd.patch.index.y == legacy_stitched.identifiers.patch.y 

827 assert cell_coadd.band == legacy_stitched.identifiers.band 

828 assert tract_bbox.contains(cell_coadd.patch.outer_bbox) 

829 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox) 

830 assert cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox) 

831 assert cell_coadd.unit == u.Unit(legacy_cell_coadd.common.units.value) 

832 assert cell_coadd.bounds.bbox.contains(cell_coadd.bbox) 

833 assert cell_coadd.grid.bbox.contains(cell_coadd.bbox) 

834 compare_sky_projection_to_legacy_wcs( 

835 cell_coadd.sky_projection, 

836 legacy_cell_coadd.wcs, 

837 TractFrame( 

838 skymap=legacy_cell_coadd.identifiers.skymap, 

839 tract=legacy_cell_coadd.identifiers.tract, 

840 bbox=tract_bbox, 

841 ), 

842 cell_coadd.bbox, 

843 is_fits=True, 

844 ) 

845 assert cell_coadd.sky_projection is cell_coadd.mask.sky_projection 

846 assert cell_coadd.sky_projection is cell_coadd.variance.sky_projection 

847 compare_psf_to_legacy( 

848 cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points 

849 ) 

850 compare_aperture_corrections_to_legacy( 

851 cell_coadd.aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox 

852 ) 

853 compare_cell_coadd_provenance_to_legacy(cell_coadd.provenance, legacy_cell_coadd) 

854 if alternates: 

855 if sky_projection := alternates.get("sky_projection"): 

856 compare_sky_projection_to_legacy_wcs( 

857 sky_projection, 

858 legacy_stitched.wcs, 

859 TractFrame( 

860 skymap=legacy_cell_coadd.identifiers.skymap, 

861 tract=legacy_cell_coadd.identifiers.tract, 

862 bbox=tract_bbox, 

863 ), 

864 cell_coadd.bbox, 

865 is_fits=True, 

866 ) 

867 if psf := alternates.get("psf"): 

868 compare_psf_to_legacy(psf, legacy_stitched.psf, points=psf_points) 

869 if aperture_corrections := alternates.get("aperture_corrections"): 

870 compare_aperture_corrections_to_legacy( 

871 aperture_corrections, legacy_stitched.ap_corr_map, cell_coadd.bbox 

872 ) 

873 if provenance := alternates.get("provenance"): 

874 compare_cell_coadd_provenance_to_legacy(provenance, legacy_cell_coadd) 

875 

876 

877def compare_cell_coadd_provenance_to_legacy( 

878 provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd 

879) -> None: 

880 """Compare a `.cells.CoaddProvenance` object to a legacy 

881 `lsst.cell_coadds.MultipleCellCoadd` object. 

882 

883 Parameters 

884 ---------- 

885 provenance 

886 New provenance object to test. 

887 legacy_cell_coadd 

888 Legacy coadd to test against. 

889 """ 

890 from lsst.cell_coadds import ObservationIdentifiers 

891 

892 for legacy_cell in legacy_cell_coadd.cells.values(): 

893 cell_index = CellIJ.from_legacy(legacy_cell.identifiers.cell) 

894 prov = provenance[cell_index] 

895 assert prov is not None, f"Legacy cell {cell_index} has inputs, so its provenance is not empty." 

896 legacy_table = astropy.table.Table( 

897 rows=[ 

898 [ 

899 ids.instrument, 

900 ids.visit, 

901 ids.detector, 

902 ids.day_obs, 

903 ids.physical_filter, 

904 legacy_input.overlaps_center, 

905 legacy_input.overlap_fraction, 

906 legacy_input.weight, 

907 legacy_input.psf_shape.getIxx(), 

908 legacy_input.psf_shape.getIyy(), 

909 legacy_input.psf_shape.getIxy(), 

910 legacy_input.psf_shape_flag, 

911 ] 

912 for ids, legacy_input in legacy_cell.inputs.items() 

913 ], 

914 dtype=[ 

915 np.object_, 

916 np.uint64, 

917 np.uint16, 

918 np.uint32, 

919 np.object_, 

920 np.bool_, 

921 np.float64, 

922 np.float64, 

923 np.float64, 

924 np.float64, 

925 np.float64, 

926 np.bool_, 

927 ], 

928 names=[ 

929 "instrument", 

930 "visit", 

931 "detector", 

932 "day_obs", 

933 "physical_filter", 

934 "overlaps_center", 

935 "overlap_fraction", 

936 "weight", 

937 "psf_shape_xx", 

938 "psf_shape_yy", 

939 "psf_shape_xy", 

940 "psf_shape_flag", 

941 ], 

942 ) 

943 # For a single cell all 'inputs' are also 'contributions'. 

944 assert len(legacy_cell.inputs) == len(prov.inputs) 

945 assert len(legacy_cell.inputs) == len(prov.contributions) 

946 prov.inputs.sort(["instrument", "visit", "detector"]) 

947 prov.contributions.sort(["instrument", "visit", "detector"]) 

948 legacy_table.sort(["instrument", "visit", "detector"]) 

949 assert_values_equal(prov.inputs["instrument"], prov.contributions["instrument"], label="inputs") 

950 assert_values_equal(prov.inputs["visit"], prov.contributions["visit"], label="inputs") 

951 assert_values_equal(prov.inputs["detector"], prov.contributions["detector"], label="inputs") 

952 assert_values_equal(prov.inputs["instrument"], legacy_table["instrument"], label="inputs") 

953 assert_values_equal(prov.inputs["visit"], legacy_table["visit"], label="inputs") 

954 assert_values_equal(prov.inputs["detector"], legacy_table["detector"], label="inputs") 

955 assert_values_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"], label="inputs") 

956 assert_values_equal(prov.inputs["day_obs"], legacy_table["day_obs"], label="inputs") 

957 assert_values_equal( 

958 prov.contributions["overlaps_center"], legacy_table["overlaps_center"], label="inputs" 

959 ) 

960 assert_values_equal( 

961 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"], label="inputs" 

962 ) 

963 assert_values_equal(prov.contributions["weight"], legacy_table["weight"], label="inputs") 

964 assert_values_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"], label="inputs") 

965 assert_values_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"], label="inputs") 

966 assert_values_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"], label="inputs") 

967 assert_values_equal( 

968 prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"], label="inputs" 

969 ) 

970 for row in prov.inputs: 

971 polygon_key = ObservationIdentifiers(**{k: row[k] for k in row.keys() if k != "polygon"}) 

972 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key] 

973 assert legacy_polygon == row["polygon"].to_legacy() 

974 

975 

976def compare_psf_to_legacy( 

977 psf: PointSpreadFunction, 

978 legacy_psf: Any, 

979 points: YX[np.ndarray] | XY[np.ndarray] | None = None, 

980 expect_legacy_raise_on_out_of_bounds: bool = False, 

981) -> int: 

982 """Compare a PSF model object to its legacy interface. 

983 

984 Parameters 

985 ---------- 

986 psf 

987 Point-spread function to test. 

988 legacy_psf 

989 Legacy `lsst.afw.detection.Psf` instance to compare with. 

990 points 

991 Points to evaluate the PSFs at. If not provided, the intersection of 

992 the PSF bounding boxes are used to generate points on a grid. 

993 expect_legacy_raise_on_out_of_bounds 

994 If `True`, expect ``legacy_psf`` to raise 

995 `lsst.afw.detection.InvalidPsfError` when evaluated at a position 

996 considered out-of-bounds by ``psf``. 

997 

998 Returns 

999 ------- 

1000 `int` 

1001 The number of points actually tested. 

1002 """ 

1003 from lsst.afw.detection import InvalidPsfError 

1004 

1005 if points is None: 

1006 points = psf.bounds.bbox.meshgrid(n=3).map(np.ravel) 

1007 legacy_points = arrays_to_legacy_points(points.x, points.y) 

1008 n_points_tested: int = 0 

1009 for p in legacy_points: 

1010 if not psf.bounds.contains(x=p.x, y=p.y): 

1011 if expect_legacy_raise_on_out_of_bounds: 

1012 with pytest.raises(InvalidPsfError): 

1013 legacy_psf.computeKernelImage(p) 

1014 continue 

1015 assert psf.kernel_bbox == Box.from_legacy(legacy_psf.computeKernelBBox(p)) 

1016 assert psf.compute_kernel_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeKernelImage(p)) 

1017 assert psf.compute_stellar_bbox(x=p.x, y=p.y) == Box.from_legacy(legacy_psf.computeImageBBox(p)) 

1018 assert psf.compute_stellar_image(x=p.x, y=p.y) == Image.from_legacy(legacy_psf.computeImage(p)) 

1019 n_points_tested += 1 

1020 return n_points_tested 

1021 

1022 

1023def compare_field_to_legacy( 

1024 field: BaseField, 

1025 legacy_field: Any, 

1026 subimage_bbox: Box, 

1027) -> None: 

1028 """Test a Field object by comparing it to an equivalent 

1029 `lsst.afw.math.BoundedField`. 

1030 

1031 Parameters 

1032 ---------- 

1033 field 

1034 Field to test. 

1035 legacy_field : ``lsst.afw.math.BoundedField`` 

1036 Equivalent legacy bounded field. 

1037 subimage_bbox 

1038 Bounding box for full-image tests. 

1039 """ 

1040 from lsst.afw.math import BoundedField as LegacyBoundedField 

1041 

1042 assert field.bounds.bbox == Box.from_legacy(legacy_field.getBBox()) 

1043 # Pixel coordinates to test the numpy array interface with. 

1044 pixel_xy = field.bounds.bbox.meshgrid(n=5).map(np.ravel) 

1045 if not isinstance(field.bounds, Box): 1045 ↛ 1046line 1045 didn't jump to line 1046 because the condition on line 1045 was never true

1046 mask = field.bounds.contains(x=pixel_xy.x, y=pixel_xy.y) 

1047 pixel_xy = pixel_xy.map(lambda v: v[mask]) 

1048 try: 

1049 assert_values_equal( 

1050 field(x=pixel_xy.x, y=pixel_xy.y), 

1051 legacy_field.evaluate(pixel_xy.x, pixel_xy.y), 

1052 # The field coefficients come from a numerical fit, so this can 

1053 # be bigger than just round-off error epsilon. 

1054 rtol=1e-5, 

1055 ) 

1056 except AssertionError as err: 

1057 err.add_note(f"evaluated at {pixel_xy}") 

1058 raise 

1059 if not isinstance(legacy_field, LegacyBoundedField): 1059 ↛ 1062line 1059 didn't jump to line 1062 because the condition on line 1059 was never true

1060 # Legacy StitchedApertureCorrection objects are not true BoundedFields 

1061 # and don't have addToImage. 

1062 return 

1063 legacy_image_1 = Image(0, bbox=subimage_bbox, dtype=np.float64).to_legacy() 

1064 legacy_field.addToImage(legacy_image_1, overlapOnly=True) 

1065 assert_images_equal( 

1066 field.render(subimage_bbox), Image.from_legacy(legacy_image_1, unit=field.unit), rtol=1e-13 

1067 ) 

1068 

1069 

1070def compare_aperture_corrections_to_legacy( 

1071 aperture_corrections: Mapping[str, BaseField], 

1072 legacy_ap_corr_map: Any, 

1073 subimage_bbox: Box, 

1074) -> None: 

1075 """Test an aperture correction `dict` by comparing it to an equivalent 

1076 `lsst.afw.image.ApCorrMap`. 

1077 

1078 Parameters 

1079 ---------- 

1080 aperture_corrections 

1081 Dictionary to test. 

1082 legacy_ap_corr_map : ``lsst.afw.image.ApCorrMap`` 

1083 Equivalent legacy aperture correction map. 

1084 subimage_bbox 

1085 Bounding box for full-image tests. 

1086 """ 

1087 assert aperture_corrections.keys() == set(legacy_ap_corr_map.keys()) 

1088 for name, field in aperture_corrections.items(): 

1089 compare_field_to_legacy(field, legacy_ap_corr_map[name], subimage_bbox) 

1090 

1091 

1092def compare_observation_summary_stats_to_legacy( 

1093 summary_stats: ObservationSummaryStats, 

1094 legacy_summary_stats: Any, 

1095) -> None: 

1096 """Test an ObservationSummaryStats object by comparing it to an equivalent 

1097 `lsst.afw.image.ExposureSummaryStats`. 

1098 

1099 Parameters 

1100 ---------- 

1101 summary_stats 

1102 Struct to test. 

1103 legacy_summary_stats : ``lsst.afw.image.ExposureSummaryStats`` 

1104 Equivalent legacy struct. 

1105 """ 

1106 for field in dataclasses.fields(legacy_summary_stats): 

1107 # Always skip version since that field has no meaning in the new 

1108 # type. 

1109 if field.name == "version": 

1110 continue 

1111 a = getattr(legacy_summary_stats, field.name) 

1112 b = getattr(summary_stats, field.name) 

1113 if isinstance(b, tuple): 

1114 for ai, bi in zip(a, b): 

1115 assert ai == bi or (math.isnan(ai) and math.isnan(bi)), f"{field.name}: {a} != {b}" 

1116 else: 

1117 assert a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}" 

1118 

1119 

1120def compare_sky_projection_to_legacy_wcs[F: Frame]( 

1121 sky_projection: SkyProjection[F], 

1122 legacy_wcs: Any, 

1123 pixel_frame: F, 

1124 subimage_bbox: Box, 

1125 is_fits: bool = False, 

1126) -> None: 

1127 """Test a Projection object by comparing it to an equivalent 

1128 `lsst.afw.geom.SkyWcs`. 

1129 

1130 Parameters 

1131 ---------- 

1132 sky_projection 

1133 Projection to test. 

1134 legacy_wcs : ``lsst.afw.geom.SkyWcs`` 

1135 Equivalent legacy WCS. 

1136 pixel_frame 

1137 Expected pixel frame for the sky_projection. 

1138 subimage_bbox 

1139 Bounding box of points to generate for tests. 

1140 is_fits 

1141 Whether this sky_projection is expected to be exactly representable as 

1142 a FITS WCS. If `False` it is assumed to have a FITS approximation 

1143 attached instead. 

1144 """ 

1145 # Pixel coordinates to test on over the subimage region of interest: 

1146 pixel_xy = subimage_bbox.meshgrid(step=50).map(np.ravel) 

1147 # Array indices of those pixel values (subtract off bbox starts): 

1148 subimage_array_xy = XY(x=pixel_xy.x - subimage_bbox.x.start, y=pixel_xy.y - subimage_bbox.y.start) 

1149 sky_coords = legacy_coords_to_astropy( 

1150 legacy_wcs.pixelToSky(arrays_to_legacy_points(pixel_xy.x, pixel_xy.y)) 

1151 ) 

1152 # Test transforming with the Projection itself, which also tests its 

1153 # nested Transform and an Astropy High-Level WCS view with no origin 

1154 # change. 

1155 check_projection(sky_projection, pixel_xy, sky_coords, pixel_frame, pixel_atol=_AST_INV_ATOL) 

1156 # Also test the Astropy High-Level WCS view with an origin change to 

1157 # array indices. 

1158 check_astropy_wcs_interface( 

1159 sky_projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5 

1160 ) 

1161 if is_fits: 1161 ↛ 1180line 1161 didn't jump to line 1180 because the condition on line 1161 was always true

1162 fits_wcs = sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=True) 

1163 assert fits_wcs is not None 

1164 check_astropy_wcs_interface(fits_wcs, subimage_array_xy, sky_coords, pixel_atol=1e-5) 

1165 # Use that FITS approximation to check that we can make a 

1166 # Projection from a FITS WCS, too. 

1167 fits_projection = SkyProjection.from_fits_wcs(fits_wcs, pixel_frame) 

1168 check_projection( 

1169 fits_projection, 

1170 subimage_array_xy, 

1171 sky_coords, 

1172 pixel_frame, 

1173 pixel_atol=1e-5, 

1174 ) 

1175 # We want Projections we create from a FITS WCS to be backed by an 

1176 # AST FrameSet so we can convert them into legacy 

1177 # `lsst.afw.geom.SkyWcs` objects if desired. 

1178 assert "Begin FrameSet" in fits_projection.show() 

1179 else: 

1180 assert sky_projection.as_fits_wcs(subimage_bbox, allow_approximation=False) is None 

1181 # The legacy SkyWcs should instead have a FITS approximation 

1182 # attached; run the same tests on that. 

1183 assert sky_projection.fits_approximation is not None 

1184 compare_sky_projection_to_legacy_wcs( 

1185 sky_projection.fits_approximation, 

1186 legacy_wcs.getFitsApproximation(), 

1187 pixel_frame, 

1188 subimage_bbox, 

1189 is_fits=True, 

1190 ) 

1191 

1192 

1193def check_transform[I: Frame, O: Frame]( 

1194 transform: Transform[I, O], 

1195 input_xy: XY[np.ndarray], 

1196 output_xy: XY[np.ndarray], 

1197 in_frame: Frame, 

1198 out_frame: Frame, 

1199 *, 

1200 check_inverted: bool = True, 

1201 in_atol: u.Quantity | float | None = None, 

1202 out_atol: u.Quantity | float | None = None, 

1203 angle_wrapped: bool = False, 

1204) -> None: 

1205 """Test Transform against known arrays of input and output points. 

1206 

1207 Parameters 

1208 ---------- 

1209 transform 

1210 Transform to test. 

1211 input_xy 

1212 Arrays of input points. 

1213 output_xy 

1214 Arrays of output points. 

1215 in_frame 

1216 Expected input frame. 

1217 out_frame 

1218 Expected output frame. 

1219 check_inverted 

1220 If `True`, recurse (once) to test the inverse transform. 

1221 in_atol 

1222 Expected absolute precision of input points. 

1223 out_atol 

1224 Expected absolute precision of output points. 

1225 angle_wrapped 

1226 If `True`, treat comparisons on sky-coordinate axes (as determined 

1227 by the frames) as angles modulo 2 pi, so differences that cross 

1228 zero or 360 degrees compare correctly. 

1229 """ 

1230 assert transform.in_frame == in_frame 

1231 assert transform.out_frame == out_frame 

1232 in_atol = 0.0 if in_atol is None else in_atol 

1233 out_atol = 0.0 if out_atol is None else out_atol 

1234 in_atol_v = in_atol.to_value(in_frame.unit) if isinstance(in_atol, u.Quantity) else float(in_atol) 

1235 out_atol_v = out_atol.to_value(out_frame.unit) if isinstance(out_atol, u.Quantity) else float(out_atol) 

1236 out_wrapped = angle_wrapped and out_frame == SkyFrame.ICRS 

1237 in_wrapped = angle_wrapped and in_frame == SkyFrame.ICRS 

1238 

1239 def assert_close( 

1240 test: np.ndarray | float | u.Quantity, 

1241 expected: np.ndarray | float | u.Quantity, 

1242 atol: float | u.Quantity, 

1243 *, 

1244 wrapped: bool, 

1245 ) -> None: 

1246 """Compare one axis, wrap-aware if that axis is a sky angle. 

1247 

1248 Parameters 

1249 ---------- 

1250 test 

1251 Values to test. 

1252 expected 

1253 Expected values. 

1254 atol 

1255 Absolute tolerance, in the units the values are expressed in. 

1256 wrapped 

1257 If `True`, treat the values as angles modulo 2 pi (radians). 

1258 """ 

1259 if not wrapped: 

1260 assert_values_equal(test, expected, atol=atol, rtol=_AST_RTOL) 

1261 return 

1262 test_values = np.asarray(getattr(test, "value", test)) 

1263 expected_values = np.asarray(getattr(expected, "value", expected)) 

1264 delta = (test_values - expected_values + np.pi) % (2.0 * np.pi) - np.pi 

1265 atol_v = atol.to_value(u.rad) if isinstance(atol, u.Quantity) else float(atol) 

1266 assert_values_equal(delta, np.zeros_like(delta), atol=atol_v, equal_nan=False) 

1267 

1268 # Test array interfaces. 

1269 test_output_xy = transform.apply_forward(x=input_xy.x, y=input_xy.y) 

1270 assert_close(test_output_xy.x, output_xy.x, out_atol_v, wrapped=out_wrapped) 

1271 assert_close(test_output_xy.y, output_xy.y, out_atol_v, wrapped=out_wrapped) 

1272 test_input_xy = transform.apply_inverse(x=output_xy.x, y=output_xy.y) 

1273 assert_close(test_input_xy.x, input_xy.x, in_atol_v, wrapped=in_wrapped) 

1274 assert_close(test_input_xy.y, input_xy.y, in_atol_v, wrapped=in_wrapped) 

1275 # Test scalar interfaces with numpy scalars. 

1276 for input_x, input_y, output_x, output_y in zip(input_xy.x, input_xy.y, output_xy.x, output_xy.y): 

1277 forward_xy = transform.apply_forward(x=input_x, y=input_y) 

1278 assert_close(forward_xy.x, output_x, out_atol_v, wrapped=out_wrapped) 

1279 assert_close(forward_xy.y, output_y, out_atol_v, wrapped=out_wrapped) 

1280 inverse_xy = transform.apply_inverse(x=output_x, y=output_y) 

1281 assert_close(inverse_xy.x, input_x, in_atol_v, wrapped=in_wrapped) 

1282 assert_close(inverse_xy.y, input_y, in_atol_v, wrapped=in_wrapped) 

1283 # Test quantity array interfaces. 

1284 input_q_xy = XY(x=input_xy.x * transform.in_frame.unit, y=input_xy.y * transform.in_frame.unit) 

1285 output_q_xy = XY(x=output_xy.x * transform.out_frame.unit, y=output_xy.y * transform.out_frame.unit) 

1286 test_output_q_xy = transform.apply_forward_q(x=input_q_xy.x, y=input_q_xy.y) 

1287 assert_close(test_output_q_xy.x, output_q_xy.x, out_atol, wrapped=out_wrapped) 

1288 assert_close(test_output_q_xy.y, output_q_xy.y, out_atol, wrapped=out_wrapped) 

1289 test_input_q_xy = transform.apply_inverse_q(x=output_q_xy.x, y=output_q_xy.y) 

1290 assert_close(test_input_q_xy.x, input_q_xy.x, in_atol, wrapped=in_wrapped) 

1291 assert_close(test_input_q_xy.y, input_q_xy.y, in_atol, wrapped=in_wrapped) 

1292 # Test quantity scalar interfaces. 

1293 for input_q_x, input_q_y, output_q_x, output_q_y in zip( 

1294 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y 

1295 ): 

1296 forward_q_xy = transform.apply_forward_q(x=input_q_x, y=input_q_y) 

1297 assert_close(forward_q_xy.x, output_q_x, out_atol, wrapped=out_wrapped) 

1298 assert_close(forward_q_xy.y, output_q_y, out_atol, wrapped=out_wrapped) 

1299 inverse_q_xy = transform.apply_inverse_q(x=output_q_x, y=output_q_y) 

1300 assert_close(inverse_q_xy.x, input_q_x, in_atol, wrapped=in_wrapped) 

1301 assert_close(inverse_q_xy.y, input_q_y, in_atol, wrapped=in_wrapped) 

1302 if check_inverted: 

1303 # Test the inverse transform. 

1304 check_transform( 

1305 transform.inverted(), 

1306 output_xy, 

1307 input_xy, 

1308 out_frame, 

1309 in_frame, 

1310 check_inverted=False, 

1311 out_atol=in_atol, 

1312 in_atol=out_atol, 

1313 angle_wrapped=angle_wrapped, 

1314 ) 

1315 

1316 

1317def check_projection[P: Frame]( 

1318 sky_projection: SkyProjection[P], 

1319 pixel_xy: XY[np.ndarray], 

1320 sky_coords: SkyCoord, 

1321 pixel_frame: Frame, 

1322 *, 

1323 pixel_atol: float | None = None, 

1324 sky_atol: u.Quantity | None = None, 

1325) -> None: 

1326 """Test a `.SkyProjection` instance against known arrays of pixel and sky 

1327 coordinates. 

1328 

1329 Parameters 

1330 ---------- 

1331 sky_projection 

1332 Projection to test. 

1333 pixel_xy 

1334 Arrays of pixel coordinates. 

1335 sky_coords 

1336 Corresponding sky coordinates. 

1337 pixel_frame 

1338 Expected pixel frame. 

1339 pixel_atol 

1340 Expected absolute precision of pixel points. 

1341 sky_atol 

1342 Expected absolute precision of sky coordinates, as a maximum 

1343 great-circle separation. 

1344 """ 

1345 assert sky_projection.pixel_frame == pixel_frame 

1346 assert sky_projection.sky_frame == SkyFrame.ICRS 

1347 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol) 

1348 pixel_atol = 0.0 if pixel_atol is None else pixel_atol 

1349 pixel_atol_q = pixel_atol * u.pix 

1350 # Test array interfaces. 

1351 test_pixel_xy = cast(XY[np.ndarray], sky_projection.sky_to_pixel(sky_coords)) 

1352 assert_values_equal(test_pixel_xy.x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL) 

1353 assert_values_equal(test_pixel_xy.y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL) 

1354 test_sky_astropy = sky_projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y) 

1355 assert_sky_coords_close(test_sky_astropy, sky_coords, sky_atol, label="pixel_to_sky") 

1356 # Test scalar interfaces. 

1357 for pixel_x, pixel_y, sky_single in zip(pixel_xy.x, pixel_xy.y, sky_coords): 

1358 assert_values_equal( 

1359 sky_projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol, rtol=_AST_RTOL 

1360 ) 

1361 assert_values_equal( 

1362 sky_projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol, rtol=_AST_RTOL 

1363 ) 

1364 test_sky_single = sky_projection.pixel_to_sky(x=pixel_x, y=pixel_y) 

1365 assert_sky_coords_close(test_sky_single, sky_single, sky_atol, label="pixel_to_sky") 

1366 # Test the underlying Transform object. 

1367 sky_xy = XY(x=sky_coords.ra.to_value(u.rad), y=sky_coords.dec.to_value(u.rad)) 

1368 check_transform( 

1369 sky_projection.pixel_to_sky_transform, 

1370 pixel_xy, 

1371 sky_xy, 

1372 pixel_frame, 

1373 SkyFrame.ICRS, 

1374 check_inverted=False, 

1375 in_atol=pixel_atol_q, 

1376 out_atol=sky_atol, 

1377 angle_wrapped=True, 

1378 ) 

1379 check_transform( 

1380 sky_projection.sky_to_pixel_transform, 

1381 sky_xy, 

1382 pixel_xy, 

1383 SkyFrame.ICRS, 

1384 pixel_frame, 

1385 check_inverted=False, 

1386 in_atol=sky_atol, 

1387 out_atol=pixel_atol_q, 

1388 angle_wrapped=True, 

1389 ) 

1390 # Test the Astropy interface adapter. 

1391 check_astropy_wcs_interface( 

1392 sky_projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol 

1393 ) 

1394 

1395 

1396def assert_sky_projections_equal( 

1397 a: SkyProjection[Any] | None, 

1398 b: SkyProjection[Any] | None, 

1399 expect_identity: bool | None = None, 

1400) -> None: 

1401 """Test that two `.SkyProjection` instances are equivalent. 

1402 

1403 Parameters 

1404 ---------- 

1405 a 

1406 First sky projection to compare. 

1407 b 

1408 Second sky projection to compare. 

1409 expect_identity 

1410 If not `None`, assert whether ``a`` and ``b`` are the same object. 

1411 """ 

1412 if a is None and b is None: 

1413 return 

1414 assert a is not None 

1415 assert b is not None 

1416 match expect_identity: 

1417 case True: 1417 ↛ 1418line 1417 didn't jump to line 1418 because the pattern on line 1417 never matched

1418 assert a is b 

1419 return 

1420 case False: 

1421 assert a is not b 

1422 case None if a is b: 

1423 return 

1424 assert a == b 

1425 

1426 

1427def check_astropy_wcs_interface( 

1428 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS, 

1429 pixel_xy: XY[np.ndarray], 

1430 sky_coords: SkyCoord, 

1431 *, 

1432 pixel_atol: float | None = None, 

1433 sky_atol: u.Quantity | None = None, 

1434) -> None: 

1435 """Test an Astropy WCS instance against known arrays of pixel and 

1436 sky coordinates. 

1437 

1438 Parameters 

1439 ---------- 

1440 wcs 

1441 Astropy WCS object to test. 

1442 pixel_xy 

1443 Arrays of pixel coordinates. 

1444 sky_coords 

1445 Corresponding sky coordinates. 

1446 pixel_atol 

1447 Expected absolute precision of pixel points. 

1448 sky_atol 

1449 Expected absolute precision of sky coordinates, as a maximum 

1450 great-circle separation. 

1451 """ 

1452 pixel_atol = 0.0 if pixel_atol is None else pixel_atol 

1453 sky_atol = _AST_SKY_ATOL if sky_atol is None else u.Quantity(sky_atol) 

1454 test_x, test_y = wcs.world_to_pixel(sky_coords) 

1455 assert_values_equal(test_x, pixel_xy.x, atol=pixel_atol, rtol=_AST_RTOL) 

1456 assert_values_equal(test_y, pixel_xy.y, atol=pixel_atol, rtol=_AST_RTOL) 

1457 test_sky_coords = wcs.pixel_to_world(pixel_xy.x, pixel_xy.y) 

1458 assert_sky_coords_close(test_sky_coords, sky_coords, sky_atol, label="pixel_to_world") 

1459 

1460 

1461def legacy_points_to_xy_array(legacy_points: list[Any]) -> XY[np.ndarray]: 

1462 """Convert a list of ``lsst.geom.Point2D`` objects to an `.XY` array. 

1463 

1464 Parameters 

1465 ---------- 

1466 legacy_points 

1467 Legacy ``lsst.geom.Point2D`` objects to convert. 

1468 """ 

1469 return XY(x=np.array([p.x for p in legacy_points]), y=np.array([p.y for p in legacy_points])) 

1470 

1471 

1472def legacy_coords_to_astropy(legacy_coords: list[Any]) -> SkyCoord: 

1473 """Convert a list of ``lsst.geom.SpherePoint`` objects to an Astropy 

1474 coordinate object. 

1475 

1476 Parameters 

1477 ---------- 

1478 legacy_coords 

1479 Legacy ``lsst.geom.SpherePoint`` objects to convert. 

1480 """ 

1481 return SkyCoord( 

1482 ra=np.array([p.getRa().asRadians() for p in legacy_coords]) * u.rad, 

1483 dec=np.array([p.getDec().asRadians() for p in legacy_coords]) * u.rad, 

1484 ) 

1485 

1486 

1487def arrays_to_legacy_points(x: np.ndarray, y: np.ndarray) -> list[Any]: 

1488 """Convert arrays of ``x`` and ``y`` to a list of ``lsst.geom.Point2D``. 

1489 

1490 Parameters 

1491 ---------- 

1492 x 

1493 X coordinates of the points. 

1494 y 

1495 Y coordinates of the points. 

1496 """ 

1497 from lsst.geom import Point2D 

1498 

1499 return [Point2D(x=xv, y=yv) for xv, yv in zip(x, y)] 

1500 

1501 

1502def compare_amplifier_to_legacy( 

1503 amplifier: Amplifier, 

1504 legacy_amplifier: Any, 

1505 *, 

1506 is_raw_assembled: bool, 

1507 expect_nominal_calibrations: bool = True, 

1508) -> None: 

1509 """Compare an `~.cameras.Amplifier` to a legacy 

1510 `lsst.afw.cameraGeom.Amplifier`. 

1511 

1512 Parameters 

1513 ---------- 

1514 amplifier 

1515 Amplifier to compare. 

1516 legacy_amplifier 

1517 Legacy `lsst.afw.cameraGeom.Amplifier` to compare against. 

1518 is_raw_assembled 

1519 Whether the raw geometry is expected to be the assembled-raw 

1520 geometry (`True`) or the unassembled-raw geometry (`False`). 

1521 expect_nominal_calibrations 

1522 Whether the amplifier is expected to carry nominal calibrations. 

1523 """ 

1524 assert legacy_amplifier.getName() == amplifier.name 

1525 assert Box.from_legacy(legacy_amplifier.getBBox()) == amplifier.bbox 

1526 if is_raw_assembled: 

1527 raw_geom = amplifier.assembled_raw_geometry 

1528 else: 

1529 raw_geom = amplifier.unassembled_raw_geometry 

1530 assert raw_geom is not None 

1531 assert ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()) == raw_geom.readout_corner 

1532 assert Box.from_legacy(legacy_amplifier.getRawBBox()) == raw_geom.bbox 

1533 assert Box.from_legacy(legacy_amplifier.getRawDataBBox()) == raw_geom.data_bbox 

1534 assert legacy_amplifier.getRawFlipX() == raw_geom.flip_x 

1535 assert legacy_amplifier.getRawFlipY() == raw_geom.flip_y 

1536 assert legacy_amplifier.getRawXYOffset().getX() == raw_geom.x_offset 

1537 assert legacy_amplifier.getRawXYOffset().getY() == raw_geom.y_offset 

1538 assert ( 

1539 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()) == raw_geom.horizontal_overscan_bbox 

1540 ) 

1541 assert Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()) == raw_geom.vertical_overscan_bbox 

1542 assert Box.from_legacy(legacy_amplifier.getRawPrescanBBox()) == raw_geom.horizontal_prescan_bbox 

1543 if expect_nominal_calibrations: 

1544 assert amplifier.nominal_calibrations is not None 

1545 assert_equal_allow_nan(legacy_amplifier.getGain(), amplifier.nominal_calibrations.gain) 

1546 assert_equal_allow_nan(legacy_amplifier.getReadNoise(), amplifier.nominal_calibrations.read_noise) 

1547 assert_equal_allow_nan(legacy_amplifier.getSaturation(), amplifier.nominal_calibrations.saturation) 

1548 assert_equal_allow_nan( 

1549 legacy_amplifier.getSuspectLevel(), amplifier.nominal_calibrations.suspect_level 

1550 ) 

1551 assert_values_equal( 

1552 legacy_amplifier.getLinearityCoeffs(), 

1553 amplifier.nominal_calibrations.linearity_coefficients, 

1554 label="linearity_coefficients", 

1555 ) 

1556 assert legacy_amplifier.getLinearityType() == amplifier.nominal_calibrations.linearity_type 

1557 

1558 

1559def compare_detector_to_legacy( 

1560 detector: Detector, 

1561 legacy_detector: Any, 

1562 *, 

1563 is_raw_assembled: bool, 

1564 expect_nominal_calibrations: bool = True, 

1565) -> None: 

1566 """Compare a `~.cameras.Detector` to a `lsst.afw.cameraGeom.Detector`. 

1567 

1568 Parameters 

1569 ---------- 

1570 detector 

1571 Detector to compare. 

1572 legacy_detector 

1573 Legacy `lsst.afw.cameraGeom.Detector` to compare against. 

1574 is_raw_assembled 

1575 Whether the raw geometry is expected to be the assembled-raw 

1576 geometry (`True`) or the unassembled-raw geometry (`False`). 

1577 expect_nominal_calibrations 

1578 Whether the detector's amplifiers are expected to carry nominal 

1579 calibrations. 

1580 """ 

1581 from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS 

1582 

1583 assert legacy_detector.getName() == detector.name 

1584 assert legacy_detector.getId() == detector.id 

1585 assert DetectorType.from_legacy(legacy_detector.getType()) == detector.type 

1586 assert Box.from_legacy(legacy_detector.getBBox()) == detector.bbox 

1587 assert legacy_detector.getSerial() == detector.serial 

1588 legacy_orientation = legacy_detector.getOrientation() 

1589 assert legacy_orientation.getFpPosition3().getX() == detector.orientation.focal_plane_x 

1590 assert legacy_orientation.getFpPosition3().getY() == detector.orientation.focal_plane_y 

1591 assert legacy_orientation.getFpPosition3().getZ() == detector.orientation.focal_plane_z 

1592 assert legacy_orientation.getReferencePoint().getX() == detector.orientation.pixel_reference_x 

1593 assert legacy_orientation.getReferencePoint().getY() == detector.orientation.pixel_reference_y 

1594 assert legacy_orientation.getYaw().asRadians() == detector.orientation.yaw.to_value(u.rad) 

1595 assert legacy_orientation.getPitch().asRadians() == detector.orientation.pitch.to_value(u.rad) 

1596 assert legacy_orientation.getRoll().asRadians() == detector.orientation.roll.to_value(u.rad) 

1597 assert legacy_detector.getPixelSize().getX() == detector.pixel_size 

1598 assert legacy_detector.getPhysicalType() == detector.physical_type 

1599 for amplifier, legacy_amplifier in zip(detector.amplifiers, legacy_detector.getAmplifiers(), strict=True): 

1600 compare_amplifier_to_legacy( 

1601 amplifier, 

1602 legacy_amplifier, 

1603 is_raw_assembled=is_raw_assembled, 

1604 expect_nominal_calibrations=expect_nominal_calibrations, 

1605 ) 

1606 pixel_xy = detector.bbox.meshgrid(n=3).map(lambda z: z.ravel().astype(np.float64)) 

1607 pixel_legacy_points = arrays_to_legacy_points(y=pixel_xy.y, x=pixel_xy.x) 

1608 fp_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FOCAL_PLANE) 

1609 check_transform( 

1610 detector.to_focal_plane, 

1611 pixel_xy, 

1612 legacy_points_to_xy_array(fp_legacy_points), 

1613 detector.frame, 

1614 detector.to_focal_plane.out_frame, 

1615 in_atol=1e-9 * u.pix, 

1616 out_atol=1e-7 * detector.to_focal_plane.out_frame.unit, 

1617 ) 

1618 fa_legacy_points = legacy_detector.transform(pixel_legacy_points, PIXELS, FIELD_ANGLE) 

1619 check_transform( 

1620 detector.to_field_angle, 

1621 pixel_xy, 

1622 legacy_points_to_xy_array(fa_legacy_points), 

1623 detector.frame, 

1624 detector.to_field_angle.out_frame, 

1625 in_atol=1e-9 * u.pix, 

1626 out_atol=1e-7 * u.arcsec, 

1627 ) 

1628 

1629 

1630def iter_concrete_archive_tree_subclasses() -> Iterator[type[ArchiveTree]]: 

1631 """Yield every importable concrete `.serialization.ArchiveTree` subclass. 

1632 

1633 Walks the ``ArchiveTree.__subclasses__()`` tree, skipping abstract 

1634 classes. Importing this module already imports every ``lsst.images`` 

1635 module that defines a subclass, so the tree is fully populated by the time 

1636 this is called. 

1637 

1638 This discovery is deliberately separate from 

1639 `check_archive_tree_class_invariants` so that the per-class check stays 

1640 usable on a single class even if this metaprogramming is removed later. 

1641 """ 

1642 seen: set[type] = set() 

1643 stack: list[type] = [ArchiveTree] 

1644 while stack: 

1645 kls = stack.pop() 

1646 for sub in kls.__subclasses__(): 

1647 if sub in seen: 1647 ↛ 1648line 1647 didn't jump to line 1648 because the condition on line 1647 was never true

1648 continue 

1649 seen.add(sub) 

1650 stack.append(sub) 

1651 if not getattr(sub, "__abstractmethods__", None): 1651 ↛ 1646line 1651 didn't jump to line 1646 because the condition on line 1651 was always true

1652 yield sub 

1653 

1654 

1655def check_bounds_contains_broadcasting(bounds: Bounds) -> None: 

1656 """Verify that `~lsst.images.Bounds.contains` accepts array-like inputs. 

1657 

1658 Uses the scalar overload as the reference and checks that 1-D arrays, 

1659 list inputs, mixed scalar-plus-array inputs, 2-D broadcast inputs, and 

1660 `.XY` / `.YX` positional-argument forms all produce results consistent 

1661 with calling the scalar overload on each ``(x, y)`` pair individually. 

1662 

1663 Parameters 

1664 ---------- 

1665 bounds 

1666 The `~lsst.images.Bounds` implementation to exercise. 

1667 """ 

1668 bbox = bounds.bbox 

1669 # One point outside each boundary, one on each boundary, one in the 

1670 # interior. The boundary points (start and stop) straddle inside/outside 

1671 # (start is inside, stop is outside) so we are guaranteed a mix of True 

1672 # and False results without hard-coding expected values. 

1673 ys = np.array( 

1674 [ 

1675 bbox.y.start - 1, 

1676 bbox.y.start, 

1677 (bbox.y.start + bbox.y.stop) // 2, 

1678 bbox.y.stop, 

1679 ] 

1680 ) 

1681 xs = np.array( 

1682 [ 

1683 bbox.x.start - 1, 

1684 bbox.x.start, 

1685 (bbox.x.start + bbox.x.stop) // 2, 

1686 bbox.x.stop, 

1687 ] 

1688 ) 

1689 # 2-D reference: expected[i, j] == contains(x=xs[j], y=ys[i]). 

1690 expected = np.array( 

1691 [[bounds.contains(x=int(xi), y=int(yi)) for xi in xs] for yi in ys], 

1692 dtype=bool, 

1693 ) 

1694 # 1-D ndarray (diagonal pairs). 

1695 assert_values_equal( 

1696 bounds.contains(x=xs, y=ys), 

1697 np.diagonal(expected), 

1698 label="bounds.contains", 

1699 ) 

1700 # list inputs (array-like). 

1701 assert_values_equal( 

1702 bounds.contains(x=xs.tolist(), y=ys.tolist()), 

1703 np.diagonal(expected), 

1704 label="bounds.contains", 

1705 ) 

1706 # Mixed: scalar y, 1-D array x — produces a 1-D result. 

1707 fixed_yi = 1 # index into ys; ys[1] == bbox.y.start (on the boundary) 

1708 assert_values_equal( 

1709 bounds.contains(x=xs, y=int(ys[fixed_yi])), 

1710 expected[fixed_yi, :], 

1711 label="bounds.contains", 

1712 ) 

1713 # Mixed: 1-D array y, scalar x. 

1714 fixed_xi = 1 # index into xs; xs[1] == bbox.x.start (on the boundary) 

1715 assert_values_equal( 

1716 bounds.contains(x=int(xs[fixed_xi]), y=ys), 

1717 expected[:, fixed_xi], 

1718 label="bounds.contains", 

1719 ) 

1720 # Float scalars: results must match the int-scalar reference. 

1721 assert bounds.contains(x=float(xs[fixed_xi]), y=float(ys[fixed_yi])) == expected[fixed_yi, fixed_xi] 

1722 # XY / YX scalar: results must match the keyword-scalar reference. 

1723 assert bounds.contains(XY(x=int(xs[fixed_xi]), y=int(ys[fixed_yi]))) == expected[fixed_yi, fixed_xi] 

1724 assert bounds.contains(YX(y=int(ys[fixed_yi]), x=int(xs[fixed_xi]))) == expected[fixed_yi, fixed_xi] 

1725 # XY / YX array: results must match the 1-D ndarray reference. 

1726 assert_values_equal(bounds.contains(XY(x=xs, y=ys)), np.diagonal(expected), label="bounds.contains") 

1727 assert_values_equal(bounds.contains(YX(y=ys, x=xs)), np.diagonal(expected), label="bounds.contains") 

1728 # 2-D broadcast: y shape (N, 1) × x shape (1, M) → (N, M). 

1729 assert_values_equal( 

1730 bounds.contains(x=xs.reshape(1, -1), y=ys.reshape(-1, 1)), 

1731 expected, 

1732 label="bounds.contains", 

1733 ) 

1734 

1735 

1736def check_archive_tree_class_invariants(tree_type: type[ArchiveTree]) -> None: 

1737 """Assert that one concrete `.serialization.ArchiveTree` subclass declares 

1738 well-formed schema-version constants and an in-memory type. 

1739 

1740 Checks that ``SCHEMA_NAME``, ``SCHEMA_VERSION``, ``MIN_READ_VERSION`` and 

1741 ``PUBLIC_TYPE`` are present and well-typed, that the version is 

1742 ``major.minor.patch``, and that ``MIN_READ_VERSION`` does not exceed the 

1743 schema major. 

1744 

1745 Parameters 

1746 ---------- 

1747 tree_type 

1748 The concrete `.serialization.ArchiveTree` subclass to check. 

1749 """ 

1750 assert hasattr(tree_type, "SCHEMA_NAME"), f"{tree_type.__name__} lacks SCHEMA_NAME" 

1751 assert hasattr(tree_type, "SCHEMA_VERSION"), f"{tree_type.__name__} lacks SCHEMA_VERSION" 

1752 assert hasattr(tree_type, "MIN_READ_VERSION"), f"{tree_type.__name__} lacks MIN_READ_VERSION" 

1753 assert hasattr(tree_type, "PUBLIC_TYPE"), f"{tree_type.__name__} lacks PUBLIC_TYPE" 

1754 assert isinstance(tree_type.SCHEMA_NAME, str) 

1755 assert len(tree_type.SCHEMA_NAME) > 0 

1756 # Allow an optional PEP 440 development-release suffix (e.g. 1.0.0.dev0) 

1757 # for schemas still in development. 

1758 assert re.fullmatch(r"^\d+\.\d+\.\d+(\.dev\d+)?$", tree_type.SCHEMA_VERSION) 

1759 assert isinstance(tree_type.MIN_READ_VERSION, int) 

1760 assert tree_type.MIN_READ_VERSION >= 1 

1761 assert isinstance(tree_type.PUBLIC_TYPE, type) 

1762 major = int(tree_type.SCHEMA_VERSION.split(".")[0]) 

1763 assert tree_type.MIN_READ_VERSION <= major