Coverage for python/lsst/images/tests/verify

1# This file is part of lsst-images.

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.

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 random

17from collections.abc import Iterable, Sequence

18from typing import TYPE_CHECKING, Any, cast

20import astropy.io.fits

21import click

22import fsspec

23import numpy as np

24from astro_metadata_translator import ObservationInfo

26from lsst.images import (

27 BackgroundMap,

28 Box,

29 DifferenceImage,

30 DifferenceImageTemplateInfo,

31 VisitImage,

32 get_legacy_difference_image_mask_planes,

33 obs_info_from_legacy,

34)

35from lsst.images.convolution_kernels import ConvolutionKernel

36from lsst.images.tests import (

37 assert_obs_metadata_fields_equal,

38 compare_masked_image_to_legacy,

39 compare_visit_image_to_legacy,

40)

42if TYPE_CHECKING:

43 import tqdm

45 from lsst.daf.butler import Butler, DataCoordinate, DimensionRecord

48# These defaults match the rewrite tasks used for DP2 in pipe_tasks. That makes

49# this a bit of a cyclic DRY violation, but not a serious problem.

50VISIT_IMAGE_BACKGROUNDS: tuple[str, ...] = ("subtracted", "skyCorr")

51DIFFERENCE_IMAGE_BACKGROUNDS: tuple[str, ...] = ()

54def _check_kernel(kernel: ConvolutionKernel | None) -> None:

55 """Sanity-check a DifferenceImage kernel."""

56 assert isinstance(kernel, ConvolutionKernel), f"kernel has type {type(kernel)}"

57 assert kernel.bounds is not None

58 assert kernel.bounds.bbox is not None

59 xy = kernel.bounds.bbox.meshgrid(3)

60 for x, y in zip(xy.x.ravel(), xy.y.ravel(), strict=True):

61 im = kernel.compute_kernel_image(x=int(x), y=int(y))

62 assert im.array.size > 0

63 assert np.isfinite(im.array).all(), "kernel has non-finite values"

64 assert np.abs(im.array).max() > 0, "kernel is identically zero"

67def _check_templates(templates: list[DifferenceImageTemplateInfo] | None, detector_bbox: Box) -> None:

68 """Sanity-check the DifferenceImage template-info list."""

69 assert isinstance(templates, list), f"templates has type {type(templates)}"

70 assert len(templates) > 0, "no template info attached"

71 for info in templates:

72 assert detector_bbox.contains(info.bounds.bbox), "template bounds outside detector"

73 if info.psf_shape_flag:

74 continue

75 assert info.psf_shape_xx * info.psf_shape_yy - info.psf_shape_xy**2 > 0, (

76 "template PSF shape is not a valid ellipse"

77 )

80def _check_backgrounds(backgrounds: BackgroundMap, bbox: Box, *, expected: Sequence[str] = ()) -> None:

81 """Sanity-check the backgrounds attached to an image.

83 ``expected`` names must be present. Each attached background's field must

84 cover the image bbox and evaluate to finite values over a coarse grid, and

85 the map's ``subtracted`` attribute must point at a present background.

86 """

87 present = set(backgrounds)

88 for name in expected:

89 assert name in present, f"expected background {name!r} not attached"

90 for name, bg in backgrounds.items():

91 field = bg.field

92 assert field.bounds.bbox.contains(bbox), f"background {name!r} does not cover the image bbox"

93 grid = field.bounds.bbox.meshgrid(2)

94 vals = field(x=grid.x, y=grid.y)

95 assert np.isfinite(vals).all(), f"background {name!r} has non-finite values"

96 if (subtracted := backgrounds.subtracted) is not None:

97 assert subtracted.name in present, "subtracted background designation not in map"

100@click.command("verify-rewrite")

101@click.argument("repo")

102@click.argument("dataset_type")

103@click.argument("collection")

104@click.option("--where", default="", help="Query string to constraint the comparison.")

105@click.option("--old-prefix", default="legacy_", help="Prefix for the old dataset type names.")

106@click.option("--new-prefix", default="", help="Prefix for the new dataset type names.")

107@click.option(

108 "--require-compressed/--no-require-compressed",

109 default=True,

110 help="Check that the new data product is lossy-compressed.",

111)

112@click.option("--check-kernel/--no-check-kernel", default=True, help="Sanity-check DifferenceImage.kernel.")

113@click.option(

114 "--check-templates/--no-check-templates", default=True, help="Sanity-check DifferenceImage.templates."

115)

116@click.option(

117 "--check-backgrounds/--no-check-backgrounds", default=True, help="Sanity-check attached backgrounds."

118)

119@click.option(

120 "--check-obs-info/--no-check-obs-info",

121 default=True,

122 help=(

123 "Check the new dataset's ObservationInfo against one re-derived from the legacy "

124 "exposure and the repo's exposure dimension record."

125 ),

126)

127def verify_rewrite(

128 *,

129 repo: str,

130 dataset_type: str,

131 collection: str,

132 where: str,

133 require_compressed: bool,

134 old_prefix: str,

135 new_prefix: str,

136 check_kernel: bool,

137 check_templates: bool,

138 check_backgrounds: bool,

139 check_obs_info: bool,

140) -> None: # numpydoc ignore=PR01

141 """Compare rewritten images in COLLECTION against the originals in

142 COLLECTION of REPO.

143 """

144 try:

145 from lsst.afw.image import Exposure # noqa: F401

146 from lsst.daf.butler import Butler

147 except ImportError as err:

148 err.add_note("verify-rewrite requires a full Rubin development environment.")

149 raise

150 with Butler.from_config(repo, collections=[collection]) as butler:

151 verifier = RewriteVerifier(butler, dataset_type, old_prefix=old_prefix, new_prefix=new_prefix)

152 # Query for and compare datasets.

153 for data_id in verifier.process(where):

154 # Test visit_image vs. future_visit_image.

155 verifier.compare_images(

156 data_id,

157 check_kernel=check_kernel,

158 check_templates=check_templates,

159 check_backgrounds=check_backgrounds,

160 check_obs_info=check_obs_info,

161 )

162 if require_compressed:

163 verifier.require_compressed(data_id)

164

165 if verifier.n_problems:

166 raise click.exceptions.Exit(1)

167

168

169class RewriteVerifier:

170 def __init__(

171 self,

172 butler: Butler,

173 base_dataset_type: str,

174 *,

175 old_prefix: str,

176 new_prefix: str,

177 ) -> None:

178 self.butler = butler

179 self.base_dataset_type = base_dataset_type

180 self.old_prefix = old_prefix

181 self.new_prefix = new_prefix

182 self._progress: tqdm.tqdm | None = None

183 self._n_problems = 0

184 self._exposure_records: dict[int, DimensionRecord] = {}

185

186 def _report_problem(self) -> None:

187 self._n_problems += 1

188

189 @property

190 def n_problems(self) -> int:

191 return self._n_problems

192

193 def process(self, where: str, *dimensions: str) -> Iterable[DataCoordinate]:

194 import tqdm

195

196 old_dataset_type = self.butler.get_dataset_type(f"{self.old_prefix}{self.base_dataset_type}")

197 dimension_group = old_dataset_type.dimensions.union(self.butler.dimensions.conform(dimensions))

198 visit_group = self.butler.dimensions["visit"].minimal_group

199 with self.butler.query() as query:

200 data_ids = list(

201 tqdm.tqdm(

202 query.where(where).join_dataset_search(old_dataset_type.name).data_ids(dimension_group),

203 desc=f"querying for {old_dataset_type.name}",

204 )

205 )

206 new_count = query.datasets(f"{self.new_prefix}{self.base_dataset_type}").where(where).count()

207 assert len(data_ids) == new_count, f"Count mismatch: new ({new_count}) != old ({len(data_ids)})."

208 # Fetch the exposure dimension records for all of these data IDs

209 # up front. Assume there's only one instrument in play.

210 visit_data_ids = {d.subset(visit_group) for d in data_ids}

211 self._exposure_records = {

212 record.id: record

213 for record in query.join_data_coordinates(visit_data_ids).dimension_records("exposure")

214 }

215 random.shuffle(data_ids)

216 self._progress = tqdm.tqdm(data_ids)

217 yield from self._progress

218

219 def compare_images(

220 self,

221 data_id: DataCoordinate,

222 *,

223 check_kernel: bool = True,

224 check_templates: bool = True,

225 check_backgrounds: bool = True,

226 check_obs_info: bool = True,

227 ) -> None:

228 old = self.butler.get(f"{self.old_prefix}{self.base_dataset_type}", data_id)

229 new = self.butler.get(f"{self.new_prefix}{self.base_dataset_type}", data_id)

230 if check_obs_info:

231 try:

232 self._check_obs_info(new, old, data_id)

233 except Exception as err:

234 self.print_error(data_id, err)

235 return

236 expected_backgrounds = None

237 plane_map = None

238 if isinstance(new, DifferenceImage):

239 expected_backgrounds = DIFFERENCE_IMAGE_BACKGROUNDS

240 plane_map = get_legacy_difference_image_mask_planes()

241 if check_kernel or check_templates:

242 try:

243 if check_kernel:

244 _check_kernel(new.kernel)

245 if check_templates:

246 _check_templates(new.templates, new.bbox)

247 except Exception as err:

248 self.print_error(data_id, err)

249 return

250 if isinstance(new, VisitImage):

251 try:

252 compare_visit_image_to_legacy(

253 new,

254 old,

255 expect_view=False,

256 plane_map=plane_map,

257 check_photometric_scaling=False,

258 instrument=cast(str, data_id["instrument"]),

259 visit=cast(int, data_id["visit"]),

260 detector=cast(int, data_id["detector"]),

261 )

262 except Exception as err:

263 self.print_error(data_id, err)

264 return

265 if expected_backgrounds is None:

266 expected_backgrounds = VISIT_IMAGE_BACKGROUNDS

267 else:

268 try:

269 compare_masked_image_to_legacy(new, old, expect_view=False, plane_map=plane_map)

270 except Exception as err:

271 self.print_error(data_id, err)

272 return

273 if check_backgrounds and expected_backgrounds is not None:

274 _check_backgrounds(new.backgrounds, new.bbox, expected=expected_backgrounds)

275

276 def _check_obs_info(self, new: Any, old: Any, data_id: DataCoordinate) -> None:

277 """Compare the new dataset's ObservationInfo against one re-derived

278 from the legacy exposure and the repo's exposure dimension record.

279

280 All mismatching fields are collected and reported together so that

281 one early failure cannot hide the rest.

282 """

283 record = self._exposure_records.get(cast(int, data_id["visit"]))

284 if record is None:

285 raise AssertionError("no exposure record for visit")

286 telescope = old.getMetadata().get("TELESCOP", None)

287 expected = obs_info_from_legacy(

288 old.info.getVisitInfo(),

289 record,

290 old.getDetector(),

291 detector_exposure_id=old.info.getId(),

292 telescope=telescope,

293 )

294 mismatches: list[str] = []

295 for field in ObservationInfo.model_fields:

296 if field == "warnings":

297 # Translator diagnostics are not part of the science

298 # metadata and need not round-trip.

299 continue

300 actual_value = getattr(new.obs_info, field, None)

301 expected_value = getattr(expected, field, None)

302 if actual_value is None and expected_value is None:

303 continue

304 if actual_value is None or expected_value is None:

305 mismatches.append(f"{field}: {actual_value!r} != {expected_value!r}")

306 continue

307 try:

308 assert_obs_metadata_fields_equal(actual_value, expected_value, label=field)

309 except Exception as err:

310 detail = str(err)

311 if not detail.startswith(f"{field}:"):

312 detail = f"{field}: {detail}"

313 if notes := getattr(err, "__notes__", ()):

314 detail += f" ({' -> '.join(notes)})"

315 mismatches.append(detail)

316 if mismatches:

317 raise AssertionError(

318 f"{len(mismatches)} ObservationInfo field(s) differ: " + "; ".join(mismatches)

319 )

320

321 def require_compressed(self, data_id: DataCoordinate) -> None:

322 ref = self.butler.find_dataset(f"{self.new_prefix}{self.base_dataset_type}", data_id)

323 assert ref is not None, f"new dataset for {data_id} is missing"

324 path = self.butler.getURI(ref)

325 fs: fsspec.AbstractFileSystem

326 fs, fp = path.to_fsspec()

327 with fs.open(fp) as stream:

328 with astropy.io.fits.open(stream, disable_image_compression=True) as hdu_list:

329 if (image_zcmptype := hdu_list["IMAGE"].header.get("ZCMPTYPE")) != "RICE_1":

330 self._report_problem()

331 self.print(data_id, f"IMAGE HDU has ZCMPTYPE={image_zcmptype!r}")

332 if (mask_zcmptype := hdu_list["MASK"].header.get("ZCMPTYPE")) != "GZIP_2":

333 self._report_problem()

334 self.print(data_id, f"MASK HDU has ZCMPTYPE={mask_zcmptype!r}")

335 if (variance_zcmptype := hdu_list["VARIANCE"].header.get("ZCMPTYPE")) != "RICE_1":

336 self._report_problem()

337 self.print(data_id, f"VARIANCE HDU has ZCMPTYPE={variance_zcmptype!r}")

338

339 def print(self, data_id: DataCoordinate, msg: str) -> None:

340 if self._progress is not None:

341 self._progress.write(f"{data_id}: {msg}")

342 else:

343 print(f"{data_id}: {msg}")

344

345 def print_error(self, data_id: DataCoordinate, err: Exception) -> None:

346 self._report_problem()

347 message = f"{err} ({type(err).__name__})"

348 notes = getattr(err, "__notes__", ())

349 if notes:

350 # The comparison utilities attach the failing component's name as

351 # an exception note (a hierarchy from outer to inner when blocks

352 # are nested, though they never are today). Prefix the message

353 # with that path so it reads as

354 # "<data-id>\n <component>: <message>"

355 # rather than an indented note trailing the message.

356 component = " -> ".join(notes)

357 block = f"{data_id}\n {component}: {message}"

358 else:

359 block = f"{data_id}: {message}"

360 if self._progress is not None: 360 ↛ 361line 360 didn't jump to line 361 because the condition on line 360 was never true

361 self._progress.write(block)

362 else:

363 print(block)

364

365

366if __name__ == "__main__":

367 verify_rewrite()

Coverage for python/lsst/images/tests/verify_rewrite.py: 33%

195 statements