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

405 statements  

« prev     ^ index     » next       coverage.py v7.14.1, created at 2026-05-29 08:40 +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. 

11 

12from __future__ import annotations 

13 

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

23 "assert_psfs_equal", 

24 "assert_visit_images_equal", 

25 "check_astropy_wcs_interface", 

26 "check_projection", 

27 "check_transform", 

28 "compare_amplifier_to_legacy", 

29 "compare_aperture_corrections_to_legacy", 

30 "compare_cell_coadd_to_legacy", 

31 "compare_detector_to_legacy", 

32 "compare_field_to_legacy", 

33 "compare_image_to_legacy", 

34 "compare_mask_to_legacy", 

35 "compare_masked_image_to_legacy", 

36 "compare_observation_summary_stats_to_legacy", 

37 "compare_photo_calib_to_legacy", 

38 "compare_projection_to_legacy_wcs", 

39 "compare_psf_to_legacy", 

40 "compare_visit_image_to_legacy", 

41 "legacy_coords_to_astropy", 

42 "legacy_points_to_xy_array", 

43) 

44 

45import dataclasses 

46import math 

47import unittest 

48from collections.abc import Mapping 

49from typing import TYPE_CHECKING, Any, Literal, cast 

50 

51import astropy.units as u 

52import astropy.wcs.wcsapi 

53import numpy as np 

54from astropy.coordinates import SkyCoord 

55 

56from .._geom import XY, YX, Box 

57from .._image import Image 

58from .._mask import Mask, MaskPlane 

59from .._masked_image import MaskedImage 

60from .._observation_summary_stats import ObservationSummaryStats 

61from .._transforms import DetectorFrame, Frame, Projection, SkyFrame, TractFrame, Transform 

62from .._visit_image import VisitImage 

63from ..aperture_corrections import ApertureCorrectionMap 

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

65from ..cells import CellCoadd, CellIJ, CoaddProvenance 

66from ..fields import BaseField, ChebyshevField 

67from ..psfs import PointSpreadFunction 

68 

69if TYPE_CHECKING: 

70 try: 

71 from lsst.cell_coadds import MultipleCellCoadd 

72 except ImportError: 

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

74 try: 

75 from lsst.afw.image import PhotoCalib as LegacyPhotoCalib 

76 except ImportError: 

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

78 

79 

80def assert_close( 

81 tc: unittest.TestCase, 

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

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

84 **kwargs: Any, 

85) -> None: 

86 """Test that two arrays, floats, or quantities are almost equal. 

87 

88 Parameters 

89 ---------- 

90 tc 

91 Test case object with assert methods to use. 

92 a 

93 Array, float, or quantity to compare. 

94 b 

95 Array, float, or quantity to compare. 

96 **kwargs 

97 Forwarded to `astropy.units.allclose`. 

98 """ 

99 tc.assertTrue(u.allclose(a, b, **kwargs), msg=f"{a} != {b}") 

100 

101 

102def assert_equal_allow_nan(tc: unittest.TestCase, a: float, b: float) -> None: 

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

104 try: 

105 tc.assertEqual(a, b) 

106 except AssertionError: 

107 if not (math.isnan(a) and math.isnan(b)): 

108 raise 

109 

110 

111def assert_images_equal( 

112 tc: unittest.TestCase, 

113 a: Image, 

114 b: Image, 

115 *, 

116 rtol: float = 0.0, 

117 atol: float = 0.0, 

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

119) -> None: 

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

121 tc.assertEqual(a.bbox, b.bbox) 

122 tc.assertEqual(a.unit, b.unit) 

123 assert_projections_equal(tc, a.projection, b.projection) 

124 if expect_view is not None: 

125 tc.assertEqual(np.may_share_memory(a.array, b.array), bool(expect_view)) 

126 if expect_view == "array": 

127 tc.assertEqual(a.metadata, b.metadata) 

128 else: 

129 tc.assertEqual(a.metadata is b.metadata, expect_view) 

130 if not expect_view: 

131 assert_close(tc, a.array, b.array, atol=atol, rtol=rtol) 

132 tc.assertEqual(a.metadata, b.metadata) 

133 

134 

135def assert_masks_equal(tc: unittest.TestCase, a: Mask, b: Mask) -> None: 

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

137 tc.assertEqual(a.bbox, b.bbox) 

138 tc.assertEqual(a.schema, b.schema) 

139 tc.assertEqual(a.metadata, b.metadata) 

140 assert_projections_equal(tc, a.projection, b.projection) 

141 np.testing.assert_array_equal(a.array, b.array) 

142 

143 

144def assert_masked_images_equal( 

145 tc: unittest.TestCase, 

146 a: MaskedImage, 

147 b: MaskedImage, 

148 *, 

149 rtol: float = 0.0, 

150 atol: float = 0.0, 

151 expect_view: bool | None = None, 

152) -> None: 

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

154 tc.assertEqual(a.metadata, b.metadata) 

155 assert_projections_equal(tc, a.projection, b.projection) 

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

157 assert_masks_equal(tc, a.mask, b.mask) 

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

159 

160 

161def assert_psfs_equal( 

162 tc: unittest.TestCase, 

163 psf1: PointSpreadFunction, 

164 psf2: PointSpreadFunction, 

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

166) -> int: 

167 """Compare two PSF objets. 

168 

169 Parameters 

170 ---------- 

171 tc 

172 Test case object with assert methods to use. 

173 psf1 

174 Point-spread function to test. 

175 psf2 

176 The other point-spread function to test. 

177 points 

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

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

180 

181 Returns 

182 ------- 

183 `int` 

184 The number of points actually tested. 

185 """ 

186 if points is None: 

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

188 

189 tc.assertEqual(psf1.kernel_bbox, psf2.kernel_bbox) 

190 

191 n_points_tested: int = 0 

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

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

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

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

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

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

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

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

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

201 tc.assertEqual( 

202 contains1, 

203 contains2, 

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

205 ) 

206 if not contains1: 

207 continue 

208 tc.assertEqual(psf1.compute_kernel_image(x=x, y=y), psf2.compute_kernel_image(x=x, y=y)) 

209 tc.assertEqual(psf1.compute_stellar_bbox(x=x, y=y), psf2.compute_stellar_bbox(x=x, y=y)) 

210 tc.assertEqual(psf1.compute_stellar_image(x=x, y=y), psf2.compute_stellar_image(x=x, y=y)) 

211 n_points_tested += 1 

212 return n_points_tested 

213 

214 

215def assert_visit_images_equal( 

216 tc: unittest.TestCase, 

217 a: VisitImage, 

218 b: VisitImage, 

219 *, 

220 expect_view: bool | None = None, 

221) -> None: 

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

223 

224 Extends `assert_masked_images_equal` with the VisitImage-specific 

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

226 corrections, photometric scaling, backgrounds, polygon bounds, 

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

228 silently miss differences in any of them. 

229 """ 

230 assert_masked_images_equal(tc, a, b, expect_view=expect_view) 

231 tc.assertEqual(a.summary_stats, b.summary_stats) 

232 tc.assertEqual(a.physical_filter, b.physical_filter) 

233 tc.assertEqual(a.band, b.band) 

234 tc.assertEqual(a.obs_info, b.obs_info) 

235 tc.assertEqual(a.detector, b.detector) 

236 tc.assertEqual(dict(a.aperture_corrections), dict(b.aperture_corrections)) 

237 tc.assertEqual(a.photometric_scaling, b.photometric_scaling) 

238 tc.assertEqual(dict(a.backgrounds), dict(b.backgrounds)) 

239 tc.assertEqual(a.backgrounds.subtracted, b.backgrounds.subtracted) 

240 tc.assertEqual(a.bounds, b.bounds) 

241 assert_psfs_equal(tc, a.psf, b.psf) 

242 

243 

244def assert_cell_coadds_equal( 

245 tc: unittest.TestCase, 

246 a: CellCoadd, 

247 b: CellCoadd, 

248 *, 

249 expect_view: bool | None = None, 

250) -> None: 

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

252 

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

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

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

256 silently miss differences in any of them. 

257 """ 

258 assert_masked_images_equal(tc, a, b, expect_view=expect_view) 

259 tc.assertEqual(a.band, b.band) 

260 tc.assertEqual(a.patch, b.patch) 

261 tc.assertEqual(a.tract, b.tract) 

262 tc.assertEqual(a.grid, b.grid) 

263 tc.assertEqual(a.bounds.missing, b.bounds.missing) 

264 tc.assertEqual(dict(a.backgrounds), dict(b.backgrounds)) 

265 tc.assertEqual(a.backgrounds.subtracted, b.backgrounds.subtracted) 

266 assert_psfs_equal(tc, a.psf, b.psf) 

267 

268 

269def compare_image_to_legacy( 

270 tc: unittest.TestCase, image: Image, legacy_image: Any, expect_view: bool | None = None 

271) -> None: 

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

273 tc.assertEqual(image.bbox, Box.from_legacy(legacy_image.getBBox())) 

274 if expect_view is not None: 

275 tc.assertEqual(np.may_share_memory(image.array, legacy_image.array), expect_view) 

276 if not expect_view: 

277 np.testing.assert_array_equal(image.array, legacy_image.array) 

278 

279 

280def compare_mask_to_legacy( 

281 tc: unittest.TestCase, mask: Mask, legacy_mask: Any, plane_map: Mapping[str, MaskPlane] | None = None 

282) -> None: 

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

284 tc.assertEqual(mask.bbox, Box.from_legacy(legacy_mask.getBBox())) 

285 if plane_map is None: 

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

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

288 np.testing.assert_array_equal( 

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

290 mask.get(new_plane.name), 

291 ) 

292 

293 

294def compare_masked_image_to_legacy( 

295 tc: unittest.TestCase, 

296 masked_image: MaskedImage, 

297 legacy_masked_image: Any, 

298 *, 

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

300 expect_view: bool | None = None, 

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

302) -> None: 

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

304 object. 

305 

306 Parameters 

307 ---------- 

308 tc 

309 Test case to use for asserts. 

310 masked_image 

311 New image to test. 

312 legacy_masked_image 

313 Legacy image to test against. 

314 plane_map 

315 Mapping between new and legacy mask planes. 

316 expect_view 

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

318 memory. 

319 alternates 

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

321 check against the legacy versions of those components. 

322 """ 

323 compare_image_to_legacy(tc, masked_image.image, legacy_masked_image.getImage(), expect_view=expect_view) 

324 compare_mask_to_legacy(tc, masked_image.mask, legacy_masked_image.getMask(), plane_map=plane_map) 

325 compare_image_to_legacy( 

326 tc, masked_image.variance, legacy_masked_image.getVariance(), expect_view=expect_view 

327 ) 

328 if alternates: 

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

330 compare_image_to_legacy(tc, image, legacy_masked_image.getImage(), expect_view=expect_view) 

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

332 compare_mask_to_legacy(tc, mask, legacy_masked_image.getMask(), plane_map=plane_map) 

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

334 compare_image_to_legacy(tc, variance, legacy_masked_image.getVariance(), expect_view=expect_view) 

335 

336 

337def compare_visit_image_to_legacy( 

338 tc: unittest.TestCase, 

339 visit_image: VisitImage, 

340 legacy_exposure: Any, 

341 *, 

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

343 expect_view: bool | None = None, 

344 instrument: str, 

345 visit: int, 

346 detector: int, 

347 applied_legacy_photo_calib: LegacyPhotoCalib | None = None, 

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

349) -> None: 

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

351 object. 

352 

353 Parameters 

354 ---------- 

355 tc 

356 Test case to use for asserts. 

357 visit_image 

358 New image to test. 

359 legacy_exposure 

360 Legacy image to test against. 

361 plane_map 

362 Mapping between new and legacy mask planes. 

363 expect_view 

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

365 memory. 

366 instrument 

367 Expected instrument name. 

368 visit 

369 Expected visit ID. 

370 detector 

371 Expected detector ID. 

372 alternates 

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

374 check against the legacy versions of those components. 

375 """ 

376 compare_masked_image_to_legacy( 

377 tc, 

378 visit_image, 

379 legacy_exposure.getMaskedImage(), 

380 plane_map=plane_map, 

381 expect_view=expect_view, 

382 alternates=alternates, 

383 ) 

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

385 compare_projection_to_legacy_wcs( 

386 tc, 

387 visit_image.projection, 

388 legacy_exposure.getWcs(), 

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

390 visit_image.bbox, 

391 ) 

392 tc.assertIs(visit_image.projection, visit_image.mask.projection) 

393 tc.assertIs(visit_image.projection, visit_image.variance.projection) 

394 compare_psf_to_legacy(tc, visit_image.psf, legacy_exposure.getPsf()) 

395 compare_observation_summary_stats_to_legacy( 

396 tc, visit_image.summary_stats, legacy_exposure.info.getSummaryStats() 

397 ) 

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

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

400 # get expensive otherwisre. 

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

402 compare_aperture_corrections_to_legacy( 

403 tc, visit_image.aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox 

404 ) 

405 compare_photo_calib_to_legacy( 

406 tc, 

407 visit_image.photometric_scaling, 

408 legacy_exposure.info.getPhotoCalib(), 

409 applied_legacy_photo_calib=applied_legacy_photo_calib, 

410 subimage_bbox=tiny_bbox, 

411 ) 

412 if alternates: 

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

414 tc.assertEqual(bbox, visit_image.bbox) 

415 if projection := alternates.get("projection"): 

416 compare_projection_to_legacy_wcs( 

417 tc, 

418 projection, 

419 legacy_exposure.getWcs(), 

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

421 visit_image.bbox, 

422 ) 

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

424 compare_psf_to_legacy(tc, psf, legacy_exposure.getPsf()) 

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

426 compare_observation_summary_stats_to_legacy( 

427 tc, summary_stats, legacy_exposure.info.getSummaryStats() 

428 ) 

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

430 compare_detector_to_legacy(tc, detector_obj, legacy_exposure.getDetector(), is_raw_assembled=True) 

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

432 visitInfo = legacy_exposure.visitInfo 

433 tc.assertEqual(obs_info.instrument, visitInfo.getInstrumentLabel()) 

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

435 compare_aperture_corrections_to_legacy( 

436 tc, aperture_corrections, legacy_exposure.info.getApCorrMap(), tiny_bbox 

437 ) 

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

439 compare_photo_calib_to_legacy( 

440 tc, 

441 photometric_scaling, 

442 legacy_exposure.info.getPhotoCalib(), 

443 applied_legacy_photo_calib=applied_legacy_photo_calib, 

444 subimage_bbox=tiny_bbox, 

445 ) 

446 

447 

448def compare_photo_calib_to_legacy( 

449 tc: unittest.TestCase, 

450 photometric_scaling: BaseField | None, 

451 legacy_photo_calib: LegacyPhotoCalib, 

452 *, 

453 applied_legacy_photo_calib: LegacyPhotoCalib | None = None, 

454 subimage_bbox: Box, 

455) -> None: 

456 if legacy_photo_calib._isConstant: 

457 if legacy_photo_calib.getCalibrationMean() == 1.0: 

458 if applied_legacy_photo_calib is None: 

459 tc.assertIsNone(photometric_scaling) 

460 return 

461 else: 

462 legacy_photo_calib = applied_legacy_photo_calib 

463 if legacy_photo_calib._isConstant: 

464 assert isinstance(photometric_scaling, ChebyshevField) 

465 assert_close( 

466 tc, photometric_scaling.coefficients, np.array([[legacy_photo_calib.getCalibrationMean()]]) 

467 ) 

468 else: 

469 assert photometric_scaling is not None 

470 compare_field_to_legacy( 

471 tc, 

472 photometric_scaling / legacy_photo_calib.getCalibrationMean(), 

473 legacy_photo_calib.computeScaledCalibration(), 

474 subimage_bbox, 

475 ) 

476 

477 

478def compare_cell_coadd_to_legacy( 

479 tc: unittest.TestCase, 

480 cell_coadd: CellCoadd, 

481 legacy_cell_coadd: MultipleCellCoadd, 

482 *, 

483 tract_bbox: Box, 

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

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

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

487) -> None: 

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

489 `lsst.cell_coadds.MultipleCellCoadd` object. 

490 

491 Parameters 

492 ---------- 

493 tc 

494 Test case to use for asserts. 

495 cell_coadd 

496 New coadd to test. 

497 legacy_cell_coadd 

498 Legacy coadd to test against. 

499 tract_bbox 

500 Bounding box of the full tract. 

501 psf_points 

502 Points to use to compare the PSFs. 

503 plane_map 

504 Mapping between new and legacy mask planes. 

505 alternates 

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

507 check against the legacy versions of those components. 

508 """ 

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

510 compare_image_to_legacy(tc, cell_coadd.image, legacy_stitched.image, expect_view=False) 

511 compare_mask_to_legacy(tc, cell_coadd.mask, legacy_stitched.mask, plane_map=plane_map) 

512 compare_image_to_legacy(tc, cell_coadd.variance, legacy_stitched.variance, expect_view=False) 

513 if legacy_stitched.mask_fractions is not None: 

514 compare_image_to_legacy( 

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

516 ) 

517 for n in range(legacy_stitched.n_noise_realizations): 

518 compare_image_to_legacy( 

519 tc, cell_coadd.noise_realizations[n], legacy_stitched.noise_realizations[n], expect_view=False 

520 ) 

521 tc.assertEqual(cell_coadd.skymap, legacy_stitched.identifiers.skymap) 

522 tc.assertEqual(cell_coadd.tract, legacy_stitched.identifiers.tract) 

523 tc.assertEqual(cell_coadd.patch.index.x, legacy_stitched.identifiers.patch.x) 

524 tc.assertEqual(cell_coadd.patch.index.y, legacy_stitched.identifiers.patch.y) 

525 tc.assertEqual(cell_coadd.band, legacy_stitched.identifiers.band) 

526 tc.assertTrue(tract_bbox.contains(cell_coadd.patch.outer_bbox)) 

527 tc.assertTrue(cell_coadd.patch.outer_bbox.contains(cell_coadd.patch.inner_bbox)) 

528 tc.assertTrue(cell_coadd.patch.outer_bbox.contains(cell_coadd.bbox)) 

529 tc.assertEqual(cell_coadd.unit, u.Unit(legacy_cell_coadd.common.units.value)) 

530 tc.assertTrue(cell_coadd.bounds.bbox.contains(cell_coadd.bbox)) 

531 tc.assertTrue(cell_coadd.grid.bbox.contains(cell_coadd.bbox)) 

532 compare_projection_to_legacy_wcs( 

533 tc, 

534 cell_coadd.projection, 

535 legacy_cell_coadd.wcs, 

536 TractFrame( 

537 skymap=legacy_cell_coadd.identifiers.skymap, 

538 tract=legacy_cell_coadd.identifiers.tract, 

539 bbox=tract_bbox, 

540 ), 

541 cell_coadd.bbox, 

542 is_fits=True, 

543 ) 

544 tc.assertIs(cell_coadd.projection, cell_coadd.mask.projection) 

545 tc.assertIs(cell_coadd.projection, cell_coadd.variance.projection) 

546 compare_psf_to_legacy( 

547 tc, cell_coadd.psf, legacy_stitched.psf, expect_legacy_raise_on_out_of_bounds=True, points=psf_points 

548 ) 

549 compare_cell_coadd_provenance_to_legacy(tc, cell_coadd.provenance, legacy_cell_coadd) 

550 if alternates: 

551 if projection := alternates.get("projection"): 

552 compare_projection_to_legacy_wcs( 

553 tc, 

554 projection, 

555 legacy_stitched.wcs, 

556 TractFrame( 

557 skymap=legacy_cell_coadd.identifiers.skymap, 

558 tract=legacy_cell_coadd.identifiers.tract, 

559 bbox=tract_bbox, 

560 ), 

561 cell_coadd.bbox, 

562 is_fits=True, 

563 ) 

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

565 compare_psf_to_legacy(tc, psf, legacy_stitched.psf, points=psf_points) 

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

567 compare_cell_coadd_provenance_to_legacy(tc, provenance, legacy_cell_coadd) 

568 

569 

570def compare_cell_coadd_provenance_to_legacy( 

571 tc: unittest.TestCase, provenance: CoaddProvenance, legacy_cell_coadd: MultipleCellCoadd 

572) -> None: 

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

574 `lsst.cell_coadds.MultipleCellCoadd` object. 

575 

576 Parameters 

577 ---------- 

578 tc 

579 Test case to use for asserts. 

580 provenance 

581 New provenance object to test. 

582 legacy_cell_coadd 

583 Legacy coadd to test against. 

584 """ 

585 from lsst.cell_coadds import ObservationIdentifiers 

586 

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

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

589 prov = provenance[cell_index] 

590 legacy_table = astropy.table.Table( 

591 rows=[ 

592 [ 

593 ids.instrument, 

594 ids.visit, 

595 ids.detector, 

596 ids.day_obs, 

597 ids.physical_filter, 

598 legacy_input.overlaps_center, 

599 legacy_input.overlap_fraction, 

600 legacy_input.weight, 

601 legacy_input.psf_shape.getIxx(), 

602 legacy_input.psf_shape.getIyy(), 

603 legacy_input.psf_shape.getIxy(), 

604 legacy_input.psf_shape_flag, 

605 ] 

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

607 ], 

608 dtype=[ 

609 np.object_, 

610 np.uint64, 

611 np.uint16, 

612 np.uint32, 

613 np.object_, 

614 np.bool_, 

615 np.float64, 

616 np.float64, 

617 np.float64, 

618 np.float64, 

619 np.float64, 

620 np.bool_, 

621 ], 

622 names=[ 

623 "instrument", 

624 "visit", 

625 "detector", 

626 "day_obs", 

627 "physical_filter", 

628 "overlaps_center", 

629 "overlap_fraction", 

630 "weight", 

631 "psf_shape_xx", 

632 "psf_shape_yy", 

633 "psf_shape_xy", 

634 "psf_shape_flag", 

635 ], 

636 ) 

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

638 tc.assertEqual(len(legacy_cell.inputs), len(prov.inputs)) 

639 tc.assertEqual(len(legacy_cell.inputs), len(prov.contributions)) 

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

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

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

643 np.testing.assert_array_equal(prov.inputs["instrument"], prov.contributions["instrument"]) 

644 np.testing.assert_array_equal(prov.inputs["visit"], prov.contributions["visit"]) 

645 np.testing.assert_array_equal(prov.inputs["detector"], prov.contributions["detector"]) 

646 np.testing.assert_array_equal(prov.inputs["instrument"], legacy_table["instrument"]) 

647 np.testing.assert_array_equal(prov.inputs["visit"], legacy_table["visit"]) 

648 np.testing.assert_array_equal(prov.inputs["detector"], legacy_table["detector"]) 

649 np.testing.assert_array_equal(prov.inputs["physical_filter"], legacy_table["physical_filter"]) 

650 np.testing.assert_array_equal(prov.inputs["day_obs"], legacy_table["day_obs"]) 

651 np.testing.assert_array_equal(prov.contributions["overlaps_center"], legacy_table["overlaps_center"]) 

652 np.testing.assert_array_equal( 

653 prov.contributions["overlap_fraction"], legacy_table["overlap_fraction"] 

654 ) 

655 np.testing.assert_array_equal(prov.contributions["weight"], legacy_table["weight"]) 

656 np.testing.assert_array_equal(prov.contributions["psf_shape_xx"], legacy_table["psf_shape_xx"]) 

657 np.testing.assert_array_equal(prov.contributions["psf_shape_yy"], legacy_table["psf_shape_yy"]) 

658 np.testing.assert_array_equal(prov.contributions["psf_shape_xy"], legacy_table["psf_shape_xy"]) 

659 np.testing.assert_array_equal(prov.contributions["psf_shape_flag"], legacy_table["psf_shape_flag"]) 

660 for row in prov.inputs: 

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

662 legacy_polygon = legacy_cell_coadd.common.visit_polygons[polygon_key] 

663 tc.assertEqual(legacy_polygon, row["polygon"].to_legacy()) 

664 

665 

666def compare_psf_to_legacy( 

667 tc: unittest.TestCase, 

668 psf: PointSpreadFunction, 

669 legacy_psf: Any, 

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

671 expect_legacy_raise_on_out_of_bounds: bool = False, 

672) -> int: 

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

674 

675 Parameters 

676 ---------- 

677 tc 

678 Test case object with assert methods to use. 

679 psf 

680 Point-spread function to test. 

681 legacy_psf 

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

683 points 

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

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

686 expect_legacy_raise_on_out_of_bounds 

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

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

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

690 

691 Returns 

692 ------- 

693 `int` 

694 The number of points actually tested. 

695 """ 

696 from lsst.afw.detection import InvalidPsfError 

697 

698 if points is None: 

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

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

701 n_points_tested: int = 0 

702 for p in legacy_points: 

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

704 if expect_legacy_raise_on_out_of_bounds: 

705 with tc.assertRaises(InvalidPsfError): 

706 legacy_psf.computeKernelImage(p) 

707 continue 

708 tc.assertEqual(psf.kernel_bbox, Box.from_legacy(legacy_psf.computeKernelBBox(p))) 

709 tc.assertEqual( 

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

711 ) 

712 tc.assertEqual( 

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

714 ) 

715 tc.assertEqual(psf.compute_stellar_image(x=p.x, y=p.y), Image.from_legacy(legacy_psf.computeImage(p))) 

716 n_points_tested += 1 

717 return n_points_tested 

718 

719 

720def compare_field_to_legacy( 

721 tc: unittest.TestCase, 

722 field: BaseField, 

723 legacy_field: Any, 

724 subimage_bbox: Box, 

725) -> None: 

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

727 `lsst.afw.math.BoundedField`. 

728 

729 Parameters 

730 ---------- 

731 tc 

732 Test case object with assert methods to use. 

733 field 

734 Field to test. 

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

736 Equivalent legacy bounded field. 

737 subimage_bbox 

738 Bounding box for full-image tests. 

739 """ 

740 tc.assertEqual(field.bounds.bbox, Box.from_legacy(legacy_field.getBBox())) 

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

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

743 assert_close(tc, field(x=pixel_xy.x, y=pixel_xy.y), legacy_field.evaluate(pixel_xy.x, pixel_xy.y)) 

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

745 legacy_field.addToImage(legacy_image_1, overlapOnly=True) 

746 assert_images_equal( 

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

748 ) 

749 

750 

751def compare_aperture_corrections_to_legacy( 

752 tc: unittest.TestCase, 

753 aperture_corrections: ApertureCorrectionMap, 

754 legacy_ap_corr_map: Any, 

755 subimage_bbox: Box, 

756) -> None: 

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

758 `lsst.afw.image.ApCorrMap`. 

759 

760 Parameters 

761 ---------- 

762 tc 

763 Test case object with assert methods to use. 

764 aperture_corrections 

765 Dictionary to test. 

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

767 Equivalent legacy aperture correction map. 

768 subimage_bbox 

769 Bounding box for full-image tests. 

770 """ 

771 tc.assertEqual(aperture_corrections.keys(), set(legacy_ap_corr_map.keys())) 

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

773 compare_field_to_legacy(tc, field, legacy_ap_corr_map[name], subimage_bbox) 

774 

775 

776def compare_observation_summary_stats_to_legacy( 

777 tc: unittest.TestCase, 

778 summary_stats: ObservationSummaryStats, 

779 legacy_summary_stats: Any, 

780) -> None: 

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

782 `lsst.afw.image.ExposureSummaryStats`. 

783 

784 Parameters 

785 ---------- 

786 tc 

787 Test case object with assert methods to use. 

788 summary_stats 

789 Struct to test. 

790 legacy : ``lsst.afw.image.ExposureSummaryStats`` 

791 Equivalent legacy struct. 

792 """ 

793 for field in dataclasses.fields(legacy_summary_stats): 

794 a = getattr(legacy_summary_stats, field.name) 

795 b = getattr(summary_stats, field.name) 

796 if isinstance(b, tuple): 

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

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

799 else: 

800 tc.assertTrue(a == b or (math.isnan(a) and math.isnan(b)), f"{field.name}: {a} != {b}") 

801 

802 

803def compare_projection_to_legacy_wcs[F: Frame]( 

804 tc: unittest.TestCase, 

805 projection: Projection[F], 

806 legacy_wcs: Any, 

807 pixel_frame: F, 

808 subimage_bbox: Box, 

809 is_fits: bool = False, 

810) -> None: 

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

812 `lsst.afw.geom.SkyWcs`. 

813 

814 Parameters 

815 ---------- 

816 tc 

817 Test case object with assert methods to use. 

818 projection 

819 Projection to test. 

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

821 Equivalent legacy WCS. 

822 pixel_frame 

823 Expected pixel frame for the projection. 

824 subimage_bbox 

825 Bounding box of points to generate for tests. 

826 is_fits 

827 Whether this projection is expected to be exactly representable as a 

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

829 attached instead. 

830 """ 

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

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

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

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

835 sky_coords = legacy_coords_to_astropy( 

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

837 ) 

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

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

840 # change. 

841 check_projection(tc, projection, pixel_xy, sky_coords, pixel_frame) 

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

843 # array indices. 

844 check_astropy_wcs_interface( 

845 tc, projection.as_astropy(subimage_bbox), subimage_array_xy, sky_coords, pixel_atol=1e-5 

846 ) 

847 if is_fits: 

848 fits_wcs = projection.as_fits_wcs(subimage_bbox, allow_approximation=True) 

849 assert fits_wcs is not None 

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

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

852 # Projection from a FITS WCS, too. 

853 fits_projection = Projection.from_fits_wcs(fits_wcs, pixel_frame) 

854 check_projection( 

855 tc, 

856 fits_projection, 

857 subimage_array_xy, 

858 sky_coords, 

859 pixel_frame, 

860 pixel_atol=1e-5, 

861 ) 

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

863 # AST FrameSet so we can convert them into legacy 

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

865 tc.assertIn("Begin FrameSet", fits_projection.show()) 

866 else: 

867 tc.assertIsNone(projection.as_fits_wcs(subimage_bbox, allow_approximation=False)) 

868 # The legacy SkyWcs should instead have a FITS approximation 

869 # attached; run the same tests on that. 

870 assert projection.fits_approximation is not None 

871 compare_projection_to_legacy_wcs( 

872 tc, 

873 projection.fits_approximation, 

874 legacy_wcs.getFitsApproximation(), 

875 pixel_frame, 

876 subimage_bbox, 

877 is_fits=True, 

878 ) 

879 

880 

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

882 tc: unittest.TestCase, 

883 transform: Transform[I, O], 

884 input_xy: XY[np.ndarray], 

885 output_xy: XY[np.ndarray], 

886 in_frame: Frame, 

887 out_frame: Frame, 

888 *, 

889 check_inverted: bool = True, 

890 in_atol: u.Quantity | None = None, 

891 out_atol: u.Quantity | None = None, 

892) -> None: 

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

894 

895 Parameters 

896 ---------- 

897 tc 

898 Test case object with assert methods to use. 

899 transform 

900 Transform to test. 

901 input_xy 

902 Arrays of input points. 

903 output_xy 

904 Arrays of output points. 

905 in_frame 

906 Expected input frame. 

907 out_frame 

908 Expected output frame. 

909 check_inverted 

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

911 in_atol 

912 Expected absolute precision of input points. 

913 out_atol 

914 Expected absolute precision of output points. 

915 """ 

916 tc.assertEqual(transform.in_frame, in_frame) 

917 tc.assertEqual(transform.out_frame, out_frame) 

918 in_atol_v = in_atol.to_value(in_frame.unit) if in_atol is not None else None 

919 out_atol_v = out_atol.to_value(out_frame.unit) if out_atol is not None else None 

920 # Test array interfaces. 

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

922 assert_close(tc, test_output_xy.x, output_xy.x, atol=out_atol_v) 

923 assert_close(tc, test_output_xy.y, output_xy.y, atol=out_atol_v) 

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

925 assert_close(tc, test_input_xy.x, input_xy.x, atol=in_atol_v) 

926 assert_close(tc, test_input_xy.y, input_xy.y, atol=in_atol_v) 

927 # Test scalar interfaces with numpy scalars. 

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

929 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).x, output_x, atol=out_atol_v) 

930 assert_close(tc, transform.apply_forward(x=input_x, y=input_y).y, output_y, atol=out_atol_v) 

931 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).x, input_x, atol=in_atol_v) 

932 assert_close(tc, transform.apply_inverse(x=output_x, y=output_y).y, input_y, atol=in_atol_v) 

933 # Test quantity array interfaces. 

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

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

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

937 assert_close(tc, test_output_q_xy.x, output_q_xy.x, atol=out_atol) 

938 assert_close(tc, test_output_q_xy.y, output_q_xy.y, atol=out_atol) 

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

940 assert_close(tc, test_input_q_xy.x, input_q_xy.x, atol=in_atol) 

941 assert_close(tc, test_input_q_xy.y, input_q_xy.y, atol=in_atol) 

942 # Test quantity scalar interfaces. 

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

944 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y 

945 ): 

946 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).x, output_q_x, atol=out_atol) 

947 assert_close(tc, transform.apply_forward_q(x=input_q_x, y=input_q_y).y, output_q_y, atol=out_atol) 

948 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).x, input_q_x, atol=in_atol) 

949 assert_close(tc, transform.apply_inverse_q(x=output_q_x, y=output_q_y).y, input_q_y, atol=in_atol) 

950 if check_inverted: 

951 # Test the inverse transform. 

952 check_transform( 

953 tc, 

954 transform.inverted(), 

955 output_xy, 

956 input_xy, 

957 out_frame, 

958 in_frame, 

959 check_inverted=False, 

960 out_atol=in_atol, 

961 in_atol=out_atol, 

962 ) 

963 

964 

965def check_projection[P: Frame]( 

966 tc: unittest.TestCase, 

967 projection: Projection[P], 

968 pixel_xy: XY[np.ndarray], 

969 sky_coords: SkyCoord, 

970 pixel_frame: Frame, 

971 *, 

972 pixel_atol: float | None = None, 

973 sky_atol: u.Quantity | None = None, 

974) -> None: 

975 """Test a `.Projection` instance against known arrays of pixel and sky 

976 coordinates. 

977 

978 Parameters 

979 ---------- 

980 tc 

981 Test case object with assert methods to use. 

982 projection 

983 Projection to test. 

984 pixel_xy 

985 Arrays of pixel coordinates. 

986 sky_coords 

987 Corresponding sky coordinates. 

988 pixel_frame 

989 Expected pixel frame. 

990 pixel_atol 

991 Expected absolute precision of pixel points. 

992 sky_atol 

993 Expected absolute precision of sky coordinates. 

994 """ 

995 tc.assertEqual(projection.pixel_frame, pixel_frame) 

996 tc.assertEqual(projection.sky_frame, SkyFrame.ICRS) 

997 sky_atol_v = sky_atol.to_value(SkyFrame.ICRS.unit) if sky_atol is not None else None 

998 pixel_atol_q = pixel_atol * u.pix if pixel_atol is not None else None 

999 # Test array interfaces. 

1000 test_pixel_xy = cast(XY[np.ndarray], projection.sky_to_pixel(sky_coords)) 

1001 assert_close(tc, test_pixel_xy.x, pixel_xy.x, atol=pixel_atol) 

1002 assert_close(tc, test_pixel_xy.y, pixel_xy.y, atol=pixel_atol) 

1003 test_sky_astropy = projection.pixel_to_sky(x=pixel_xy.x, y=pixel_xy.y) 

1004 assert_close(tc, test_sky_astropy.ra, sky_coords.ra, atol=sky_atol_v) 

1005 assert_close(tc, test_sky_astropy.dec, sky_coords.dec, atol=sky_atol_v) 

1006 # Test scalar interfaces. 

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

1008 assert_close(tc, projection.sky_to_pixel(sky_single).x, pixel_x, atol=pixel_atol) 

1009 assert_close(tc, projection.sky_to_pixel(sky_single).y, pixel_y, atol=pixel_atol) 

1010 assert_close(tc, projection.pixel_to_sky(x=pixel_x, y=pixel_y).ra, sky_single.ra, atol=sky_atol_v) 

1011 assert_close(tc, projection.pixel_to_sky(x=pixel_x, y=pixel_y).dec, sky_single.dec, atol=sky_atol_v) 

1012 # Test the underlying Transform object. 

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

1014 check_transform( 

1015 tc, 

1016 projection.pixel_to_sky_transform, 

1017 pixel_xy, 

1018 sky_xy, 

1019 pixel_frame, 

1020 SkyFrame.ICRS, 

1021 check_inverted=False, 

1022 in_atol=pixel_atol_q, 

1023 out_atol=sky_atol, 

1024 ) 

1025 check_transform( 

1026 tc, 

1027 projection.sky_to_pixel_transform, 

1028 sky_xy, 

1029 pixel_xy, 

1030 SkyFrame.ICRS, 

1031 pixel_frame, 

1032 check_inverted=False, 

1033 in_atol=sky_atol, 

1034 out_atol=pixel_atol_q, 

1035 ) 

1036 # Test the Astropy interface adapter. 

1037 check_astropy_wcs_interface( 

1038 tc, projection.as_astropy(), pixel_xy, sky_coords, pixel_atol=pixel_atol, sky_atol=sky_atol 

1039 ) 

1040 

1041 

1042def assert_projections_equal( 

1043 tc: unittest.TestCase, 

1044 a: Projection[Any] | None, 

1045 b: Projection[Any] | None, 

1046 expect_identity: bool | None = None, 

1047) -> None: 

1048 """Test that two `.Projection` instances are equivalent.""" 

1049 if a is None and b is None: 

1050 return 

1051 assert a is not None and b is not None 

1052 match expect_identity: 

1053 case True: 

1054 tc.assertIs(a, b) 

1055 return 

1056 case False: 

1057 tc.assertIsNot(a, b) 

1058 case None if a is b: 

1059 return 

1060 tc.assertEqual(a.pixel_frame, b.pixel_frame) 

1061 tc.assertEqual(a.show(simplified=True), b.show(simplified=True)) 

1062 assert_projections_equal( 

1063 tc, a.fits_approximation, cast(Projection[Any], b.fits_approximation), expect_identity=False 

1064 ) 

1065 

1066 

1067def check_astropy_wcs_interface( 

1068 tc: unittest.TestCase, 

1069 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS, 

1070 pixel_xy: XY[np.ndarray], 

1071 sky_coords: SkyCoord, 

1072 *, 

1073 pixel_atol: float | None = None, 

1074 sky_atol: u.Quantity | None = None, 

1075) -> None: 

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

1077 sky coordinates. 

1078 

1079 Parameters 

1080 ---------- 

1081 tc 

1082 Test case object with assert methods to use. 

1083 wcs 

1084 Astropy WCS object to test. 

1085 pixel_xy 

1086 Arrays of pixel coordinates. 

1087 sky_coords 

1088 Corresponding sky coordinates. 

1089 pixel_atol 

1090 Expected absolute precision of pixel points. 

1091 sky_atol 

1092 Expected absolute precision of sky coordinates. 

1093 """ 

1094 test_x, test_y = wcs.world_to_pixel(sky_coords) 

1095 assert_close(tc, test_x, pixel_xy.x, atol=pixel_atol) 

1096 assert_close(tc, test_y, pixel_xy.y, atol=pixel_atol) 

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

1098 assert_close(tc, test_sky_coords.ra, sky_coords.ra, atol=sky_atol) 

1099 assert_close(tc, test_sky_coords.dec, sky_coords.dec, atol=sky_atol) 

1100 

1101 

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

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

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

1105 

1106 

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

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

1109 coordinate object. 

1110 """ 

1111 return SkyCoord( 

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

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

1114 ) 

1115 

1116 

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

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

1119 from lsst.geom import Point2D 

1120 

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

1122 

1123 

1124def compare_amplifier_to_legacy( 

1125 tc: unittest.TestCase, 

1126 amplifier: Amplifier, 

1127 legacy_amplifier: Any, 

1128 *, 

1129 is_raw_assembled: bool, 

1130 expect_nominal_calibrations: bool = True, 

1131) -> None: 

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

1133 `lsst.afw.cameraGeom.Amplifier`. 

1134 """ 

1135 tc.assertEqual(legacy_amplifier.getName(), amplifier.name) 

1136 tc.assertEqual(Box.from_legacy(legacy_amplifier.getBBox()), amplifier.bbox) 

1137 if is_raw_assembled: 

1138 raw_geom = amplifier.assembled_raw_geometry 

1139 else: 

1140 raw_geom = amplifier.unassembled_raw_geometry 

1141 assert raw_geom is not None 

1142 tc.assertEqual(ReadoutCorner.from_legacy(legacy_amplifier.getReadoutCorner()), raw_geom.readout_corner) 

1143 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawBBox()), raw_geom.bbox) 

1144 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawDataBBox()), raw_geom.data_bbox) 

1145 tc.assertEqual(legacy_amplifier.getRawFlipX(), raw_geom.flip_x) 

1146 tc.assertEqual(legacy_amplifier.getRawFlipY(), raw_geom.flip_y) 

1147 tc.assertEqual(legacy_amplifier.getRawXYOffset().getX(), raw_geom.x_offset) 

1148 tc.assertEqual(legacy_amplifier.getRawXYOffset().getY(), raw_geom.y_offset) 

1149 tc.assertEqual( 

1150 Box.from_legacy(legacy_amplifier.getRawHorizontalOverscanBBox()), raw_geom.horizontal_overscan_bbox 

1151 ) 

1152 tc.assertEqual( 

1153 Box.from_legacy(legacy_amplifier.getRawVerticalOverscanBBox()), raw_geom.vertical_overscan_bbox 

1154 ) 

1155 tc.assertEqual(Box.from_legacy(legacy_amplifier.getRawPrescanBBox()), raw_geom.horizontal_prescan_bbox) 

1156 if expect_nominal_calibrations: 

1157 assert amplifier.nominal_calibrations is not None 

1158 assert_equal_allow_nan(tc, legacy_amplifier.getGain(), amplifier.nominal_calibrations.gain) 

1159 assert_equal_allow_nan(tc, legacy_amplifier.getReadNoise(), amplifier.nominal_calibrations.read_noise) 

1160 assert_equal_allow_nan( 

1161 tc, legacy_amplifier.getSaturation(), amplifier.nominal_calibrations.saturation 

1162 ) 

1163 assert_equal_allow_nan( 

1164 tc, legacy_amplifier.getSuspectLevel(), amplifier.nominal_calibrations.suspect_level 

1165 ) 

1166 np.testing.assert_array_equal( 

1167 legacy_amplifier.getLinearityCoeffs(), amplifier.nominal_calibrations.linearity_coefficients 

1168 ) 

1169 tc.assertEqual(legacy_amplifier.getLinearityType(), amplifier.nominal_calibrations.linearity_type) 

1170 

1171 

1172def compare_detector_to_legacy( 

1173 tc: unittest.TestCase, 

1174 detector: Detector, 

1175 legacy_detector: Any, 

1176 *, 

1177 is_raw_assembled: bool, 

1178 expect_nominal_calibrations: bool = True, 

1179) -> None: 

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

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

1182 

1183 tc.assertEqual(legacy_detector.getName(), detector.name) 

1184 tc.assertEqual(legacy_detector.getId(), detector.id) 

1185 tc.assertEqual(DetectorType.from_legacy(legacy_detector.getType()), detector.type) 

1186 tc.assertEqual(Box.from_legacy(legacy_detector.getBBox()), detector.bbox) 

1187 tc.assertEqual(legacy_detector.getSerial(), detector.serial) 

1188 legacy_orientation = legacy_detector.getOrientation() 

1189 tc.assertEqual(legacy_orientation.getFpPosition3().getX(), detector.orientation.focal_plane_x) 

1190 tc.assertEqual(legacy_orientation.getFpPosition3().getY(), detector.orientation.focal_plane_y) 

1191 tc.assertEqual(legacy_orientation.getFpPosition3().getZ(), detector.orientation.focal_plane_z) 

1192 tc.assertEqual(legacy_orientation.getReferencePoint().getX(), detector.orientation.pixel_reference_x) 

1193 tc.assertEqual(legacy_orientation.getReferencePoint().getY(), detector.orientation.pixel_reference_y) 

1194 tc.assertEqual(legacy_orientation.getYaw().asRadians(), detector.orientation.yaw.to_value(u.rad)) 

1195 tc.assertEqual(legacy_orientation.getPitch().asRadians(), detector.orientation.pitch.to_value(u.rad)) 

1196 tc.assertEqual(legacy_orientation.getRoll().asRadians(), detector.orientation.roll.to_value(u.rad)) 

1197 tc.assertEqual(legacy_detector.getPixelSize().getX(), detector.pixel_size) 

1198 tc.assertEqual(legacy_detector.getPhysicalType(), detector.physical_type) 

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

1200 compare_amplifier_to_legacy( 

1201 tc, 

1202 amplifier, 

1203 legacy_amplifier, 

1204 is_raw_assembled=is_raw_assembled, 

1205 expect_nominal_calibrations=expect_nominal_calibrations, 

1206 ) 

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

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

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

1210 check_transform( 

1211 tc, 

1212 detector.to_focal_plane, 

1213 pixel_xy, 

1214 legacy_points_to_xy_array(fp_legacy_points), 

1215 detector.frame, 

1216 detector.to_focal_plane.out_frame, 

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

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

1219 ) 

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

1221 check_transform( 

1222 tc, 

1223 detector.to_field_angle, 

1224 pixel_xy, 

1225 legacy_points_to_xy_array(fa_legacy_points), 

1226 detector.frame, 

1227 detector.to_field_angle.out_frame, 

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

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

1230 )