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 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 BackgroundMap,

29 Box,

30 DifferenceImage,

31 DifferenceImageTemplateInfo,

32 VisitImage,

33 get_legacy_difference_image_mask_planes,

34 obs_info_from_legacy,

35)

36from lsst.images.convolution_kernels import ConvolutionKernel

37from lsst.images.tests import (

38 assert_obs_metadata_fields_equal,

39 compare_masked_image_to_legacy,

40 compare_visit_image_to_legacy,

41)

43if TYPE_CHECKING:

44 import tqdm

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

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

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

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

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

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

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

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

58 assert kernel.bounds is not None

59 assert kernel.bounds.bbox is not None

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

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

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

63 assert im.array.size > 0

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

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

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

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

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

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

72 for info in templates:

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

74 if info.psf_shape_flag:

75 continue

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

77 "template PSF shape is not a valid ellipse"

78 )

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

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

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

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

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

87 """

88 present = set(backgrounds)

89 for name in expected:

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

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

92 field = bg.field

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

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

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

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

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

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

100

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

102@click.argument("repo")

103@click.argument("dataset_type")

104@click.argument("collection")

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

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

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

108@click.option(

109 "--require-compressed/--no-require-compressed",

110 default=True,

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

112)

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

114@click.option(

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

116)

117@click.option(

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

119)

120@click.option(

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

122 default=True,

123 help=(

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

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

126 ),

127)

128def verify_rewrite(

129 *,

130 repo: str,

131 dataset_type: str,

132 collection: str,

133 where: str,

134 require_compressed: bool,

135 old_prefix: str,

136 new_prefix: str,

137 check_kernel: bool,

138 check_templates: bool,

139 check_backgrounds: bool,

140 check_obs_info: bool,

141) -> None: # numpydoc ignore=PR01

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

143 COLLECTION of REPO.

144 """

145 try:

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

147 from lsst.daf.butler import Butler

148 except ImportError as err:

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

150 raise

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

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

153 # Query for and compare datasets.

154 for data_id in verifier.process(where):

155 # Test visit_image vs. future_visit_image.

156 verifier.compare_images(

157 data_id,

158 check_kernel=check_kernel,

159 check_templates=check_templates,

160 check_backgrounds=check_backgrounds,

161 check_obs_info=check_obs_info,

162 )

163 if require_compressed:

164 verifier.require_compressed(data_id)

165

166 if verifier.n_problems:

167 raise click.exceptions.Exit(1)

168

169

170class RewriteVerifier:

171 def __init__(

172 self,

173 butler: Butler,

174 base_dataset_type: str,

175 *,

176 old_prefix: str,

177 new_prefix: str,

178 ) -> None:

179 self.butler = butler

180 self.base_dataset_type = base_dataset_type

181 self.old_prefix = old_prefix

182 self.new_prefix = new_prefix

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

184 self._n_problems = 0

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

186

187 def _report_problem(self) -> None:

188 self._n_problems += 1

189

190 @property

191 def n_problems(self) -> int:

192 return self._n_problems

193

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

195 import tqdm

196

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

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

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

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

201 data_ids = list(

202 tqdm.tqdm(

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

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

205 )

206 )

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

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

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

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

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

212 self._exposure_records = {

213 record.id: record

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

215 }

216 random.shuffle(data_ids)

217 self._progress = tqdm.tqdm(data_ids)

218 yield from self._progress

219

220 def compare_images(

221 self,

222 data_id: DataCoordinate,

223 *,

224 check_kernel: bool = True,

225 check_templates: bool = True,

226 check_backgrounds: bool = True,

227 check_obs_info: bool = True,

228 ) -> None:

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

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

231 if check_obs_info:

232 try:

233 self._check_obs_info(new, old, data_id)

234 except Exception as err:

235 self.print_error(data_id, err)

236 return

237 expected_backgrounds = None

238 plane_map = None

239 if isinstance(new, DifferenceImage):

240 expected_backgrounds = DIFFERENCE_IMAGE_BACKGROUNDS

241 plane_map = get_legacy_difference_image_mask_planes()

242 if check_kernel or check_templates:

243 try:

244 if check_kernel:

245 _check_kernel(new.kernel)

246 if check_templates:

247 _check_templates(new.templates, new.bbox)

248 except Exception as err:

249 self.print_error(data_id, err)

250 return

251 if isinstance(new, VisitImage):

252 try:

253 compare_visit_image_to_legacy(

254 new,

255 old,

256 expect_view=False,

257 plane_map=plane_map,

258 check_photometric_scaling=False,

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

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

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

262 )

263 except Exception as err:

264 self.print_error(data_id, err)

265 return

266 if expected_backgrounds is None:

267 expected_backgrounds = VISIT_IMAGE_BACKGROUNDS

268 else:

269 try:

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

271 except Exception as err:

272 self.print_error(data_id, err)

273 return

274 if check_backgrounds and expected_backgrounds is not None:

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

276

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

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

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

280

281 All mismatching fields are collected and reported together so that

282 one early failure cannot hide the rest.

283 """

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

285 if record is None:

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

287 # Work on a copy, since embedded ObservationInfo cards are stripped

288 # from the metadata that is passed in.

289 metadata = copy.deepcopy(old.getMetadata())

290 expected = obs_info_from_legacy(

291 old.info.getVisitInfo(),

292 record,

293 old.getDetector(),

294 detector_exposure_id=old.info.getId(),

295 telescope=metadata.get("TELESCOP", None),

296 header=metadata,

297 )

298 mismatches: list[str] = []

299 for field in ObservationInfo.model_fields:

300 if field == "warnings":

301 # Translator diagnostics are not part of the science

302 # metadata and need not round-trip.

303 continue

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

305 expected_value = getattr(expected, field, None)

306 if actual_value is None and expected_value is None:

307 continue

308 if actual_value is None or expected_value is None:

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

310 continue

311 try:

312 assert_obs_metadata_fields_equal(actual_value, expected_value, label=field)

313 except Exception as err:

314 detail = str(err)

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

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

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

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

319 mismatches.append(detail)

320 if mismatches:

321 raise AssertionError(

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

323 )

324

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

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

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

328 path = self.butler.getURI(ref)

329 fs: fsspec.AbstractFileSystem

330 fs, fp = path.to_fsspec()

331 with fs.open(fp) as stream:

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

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

334 self._report_problem()

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

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

337 self._report_problem()

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

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

340 self._report_problem()

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

342

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

344 if self._progress is not None:

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

346 else:

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

348

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

350 self._report_problem()

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

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

353 if notes:

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

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

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

357 # with that path so it reads as

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

359 # rather than an indented note trailing the message.

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

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

362 else:

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

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

365 self._progress.write(block)

366 else:

367 print(block)

368

369

370if __name__ == "__main__":

371 verify_rewrite()

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

196 statements