Coverage for tests/test_visit_image.py: 55%

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 

14import dataclasses 

15import os 

16import warnings 

17from typing import Any, Literal 

18 

19import astropy.io.fits 

20import astropy.units as u 

21import astropy.wcs 

22import numpy as np 

23import pytest 

24from astro_metadata_translator import ObservationInfo 

25 

26from lsst.images import ( 

27 BackgroundMap, 

28 Box, 

29 DetectorFrame, 

30 DifferenceImage, 

31 Image, 

32 MaskPlane, 

33 MaskSchema, 

34 ObservationSummaryStats, 

35 Polygon, 

36 SkyProjectionAstropyView, 

37 TractFrame, 

38 VisitImage, 

39 get_legacy_difference_image_mask_planes, 

40 get_legacy_visit_image_mask_planes, 

41) 

42from lsst.images.aperture_corrections import ApertureCorrectionMap, aperture_corrections_to_legacy 

43from lsst.images.cameras import Detector 

44from lsst.images.fields import ChebyshevField, SplineField, SumField, field_from_legacy_photo_calib 

45from lsst.images.fits import ExtensionKey, FitsOpaqueMetadata 

46from lsst.images.psfs import GaussianPointSpreadFunction, PointSpreadFunction 

47from lsst.images.serialization import read_archive 

48from lsst.images.tests import ( 

49 DP2_VISIT_DETECTOR_DATA_ID, 

50 RoundtripFits, 

51 RoundtripJson, 

52 RoundtripNdf, 

53 TemporaryButler, 

54 assert_close, 

55 assert_masked_images_equal, 

56 assert_sky_projections_equal, 

57 assert_visit_images_equal, 

58 compare_aperture_corrections_to_legacy, 

59 compare_detector_to_legacy, 

60 compare_photo_calib_to_legacy, 

61 compare_visit_image_to_legacy, 

62 make_random_sky_projection, 

63) 

64 

65try: 

66 import h5py # noqa: F401 

67 

68 HAVE_H5PY = True 

69except ImportError: 

70 HAVE_H5PY = False 

71 

72try: 

73 from lsst.afw.image import Exposure as LegacyExposure 

74 from lsst.afw.image import VisitInfo as LegacyVisitInfo 

75except ImportError: 

76 type LegacyExposure = Any # type: ignore[no-redef] 

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

78 

79EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None) 

80LOCAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "data") 

81 

82skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed") 

83 

84 

85@pytest.fixture(scope="session") 

86def visit_image_components() -> dict[str, Any]: 

87 """Return a dictionary of VisitImage components.""" 

88 rng = np.random.default_rng(500) 

89 det_frame = DetectorFrame(instrument="Inst", visit=1234, detector=1, bbox=Box.factory[1:4096, 1:4096]) 

90 mask_schema = MaskSchema([MaskPlane("M1", "D1")]) 

91 obs_info = ObservationInfo(instrument="LSSTCam", detector_num=4, physical_filter="r1") 

92 summary_stats = ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4) 

93 gaussian_psf = GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]) 

94 aperture_corrections: ApertureCorrectionMap = { 

95 "flux1": ChebyshevField(det_frame.bbox, np.array([0.75])), 

96 "flux2": ChebyshevField(det_frame.bbox, np.array([0.625])), 

97 } 

98 detector = read_archive(os.path.join(LOCAL_DATA_DIR, "detector.json"), Detector) 

99 image = Image(42, shape=(1024, 1024), unit=u.nJy) 

100 variance = Image(5.0, shape=(1024, 1024), unit=u.nJy * u.nJy) 

101 # polygon is the lower triangle of the image. 

102 polygon = Polygon(x_vertices=[-0.5, 1023.5, -0.5], y_vertices=[-0.5, -0.5, 1023.5]) 

103 sky_projection = make_random_sky_projection(rng, det_frame, det_frame.bbox) 

104 return { 

105 "mask_schema": mask_schema, 

106 "obs_info": obs_info, 

107 "summary_stats": summary_stats, 

108 "gaussian_psf": gaussian_psf, 

109 "aperture_corrections": aperture_corrections, 

110 "detector": detector, 

111 "image": image, 

112 "variance": variance, 

113 "polygon": polygon, 

114 "sky_projection": sky_projection, 

115 } 

116 

117 

118def make_visit_image(components: dict[str, Any]) -> VisitImage: 

119 """Construct a new VisitImage with most components populated.""" 

120 det_frame = components["sky_projection"].pixel_frame 

121 opaque = FitsOpaqueMetadata() 

122 hdr = astropy.io.fits.Header() 

123 with warnings.catch_warnings(): 

124 # Silence warnings about long keys becoming HIERARCH. 

125 warnings.simplefilter("ignore", category=astropy.io.fits.verify.VerifyWarning) 

126 hdr.update({"PLATFORM": "lsstcam", "LSST BUTLER ID": "123456789"}) 

127 opaque.extract_legacy_primary_header(hdr) 

128 # API signature suggests sky_projection and obs_info can be None but 

129 # they are required (unless you pass them in via the image plane). 

130 vi = VisitImage( 

131 components["image"], 

132 variance=components["variance"], 

133 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]), 

134 mask_schema=components["mask_schema"], 

135 sky_projection=components["sky_projection"], 

136 obs_info=components["obs_info"], 

137 summary_stats=components["summary_stats"], 

138 detector=components["detector"], 

139 bounds=components["polygon"], 

140 aperture_corrections=components["aperture_corrections"], 

141 band="r", 

142 ) 

143 vi.backgrounds.add( 

144 "standard", 

145 ChebyshevField(det_frame.bbox, np.array([[2.0]])), 

146 description="Background subtracted from the image.", 

147 is_subtracted=True, 

148 ) 

149 vi._opaque_metadata = opaque 

150 return vi 

151 

152 

153def make_simplest_visit_image(components: dict[str, Any]) -> VisitImage: 

154 """Construct a VisitImage with the minimal set of components populated.""" 

155 return VisitImage( 

156 components["image"], 

157 psf=GaussianPointSpreadFunction(2.5, stamp_size=33, bounds=Box.factory[-10:10, -12:13]), 

158 mask_schema=components["mask_schema"], 

159 sky_projection=components["sky_projection"], 

160 detector=components["detector"], 

161 obs_info=components["obs_info"], 

162 band="r", 

163 ) 

164 

165 

166def _make_sum_background_visit_image(components: dict[str, Any], visit_image: VisitImage) -> VisitImage: 

167 """Return a VisitImage whose subtracted background is a SumField.""" 

168 rng = np.random.default_rng(42) 

169 bbox = visit_image.image.sky_projection.pixel_frame.bbox 

170 bin_y = bbox.y.linspace(6) 

171 bin_x = bbox.x.linspace(7) 

172 spline_a = SplineField( 

173 bbox, 

174 rng.standard_normal(size=(bin_y.size, bin_x.size)), 

175 y=bin_y, 

176 x=bin_x, 

177 ) 

178 spline_b = SplineField( 

179 bbox, 

180 rng.standard_normal(size=(bin_y.size, bin_x.size)), 

181 y=bin_y, 

182 x=bin_x, 

183 ) 

184 sum_field = SumField([spline_a, spline_b]) 

185 bg_map = BackgroundMap() 

186 bg_map.add( 

187 "stacked", 

188 sum_field, 

189 description="Two-operand SumField subtracted background.", 

190 is_subtracted=True, 

191 ) 

192 return VisitImage( 

193 components["image"], 

194 variance=components["variance"], 

195 psf=components["gaussian_psf"], 

196 mask_schema=components["mask_schema"], 

197 sky_projection=components["sky_projection"], 

198 obs_info=components["obs_info"], 

199 summary_stats=components["summary_stats"], 

200 detector=components["detector"], 

201 band="r", 

202 backgrounds=bg_map, 

203 ) 

204 

205 

206def _check_sum_background_round_trip(result: VisitImage, original: VisitImage) -> None: 

207 """Assert that a round-tripped SumField background matches the original.""" 

208 subtracted = result.backgrounds.subtracted 

209 assert subtracted is not None 

210 assert isinstance(subtracted.field, SumField) 

211 original_subtracted = original.backgrounds.subtracted 

212 assert original_subtracted is not None 

213 original_field = original_subtracted.field 

214 assert isinstance(original_field, SumField) 

215 round_field = subtracted.field 

216 assert isinstance(round_field, SumField) 

217 assert len(round_field.operands) == len(original_field.operands) 

218 for round_op, orig_op in zip(round_field.operands, original_field.operands, strict=True): 

219 assert round_op == orig_op 

220 

221 

222def test_basics(visit_image_components: dict[str, Any]) -> None: 

223 """Verify VisitImage constructor patterns and required-argument checks.""" 

224 c = visit_image_components 

225 # Test default fill of variance. 

226 visit = make_simplest_visit_image(c) 

227 assert visit.variance.array[0, 0] == 1.0 

228 assert visit[...] is not visit 

229 assert str(visit) == "VisitImage(Image([y=0:1024, x=0:1024], int64), ['M1'])" 

230 assert repr(visit) == ( 

231 "VisitImage(Image(..., bbox=Box(y=Interval(start=0, stop=1024), x=Interval(start=0, stop=1024))," 

232 " dtype=dtype('int64')), mask_schema=MaskSchema([MaskPlane(name='M1', description='D1')]," 

233 " dtype=dtype('uint8')))" 

234 ) 

235 

236 astropy_wcs = visit.astropy_wcs 

237 assert isinstance(astropy_wcs, SkyProjectionAstropyView) 

238 approx_wcs = visit.fits_wcs 

239 assert isinstance(approx_wcs, astropy.wcs.WCS) 

240 

241 with pytest.raises(TypeError): 

242 # Requires a PSF. 

243 VisitImage( 

244 c["image"], 

245 mask_schema=c["mask_schema"], 

246 sky_projection=c["sky_projection"], 

247 obs_info=c["obs_info"], 

248 detector=c["detector"], 

249 band="r", 

250 ) 

251 

252 with pytest.raises(TypeError): 

253 # Requires ObservationInfo. 

254 VisitImage( 

255 c["image"], 

256 psf=c["gaussian_psf"], 

257 mask_schema=c["mask_schema"], 

258 sky_projection=c["sky_projection"], 

259 detector=c["detector"], 

260 band="r", 

261 ) 

262 

263 with pytest.raises(TypeError): 

264 # Requires a sky_projection. 

265 VisitImage( 

266 c["image"], 

267 psf=c["gaussian_psf"], 

268 mask_schema=c["mask_schema"], 

269 obs_info=c["obs_info"], 

270 detector=c["detector"], 

271 band="r", 

272 ) 

273 

274 with pytest.raises(TypeError): 

275 # Requires a detector. 

276 VisitImage( 

277 c["image"], 

278 psf=c["gaussian_psf"], 

279 mask_schema=c["mask_schema"], 

280 sky_projection=c["sky_projection"], 

281 obs_info=c["obs_info"], 

282 band="r", 

283 ) 

284 

285 with pytest.raises(TypeError): 

286 # Requires some form of mask. 

287 VisitImage( 

288 c["image"], 

289 psf=c["gaussian_psf"], 

290 sky_projection=c["sky_projection"], 

291 obs_info=c["obs_info"], 

292 detector=c["detector"], 

293 band="r", 

294 ) 

295 

296 with pytest.raises(TypeError): 

297 VisitImage( 

298 Image(42, shape=(5, 5)), 

299 psf=c["gaussian_psf"], 

300 mask_schema=c["mask_schema"], 

301 sky_projection=c["sky_projection"], 

302 obs_info=c["obs_info"], 

303 detector=c["detector"], 

304 band="r", 

305 ) 

306 

307 # Requires a DetectorFrame. 

308 rng = np.random.default_rng(501) 

309 tract_frame = TractFrame(skymap="Skymap", tract=1, bbox=Box.factory[1:10, 1:10]) 

310 tract_proj = make_random_sky_projection(rng, tract_frame, Box.factory[1:4096, 1:4096]) 

311 with pytest.raises(TypeError): 

312 VisitImage( 

313 c["image"], 

314 sky_projection=tract_proj, 

315 psf=c["gaussian_psf"], 

316 mask_schema=c["mask_schema"], 

317 obs_info=c["obs_info"], 

318 detector=c["detector"], 

319 band="r", 

320 ) 

321 

322 # Variance unit mismatch. 

323 with pytest.raises(ValueError): 

324 VisitImage( 

325 c["image"], 

326 variance=c["image"], 

327 psf=c["gaussian_psf"], 

328 mask_schema=c["mask_schema"], 

329 sky_projection=c["sky_projection"], 

330 obs_info=c["obs_info"], 

331 detector=c["detector"], 

332 band="r", 

333 ) 

334 

335 

336def test_copy_and_slice(visit_image_components: dict[str, Any]) -> None: 

337 """Verify that copy deep-copies arrays and components while slice shares 

338 them. 

339 """ 

340 c = visit_image_components 

341 visit_image = make_visit_image(c) 

342 copy = visit_image.copy() 

343 copy.image.array[0, 0] = 30.0 

344 assert visit_image.image.array[0, 0] == 42.0 

345 assert copy.image.array[0, 0] == 30.0 

346 subvisit = visit_image[Box.factory[0:5, 0:5]] 

347 # Check summary stats. 

348 assert copy.summary_stats == visit_image.summary_stats 

349 assert copy.summary_stats is not visit_image.summary_stats 

350 assert subvisit.summary_stats == visit_image.summary_stats 

351 assert subvisit.summary_stats is visit_image.summary_stats 

352 # Check aperture corrections. 

353 assert copy.aperture_corrections.keys() == visit_image.aperture_corrections.keys() 

354 assert copy.aperture_corrections is not visit_image.aperture_corrections 

355 assert subvisit.aperture_corrections.keys() == visit_image.aperture_corrections.keys() 

356 assert subvisit.aperture_corrections is visit_image.aperture_corrections 

357 # Check backgrounds. 

358 assert copy.backgrounds.keys() == visit_image.backgrounds.keys() 

359 assert copy.backgrounds is not visit_image.backgrounds 

360 assert subvisit.backgrounds.keys() == visit_image.backgrounds.keys() 

361 assert subvisit.backgrounds is visit_image.backgrounds 

362 # Check bounds. 

363 assert copy.bounds is c["polygon"] 

364 assert subvisit.bounds == subvisit.bbox # original polygon wholly encloses subvisit.bbox 

365 

366 

367def test_obs_info(visit_image_components: dict[str, Any]) -> None: 

368 """Verify that ObservationInfo is present and carries the expected 

369 instrument. 

370 """ 

371 visit_image = make_visit_image(visit_image_components) 

372 assert visit_image.obs_info is not None 

373 assert visit_image.obs_info.instrument == "LSSTCam" 

374 

375 

376def test_summary_stats(visit_image_components: dict[str, Any]) -> None: 

377 """Verify ObservationSummaryStats equality and inequality comparisons.""" 

378 summary_stats = visit_image_components["summary_stats"] 

379 assert summary_stats == ObservationSummaryStats(psfSigma=2.5, zeroPoint=31.4) 

380 assert summary_stats != ObservationSummaryStats(psfSigma=2.5) 

381 assert summary_stats != ObservationSummaryStats(psfSigma=2.5, raCorners=(5.2, 5.4, 5.4, 5.2)) 

382 

383 

384@skip_no_h5py 

385def test_round_trip_ndf(visit_image_components: dict[str, Any]) -> None: 

386 """Verify NDF round-trip produces a VisitImage equal to the original.""" 

387 visit_image = make_visit_image(visit_image_components) 

388 with RoundtripNdf(visit_image, "VisitImage") as roundtrip: 

389 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False) 

390 

391 

392@skip_no_h5py 

393def test_fits_ndf_consistency(visit_image_components: dict[str, Any]) -> None: 

394 """Verify FITS and NDF backends produce equal VisitImages on round-trip.""" 

395 visit_image = make_visit_image(visit_image_components) 

396 with RoundtripFits(visit_image) as fits_rt, RoundtripNdf(visit_image) as ndf_rt: 

397 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False) 

398 assert_visit_images_equal(visit_image, ndf_rt.result, expect_view=False) 

399 assert_visit_images_equal(fits_rt.result, ndf_rt.result, expect_view=False) 

400 

401 

402def test_fits_json_consistency(visit_image_components: dict[str, Any]) -> None: 

403 """Verify FITS and JSON backends produce equal VisitImages.""" 

404 visit_image = make_visit_image(visit_image_components) 

405 with ( 

406 RoundtripFits(visit_image) as fits_rt, 

407 RoundtripJson(visit_image) as json_rt, 

408 ): 

409 assert_visit_images_equal(visit_image, fits_rt.result, expect_view=False) 

410 assert_visit_images_equal(visit_image, json_rt.result, expect_view=False) 

411 assert_visit_images_equal(fits_rt.result, json_rt.result, expect_view=False) 

412 

413 

414def test_read_write(visit_image_components: dict[str, Any]) -> None: 

415 """Verify a VisitImage round-trips through FITS with correct compression. 

416 

417 Checks compression headers, subimage reads, equality, and opaque 

418 metadata. Contains only butler-free assertions; component reads live 

419 in `test_read_write_components`. 

420 """ 

421 visit_image = make_visit_image(visit_image_components) 

422 with RoundtripFits(visit_image, "VisitImage") as roundtrip: 

423 # Check that we're still using the right compression, and that we 

424 # wrote WCSs. 

425 fits = roundtrip.inspect() 

426 assert fits[1].header["ZCMPTYPE"] == "GZIP_2" 

427 assert fits[1].header["CTYPE1"] == "RA---TAN" 

428 assert fits[2].header["ZCMPTYPE"] == "GZIP_2" 

429 assert fits[2].header["CTYPE1"] == "RA---TAN" 

430 assert fits[3].header["ZCMPTYPE"] == "GZIP_2" 

431 assert fits[3].header["CTYPE1"] == "RA---TAN" 

432 # Check a subimage read (no component arg — does not trigger a skip). 

433 subbox = Box.factory[8:13, 9:30] 

434 subimage = roundtrip.get(bbox=subbox) 

435 assert_masked_images_equal(subimage, visit_image[subbox], expect_view=False) 

436 

437 assert_visit_images_equal(roundtrip.result, visit_image, expect_view=False) 

438 # Check that the round-tripped headers are the same (up to card order). 

439 assert len(roundtrip.result._opaque_metadata.headers[ExtensionKey()]) == 1 

440 assert dict(visit_image._opaque_metadata.headers[ExtensionKey()]) == dict( 

441 roundtrip.result._opaque_metadata.headers[ExtensionKey()] 

442 ) 

443 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("IMAGE")] 

444 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("MASK")] 

445 assert not roundtrip.result._opaque_metadata.headers[ExtensionKey("VARIANCE")] 

446 # Spot-check the concrete background contents (names, field types, 

447 # subtracted entry) against the known fixture, so the equality check 

448 # above is not vacuously satisfied by empty background maps. 

449 assert isinstance(roundtrip.result.backgrounds, BackgroundMap) 

450 assert roundtrip.result.backgrounds.keys() == {"standard"} 

451 assert isinstance(roundtrip.result.backgrounds["standard"].field, ChebyshevField) 

452 assert roundtrip.result.backgrounds.subtracted.name == "standard" 

453 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image." 

454 

455 

456def test_read_write_components(visit_image_components: dict[str, Any]) -> None: 

457 """Verify component reads and storage-class overrides round-trip correctly. 

458 

459 Requires a butler; skips when `lsst.daf.butler` is absent. 

460 Butler-free assertions live in `test_read_write`. 

461 """ 

462 c = visit_image_components 

463 visit_image = make_visit_image(c) 

464 with RoundtripFits(visit_image, "VisitImage") as roundtrip: 

465 subbox = Box.factory[8:13, 9:30] 

466 subimage = roundtrip.get(bbox=subbox) 

467 

468 # Get an explicit masked image to compare with the subimage. 

469 subimage_masked = roundtrip.get("masked_image", bbox=subbox) 

470 assert_masked_images_equal(subimage_masked, subimage, expect_view=False) 

471 

472 # Get the same masked image in a multi-component get and ensure 

473 # it is the same thing. 

474 components = roundtrip.get("components", components=["masked_image", "psf"], bbox=subbox) 

475 assert set(components) == {"masked_image", "psf"} 

476 assert_masked_images_equal(components["masked_image"], subimage_masked, expect_view=False) 

477 

478 assert roundtrip.get("bbox") == visit_image.bbox 

479 

480 obs_info = roundtrip.get("obs_info") 

481 assert isinstance(obs_info, ObservationInfo) 

482 assert obs_info == visit_image.obs_info 

483 

484 summary_stats = roundtrip.get("summary_stats") 

485 assert isinstance(summary_stats, ObservationSummaryStats) 

486 assert summary_stats == visit_image.summary_stats 

487 

488 psf = roundtrip.get("psf") 

489 assert isinstance(psf, GaussianPointSpreadFunction) 

490 assert psf.kernel_bbox == c["gaussian_psf"].kernel_bbox 

491 

492 backgrounds = roundtrip.get("backgrounds") 

493 assert isinstance(backgrounds, BackgroundMap) 

494 assert backgrounds.keys() == {"standard"} 

495 assert isinstance(backgrounds["standard"].field, ChebyshevField) 

496 assert backgrounds.subtracted.name == "standard" 

497 assert roundtrip.result.backgrounds.subtracted.description == "Background subtracted from the image." 

498 

499 # Test some components get edge cases. 

500 components = roundtrip.get("components", components="image") 

501 assert isinstance(components["image"], Image) 

502 

503 components = roundtrip.get("components") 

504 assert set(components) == { 

505 "image", 

506 "variance", 

507 "psf", 

508 "bbox", 

509 "mask", 

510 "obs_info", 

511 "backgrounds", 

512 "detector", 

513 "aperture_corrections", 

514 "sky_projection", 

515 "summary_stats", 

516 "photometric_scaling", 

517 } 

518 

519 # Butler morphs RuntimeError to ValueError. 

520 with pytest.raises(ValueError): 

521 roundtrip.get("components", components=["image", "nonexistent"]) 

522 

523 with pytest.raises(ValueError): 

524 roundtrip.get("components", components=["image", "components"]) 

525 

526 with pytest.raises(ValueError): 

527 roundtrip.get("components", components=[]) 

528 

529 with pytest.raises(ValueError): 

530 # PSF does not know how to use bbox so this fails. 

531 roundtrip.get("components", components="psf", bbox=subbox) 

532 

533 

534def test_sum_background_round_trip_fits(visit_image_components: dict[str, Any]) -> None: 

535 """Verify FITS backend keeps two same-named SumField operands as distinct 

536 EXTVERs. 

537 """ 

538 visit_image = make_visit_image(visit_image_components) 

539 visit = _make_sum_background_visit_image(visit_image_components, visit_image) 

540 with RoundtripFits(visit) as roundtrip: 

541 _check_sum_background_round_trip(roundtrip.result, visit) 

542 

543 

544@skip_no_h5py 

545def test_sum_background_round_trip_ndf(visit_image_components: dict[str, Any]) -> None: 

546 """Verify NDF backend disambiguates the repeated ``data`` leaf, just as 

547 the FITS backend does. 

548 """ 

549 visit_image = make_visit_image(visit_image_components) 

550 visit = _make_sum_background_visit_image(visit_image_components, visit_image) 

551 with RoundtripNdf(visit) as roundtrip: 

552 _check_sum_background_round_trip(roundtrip.result, visit) 

553 

554 

555@dataclasses.dataclass 

556class _LegacyTestData: 

557 filename: str 

558 plane_map: dict[str, MaskPlane] = dataclasses.field(default_factory=get_legacy_visit_image_mask_planes) 

559 unit: u.Unit = u.nJy 

560 storage_class: str = "VisitImage" 

561 read_cls: type[VisitImage] = VisitImage 

562 legacy_exposure: LegacyExposure = dataclasses.field(init=False) 

563 

564 @classmethod 

565 def get( 

566 cls, which: Literal["visit_image", "preliminary_visit_image", "difference_image"] 

567 ) -> _LegacyTestData: 

568 if EXTERNAL_DATA_DIR is None: 568 ↛ 570line 568 didn't jump to line 570 because the condition on line 568 was always true

569 pytest.skip("TESTDATA_IMAGES is not set up.") 

570 result = cls( 

571 os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", f"{which}.fits"), 

572 ) 

573 match which: 

574 case "preliminary_visit_image": 

575 result.unit = u.electron 

576 case "difference_image": 

577 result.storage_class = "DifferenceImage" 

578 result.read_cls = DifferenceImage 

579 result.plane_map = get_legacy_difference_image_mask_planes() 

580 case "visit_image": 

581 pass 

582 try: 

583 from lsst.afw.image import ExposureFitsReader 

584 

585 result.legacy_exposure = ExposureFitsReader(result.filename).read() 

586 except ImportError: 

587 pytest.skip("lsst.afw.image is not available; cannot read legacy exposures") 

588 result.visit_image = result.read_cls.read_legacy( 

589 result.filename, preserve_quantization=True, plane_map=result.plane_map 

590 ) 

591 return result 

592 

593 

594@pytest.fixture(scope="session", params=["visit_image", "preliminary_visit_image", "difference_image"]) 

595def legacy_test_data(request: pytest.FixtureRequest) -> _LegacyTestData: 

596 """Return legacy test data. 

597 

598 Tests that depend on this parameterized fixture run on all of the legacy 

599 test images. 

600 """ 

601 return _LegacyTestData.get(request.param) 

602 

603 

604@pytest.fixture(scope="session", params=["visit_image", "difference_image"]) 

605def legacy_test_data_calibrated(request: pytest.FixtureRequest) -> _LegacyTestData: 

606 """Return legacy test data for calibrated images only. 

607 

608 Tests that depend on this parameterized fixture do not run on 

609 preliminary_visit_image, since that has 'electron' pixel units 

610 """ 

611 return _LegacyTestData.get(request.param) 

612 

613 

614def _check_legacy_obs_info(obs_info: ObservationInfo | None) -> None: 

615 """Assert obs_info carries expected LSSTCam/DP2 field values.""" 

616 assert isinstance(obs_info, ObservationInfo) 

617 assert obs_info.instrument == "LSSTCam" 

618 assert obs_info.detector_num == 85, obs_info 

619 assert obs_info.detector_unique_name == "R21_S11", obs_info 

620 assert obs_info.physical_filter == "r_57", obs_info 

621 

622 

623def test_legacy_errors(legacy_test_data: _LegacyTestData) -> None: 

624 """Verify that from_legacy and read_legacy raise ValueError on 

625 conflicting arguments. 

626 """ 

627 with pytest.raises(ValueError): 

628 VisitImage.from_legacy(legacy_test_data.legacy_exposure, instrument="HSC") 

629 with pytest.raises(ValueError): 

630 VisitImage.from_legacy(legacy_test_data.legacy_exposure, visit=123456) 

631 with pytest.raises(ValueError): 

632 VisitImage.from_legacy(legacy_test_data.legacy_exposure, unit=u.mJy) 

633 visit = VisitImage.from_legacy( 

634 legacy_test_data.legacy_exposure, 

635 instrument="LSSTCam", 

636 unit=legacy_test_data.unit, 

637 visit=2025052000177, 

638 ) 

639 assert visit.unit == legacy_test_data.unit 

640 

641 with pytest.raises(ValueError): 

642 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, instrument="HSC") 

643 with pytest.raises(ValueError): 

644 legacy_test_data.read_cls.read_legacy(legacy_test_data.filename, visit=123456) 

645 

646 

647def test_component_reads(legacy_test_data: _LegacyTestData) -> None: 

648 """Verify that individual components can be read from a legacy FITS 

649 file. 

650 """ 

651 visit = VisitImage.read_legacy(legacy_test_data.filename) 

652 proj = VisitImage.read_legacy(legacy_test_data.filename, component="sky_projection") 

653 assert_sky_projections_equal(proj, visit.sky_projection, expect_identity=False) 

654 image = VisitImage.read_legacy(legacy_test_data.filename, component="image") 

655 assert image == visit.image 

656 assert_sky_projections_equal(proj, image.sky_projection, expect_identity=False) 

657 variance = VisitImage.read_legacy(legacy_test_data.filename, component="variance") 

658 assert variance == visit.variance 

659 assert_sky_projections_equal(proj, variance.sky_projection, expect_identity=False) 

660 mask = VisitImage.read_legacy(legacy_test_data.filename, component="mask") 

661 assert mask == visit.mask 

662 assert_sky_projections_equal(proj, mask.sky_projection, expect_identity=False) 

663 psf = VisitImage.read_legacy(legacy_test_data.filename, component="psf") 

664 assert isinstance(psf, PointSpreadFunction) 

665 obs_info = VisitImage.read_legacy(legacy_test_data.filename, component="obs_info") 

666 _check_legacy_obs_info(obs_info) 

667 summary_stats = VisitImage.read_legacy(legacy_test_data.filename, component="summary_stats") 

668 assert isinstance(summary_stats, ObservationSummaryStats) 

669 assert summary_stats.nPsfStar == legacy_test_data.legacy_exposure.info.getSummaryStats().nPsfStar 

670 compare_aperture_corrections_to_legacy( 

671 VisitImage.read_legacy(legacy_test_data.filename, component="aperture_corrections"), 

672 legacy_test_data.legacy_exposure.info.getApCorrMap(), 

673 visit.bbox, 

674 ) 

675 detector = VisitImage.read_legacy(legacy_test_data.filename, component="detector") 

676 compare_detector_to_legacy( 

677 detector, legacy_test_data.legacy_exposure.getDetector(), is_raw_assembled=True 

678 ) 

679 photometric_scaling = VisitImage.read_legacy(legacy_test_data.filename, component="photometric_scaling") 

680 compare_photo_calib_to_legacy( 

681 photometric_scaling, 

682 legacy_test_data.legacy_exposure.getPhotoCalib(), 

683 subimage_bbox=visit.bbox, 

684 ) 

685 

686 

687def test_legacy_obs_info(legacy_test_data: _LegacyTestData) -> None: 

688 """Verify that ObservationInfo is constructed correctly from a legacy 

689 exposure. 

690 """ 

691 legacy = VisitImage.from_legacy(legacy_test_data.legacy_exposure, plane_map=legacy_test_data.plane_map) 

692 assert legacy.obs_info is not None 

693 assert legacy.obs_info == legacy_test_data.visit_image.obs_info 

694 assert legacy.obs_info is not None # for mypy 

695 assert legacy.obs_info.instrument == "LSSTCam" 

696 assert legacy.obs_info.detector_num == 85, legacy.obs_info 

697 assert legacy.obs_info.detector_unique_name == "R21_S11", legacy.obs_info 

698 assert legacy.obs_info.physical_filter == "r_57", legacy.obs_info 

699 

700 

701def test_aperture_corrections_to_legacy(legacy_test_data: _LegacyTestData) -> None: 

702 """Verify that aperture corrections round-trip through a legacy 

703 ApCorrMap. 

704 """ 

705 ap_corrections = legacy_test_data.visit_image.aperture_corrections 

706 legacy_ap_corr_map = aperture_corrections_to_legacy(ap_corrections) 

707 compare_aperture_corrections_to_legacy( 

708 ap_corrections, 

709 legacy_ap_corr_map, 

710 legacy_test_data.visit_image.bbox, 

711 ) 

712 

713 

714def _check_legacy_headers(visit_image: VisitImage) -> None: 

715 """Assert that primary and extension headers are stripped correctly.""" 

716 header = visit_image._opaque_metadata.headers[ExtensionKey()] 

717 assert "EXPTIME" in header 

718 assert header["PLATFORM"] == "lsstcam" 

719 assert "LSST BUTLER ID" not in header 

720 assert "AR HDU" not in header 

721 assert "A_ORDER" not in header 

722 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("IMAGE"), astropy.io.fits.Header()) 

723 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("MASK"), astropy.io.fits.Header()) 

724 assert not visit_image._opaque_metadata.headers.get(ExtensionKey("VARIANCE"), astropy.io.fits.Header()) 

725 

726 

727def test_read_legacy_headers(legacy_test_data: _LegacyTestData) -> None: 

728 """Verify that headers were stripped and interpreted correctly in 

729 read_legacy. 

730 """ 

731 assert legacy_test_data.visit_image.unit == legacy_test_data.unit 

732 _check_legacy_headers(legacy_test_data.visit_image) 

733 

734 

735def test_from_legacy_headers(legacy_test_data: _LegacyTestData) -> None: