Coverage for python/lsst/images/tests/verify_rewrite.py: 34%
214 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-08 01:42 -0800
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
14__all__ = ()
16import copy
17import random
18from collections.abc import Iterable, Sequence
19from typing import TYPE_CHECKING, Any, cast
21import astropy.io.fits
22import click
23import fsspec
24import numpy as np
25from astro_metadata_translator import ObservationInfo
27from lsst.images import (
28 YX,
29 BackgroundMap,
30 Box,
31 DifferenceImage,
32 DifferenceImageTemplateInfo,
33 VisitImage,
34 get_legacy_difference_image_mask_planes,
35 obs_info_from_legacy,
36)
37from lsst.images.cells import CellCoadd
38from lsst.images.convolution_kernels import ConvolutionKernel
39from lsst.images.tests import (
40 assert_obs_metadata_fields_equal,
41 cell_center_points,
42 compare_cell_coadd_to_legacy_lazy,
43 compare_masked_image_to_legacy,
44 compare_visit_image_to_legacy,
45)
47if TYPE_CHECKING:
48 import tqdm
50 from lsst.daf.butler import Butler, DataCoordinate, DimensionRecord
53# These defaults match the rewrite tasks used for DP2 in pipe_tasks. That makes
54# this a bit of a cyclic DRY violation, but not a serious problem.
55VISIT_IMAGE_BACKGROUNDS: tuple[str, ...] = ("subtracted", "skyCorr")
56DIFFERENCE_IMAGE_BACKGROUNDS: tuple[str, ...] = ()
59def _check_kernel(kernel: ConvolutionKernel | None) -> None:
60 """Sanity-check a DifferenceImage kernel."""
61 assert isinstance(kernel, ConvolutionKernel), f"kernel has type {type(kernel)}"
62 assert kernel.bounds is not None
63 assert kernel.bounds.bbox is not None
64 xy = kernel.bounds.bbox.meshgrid(3)
65 for x, y in zip(xy.x.ravel(), xy.y.ravel(), strict=True):
66 im = kernel.compute_kernel_image(x=int(x), y=int(y))
67 assert im.array.size > 0
68 assert np.isfinite(im.array).all(), "kernel has non-finite values"
69 assert np.abs(im.array).max() > 0, "kernel is identically zero"
72def _check_templates(templates: list[DifferenceImageTemplateInfo] | None, detector_bbox: Box) -> None:
73 """Sanity-check the DifferenceImage template-info list."""
74 assert isinstance(templates, list), f"templates has type {type(templates)}"
75 assert len(templates) > 0, "no template info attached"
76 for info in templates:
77 assert detector_bbox.contains(info.bounds.bbox), "template bounds outside detector"
78 if info.psf_shape_flag:
79 continue
80 assert info.psf_shape_xx * info.psf_shape_yy - info.psf_shape_xy**2 > 0, (
81 "template PSF shape is not a valid ellipse"
82 )
85def _check_backgrounds(backgrounds: BackgroundMap, bbox: Box, *, expected: Sequence[str] = ()) -> None:
86 """Sanity-check the backgrounds attached to an image.
88 ``expected`` names must be present. Each attached background's field must
89 cover the image bbox and evaluate to finite values over a coarse grid, and
90 the map's ``subtracted`` attribute must point at a present background.
91 """
92 present = set(backgrounds)
93 for name in expected:
94 assert name in present, f"expected background {name!r} not attached"
95 for name, bg in backgrounds.items():
96 field = bg.field
97 assert field.bounds.bbox.contains(bbox), f"background {name!r} does not cover the image bbox"
98 grid = field.bounds.bbox.meshgrid(2)
99 vals = field(x=grid.x, y=grid.y)
100 assert np.isfinite(vals).all(), f"background {name!r} has non-finite values"
101 if (subtracted := backgrounds.subtracted) is not None:
102 assert subtracted.name in present, "subtracted background designation not in map"
105def _spot_check_points(points: YX[np.ndarray], n: int, rng: random.Random) -> YX[np.ndarray]:
106 """Sample up to ``n`` of the given points without replacement, keeping
107 x/y pairs together.
109 ``rng`` makes the sampling reproducible across runs.
110 """
111 n = min(n, len(points.y))
112 indices = rng.sample(range(len(points.y)), n)
113 return YX(y=points.y[indices], x=points.x[indices])
116@click.command("verify-rewrite")
117@click.argument("repo")
118@click.argument("dataset_type")
119@click.argument("collection")
120@click.option("--where", default="", help="Query string to constraint the comparison.")
121@click.option("--old-prefix", default="legacy_", help="Prefix for the old dataset type names.")
122@click.option("--new-prefix", default="", help="Prefix for the new dataset type names.")
123@click.option(
124 "--require-compressed/--no-require-compressed",
125 default=True,
126 help="Check that the new data product is lossy-compressed.",
127)
128@click.option("--check-kernel/--no-check-kernel", default=True, help="Sanity-check DifferenceImage.kernel.")
129@click.option(
130 "--check-templates/--no-check-templates", default=True, help="Sanity-check DifferenceImage.templates."
131)
132@click.option(
133 "--check-backgrounds/--no-check-backgrounds", default=True, help="Sanity-check attached backgrounds."
134)
135@click.option(
136 "--check-obs-info/--no-check-obs-info",
137 default=True,
138 help=(
139 "Check the new dataset's ObservationInfo against one re-derived from the legacy "
140 "exposure and the repo's exposure dimension record."
141 ),
142)
143@click.option(
144 "--pixel-atol-factor",
145 default=0.5,
146 help="CellCoadd pixel atol is this times std() of the legacy plane.",
147)
148@click.option(
149 "--coadd-check-points",
150 default=4,
151 help=(
152 "Cell centers per coadd used to spot-check PSF and aperture corrections."
153 " Ignored for other dataset types."
154 ),
155)
156def verify_rewrite(
157 *,
158 repo: str,
159 dataset_type: str,
160 collection: str,
161 where: str,
162 require_compressed: bool,
163 old_prefix: str,
164 new_prefix: str,
165 check_kernel: bool,
166 check_templates: bool,
167 check_backgrounds: bool,
168 check_obs_info: bool,
169 pixel_atol_factor: float,
170 coadd_check_points: int,
171) -> None: # numpydoc ignore=PR01
172 """Compare rewritten images in COLLECTION against the originals in
173 COLLECTION of REPO.
174 """
175 try:
176 from lsst.afw.image import Exposure # noqa: F401
177 from lsst.daf.butler import Butler
178 except ImportError as err:
179 err.add_note("verify-rewrite requires a full Rubin development environment.")
180 raise
181 with Butler.from_config(repo, collections=[collection]) as butler:
182 verifier = RewriteVerifier(
183 butler,
184 dataset_type,
185 old_prefix=old_prefix,
186 new_prefix=new_prefix,
187 pixel_atol_factor=pixel_atol_factor,
188 coadd_check_points=coadd_check_points,
189 )
190 # Query for and compare datasets.
191 for data_id in verifier.process(where):
192 # Test visit_image vs. future_visit_image.
193 verifier.compare_images(
194 data_id,
195 check_kernel=check_kernel,
196 check_templates=check_templates,
197 check_backgrounds=check_backgrounds,
198 check_obs_info=check_obs_info,
199 )
200 if require_compressed:
201 verifier.require_compressed(data_id)
203 if verifier.n_problems:
204 raise click.exceptions.Exit(1)
207class RewriteVerifier:
208 def __init__(
209 self,
210 butler: Butler,
211 base_dataset_type: str,
212 *,
213 old_prefix: str,
214 new_prefix: str,
215 pixel_atol_factor: float = 0.5,
216 coadd_check_points: int = 4,
217 ) -> None:
218 self.butler = butler
219 self.base_dataset_type = base_dataset_type
220 self.old_prefix = old_prefix
221 self.new_prefix = new_prefix
222 self.pixel_atol_factor = pixel_atol_factor
223 self.coadd_check_points = coadd_check_points
224 self._progress: tqdm.tqdm | None = None
225 self._n_problems = 0
226 self._exposure_records: dict[int, DimensionRecord] = {}
228 def _report_problem(self) -> None:
229 self._n_problems += 1
231 @property
232 def n_problems(self) -> int:
233 return self._n_problems
235 def process(self, where: str, *dimensions: str) -> Iterable[DataCoordinate]:
236 import tqdm
238 old_dataset_type = self.butler.get_dataset_type(f"{self.old_prefix}{self.base_dataset_type}")
239 dimension_group = old_dataset_type.dimensions.union(self.butler.dimensions.conform(dimensions))
240 visit_group = self.butler.dimensions["visit"].minimal_group
241 with self.butler.query() as query:
242 data_ids = list(
243 tqdm.tqdm(
244 query.where(where).join_dataset_search(old_dataset_type.name).data_ids(dimension_group),
245 desc=f"querying for {old_dataset_type.name}",
246 )
247 )
248 new_count = query.datasets(f"{self.new_prefix}{self.base_dataset_type}").where(where).count()
249 assert len(data_ids) == new_count, f"Count mismatch: new ({new_count}) != old ({len(data_ids)})."
250 # Fetch the exposure dimension records for all of these data IDs
251 # up front. Assume there's only one instrument in play. Coadd
252 # (patch/band) data IDs have no visits, and nothing to fetch.
253 if "visit" in dimension_group:
254 visit_data_ids = {d.subset(visit_group) for d in data_ids}
255 self._exposure_records = {
256 record.id: record
257 for record in query.join_data_coordinates(visit_data_ids).dimension_records("exposure")
258 }
259 random.shuffle(data_ids)
260 self._progress = tqdm.tqdm(data_ids)
261 yield from self._progress
263 def compare_images(
264 self,
265 data_id: DataCoordinate,
266 *,
267 check_kernel: bool = True,
268 check_templates: bool = True,
269 check_backgrounds: bool = True,
270 check_obs_info: bool = True,
271 ) -> None:
272 old = self.butler.get(f"{self.old_prefix}{self.base_dataset_type}", data_id)
273 new = self.butler.get(f"{self.new_prefix}{self.base_dataset_type}", data_id)
274 if check_obs_info and "visit" in data_id.dimensions:
275 try:
276 self._check_obs_info(new, old, data_id)
277 except Exception as err:
278 self.print_error(data_id, err)
279 return
280 expected_backgrounds = None
281 plane_map = None
282 if isinstance(new, CellCoadd):
283 rng = random.Random(f"{self.base_dataset_type} {data_id}")
284 points = _spot_check_points(cell_center_points(new.bounds), self.coadd_check_points, rng)
285 try:
286 compare_cell_coadd_to_legacy_lazy(
287 new,
288 old,
289 skymap=str(data_id["skymap"]),
290 tract=int(data_id["tract"]),
291 patch=int(data_id["patch"]),
292 band=str(data_id["band"]),
293 tract_bbox=new.sky_projection.pixel_frame.bbox,
294 image_atol=self.pixel_atol_factor * float(np.nanstd(old.getImage().array)),
295 variance_atol=self.pixel_atol_factor * float(np.nanstd(old.getVariance().array)),
296 psf_points=points,
297 )
298 except Exception as err:
299 self.print_error(data_id, err)
300 return
301 if isinstance(new, DifferenceImage):
302 expected_backgrounds = DIFFERENCE_IMAGE_BACKGROUNDS
303 plane_map = get_legacy_difference_image_mask_planes()
304 if check_kernel or check_templates:
305 try:
306 if check_kernel:
307 _check_kernel(new.kernel)
308 if check_templates:
309 _check_templates(new.templates, new.bbox)
310 except Exception as err:
311 self.print_error(data_id, err)
312 return
313 if isinstance(new, VisitImage):
314 try:
315 compare_visit_image_to_legacy(
316 new,
317 old,
318 expect_view=False,
319 plane_map=plane_map,
320 check_photometric_scaling=False,
321 instrument=cast(str, data_id["instrument"]),
322 visit=cast(int, data_id["visit"]),
323 detector=cast(int, data_id["detector"]),
324 )
325 except Exception as err:
326 self.print_error(data_id, err)
327 return
328 if expected_backgrounds is None:
329 expected_backgrounds = VISIT_IMAGE_BACKGROUNDS
330 else:
331 try:
332 compare_masked_image_to_legacy(new, old, expect_view=False, plane_map=plane_map)
333 except Exception as err:
334 self.print_error(data_id, err)
335 return
336 if check_backgrounds and expected_backgrounds is not None:
337 _check_backgrounds(new.backgrounds, new.bbox, expected=expected_backgrounds)
339 def _check_obs_info(self, new: Any, old: Any, data_id: DataCoordinate) -> None:
340 """Compare the new dataset's ObservationInfo against one re-derived
341 from the legacy exposure and the repo's exposure dimension record.
343 All mismatching fields are collected and reported together so that
344 one early failure cannot hide the rest.
345 """
346 record = self._exposure_records.get(cast(int, data_id["visit"]))
347 if record is None:
348 raise AssertionError("no exposure record for visit")
349 # Work on a copy, since embedded ObservationInfo cards are stripped
350 # from the metadata that is passed in.
351 metadata = copy.deepcopy(old.getMetadata())
352 expected = obs_info_from_legacy(
353 old.info.getVisitInfo(),
354 record,
355 old.getDetector(),
356 detector_exposure_id=old.info.getId(),
357 telescope=metadata.get("TELESCOP", None),
358 header=metadata,
359 )
360 mismatches: list[str] = []
361 for field in ObservationInfo.model_fields:
362 if field == "warnings":
363 # Translator diagnostics are not part of the science
364 # metadata and need not round-trip.
365 continue
366 actual_value = getattr(new.obs_info, field, None)
367 expected_value = getattr(expected, field, None)
368 if actual_value is None and expected_value is None:
369 continue
370 if actual_value is None or expected_value is None:
371 mismatches.append(f"{field}: {actual_value!r} != {expected_value!r}")
372 continue
373 try:
374 assert_obs_metadata_fields_equal(actual_value, expected_value, label=field)
375 except Exception as err:
376 detail = str(err)
377 if not detail.startswith(f"{field}:"):
378 detail = f"{field}: {detail}"
379 if notes := getattr(err, "__notes__", ()):
380 detail += f" ({' -> '.join(notes)})"
381 mismatches.append(detail)
382 if mismatches:
383 raise AssertionError(
384 f"{len(mismatches)} ObservationInfo field(s) differ: " + "; ".join(mismatches)
385 )
387 def require_compressed(self, data_id: DataCoordinate) -> None:
388 ref = self.butler.find_dataset(f"{self.new_prefix}{self.base_dataset_type}", data_id)
389 assert ref is not None, f"new dataset for {data_id} is missing"
390 path = self.butler.getURI(ref)
391 fs: fsspec.AbstractFileSystem
392 fs, fp = path.to_fsspec()
393 with fs.open(fp) as stream:
394 with astropy.io.fits.open(stream, disable_image_compression=True) as hdu_list:
395 if (image_zcmptype := hdu_list["IMAGE"].header.get("ZCMPTYPE")) != "RICE_1":
396 self._report_problem()
397 self.print(data_id, f"IMAGE HDU has ZCMPTYPE={image_zcmptype!r}")
398 if (mask_zcmptype := hdu_list["MASK"].header.get("ZCMPTYPE")) != "GZIP_2":
399 self._report_problem()
400 self.print(data_id, f"MASK HDU has ZCMPTYPE={mask_zcmptype!r}")
401 if (variance_zcmptype := hdu_list["VARIANCE"].header.get("ZCMPTYPE")) != "RICE_1":
402 self._report_problem()
403 self.print(data_id, f"VARIANCE HDU has ZCMPTYPE={variance_zcmptype!r}")
405 def print(self, data_id: DataCoordinate, msg: str) -> None:
406 if self._progress is not None:
407 self._progress.write(f"{data_id}: {msg}")
408 else:
409 print(f"{data_id}: {msg}")
411 def print_error(self, data_id: DataCoordinate, err: Exception) -> None:
412 self._report_problem()
413 message = f"{err} ({type(err).__name__})"
414 notes = getattr(err, "__notes__", ())
415 if notes:
416 # The comparison utilities attach the failing component's name as
417 # an exception note (a hierarchy from outer to inner when blocks
418 # are nested, though they never are today). Prefix the message
419 # with that path so it reads as
420 # "<data-id>\n <component>: <message>"
421 # rather than an indented note trailing the message.
422 component = " -> ".join(notes)
423 block = f"{data_id}\n {component}: {message}"
424 else:
425 block = f"{data_id}: {message}"
426 if self._progress is not None: 426 ↛ 427line 426 didn't jump to line 427 because the condition on line 426 was never true
427 self._progress.write(block)
428 else:
429 print(block)
432if __name__ == "__main__":
433 verify_rewrite()