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

191 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-04-02 01:01 -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 "arrays_to_legacy_points", 

16 "assert_close", 

17 "assert_images_equal", 

18 "assert_masked_images_equal", 

19 "assert_masks_equal", 

20 "assert_projections_equal", 

21 "check_astropy_wcs_interface", 

22 "check_projection", 

23 "check_transform", 

24 "compare_image_to_legacy", 

25 "compare_mask_to_legacy", 

26 "compare_masked_image_to_legacy", 

27 "compare_projection_to_legacy_wcs", 

28 "compare_psf_to_legacy", 

29 "compare_visit_image_to_legacy", 

30 "legacy_coords_to_astropy", 

31 "legacy_points_to_xy_array", 

32) 

33 

34import unittest 

35from collections.abc import Mapping 

36from typing import Any, Literal, cast 

37 

38import astropy.units as u 

39import astropy.wcs.wcsapi 

40import numpy as np 

41from astropy.coordinates import SkyCoord 

42 

43from .._geom import XY, Box 

44from .._image import Image 

45from .._mask import Mask, MaskPlane 

46from .._masked_image import MaskedImage 

47from .._transforms import DetectorFrame, Frame, Projection, SkyFrame, Transform 

48from .._visit_image import VisitImage 

49from ..psfs import PointSpreadFunction 

50 

51 

52def assert_close( 

53 tc: unittest.TestCase, 

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

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

56 **kwargs: Any, 

57) -> None: 

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

59 

60 Parameters 

61 ---------- 

62 tc 

63 Test case object with assert methods to use. 

64 a 

65 Array, float, or quantity to compare. 

66 b 

67 Array, float, or quantity to compare. 

68 **kwargs 

69 Forwarded to `astropy.units.allclose`. 

70 """ 

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

72 

73 

74def assert_images_equal( 

75 tc: unittest.TestCase, 

76 a: Image, 

77 b: Image, 

78 *, 

79 rtol: float = 0.0, 

80 atol: float = 0.0, 

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

82) -> None: 

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

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

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

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

87 if expect_view is not None: 

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

89 if expect_view == "array": 

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

91 else: 

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

93 if not expect_view: 

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

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

96 

97 

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

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

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

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

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

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

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

105 

106 

107def assert_masked_images_equal( 

108 tc: unittest.TestCase, 

109 a: MaskedImage, 

110 b: MaskedImage, 

111 *, 

112 rtol: float = 0.0, 

113 atol: float = 0.0, 

114 expect_view: bool | None = None, 

115) -> None: 

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

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

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

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

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

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

122 

123 

124def compare_image_to_legacy( 

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

126) -> None: 

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

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

129 if expect_view is not None: 

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

131 if not expect_view: 

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

133 

134 

135def compare_mask_to_legacy( 

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

137) -> None: 

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

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

140 if plane_map is None: 

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

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

143 np.testing.assert_array_equal( 

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

145 mask.get(new_plane.name), 

146 ) 

147 

148 

149def compare_masked_image_to_legacy( 

150 tc: unittest.TestCase, 

151 masked_image: MaskedImage, 

152 legacy_masked_image: Any, 

153 *, 

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

155 expect_view: bool | None = None, 

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

157) -> None: 

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

159 object. 

160 

161 Parameters 

162 ---------- 

163 tc 

164 Test case to use for asserts. 

165 masked_image 

166 New image to test. 

167 legacy_masked_image 

168 Legacy image to test against. 

169 plane_map 

170 Mapping between new and legacy mask planes. 

171 expect_view 

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

173 memory. 

174 alternates 

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

176 check against the legacy versions of those components. 

177 """ 

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

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

180 compare_image_to_legacy( 

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

182 ) 

183 if alternates: 

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

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

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

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

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

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

190 

191 

192def compare_visit_image_to_legacy( 

193 tc: unittest.TestCase, 

194 visit_image: VisitImage, 

195 legacy_exposure: Any, 

196 *, 

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

198 expect_view: bool | None = None, 

199 instrument: str, 

200 visit: int, 

201 detector: int, 

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

203) -> None: 

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

205 object. 

206 

207 Parameters 

208 ---------- 

209 tc 

210 Test case to use for asserts. 

211 visit_image 

212 New image to test. 

213 legacy_exposure 

214 Legacy image to test against. 

215 plane_map 

216 Mapping between new and legacy mask planes. 

217 expect_view 

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

219 memory. 

220 instrument 

221 Expected instrument name. 

222 visit 

223 Expected visit ID. 

224 detetector 

225 Expected detector ID. 

226 alternates 

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

228 check against the legacy versions of those components. 

229 """ 

230 compare_masked_image_to_legacy( 

231 tc, 

232 visit_image, 

233 legacy_exposure.getMaskedImage(), 

234 plane_map=plane_map, 

235 expect_view=expect_view, 

236 alternates=alternates, 

237 ) 

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

239 compare_projection_to_legacy_wcs( 

240 tc, 

241 visit_image.projection, 

242 legacy_exposure.getWcs(), 

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

244 visit_image.bbox, 

245 ) 

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

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

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

249 if alternates: 

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

251 compare_projection_to_legacy_wcs( 

252 tc, 

253 projection, 

254 legacy_exposure.getWcs(), 

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

256 visit_image.bbox, 

257 ) 

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

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

260 

261 

262def compare_psf_to_legacy( 

263 tc: unittest.TestCase, psf: PointSpreadFunction, legacy_psf: Any, subimage_bbox: Box | None = None 

264) -> int: 

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

266 

267 Parameters 

268 ---------- 

269 tc 

270 Test case object with assert methods to use. 

271 psf 

272 Point-spread function to test. 

273 legacy_psf 

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

275 subimage_bbox 

276 Bounding box to draw test points from. Defaults to ``psf.bounds``. 

277 

278 Returns 

279 ------- 

280 `int` 

281 The number of points actually tested. 

282 """ 

283 if subimage_bbox is None: 

284 if isinstance(psf.bounds, Box): 

285 subimage_bbox = psf.bounds 

286 elif hasattr(psf.bounds, "bbox"): 

287 subimage_bbox = cast(Box, psf.bounds.bbox) 

288 else: 

289 raise TypeError("No usable PSF bounds found; 'subimage_bbox' must be provided.") 

290 

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

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

293 legacy_points = arrays_to_legacy_points(pixel_xy.x, pixel_xy.y) 

294 

295 n_points_tested: int = 0 

296 for p in legacy_points: 

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

298 continue 

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

300 tc.assertEqual( 

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

302 ) 

303 tc.assertEqual( 

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

305 ) 

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

307 n_points_tested += 1 

308 return n_points_tested 

309 

310 

311def compare_projection_to_legacy_wcs[F: Frame]( 

312 tc: unittest.TestCase, 

313 projection: Projection[F], 

314 legacy_wcs: Any, 

315 pixel_frame: F, 

316 subimage_bbox: Box, 

317 is_fits: bool = False, 

318) -> None: 

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

320 `lsst.afw.geom.SkyWcs`. 

321 

322 Parameters 

323 ---------- 

324 tc 

325 Test case object with assert methods to use. 

326 projection 

327 Projection to test. 

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

329 Equivalent legacy WCS. 

330 pixel_frame 

331 Expected pixel frame for the projection. 

332 subimage_bbox 

333 Bounding box of points to generate for tests. 

334 is_fits 

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

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

337 attached instead. 

338 """ 

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

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

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

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

343 sky_coords = legacy_coords_to_astropy( 

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

345 ) 

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

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

348 # change. 

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

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

351 # array indices. 

352 check_astropy_wcs_interface( 

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

354 ) 

355 if is_fits: 

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

357 assert fits_wcs is not None 

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

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

360 # Projection from a FITS WCS, too. 

361 fits_projection = Projection.from_fits_wcs(fits_wcs, pixel_frame) 

362 check_projection( 

363 tc, 

364 fits_projection, 

365 subimage_array_xy, 

366 sky_coords, 

367 pixel_frame, 

368 pixel_atol=1e-5, 

369 ) 

370 else: 

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

372 # The legacy SkyWcs should instead have a FITS approximation 

373 # attached; run the same tests on that. 

374 assert projection.fits_approximation is not None 

375 compare_projection_to_legacy_wcs( 

376 tc, 

377 projection.fits_approximation, 

378 legacy_wcs.getFitsApproximation(), 

379 pixel_frame, 

380 subimage_bbox, 

381 is_fits=True, 

382 ) 

383 

384 

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

386 tc: unittest.TestCase, 

387 transform: Transform[I, O], 

388 input_xy: XY[np.ndarray], 

389 output_xy: XY[np.ndarray], 

390 in_frame: Frame, 

391 out_frame: Frame, 

392 *, 

393 check_inverted: bool = True, 

394 in_atol: u.Quantity | None = None, 

395 out_atol: u.Quantity | None = None, 

396) -> None: 

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

398 

399 Parameters 

400 ---------- 

401 tc 

402 Test case object with assert methods to use. 

403 transform 

404 Transform to test. 

405 input_xy 

406 Arrays of input points. 

407 outout_xy 

408 Arrays of output points. 

409 in_frame 

410 Expected input frame. 

411 out_frame 

412 Expected output frame. 

413 check_inverted 

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

415 in_atol 

416 Expected absolute precision of input points. 

417 out_atol 

418 Expected absolute precision of output points. 

419 """ 

420 tc.assertEqual(transform.in_frame, in_frame) 

421 tc.assertEqual(transform.out_frame, out_frame) 

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

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

424 # Test array interfaces. 

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

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

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

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

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

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

431 # Test scalar interfaces with numpy scalars. 

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

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

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

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

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

437 # Test quantity array interfaces. 

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

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

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

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

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

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

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

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

446 # Test quantity scalar interfaces. 

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

448 input_q_xy.x, input_q_xy.y, output_q_xy.x, output_q_xy.y 

449 ): 

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

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

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

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

454 if check_inverted: 

455 # Test the inverse transform. 

456 check_transform( 

457 tc, 

458 transform.inverted(), 

459 output_xy, 

460 input_xy, 

461 out_frame, 

462 in_frame, 

463 check_inverted=False, 

464 out_atol=in_atol, 

465 in_atol=out_atol, 

466 ) 

467 

468 

469def check_projection[P: Frame]( 

470 tc: unittest.TestCase, 

471 projection: Projection[P], 

472 pixel_xy: XY[np.ndarray], 

473 sky_coords: SkyCoord, 

474 pixel_frame: Frame, 

475 *, 

476 pixel_atol: float | None = None, 

477 sky_atol: u.Quantity | None = None, 

478) -> None: 

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

480 coordinates. 

481 

482 Parameters 

483 ---------- 

484 tc 

485 Test case object with assert methods to use. 

486 projection 

487 Projection to test. 

488 pixel_xy 

489 Arrays of pixel coordinates. 

490 sky_coords 

491 Corresponding sky coordinates. 

492 pixel_frame 

493 Expected pixel frame. 

494 pixel_atol 

495 Expected absolute precision of pixel points. 

496 sky_atol 

497 Expected absolute precision of sky coordinates. 

498 """ 

499 tc.assertEqual(projection.pixel_frame, pixel_frame) 

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

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

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

503 # Test array interfaces. 

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

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

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

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

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

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

510 # Test scalar interfaces. 

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

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

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

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

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

516 # Test the underlying Transform object. 

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

518 check_transform( 

519 tc, 

520 projection.pixel_to_sky_transform, 

521 pixel_xy, 

522 sky_xy, 

523 pixel_frame, 

524 SkyFrame.ICRS, 

525 check_inverted=False, 

526 in_atol=pixel_atol_q, 

527 out_atol=sky_atol, 

528 ) 

529 check_transform( 

530 tc, 

531 projection.sky_to_pixel_transform, 

532 sky_xy, 

533 pixel_xy, 

534 SkyFrame.ICRS, 

535 pixel_frame, 

536 check_inverted=False, 

537 in_atol=sky_atol, 

538 out_atol=pixel_atol_q, 

539 ) 

540 # Test the Astropy interface adapter. 

541 check_astropy_wcs_interface( 

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

543 ) 

544 

545 

546def assert_projections_equal( 

547 tc: unittest.TestCase, 

548 a: Projection[Any] | None, 

549 b: Projection[Any] | None, 

550 expect_identity: bool | None = None, 

551) -> None: 

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

553 if a is None and b is None: 

554 return 

555 assert a is not None and b is not None 

556 match expect_identity: 

557 case True: 

558 tc.assertIs(a, b) 

559 return 

560 case False: 

561 tc.assertIsNot(a, b) 

562 case None if a is b: 

563 return 

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

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

566 assert_projections_equal( 

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

568 ) 

569 

570 

571def check_astropy_wcs_interface( 

572 tc: unittest.TestCase, 

573 wcs: astropy.wcs.wcsapi.BaseHighLevelWCS, 

574 pixel_xy: XY[np.ndarray], 

575 sky_coords: SkyCoord, 

576 *, 

577 pixel_atol: float | None = None, 

578 sky_atol: u.Quantity | None = None, 

579) -> None: 

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

581 sky coordinates. 

582 

583 Parameters 

584 ---------- 

585 tc 

586 Test case object with assert methods to use. 

587 wcs 

588 Astropy WCS object to test. 

589 pixel_xy 

590 Arrays of pixel coordinates. 

591 sky_coords 

592 Corresponding sky coordinates. 

593 pixel_atol 

594 Expected absolute precision of pixel points. 

595 sky_atol 

596 Expected absolute precision of sky coordinates. 

597 """ 

598 test_x, test_y = wcs.world_to_pixel(sky_coords) 

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

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

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

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

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

604 

605 

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

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

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

609 

610 

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

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

613 coordinate object. 

614 """ 

615 return SkyCoord( 

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

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

618 ) 

619 

620 

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

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

623 from lsst.geom import Point2D 

624 

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