Coverage for tests/test_cell_coadd.py: 75%
342 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-20 09:29 +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.
12from __future__ import annotations
14import copy
15import dataclasses
16import json
17import os
18import pickle
19from pathlib import Path
20from typing import Any
22import numpy as np
23import pytest
25from lsst.images import YX, BoundsError, 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 (
39 JsonRef,
40 class_for_schema,
41 parameterize_tree,
42 read_archive,
43)
44from lsst.images.tests import (
45 DP2_COADD_DATA_ID,
46 DP2_COADD_MISSING_CELL,
47 RoundtripFits,
48 RoundtripJson,
49 RoundtripNdf,
50 assert_cell_coadds_equal,
51 assert_images_equal,
52 assert_masked_images_equal,
53 assert_psfs_equal,
54 check_bounds_contains_broadcasting,
55 compare_cell_coadd_to_legacy,
56 compare_masked_image_to_legacy,
57 compare_psf_to_legacy,
58 compare_sky_projection_to_legacy_wcs,
59 current_fixture_path,
60 reset_afw_mask_planes, # noqa: F401
61)
63try:
64 import h5py # noqa: F401
66 HAVE_H5PY = True
67except ImportError:
68 HAVE_H5PY = False
70try:
71 import lsst.afw.image # noqa: F401
72 from lsst.cell_coadds import MultipleCellCoadd as LegacyMultipleCellCoadd
74 HAVE_LEGACY = True
75except ImportError:
76 HAVE_LEGACY = False
77 type LegacyMultipleCellCoadd = Any # type: ignore[no-redef]
79EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None)
80FIXTURE_DIR = Path(__file__).parent / "data" / "schemas"
82skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed")
83skip_no_legacy = pytest.mark.skipif(not HAVE_LEGACY, reason="lsst.afw (etc) could not be imported.")
86@dataclasses.dataclass
87class _LegacyTestData:
88 """A struct holding test data loaded from EXTERNAL_DATA_DIR."""
90 filename: str
91 tract_bbox: Box
92 legacy_cell_coadd: LegacyMultipleCellCoadd
93 cell_coadd: CellCoadd
94 plane_map: dict[str, MaskPlane] = dataclasses.field(default_factory=get_legacy_deep_coadd_mask_planes)
96 def make_psf_points(self, bbox: Box | None = None) -> YX[np.ndarray]:
97 """Create random PSF sample points within the given bbox."""
98 if bbox is None:
99 bbox = self.cell_coadd.bbox
100 rng = np.random.default_rng(44)
101 xc, yc = np.meshgrid(
102 np.arange(
103 bbox.x.start + self.cell_coadd.grid.cell_shape.x * 0.5,
104 bbox.x.stop,
105 self.cell_coadd.grid.cell_shape.x,
106 ),
107 np.arange(
108 bbox.y.start + self.cell_coadd.grid.cell_shape.y * 0.5,
109 bbox.y.stop,
110 self.cell_coadd.grid.cell_shape.y,
111 ),
112 )
113 return YX(
114 y=yc.ravel() + rng.uniform(-0.4, 0.4, size=yc.size),
115 x=xc.ravel() + rng.uniform(-0.4, 0.4, size=xc.size),
116 )
119@pytest.fixture
120def legacy_test_data(reset_afw_mask_planes: None) -> _LegacyTestData: # noqa: F811
121 """Return a struct of CellCoadd loaded from legacy test data.
123 Skips if ``TESTDATA_IMAGES_DIR`` is not set or if ``lsst.cell_coadds``
124 cannot be imported.
125 """
126 if EXTERNAL_DATA_DIR is None: 126 ↛ 128line 126 didn't jump to line 128 because the condition on line 126 was always true
127 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.")
128 try:
129 from lsst.cell_coadds import MultipleCellCoadd
130 except ImportError:
131 pytest.skip("lsst.cell_coadds could not be imported.")
132 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "deep_coadd_cell_predetection.fits")
133 plane_map = get_legacy_deep_coadd_mask_planes()
134 legacy_cell_coadd = MultipleCellCoadd.read_fits(filename)
135 with open(os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "skyMap.pickle"), "rb") as stream:
136 skymap = pickle.load(stream)
137 cell_coadd = CellCoadd.from_legacy_cell_coadd(
138 legacy_cell_coadd,
139 plane_map=plane_map,
140 tract_info=skymap[DP2_COADD_DATA_ID["tract"]],
141 )
142 return _LegacyTestData(
143 filename=filename,
144 tract_bbox=Box.from_legacy(skymap[DP2_COADD_DATA_ID["tract"]].getBBox()),
145 legacy_cell_coadd=legacy_cell_coadd,
146 cell_coadd=cell_coadd,
147 )
150@pytest.fixture
151def minified_cell_coadd() -> CellCoadd:
152 """Return a tiny CellCoadd from JSON data stored in this package."""
153 path = current_fixture_path(FIXTURE_DIR, "cell_coadd", variant="as_shipped")
154 return read_archive(path, CellCoadd)
157def make_subbox(full_bbox: Box) -> Box:
158 """Make a box that's useful for nontrivial subimage tests.
160 This box only overlaps (but does not fully cover) the middle 2 (of 4)
161 cells in y, while covering exactly the last column of cells in x. It does
162 not cover the missing cell.
163 """
164 return Box.factory[
165 full_bbox.y.start + 252 : full_bbox.y.stop - 175,
166 full_bbox.x.stop - 150 : full_bbox.x.stop,
167 ]
170def test_cell_coadd_repr_str_pinned(minified_cell_coadd: CellCoadd) -> None:
171 """Pin the exact str and repr output of a CellCoadd.
173 A coadd takes component types no string can express, so repr is the
174 descriptive form rather than a constructor call. Both come from the
175 report, as they do for the sibling `VisitImage`.
176 """
177 assert str(minified_cell_coadd) == "CellCoadd([y=48:60, x=36:48], tract=9813)"
178 assert repr(minified_cell_coadd) == "<CellCoadd([y=48:60, x=36:48], tract=9813)>"
181def report_fields(report: Report) -> dict[str, Any]:
182 """Return a mapping from report field label to field value."""
183 return {field.label: field.value for field in report.fields}
186def make_test_provenance() -> CoaddProvenance:
187 """Return a provenance with three input images taken over two nights,
188 whose contributions cover three cells unevenly.
189 """
190 inputs = CoaddProvenance.make_empty_input_table(3)
191 inputs["instrument"] = "LSSTCam"
192 inputs["physical_filter"] = "r_57"
193 inputs["visit"] = [101, 102, 103]
194 inputs["detector"] = [1, 2, 3]
195 inputs["day_obs"] = [20250520, 20250520, 20250521]
196 contributions = CoaddProvenance.make_empty_contribution_table(6)
197 contributions["cell_i"] = [0, 0, 0, 1, 1, 2]
198 contributions["cell_j"] = 0
199 return CoaddProvenance(inputs=inputs, contributions=contributions)
202def test_coadd_provenance_repr_str_pinned(minified_cell_coadd: CellCoadd) -> None:
203 """Pin the str and repr of a CoaddProvenance.
205 Neither table can be expressed in an eval-able string, so repr is the
206 descriptive form, and both report only what len() can answer.
207 """
208 provenance = minified_cell_coadd.provenance
209 assert str(provenance) == "CoaddProvenance(6 input images)"
210 assert repr(provenance) == "<CoaddProvenance(6 input images)>"
213def test_coadd_provenance_brief_report_skips_column_scans(minified_cell_coadd: CellCoadd) -> None:
214 """The brief report carries no fields, so repr and str never scan a
215 column of either table.
216 """
217 report = minified_cell_coadd.provenance._describe(DescribeOptions(brief=True))
218 assert not report.fields
219 assert report.to_repr() == "<CoaddProvenance(6 input images)>"
222def test_coadd_provenance_report_summarizes_its_tables(minified_cell_coadd: CellCoadd) -> None:
223 """The expanded report summarizes both tables instead of rendering them."""
224 report = minified_cell_coadd.provenance.describe()
225 assert report.type_name == "CoaddProvenance"
226 assert report_fields(report) == {
227 "instrument": "LSSTCam",
228 "physical_filter": "r_57",
229 "input images": "6 from 2 visits",
230 "day_obs": "20250520",
231 "cells": "3 with contributions",
232 "per cell": "2 input images",
233 }
234 assert list(report_fields(report)) == [
235 "instrument",
236 "physical_filter",
237 "input images",
238 "day_obs",
239 "cells",
240 "per cell",
241 ]
242 # Rich renders an astropy table as plain text and never consults its
243 # _repr_html_, so the report tabulates nothing.
244 assert not report.tables
245 assert "input images" in report._repr_html_()
248def test_coadd_provenance_report_counts_cells_against_bounds(minified_cell_coadd: CellCoadd) -> None:
249 """Given the image's cells, the report states coverage as a fraction."""
250 report = minified_cell_coadd.provenance._describe(bounds=minified_cell_coadd.bounds)
251 assert report_fields(report)["cells"] == "3 of 3 with contributions"
254def test_coadd_provenance_report_shows_partial_cell_coverage() -> None:
255 """Cells with image data but no contribution rows show up in the ratio."""
256 grid = CellGrid(bbox=Box.from_shape((30, 30)), cell_shape=YX(10, 10))
257 bounds = CellGridBounds(grid=grid, bbox=Box.factory[0:30, 0:30])
258 report = make_test_provenance()._describe(bounds=bounds)
259 assert report_fields(report)["cells"] == "3 of 9 with contributions"
262def test_coadd_provenance_report_ranges_over_nights_and_cells() -> None:
263 """Several visits, several nights and uneven cell coverage give the
264 ranged forms of each field.
265 """
266 fields = report_fields(make_test_provenance().describe())
267 assert fields["input images"] == "3 from 3 visits"
268 assert fields["day_obs"] == "20250520 - 20250521"
269 assert fields["cells"] == "3 with contributions"
270 assert fields["per cell"] == "1 - 3 input images (median 2)"
273def test_empty_coadd_provenance_cannot_be_constructed() -> None:
274 """Provenance is only meaningful when it records an input image."""
275 with pytest.raises(ValueError, match="at least one input"):
276 CoaddProvenance(
277 inputs=CoaddProvenance.make_empty_input_table(0),
278 contributions=CoaddProvenance.make_empty_contribution_table(0),
279 )
282def test_coadd_provenance_subset_without_contributions_is_none() -> None:
283 """Subsetting to cells that no image contributed to yields no provenance
284 rather than an empty one.
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 contributions = CoaddProvenance.make_empty_contribution_table(1)
293 contributions["cell_i"] = [0]
294 contributions["cell_j"] = [0]
295 contributions["instrument"] = "LSSTCam"
296 contributions["visit"] = [101]
297 contributions["detector"] = [1]
298 provenance = CoaddProvenance(inputs=inputs, contributions=contributions)
300 populated = provenance[CellIJ(0, 0)]
301 assert populated is not None
302 assert len(populated.inputs) == 1
303 assert len(populated.contributions) == 1
305 assert provenance[CellIJ(1, 1)] is None
308def test_coadd_provenance_report_counts_one_input_image_in_the_singular() -> None:
309 """Pin the singular forms of the counted phrases.
311 One input image is what the standalone ``coadd_provenance`` file holds,
312 and what a single-cell subset of a coadd holds, so this is the case the
313 ``describe`` command line subcommand shows most often.
314 """
315 inputs = CoaddProvenance.make_empty_input_table(1)
316 inputs["instrument"] = "LSSTCam"
317 inputs["physical_filter"] = "r_57"
318 inputs["visit"] = [101]
319 inputs["detector"] = [1]
320 inputs["day_obs"] = [20250520]
321 provenance = CoaddProvenance(
322 inputs=inputs, contributions=CoaddProvenance.make_empty_contribution_table(1)
323 )
324 assert str(provenance) == "CoaddProvenance(1 input image)"
325 assert repr(provenance) == "<CoaddProvenance(1 input image)>"
326 fields = report_fields(provenance.describe())
327 assert fields["input images"] == "1 from 1 visit"
328 assert fields["per cell"] == "1 input image"
331def test_cell_coadd_report_accounts_for_every_component(minified_cell_coadd: CellCoadd) -> None:
332 """Nothing the coadd carries is absent from its report."""
333 report = minified_cell_coadd.describe()
334 fields = report_fields(report)
335 assert fields["mask_fractions"] == "rejected"
336 assert fields["noise_realizations"] == "1 image (dtype float32)"
337 assert fields["aperture_corrections"] == "3 fields"
338 # Provenance is a child, so it needs no field line of its own.
339 assert "provenance" not in fields
340 assert list(report.children) == [
341 "image",
342 "mask",
343 "variance",
344 "sky_projection",
345 "psf",
346 "provenance",
347 "backgrounds",
348 ]
349 provenance = report.children["provenance"]
350 assert provenance.type_name == "CoaddProvenance"
351 # The coadd passed its cells down, so coverage is stated as a fraction.
352 assert report_fields(provenance)["cells"] == "3 of 3 with contributions"
355def test_cell_coadd_report_counts_aperture_corrections_only(minified_cell_coadd: CellCoadd) -> None:
356 """A coadd can carry dozens of aperture corrections with long names, so
357 the report counts them and never names them, detail or not.
358 """
359 plain = report_fields(minified_cell_coadd.describe())["aperture_corrections"]
360 assert plain == "3 fields"
361 detailed = report_fields(minified_cell_coadd.describe(detail=True))["aperture_corrections"]
362 assert detailed == "3 fields" # No additional detail
365def test_cell_coadd_report_states_absent_provenance(minified_cell_coadd: CellCoadd) -> None:
366 """A coadd with no provenance says so, while components that are merely
367 empty stay out of the report.
368 """
369 bare = CellCoadd(
370 minified_cell_coadd.image,
371 mask=minified_cell_coadd.mask,
372 variance=minified_cell_coadd.variance,
373 sky_projection=minified_cell_coadd.sky_projection,
374 band=minified_cell_coadd.band,
375 psf=minified_cell_coadd.psf,
376 patch=minified_cell_coadd.patch,
377 )
378 report = bare.describe()
379 fields = report_fields(report)
380 assert fields["provenance"] == "none"
381 assert "provenance" not in report.children
382 assert "mask_fractions" not in fields
383 assert "noise_realizations" not in fields
384 assert "aperture_corrections" not in fields
387def test_cell_coadd_copy_preserves_absent_optional_state(minified_cell_coadd: CellCoadd) -> None:
388 """Copying must support coadds without patch or provenance records."""
389 minified_cell_coadd._patch = None
390 minified_cell_coadd._provenance = None
392 copied = minified_cell_coadd.copy()
394 assert copied._patch is None
395 assert copied._provenance is None
396 assert_masked_images_equal(copied, minified_cell_coadd)
397 with pytest.raises(AttributeError, match="no patch"):
398 copied.patch
399 with pytest.raises(AttributeError, match="no provenance"):
400 copied.provenance
403def test_cell_grid_patch_str_uses_clean_geometry() -> None:
404 """CellGrid, PatchDefinition and CellGridBounds str drop the
405 Interval/YX/Box/CellIJ wrappers.
407 The report renders field values with str, so these must use the compact
408 geometry forms rather than pydantic's default field-by-field repr.
409 """
410 grid = CellGrid(bbox=Box.from_shape((100, 200)), cell_shape=YX(10, 20))
411 assert str(grid) == "[y=0:100, x=0:200], cell_shape=(y=10, x=20)"
413 patch = PatchDefinition(id=73, index=YX(7, 3), inner_bbox=Box.factory[1:3, 2:4], cells=grid)
414 assert str(patch) == (
415 "id=73, index=(y=7, x=3), inner_bbox=[y=1:3, x=2:4], "
416 "cells=([y=0:100, x=0:200], cell_shape=(y=10, x=20))"
417 )
418 # repr stays as the pydantic default so it remains eval-ish and distinct.
419 assert "Interval(" in repr(patch)
420 assert "YX(" in repr(patch)
422 bounds = CellGridBounds(grid=grid, bbox=Box.factory[0:40, 0:60])
423 assert str(bounds) == "[y=0:40, x=0:60] in grid ([y=0:100, x=0:200], cell_shape=(y=10, x=20))"
424 bounds_missing = CellGridBounds(
425 grid=grid, bbox=Box.factory[0:40, 0:60], missing=frozenset({CellIJ(i=1, j=1), CellIJ(i=0, j=2)})
426 )
427 assert str(bounds_missing) == (
428 "[y=0:40, x=0:60] in grid ([y=0:100, x=0:200], cell_shape=(y=10, x=20)), "
429 "missing={(i=0, j=2), (i=1, j=1)}"
430 )
431 assert "Interval(" in repr(bounds_missing)
434def test_cell_shape_accepts_both_spellings() -> None:
435 """Verify both on-disk spellings of an XY pair still validate.
437 Shipped cell_coadd 1.0.0 files carry the array spelling, so this stays
438 readable until a major bump retires the as_shipped fixture. A focused
439 check here gives a two-line failure instead of a 52 KB fixture failing
440 to read.
441 """
442 tree_cls = class_for_schema("cell_psf")
443 assert tree_cls is not None
444 model = parameterize_tree(tree_cls, JsonRef)
445 fixture = json.loads(current_fixture_path(FIXTURE_DIR, "cell_psf").read_text())
446 assert fixture["bounds"]["grid"]["cell_shape"] == {"y": 4, "x": 4}
448 legacy = copy.deepcopy(fixture)
449 legacy["bounds"]["grid"]["cell_shape"] = [4, 4]
450 tree = model.model_validate(legacy)
451 assert isinstance(tree, CellPointSpreadFunctionSerializationModel)
452 assert tree.bounds.grid.cell_shape.y == 4
453 assert tree.bounds.grid.cell_shape.x == 4
456def test_cell_psf_rejects_index_below_bounds(minified_cell_coadd: CellCoadd) -> None:
457 """An index below the PSF bounds must not wrap around its array."""
458 start = minified_cell_coadd.psf.bounds.subgrid_start
459 with pytest.raises(BoundsError, match="out of bounds"):
460 minified_cell_coadd.psf[CellIJ(i=start.i - 1, j=start.j)]
463def test_from_legacy(legacy_test_data: _LegacyTestData) -> None:
464 """Test constructing a CellCoadd by converting a legacy
465 ``MultipleCellCoadd``.
466 """
467 assert legacy_test_data.cell_coadd.bounds.missing == {CellIJ(**DP2_COADD_MISSING_CELL)}
468 assert legacy_test_data.cell_coadd.bbox == Box.factory[12900:13500, 9600:10050]
469 compare_cell_coadd_to_legacy(
470 legacy_test_data.cell_coadd,
471 legacy_test_data.legacy_cell_coadd,
472 tract_bbox=legacy_test_data.tract_bbox,
473 plane_map=legacy_test_data.plane_map,
474 psf_points=legacy_test_data.make_psf_points(),
475 )
478def test_roundtrip(legacy_test_data: _LegacyTestData) -> None:
479 """Test that a CellCoadd roundtrips through FITS."""
480 with RoundtripFits(legacy_test_data.cell_coadd, "CellCoadd") as roundtrip:
481 # Check a subimage read (no component arg — does not trigger a skip).
482 subbox = Box.factory[
483 legacy_test_data.cell_coadd.bbox.y.start + 252 : legacy_test_data.cell_coadd.bbox.y.stop - 175,
484 legacy_test_data.cell_coadd.bbox.x.stop - 150 : legacy_test_data.cell_coadd.bbox.x.stop,
485 ]
486 subimage = roundtrip.get(bbox=subbox)
487 assert_masked_images_equal(subimage, legacy_test_data.cell_coadd[subbox], expect_view=False)
488 with roundtrip.inspect() as fits:
489 for extname in ["IMAGE", "MASK", "VARIANCE", "MASK_FRACTIONS/REJECTED"] + [
490 f"NOISE_REALIZATIONS/{n}" for n in range(len(legacy_test_data.cell_coadd.noise_realizations))
491 ]:
492 assert fits[extname].header["ZTILE1"] == legacy_test_data.cell_coadd.grid.cell_shape.x
493 assert fits[extname].header["ZTILE2"] == legacy_test_data.cell_coadd.grid.cell_shape.y
494 # Fixture self-consistency: bbox and missing-cell set are as expected.
495 assert legacy_test_data.cell_coadd.bounds.missing == {CellIJ(**DP2_COADD_MISSING_CELL)}
496 assert legacy_test_data.cell_coadd.bbox == Box.factory[12900:13500, 9600:10050]
497 # Full round-trip fidelity.
498 assert_cell_coadds_equal(roundtrip.result, legacy_test_data.cell_coadd, expect_view=False)
499 compare_cell_coadd_to_legacy(
500 roundtrip.result,
501 legacy_test_data.legacy_cell_coadd,
502 tract_bbox=legacy_test_data.tract_bbox,
503 plane_map=legacy_test_data.plane_map,
504 psf_points=legacy_test_data.make_psf_points(),
505 )
508def test_roundtrip_components(legacy_test_data: _LegacyTestData) -> None:
509 """Test component and subimage reads.
511 This test will be skipped if `lsst.daf.butler` is not available instead of
512 falling back to non-butler I/O, which is why we don't want to merge it
513 with `test_roundtrip`.
514 """
515 with RoundtripFits(legacy_test_data.cell_coadd, "CellCoadd") as roundtrip:
516 subbox = make_subbox(legacy_test_data.cell_coadd.bbox)
517 subpsf = roundtrip.get("psf", bbox=subbox)
518 assert subpsf.bounds.bbox == Box(
519 y=Interval.factory[
520 legacy_test_data.cell_coadd.bbox.y.start + 150 : legacy_test_data.cell_coadd.bbox.y.stop - 150
521 ],
522 x=subbox.x,
523 )
524 assert_psfs_equal(
525 subpsf,
526 legacy_test_data.cell_coadd.psf,
527 points=legacy_test_data.make_psf_points(subbox),
528 )
529 assert roundtrip.get("bbox") == legacy_test_data.cell_coadd.bbox
530 alternates = {
531 k: roundtrip.get(k)
532 for k in [
533 "sky_projection",
534 "image",
535 "mask",
536 "variance",
537 "masked_image",
538 "psf",
539 "aperture_corrections",
540 "provenance",
541 "backgrounds",
542 "bbox",
543 "unit",
544 ]
545 }
546 # Read all the components at once.
547 all_components = roundtrip.get("components")
548 assert set(all_components) == set(alternates) - {"masked_image"}
549 assert all_components["bbox"] == alternates["bbox"]
550 assert all_components["unit"] == legacy_test_data.cell_coadd.unit
551 assert_psfs_equal(all_components["psf"], alternates["psf"])
552 assert_images_equal(all_components["image"], alternates["image"])
554 backgrounds = roundtrip.get("backgrounds")
555 assert backgrounds.keys() == set()
556 assert backgrounds.subtracted is None
558 compare_cell_coadd_to_legacy(
559 roundtrip.result,
560 legacy_test_data.legacy_cell_coadd,
561 tract_bbox=legacy_test_data.tract_bbox,
562 plane_map=legacy_test_data.plane_map,
563 alternates=alternates,
564 psf_points=legacy_test_data.make_psf_points(),
565 )
568def test_fits_compression(legacy_test_data: _LegacyTestData) -> None:
569 """Test lossy FITS compression produces the expected headers."""
570 with RoundtripFits(
571 legacy_test_data.cell_coadd,
572 storage_class="CellCoadd",
573 recipe="lossy16",
574 compression_options={
575 "image": FitsCompressionOptions.LOSSY,
576 "variance": FitsCompressionOptions.LOSSY,
577 },
578 ) as roundtrip:
579 with roundtrip.inspect() as fits:
580 for extname in ["IMAGE", "MASK", "VARIANCE", "MASK_FRACTIONS/REJECTED"] + [
581 f"NOISE_REALIZATIONS/{n}" for n in range(len(legacy_test_data.cell_coadd.noise_realizations))
582 ]:
583 assert fits[extname].header["ZTILE1"] == legacy_test_data.cell_coadd.grid.cell_shape.x
584 assert fits[extname].header["ZTILE2"] == legacy_test_data.cell_coadd.grid.cell_shape.y
585 if extname == "MASK" or extname.startswith("MASK_FRACTIONS"):
586 assert fits[extname].header["ZCMPTYPE"] == "GZIP_2"
587 else:
588 assert fits[extname].header["ZCMPTYPE"] == "RICE_1"
589 assert fits[extname].header["ZQUANTIZ"] == "SUBTRACTIVE_DITHER_2"
592def test_json_roundtrip(legacy_test_data: _LegacyTestData) -> None:
593 """Verify a CellCoadd round-trips correctly through the JSON archive."""
594 with RoundtripJson(legacy_test_data.cell_coadd) as roundtrip:
595 pass
596 assert_cell_coadds_equal(roundtrip.result, legacy_test_data.cell_coadd, expect_view=False)
599def test_to_legacy_cell_coadd(legacy_test_data: _LegacyTestData) -> None:
600 """Verify converting a CellCoadd back into a legacy MultipleCellCoadd."""
601 legacy_cell_coadd = legacy_test_data.cell_coadd.to_legacy_cell_coadd()
602 compare_cell_coadd_to_legacy(
603 legacy_test_data.cell_coadd,
604 legacy_cell_coadd,
605 tract_bbox=legacy_test_data.tract_bbox,
606 plane_map=legacy_test_data.plane_map,
607 psf_points=legacy_test_data.make_psf_points(),
608 )
609 with pytest.raises(
610 ValueError, match="MultipleCellCoadd requires its bounding box to lie on the cell grid."
611 ):
612 legacy_test_data.cell_coadd[make_subbox(legacy_test_data.cell_coadd.bbox)].to_legacy_cell_coadd()
615@skip_no_legacy
616def test_to_legacy(legacy_test_data: _LegacyTestData) -> None:
617 """Test converting a CellCoadd back into a legacy Exposure."""
618 legacy_exposure = legacy_test_data.cell_coadd.to_legacy()
619 assert legacy_exposure.getFilter().bandLabel == legacy_test_data.cell_coadd.band
620 assert Box.from_legacy(legacy_exposure.getBBox()) == legacy_test_data.cell_coadd.bbox
621 compare_masked_image_to_legacy(
622 legacy_test_data.cell_coadd,
623 legacy_exposure.maskedImage,
624 plane_map=legacy_test_data.plane_map,
625 expect_view=True,
626 )
627 compare_psf_to_legacy(
628 legacy_test_data.cell_coadd.psf,
629 legacy_exposure.getPsf(),
630 points=legacy_test_data.make_psf_points(),
631 expect_legacy_raise_on_out_of_bounds=True,
632 )
633 compare_sky_projection_to_legacy_wcs(
634 legacy_test_data.cell_coadd.sky_projection,
635 legacy_exposure.getWcs(),
636 legacy_test_data.cell_coadd.sky_projection.pixel_frame,
637 subimage_bbox=legacy_test_data.cell_coadd.bbox,
638 is_fits=True,
639 )
640 subbox = make_subbox(legacy_test_data.cell_coadd.bbox)
641 compare_masked_image_to_legacy(
642 legacy_test_data.cell_coadd[subbox],
643 legacy_test_data.cell_coadd[subbox].to_legacy().maskedImage,
644 plane_map=legacy_test_data.plane_map,
645 expect_view=True,
646 )
649@skip_no_h5py
650def test_ndf_roundtrip(legacy_test_data: _LegacyTestData) -> None:
651 """Test that CellCoadd round-trips through NDF."""
652 with RoundtripNdf(legacy_test_data.cell_coadd, "CellCoadd") as roundtrip:
653 assert_cell_coadds_equal(roundtrip.result, legacy_test_data.cell_coadd, expect_view=False)
656# The float32 image pixels are of order 10 nJy, as is the test background
657# below, so background arithmetic is only reproducible to a few ULPs at that
658# scale.
659BACKGROUND_ATOL = 1e-5
662def _add_gradient_background(coadd: CellCoadd, name: str = "pretty") -> ChebyshevField:
663 """Add a non-constant background over the coadd's full bbox and return the
664 field that was added.
665 """
666 field = ChebyshevField(
667 coadd.bbox,
668 np.array([[10.0, 2.0, 0.5], [1.0, 0.75, 0.0], [0.25, 0.0, 0.0]]),
669 unit=coadd.unit,
670 )
671 coadd.backgrounds.add(name, field, description="Gradient background for tests.")
672 return field
675def _make_background_subbox(coadd: CellCoadd) -> Box:
676 """Make a box that trims a different number of pixels from each side of the
677 coadd, so a background rendered over the wrong box cannot match by chance.
678 """
679 return Box.factory[
680 coadd.bbox.y.start + 2 : coadd.bbox.y.stop - 3,
681 coadd.bbox.x.start + 1 : coadd.bbox.x.stop - 2,
682 ]
685def test_apply_background_subimage(minified_cell_coadd: CellCoadd) -> None:
686 """Test that applying a background to a subimage subtracts only the
687 portion of the background model that overlaps the subimage.
688 """
689 field = _add_gradient_background(minified_cell_coadd)
690 subbox = _make_background_subbox(minified_cell_coadd)
691 subimage = minified_cell_coadd[subbox]
692 original = subimage.image.array.copy()
694 subimage.apply_background("pretty")
696 assert subimage.backgrounds.subtracted is not None
697 assert subimage.backgrounds.subtracted.name == "pretty"
698 assert subimage.image.bbox == subbox
699 expected = field.render(subbox, dtype=subimage.image.array.dtype).array
700 np.testing.assert_allclose(original - subimage.image.array, expected, atol=BACKGROUND_ATOL)
703def test_restore_background_subimage(minified_cell_coadd: CellCoadd) -> None:
704 """Test that restoring the original background on a subimage adds back
705 only the portion of the background model that overlaps the subimage.
706 """
707 _add_gradient_background(minified_cell_coadd)
708 subbox = _make_background_subbox(minified_cell_coadd)
709 subimage = minified_cell_coadd[subbox]
710 original = subimage.image.array.copy()
711 subimage.apply_background("pretty")
713 subimage.apply_background(None)
715 assert subimage.backgrounds.subtracted is None
716 np.testing.assert_allclose(subimage.image.array, original, atol=BACKGROUND_ATOL)
719def test_failed_background_switch_preserves_pixels_and_state(
720 minified_cell_coadd: CellCoadd,
721) -> None:
722 """An invalid replacement must not restore the current background."""
723 _add_gradient_background(minified_cell_coadd)
724 minified_cell_coadd.apply_background("pretty")
725 before = minified_cell_coadd.image.array.copy()
727 with pytest.raises(KeyError):
728 minified_cell_coadd.apply_background("missing")
730 np.testing.assert_array_equal(minified_cell_coadd.image.array, before)
731 assert minified_cell_coadd.backgrounds.subtracted is not None
732 assert minified_cell_coadd.backgrounds.subtracted.name == "pretty"
735def test_switch_background_applies_replacement(minified_cell_coadd: CellCoadd) -> None:
736 """Switching models restores the old background and subtracts the new."""
737 field = _add_gradient_background(minified_cell_coadd)
738 minified_cell_coadd.backgrounds.add("other", field * 2.0)
739 original = minified_cell_coadd.image.array.copy()
740 minified_cell_coadd.apply_background("pretty")
742 minified_cell_coadd.apply_background("other")
744 expected = (
745 original
746 - (field * 2.0).render(minified_cell_coadd.bbox, dtype=minified_cell_coadd.image.array.dtype).array
747 )
748 np.testing.assert_allclose(minified_cell_coadd.image.array, expected, atol=BACKGROUND_ATOL)
749 assert minified_cell_coadd.backgrounds.subtracted is not None
750 assert minified_cell_coadd.backgrounds.subtracted.name == "other"
753def test_apply_background_after_bounded_read(minified_cell_coadd: CellCoadd) -> None:
754 """Test that a background can be applied to a coadd read with a bbox
755 parameter.
756 """
757 field = _add_gradient_background(minified_cell_coadd)
758 subbox = _make_background_subbox(minified_cell_coadd)
759 with RoundtripJson(minified_cell_coadd, "CellCoadd") as roundtrip:
760 subimage = roundtrip.get(bbox=subbox)
761 original = subimage.image.array.copy()
763 subimage.apply_background("pretty")
765 assert subimage.image.bbox == subbox
766 expected = field.render(subbox, dtype=subimage.image.array.dtype).array
767 np.testing.assert_allclose(original - subimage.image.array, expected, atol=BACKGROUND_ATOL)
770def test_cell_grid_bounds_contains_broadcasting(minified_cell_coadd: CellCoadd) -> None:
771 """Test that CellGridBounds.contains broadcasts like a numpy ufunc."""
772 assert minified_cell_coadd.bounds.missing, "fixture should retain a missing cell"
773 check_bounds_contains_broadcasting(minified_cell_coadd.bounds)
776def test_intersection_bounds_contains_broadcasting(minified_cell_coadd: CellCoadd) -> None:
777 """Test that IntersectionBounds.contains broadcasts like a numpy ufunc."""
778 # Clip the CellGridBounds with a Box offset by 1 pixel on each side so it
779 # does not snap to any cell boundary, forcing a lazy IntersectionBounds.
780 bounds = minified_cell_coadd.bounds
781 clip = Box.factory[
782 bounds.bbox.y.start + 1 : bounds.bbox.y.stop - 1,
783 bounds.bbox.x.start + 1 : bounds.bbox.x.stop - 1,
784 ]
785 check_bounds_contains_broadcasting(bounds.intersection(clip))