Coverage for tests/test_cell_coadd.py: 71%

291 statements  

« prev     ^ index     » next       coverage.py v7.15.2, created at 2026-08-08 04:08 -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 

14import copy 

15import dataclasses 

16import json 

17import os 

18import pickle 

19from pathlib import Path 

20from typing import Any 

21 

22import numpy as np 

23import pytest 

24 

25from lsst.images import YX, Box, Interval, MaskPlane, get_legacy_deep_coadd_mask_planes 

26from lsst.images.cells import ( 

27 CellCoadd, 

28 CellGrid, 

29 CellGridBounds, 

30 CellIJ, 

31 CellPointSpreadFunctionSerializationModel, 

32 CoaddProvenance, 

33 PatchDefinition, 

34) 

35from lsst.images.describe import DescribeOptions, Report 

36from lsst.images.fields import ChebyshevField 

37from lsst.images.fits import FitsCompressionOptions 

38from lsst.images.serialization import JsonRef, class_for_schema, parameterize_tree, read_archive 

39from lsst.images.tests import ( 

40 DP2_COADD_DATA_ID, 

41 DP2_COADD_MISSING_CELL, 

42 RoundtripFits, 

43 RoundtripJson, 

44 RoundtripNdf, 

45 assert_cell_coadds_equal, 

46 assert_images_equal, 

47 assert_masked_images_equal, 

48 assert_psfs_equal, 

49 check_bounds_contains_broadcasting, 

50 compare_cell_coadd_to_legacy, 

51 compare_masked_image_to_legacy, 

52 compare_psf_to_legacy, 

53 compare_sky_projection_to_legacy_wcs, 

54 current_fixture_path, 

55) 

56 

57try: 

58 import h5py # noqa: F401 

59 

60 HAVE_H5PY = True 

61except ImportError: 

62 HAVE_H5PY = False 

63 

64try: 

65 import lsst.afw.image # noqa: F401 

66 from lsst.cell_coadds import MultipleCellCoadd as LegacyMultipleCellCoadd 

67 

68 HAVE_LEGACY = True 

69except ImportError: 

70 HAVE_LEGACY = False 

71 type LegacyMultipleCellCoadd = Any # type: ignore[no-redef] 

72 

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

74FIXTURE_DIR = Path(__file__).parent / "data" / "schemas" 

75 

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

77skip_no_legacy = pytest.mark.skipif(not HAVE_LEGACY, reason="lsst.afw (etc) could not be imported.") 

78 

79 

80@dataclasses.dataclass 

81class _LegacyTestData: 

82 """A struct holding test data loaded from EXTERNAL_DATA_DIR.""" 

83 

84 filename: str 

85 tract_bbox: Box 

86 legacy_cell_coadd: LegacyMultipleCellCoadd 

87 cell_coadd: CellCoadd 

88 plane_map: dict[str, MaskPlane] = dataclasses.field(default_factory=get_legacy_deep_coadd_mask_planes) 

89 

90 def make_psf_points(self, bbox: Box | None = None) -> YX[np.ndarray]: 

91 """Create random PSF sample points within the given bbox.""" 

92 if bbox is None: 

93 bbox = self.cell_coadd.bbox 

94 rng = np.random.default_rng(44) 

95 xc, yc = np.meshgrid( 

96 np.arange( 

97 bbox.x.start + self.cell_coadd.grid.cell_shape.x * 0.5, 

98 bbox.x.stop, 

99 self.cell_coadd.grid.cell_shape.x, 

100 ), 

101 np.arange( 

102 bbox.y.start + self.cell_coadd.grid.cell_shape.y * 0.5, 

103 bbox.y.stop, 

104 self.cell_coadd.grid.cell_shape.y, 

105 ), 

106 ) 

107 return YX( 

108 y=yc.ravel() + rng.uniform(-0.4, 0.4, size=yc.size), 

109 x=xc.ravel() + rng.uniform(-0.4, 0.4, size=xc.size), 

110 ) 

111 

112 

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

114def legacy_test_data() -> _LegacyTestData: 

115 """Return a struct of CellCoadd loaded from legacy test data. 

116 

117 Skips if ``TESTDATA_IMAGES_DIR`` is not set or if ``lsst.cell_coadds`` 

118 cannot be imported. 

119 """ 

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

121 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.") 

122 try: 

123 from lsst.cell_coadds import MultipleCellCoadd 

124 except ImportError: 

125 pytest.skip("lsst.cell_coadds could not be imported.") 

126 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "deep_coadd_cell_predetection.fits") 

127 plane_map = get_legacy_deep_coadd_mask_planes() 

128 legacy_cell_coadd = MultipleCellCoadd.read_fits(filename) 

129 with open(os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "skyMap.pickle"), "rb") as stream: 

130 skymap = pickle.load(stream) 

131 cell_coadd = CellCoadd.from_legacy_cell_coadd( 

132 legacy_cell_coadd, 

133 plane_map=plane_map, 

134 tract_info=skymap[DP2_COADD_DATA_ID["tract"]], 

135 ) 

136 return _LegacyTestData( 

137 filename=filename, 

138 tract_bbox=Box.from_legacy(skymap[DP2_COADD_DATA_ID["tract"]].getBBox()), 

139 legacy_cell_coadd=legacy_cell_coadd, 

140 cell_coadd=cell_coadd, 

141 ) 

142 

143 

144@pytest.fixture 

145def minified_cell_coadd() -> CellCoadd: 

146 """Return a tiny CellCoadd from JSON data stored in this package.""" 

147 path = current_fixture_path(FIXTURE_DIR, "cell_coadd", variant="as_shipped") 

148 return read_archive(path, CellCoadd) 

149 

150 

151def make_subbox(full_bbox: Box) -> Box: 

152 """Make a box that's useful for nontrivial subimage tests. 

153 

154 This box only overlaps (but does not fully cover) the middle 2 (of 4) 

155 cells in y, while covering exactly the last column of cells in x. It does 

156 not cover the missing cell. 

157 """ 

158 return Box.factory[ 

159 full_bbox.y.start + 252 : full_bbox.y.stop - 175, 

160 full_bbox.x.stop - 150 : full_bbox.x.stop, 

161 ] 

162 

163 

164def test_cell_coadd_repr_str_pinned(minified_cell_coadd: CellCoadd) -> None: 

165 """Pin the exact str and repr output of a CellCoadd. 

166 

167 A coadd takes component types no string can express, so repr is the 

168 descriptive form rather than a constructor call. Both come from the 

169 report, as they do for the sibling `VisitImage`. 

170 """ 

171 assert str(minified_cell_coadd) == "CellCoadd([y=48:60, x=36:48], tract=9813)" 

172 assert repr(minified_cell_coadd) == "<CellCoadd([y=48:60, x=36:48], tract=9813)>" 

173 

174 

175def report_fields(report: Report) -> dict[str, Any]: 

176 """Return a mapping from report field label to field value.""" 

177 return {field.label: field.value for field in report.fields} 

178 

179 

180def make_test_provenance() -> CoaddProvenance: 

181 """Return a provenance with three input images taken over two nights, 

182 whose contributions cover three cells unevenly. 

183 """ 

184 inputs = CoaddProvenance.make_empty_input_table(3) 

185 inputs["instrument"] = "LSSTCam" 

186 inputs["physical_filter"] = "r_57" 

187 inputs["visit"] = [101, 102, 103] 

188 inputs["detector"] = [1, 2, 3] 

189 inputs["day_obs"] = [20250520, 20250520, 20250521] 

190 contributions = CoaddProvenance.make_empty_contribution_table(6) 

191 contributions["cell_i"] = [0, 0, 0, 1, 1, 2] 

192 contributions["cell_j"] = 0 

193 return CoaddProvenance(inputs=inputs, contributions=contributions) 

194 

195 

196def test_coadd_provenance_repr_str_pinned(minified_cell_coadd: CellCoadd) -> None: 

197 """Pin the str and repr of a CoaddProvenance. 

198 

199 Neither table can be expressed in an eval-able string, so repr is the 

200 descriptive form, and both report only what len() can answer. 

201 """ 

202 provenance = minified_cell_coadd.provenance 

203 assert str(provenance) == "CoaddProvenance(6 input images)" 

204 assert repr(provenance) == "<CoaddProvenance(6 input images)>" 

205 

206 

207def test_coadd_provenance_brief_report_skips_column_scans(minified_cell_coadd: CellCoadd) -> None: 

208 """The brief report carries no fields, so repr and str never scan a 

209 column of either table. 

210 """ 

211 report = minified_cell_coadd.provenance._describe(DescribeOptions(brief=True)) 

212 assert not report.fields 

213 assert report.to_repr() == "<CoaddProvenance(6 input images)>" 

214 

215 

216def test_coadd_provenance_report_summarizes_its_tables(minified_cell_coadd: CellCoadd) -> None: 

217 """The expanded report summarizes both tables instead of rendering them.""" 

218 report = minified_cell_coadd.provenance.describe() 

219 assert report.type_name == "CoaddProvenance" 

220 assert report_fields(report) == { 

221 "instrument": "LSSTCam", 

222 "physical_filter": "r_57", 

223 "input images": "6 from 2 visits", 

224 "day_obs": "20250520", 

225 "cells": "3 with contributions", 

226 "per cell": "2 input images", 

227 } 

228 assert list(report_fields(report)) == [ 

229 "instrument", 

230 "physical_filter", 

231 "input images", 

232 "day_obs", 

233 "cells", 

234 "per cell", 

235 ] 

236 # Rich renders an astropy table as plain text and never consults its 

237 # _repr_html_, so the report tabulates nothing. 

238 assert not report.tables 

239 assert "input images" in report._repr_html_() 

240 

241 

242def test_coadd_provenance_report_counts_cells_against_bounds(minified_cell_coadd: CellCoadd) -> None: 

243 """Given the image's cells, the report states coverage as a fraction.""" 

244 report = minified_cell_coadd.provenance._describe(bounds=minified_cell_coadd.bounds) 

245 assert report_fields(report)["cells"] == "3 of 3 with contributions" 

246 

247 

248def test_coadd_provenance_report_shows_partial_cell_coverage() -> None: 

249 """Cells with image data but no contribution rows show up in the ratio.""" 

250 grid = CellGrid(bbox=Box.from_shape((30, 30)), cell_shape=YX(10, 10)) 

251 bounds = CellGridBounds(grid=grid, bbox=Box.factory[0:30, 0:30]) 

252 report = make_test_provenance()._describe(bounds=bounds) 

253 assert report_fields(report)["cells"] == "3 of 9 with contributions" 

254 

255 

256def test_coadd_provenance_report_ranges_over_nights_and_cells() -> None: 

257 """Several visits, several nights and uneven cell coverage give the 

258 ranged forms of each field. 

259 """ 

260 fields = report_fields(make_test_provenance().describe()) 

261 assert fields["input images"] == "3 from 3 visits" 

262 assert fields["day_obs"] == "20250520 - 20250521" 

263 assert fields["cells"] == "3 with contributions" 

264 assert fields["per cell"] == "1 - 3 input images (median 2)" 

265 

266 

267def test_coadd_provenance_report_handles_empty_tables() -> None: 

268 """A provenance with no rows describes without raising.""" 

269 provenance = CoaddProvenance( 

270 inputs=CoaddProvenance.make_empty_input_table(0), 

271 contributions=CoaddProvenance.make_empty_contribution_table(0), 

272 ) 

273 assert str(provenance) == "CoaddProvenance(no input images)" 

274 report = provenance.describe() 

275 assert report_fields(report) == {"input images": "none", "cells": "none"} 

276 report.__rich__() 

277 

278 

279def test_coadd_provenance_report_counts_one_input_image_in_the_singular() -> None: 

280 """Pin the singular forms of the counted phrases. 

281 

282 One input image is what the standalone ``coadd_provenance`` file holds, 

283 and what a single-cell subset of a coadd holds, so this is the case the 

284 ``describe`` command line subcommand shows most often. 

285 """ 

286 inputs = CoaddProvenance.make_empty_input_table(1) 

287 inputs["instrument"] = "LSSTCam" 

288 inputs["physical_filter"] = "r_57" 

289 inputs["visit"] = [101] 

290 inputs["detector"] = [1] 

291 inputs["day_obs"] = [20250520] 

292 provenance = CoaddProvenance( 

293 inputs=inputs, contributions=CoaddProvenance.make_empty_contribution_table(1) 

294 ) 

295 assert str(provenance) == "CoaddProvenance(1 input image)" 

296 assert repr(provenance) == "<CoaddProvenance(1 input image)>" 

297 fields = report_fields(provenance.describe()) 

298 assert fields["input images"] == "1 from 1 visit" 

299 assert fields["per cell"] == "1 input image" 

300 

301 

302def test_cell_coadd_report_accounts_for_every_component(minified_cell_coadd: CellCoadd) -> None: 

303 """Nothing the coadd carries is absent from its report.""" 

304 report = minified_cell_coadd.describe() 

305 fields = report_fields(report) 

306 assert fields["mask_fractions"] == "rejected" 

307 assert fields["noise_realizations"] == "1 image (dtype float32)" 

308 assert fields["aperture_corrections"] == "3 fields" 

309 # Provenance is a child, so it needs no field line of its own. 

310 assert "provenance" not in fields 

311 assert list(report.children) == [ 

312 "image", 

313 "mask", 

314 "variance", 

315 "sky_projection", 

316 "psf", 

317 "provenance", 

318 "backgrounds", 

319 ] 

320 provenance = report.children["provenance"] 

321 assert provenance.type_name == "CoaddProvenance" 

322 # The coadd passed its cells down, so coverage is stated as a fraction. 

323 assert report_fields(provenance)["cells"] == "3 of 3 with contributions" 

324 

325 

326def test_cell_coadd_report_counts_aperture_corrections_only(minified_cell_coadd: CellCoadd) -> None: 

327 """A coadd can carry dozens of aperture corrections with long names, so 

328 the report counts them and never names them, detail or not. 

329 """ 

330 plain = report_fields(minified_cell_coadd.describe())["aperture_corrections"] 

331 assert plain == "3 fields" 

332 detailed = report_fields(minified_cell_coadd.describe(detail=True))["aperture_corrections"] 

333 assert detailed == "3 fields" # No additional detail 

334 

335 

336def test_cell_coadd_report_states_absent_provenance(minified_cell_coadd: CellCoadd) -> None: 

337 """A coadd with no provenance says so, while components that are merely 

338 empty stay out of the report. 

339 """ 

340 bare = CellCoadd( 

341 minified_cell_coadd.image, 

342 mask=minified_cell_coadd.mask, 

343 variance=minified_cell_coadd.variance, 

344 sky_projection=minified_cell_coadd.sky_projection, 

345 band=minified_cell_coadd.band, 

346 psf=minified_cell_coadd.psf, 

347 patch=minified_cell_coadd.patch, 

348 ) 

349 report = bare.describe() 

350 fields = report_fields(report) 

351 assert fields["provenance"] == "none" 

352 assert "provenance" not in report.children 

353 assert "mask_fractions" not in fields 

354 assert "noise_realizations" not in fields 

355 assert "aperture_corrections" not in fields 

356 

357 

358def test_cell_grid_patch_str_uses_clean_geometry() -> None: 

359 """CellGrid, PatchDefinition and CellGridBounds str drop the 

360 Interval/YX/Box/CellIJ wrappers. 

361 

362 The report renders field values with str, so these must use the compact 

363 geometry forms rather than pydantic's default field-by-field repr. 

364 """ 

365 grid = CellGrid(bbox=Box.from_shape((100, 200)), cell_shape=YX(10, 20)) 

366 assert str(grid) == "[y=0:100, x=0:200], cell_shape=(y=10, x=20)" 

367 

368 patch = PatchDefinition(id=73, index=YX(7, 3), inner_bbox=Box.factory[1:3, 2:4], cells=grid) 

369 assert str(patch) == ( 

370 "id=73, index=(y=7, x=3), inner_bbox=[y=1:3, x=2:4], " 

371 "cells=([y=0:100, x=0:200], cell_shape=(y=10, x=20))" 

372 ) 

373 # repr stays as the pydantic default so it remains eval-ish and distinct. 

374 assert "Interval(" in repr(patch) 

375 assert "YX(" in repr(patch) 

376 

377 bounds = CellGridBounds(grid=grid, bbox=Box.factory[0:40, 0:60]) 

378 assert str(bounds) == "[y=0:40, x=0:60] in grid ([y=0:100, x=0:200], cell_shape=(y=10, x=20))" 

379 bounds_missing = CellGridBounds( 

380 grid=grid, bbox=Box.factory[0:40, 0:60], missing=frozenset({CellIJ(i=1, j=1), CellIJ(i=0, j=2)}) 

381 ) 

382 assert str(bounds_missing) == ( 

383 "[y=0:40, x=0:60] in grid ([y=0:100, x=0:200], cell_shape=(y=10, x=20)), " 

384 "missing={(i=0, j=2), (i=1, j=1)}" 

385 ) 

386 assert "Interval(" in repr(bounds_missing) 

387 

388 

389def test_cell_shape_accepts_both_spellings() -> None: 

390 """Verify both on-disk spellings of an XY pair still validate. 

391 

392 Shipped cell_coadd 1.0.0 files carry the array spelling, so this stays 

393 readable until a major bump retires the as_shipped fixture. A focused 

394 check here gives a two-line failure instead of a 52 KB fixture failing 

395 to read. 

396 """ 

397 tree_cls = class_for_schema("cell_psf") 

398 assert tree_cls is not None 

399 model = parameterize_tree(tree_cls, JsonRef) 

400 fixture = json.loads(current_fixture_path(FIXTURE_DIR, "cell_psf").read_text()) 

401 assert fixture["bounds"]["grid"]["cell_shape"] == {"y": 4, "x": 4} 

402 

403 legacy = copy.deepcopy(fixture) 

404 legacy["bounds"]["grid"]["cell_shape"] = [4, 4] 

405 tree = model.model_validate(legacy) 

406 assert isinstance(tree, CellPointSpreadFunctionSerializationModel) 

407 assert tree.bounds.grid.cell_shape.y == 4 

408 assert tree.bounds.grid.cell_shape.x == 4 

409 

410 

411def test_from_legacy(legacy_test_data: _LegacyTestData) -> None: 

412 """Test constructing a CellCoadd by converting a legacy 

413 ``MultipleCellCoadd``. 

414 """ 

415 assert legacy_test_data.cell_coadd.bounds.missing == {CellIJ(**DP2_COADD_MISSING_CELL)} 

416 assert legacy_test_data.cell_coadd.bbox == Box.factory[12900:13500, 9600:10050] 

417 compare_cell_coadd_to_legacy( 

418 legacy_test_data.cell_coadd, 

419 legacy_test_data.legacy_cell_coadd, 

420 tract_bbox=legacy_test_data.tract_bbox, 

421 plane_map=legacy_test_data.plane_map, 

422 psf_points=legacy_test_data.make_psf_points(), 

423 ) 

424 

425 

426def test_roundtrip(legacy_test_data: _LegacyTestData) -> None: 

427 """Test that a CellCoadd roundtrips through FITS.""" 

428 with RoundtripFits(legacy_test_data.cell_coadd, "CellCoadd") as roundtrip: 

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

430 subbox = Box.factory[ 

431 legacy_test_data.cell_coadd.bbox.y.start + 252 : legacy_test_data.cell_coadd.bbox.y.stop - 175, 

432 legacy_test_data.cell_coadd.bbox.x.stop - 150 : legacy_test_data.cell_coadd.bbox.x.stop, 

433 ] 

434 subimage = roundtrip.get(bbox=subbox) 

435 assert_masked_images_equal(subimage, legacy_test_data.cell_coadd[subbox], expect_view=False) 

436 with roundtrip.inspect() as fits: 

437 for extname in ["IMAGE", "MASK", "VARIANCE", "MASK_FRACTIONS/REJECTED"] + [ 

438 f"NOISE_REALIZATIONS/{n}" for n in range(len(legacy_test_data.cell_coadd.noise_realizations)) 

439 ]: 

440 assert fits[extname].header["ZTILE1"] == legacy_test_data.cell_coadd.grid.cell_shape.x 

441 assert fits[extname].header["ZTILE2"] == legacy_test_data.cell_coadd.grid.cell_shape.y 

442 # Fixture self-consistency: bbox and missing-cell set are as expected. 

443 assert legacy_test_data.cell_coadd.bounds.missing == {CellIJ(**DP2_COADD_MISSING_CELL)} 

444 assert legacy_test_data.cell_coadd.bbox == Box.factory[12900:13500, 9600:10050] 

445 # Full round-trip fidelity. 

446 assert_cell_coadds_equal(roundtrip.result, legacy_test_data.cell_coadd, expect_view=False) 

447 compare_cell_coadd_to_legacy( 

448 roundtrip.result, 

449 legacy_test_data.legacy_cell_coadd, 

450 tract_bbox=legacy_test_data.tract_bbox, 

451 plane_map=legacy_test_data.plane_map, 

452 psf_points=legacy_test_data.make_psf_points(), 

453 ) 

454 

455 

456def test_roundtrip_components(legacy_test_data: _LegacyTestData) -> None: 

457 """Test component and subimage reads. 

458 

459 This test will be skipped if `lsst.daf.butler` is not available instead of 

460 falling back to non-butler I/O, which is why we don't want to merge it 

461 with `test_roundtrip`. 

462 """ 

463 with RoundtripFits(legacy_test_data.cell_coadd, "CellCoadd") as roundtrip: 

464 subbox = make_subbox(legacy_test_data.cell_coadd.bbox) 

465 subpsf = roundtrip.get("psf", bbox=subbox) 

466 assert subpsf.bounds.bbox == Box( 

467 y=Interval.factory[ 

468 legacy_test_data.cell_coadd.bbox.y.start + 150 : legacy_test_data.cell_coadd.bbox.y.stop - 150 

469 ], 

470 x=subbox.x, 

471 ) 

472 assert_psfs_equal( 

473 subpsf, 

474 legacy_test_data.cell_coadd.psf, 

475 points=legacy_test_data.make_psf_points(subbox), 

476 ) 

477 assert roundtrip.get("bbox") == legacy_test_data.cell_coadd.bbox 

478 alternates = { 

479 k: roundtrip.get(k) 

480 for k in [ 

481 "sky_projection", 

482 "image", 

483 "mask", 

484 "variance", 

485 "masked_image", 

486 "psf", 

487 "aperture_corrections", 

488 "provenance", 

489 "backgrounds", 

490 "bbox", 

491 ] 

492 } 

493 # Read all the components at once. 

494 all_components = roundtrip.get("components") 

495 assert set(all_components) == set(alternates) - {"masked_image"} 

496 assert all_components["bbox"] == alternates["bbox"] 

497 assert_psfs_equal(all_components["psf"], alternates["psf"]) 

498 assert_images_equal(all_components["image"], alternates["image"]) 

499 

500 backgrounds = roundtrip.get("backgrounds") 

501 assert backgrounds.keys() == set() 

502 assert backgrounds.subtracted is None 

503 

504 compare_cell_coadd_to_legacy( 

505 roundtrip.result, 

506 legacy_test_data.legacy_cell_coadd, 

507 tract_bbox=legacy_test_data.tract_bbox, 

508 plane_map=legacy_test_data.plane_map, 

509 alternates=alternates, 

510 psf_points=legacy_test_data.make_psf_points(), 

511 ) 

512 

513 

514def test_fits_compression(legacy_test_data: _LegacyTestData) -> None: 

515 """Test lossy FITS compression produces the expected headers.""" 

516 with RoundtripFits( 

517 legacy_test_data.cell_coadd, 

518 storage_class="CellCoadd", 

519 recipe="lossy16", 

520 compression_options={ 

521 "image": FitsCompressionOptions.LOSSY, 

522 "variance": FitsCompressionOptions.LOSSY, 

523 }, 

524 ) as roundtrip: 

525 with roundtrip.inspect() as fits: 

526 for extname in ["IMAGE", "MASK", "VARIANCE", "MASK_FRACTIONS/REJECTED"] + [ 

527 f"NOISE_REALIZATIONS/{n}" for n in range(len(legacy_test_data.cell_coadd.noise_realizations)) 

528 ]: 

529 assert fits[extname].header["ZTILE1"] == legacy_test_data.cell_coadd.grid.cell_shape.x 

530 assert fits[extname].header["ZTILE2"] == legacy_test_data.cell_coadd.grid.cell_shape.y 

531 if extname == "MASK" or extname.startswith("MASK_FRACTIONS"): 

532 assert fits[extname].header["ZCMPTYPE"] == "GZIP_2" 

533 else: 

534 assert fits[extname].header["ZCMPTYPE"] == "RICE_1" 

535 assert fits[extname].header["ZQUANTIZ"] == "SUBTRACTIVE_DITHER_2" 

536 

537 

538def test_json_roundtrip(legacy_test_data: _LegacyTestData) -> None: 

539 """Verify a CellCoadd round-trips correctly through the JSON archive.""" 

540 with RoundtripJson(legacy_test_data.cell_coadd) as roundtrip: 

541 pass 

542 assert_cell_coadds_equal(roundtrip.result, legacy_test_data.cell_coadd, expect_view=False) 

543 

544 

545def test_to_legacy_cell_coadd(legacy_test_data: _LegacyTestData) -> None: 

546 """Verify converting a CellCoadd back into a legacy MultipleCellCoadd.""" 

547 legacy_cell_coadd = legacy_test_data.cell_coadd.to_legacy_cell_coadd() 

548 compare_cell_coadd_to_legacy( 

549 legacy_test_data.cell_coadd, 

550 legacy_cell_coadd, 

551 tract_bbox=legacy_test_data.tract_bbox, 

552 plane_map=legacy_test_data.plane_map, 

553 psf_points=legacy_test_data.make_psf_points(), 

554 ) 

555 with pytest.raises( 

556 ValueError, match="MultipleCellCoadd requires its bounding box to lie on the cell grid." 

557 ): 

558 legacy_test_data.cell_coadd[make_subbox(legacy_test_data.cell_coadd.bbox)].to_legacy_cell_coadd() 

559 

560 

561@skip_no_legacy 

562def test_to_legacy(legacy_test_data: _LegacyTestData) -> None: 

563 """Test converting a CellCoadd back into a legacy Exposure.""" 

564 legacy_exposure = legacy_test_data.cell_coadd.to_legacy() 

565 assert legacy_exposure.getFilter().bandLabel == legacy_test_data.cell_coadd.band 

566 assert Box.from_legacy(legacy_exposure.getBBox()) == legacy_test_data.cell_coadd.bbox 

567 compare_masked_image_to_legacy( 

568 legacy_test_data.cell_coadd, 

569 legacy_exposure.maskedImage, 

570 plane_map=legacy_test_data.plane_map, 

571 expect_view=True, 

572 ) 

573 compare_psf_to_legacy( 

574 legacy_test_data.cell_coadd.psf, 

575 legacy_exposure.getPsf(), 

576 points=legacy_test_data.make_psf_points(), 

577 expect_legacy_raise_on_out_of_bounds=True, 

578 ) 

579 compare_sky_projection_to_legacy_wcs( 

580 legacy_test_data.cell_coadd.sky_projection, 

581 legacy_exposure.getWcs(), 

582 legacy_test_data.cell_coadd.sky_projection.pixel_frame, 

583 subimage_bbox=legacy_test_data.cell_coadd.bbox, 

584 is_fits=True, 

585 ) 

586 subbox = make_subbox(legacy_test_data.cell_coadd.bbox) 

587 compare_masked_image_to_legacy( 

588 legacy_test_data.cell_coadd[subbox], 

589 legacy_test_data.cell_coadd[subbox].to_legacy().maskedImage, 

590 plane_map=legacy_test_data.plane_map, 

591 expect_view=True, 

592 ) 

593 

594 

595@skip_no_h5py 

596def test_ndf_roundtrip(legacy_test_data: _LegacyTestData) -> None: 

597 """Test that CellCoadd round-trips through NDF.""" 

598 with RoundtripNdf(legacy_test_data.cell_coadd, "CellCoadd") as roundtrip: 

599 assert_cell_coadds_equal(roundtrip.result, legacy_test_data.cell_coadd, expect_view=False) 

600 

601 

602# The float32 image pixels are of order 10 nJy, as is the test background 

603# below, so background arithmetic is only reproducible to a few ULPs at that 

604# scale. 

605BACKGROUND_ATOL = 1e-5 

606 

607 

608def _add_gradient_background(coadd: CellCoadd, name: str = "pretty") -> ChebyshevField: 

609 """Add a non-constant background over the coadd's full bbox and return the 

610 field that was added. 

611 """ 

612 field = ChebyshevField( 

613 coadd.bbox, 

614 np.array([[10.0, 2.0, 0.5], [1.0, 0.75, 0.0], [0.25, 0.0, 0.0]]), 

615 unit=coadd.unit, 

616 ) 

617 coadd.backgrounds.add(name, field, description="Gradient background for tests.") 

618 return field 

619 

620 

621def _make_background_subbox(coadd: CellCoadd) -> Box: 

622 """Make a box that trims a different number of pixels from each side of the 

623 coadd, so a background rendered over the wrong box cannot match by chance. 

624 """ 

625 return Box.factory[ 

626 coadd.bbox.y.start + 2 : coadd.bbox.y.stop - 3, 

627 coadd.bbox.x.start + 1 : coadd.bbox.x.stop - 2, 

628 ] 

629 

630 

631def test_apply_background_subimage(minified_cell_coadd: CellCoadd) -> None: 

632 """Test that applying a background to a subimage subtracts only the 

633 portion of the background model that overlaps the subimage. 

634 """ 

635 field = _add_gradient_background(minified_cell_coadd) 

636 subbox = _make_background_subbox(minified_cell_coadd) 

637 subimage = minified_cell_coadd[subbox] 

638 original = subimage.image.array.copy() 

639 

640 subimage.apply_background("pretty") 

641 

642 assert subimage.backgrounds.subtracted is not None 

643 assert subimage.backgrounds.subtracted.name == "pretty" 

644 assert subimage.image.bbox == subbox 

645 expected = field.render(subbox, dtype=subimage.image.array.dtype).array 

646 np.testing.assert_allclose(original - subimage.image.array, expected, atol=BACKGROUND_ATOL) 

647 

648 

649def test_restore_background_subimage(minified_cell_coadd: CellCoadd) -> None: 

650 """Test that restoring the original background on a subimage adds back 

651 only the portion of the background model that overlaps the subimage. 

652 """ 

653 _add_gradient_background(minified_cell_coadd) 

654 subbox = _make_background_subbox(minified_cell_coadd) 

655 subimage = minified_cell_coadd[subbox] 

656 original = subimage.image.array.copy() 

657 subimage.apply_background("pretty") 

658 

659 subimage.apply_background(None) 

660 

661 assert subimage.backgrounds.subtracted is None 

662 np.testing.assert_allclose(subimage.image.array, original, atol=BACKGROUND_ATOL) 

663 

664 

665def test_apply_background_after_bounded_read(minified_cell_coadd: CellCoadd) -> None: 

666 """Test that a background can be applied to a coadd read with a bbox 

667 parameter. 

668 """ 

669 field = _add_gradient_background(minified_cell_coadd) 

670 subbox = _make_background_subbox(minified_cell_coadd) 

671 with RoundtripJson(minified_cell_coadd, "CellCoadd") as roundtrip: 

672 subimage = roundtrip.get(bbox=subbox) 

673 original = subimage.image.array.copy() 

674 

675 subimage.apply_background("pretty") 

676 

677 assert subimage.image.bbox == subbox 

678 expected = field.render(subbox, dtype=subimage.image.array.dtype).array 

679 np.testing.assert_allclose(original - subimage.image.array, expected, atol=BACKGROUND_ATOL) 

680 

681 

682def test_cell_grid_bounds_contains_broadcasting(minified_cell_coadd: CellCoadd) -> None: 

683 """Test that CellGridBounds.contains broadcasts like a numpy ufunc.""" 

684 assert minified_cell_coadd.bounds.missing, "fixture should retain a missing cell" 

685 check_bounds_contains_broadcasting(minified_cell_coadd.bounds) 

686 

687 

688def test_intersection_bounds_contains_broadcasting(minified_cell_coadd: CellCoadd) -> None: 

689 """Test that IntersectionBounds.contains broadcasts like a numpy ufunc.""" 

690 # Clip the CellGridBounds with a Box offset by 1 pixel on each side so it 

691 # does not snap to any cell boundary, forcing a lazy IntersectionBounds. 

692 bounds = minified_cell_coadd.bounds 

693 clip = Box.factory[ 

694 bounds.bbox.y.start + 1 : bounds.bbox.y.stop - 1, 

695 bounds.bbox.x.start + 1 : bounds.bbox.x.stop - 1, 

696 ] 

697 check_bounds_contains_broadcasting(bounds.intersection(clip))