Coverage for python/lsst/images/tests/_minify_for

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.

12"""Minify a real on-disk archive into a small JSON test fixture.

14Reads a FITS or NDF file via the appropriate input archive, takes a

15small subset of the in-memory object, and writes JSON via

16``JsonOutputArchive``. Used to populate ``tests/data/schema_v1/legacy/``

17with derived-from-real test data that exercises the full read path

18including the absence-of-stamp legacy default.

20Per top-level type the subset rule is:

22 VisitImage Crop the image/mask/variance planes to a small (~16x16)

23 corner, keeping the real single-instance structures (PSF

24 such as Piff, detector frames) that synthetic fixtures

25 cannot reproduce -- the whole point of deriving a fixture

26 from real data. Homogeneous repeated collections (detector

27 amplifiers, aperture-correction entries) are trimmed to a

28 representative few, since one entry exercises the schema as

29 well as sixteen. The projection's pixel->sky mapping is

30 replaced by its linear (affine) approximation over the kept

31 box: a real TAN-SIP WCS serializes as a ~100 KB AST

32 polynomial dump, but over a 16x16 box it is linear to far

33 below a pixel, so the affine form is schema-identical and

34 orders of magnitude smaller. A Piff PSF's field

35 interpolation is truncated to a low order (the order-4

36 solution table is ~225 KB; order 0, the field-averaged PSF,

37 is schema-identical and ~13x smaller).

39 CellCoadd Crop to a small block of cells (preferring a block that

40 includes a missing cell so the sparse-grid path is

41 exercised) and then *morph* that block onto a tiny cell

42 grid: each cell's planes are decimated from the native

43 cell size down to a few pixels and re-stitched, and the

44 PSF kernels are cropped to a small odd window. The grid

45 topology (number of cells, the missing-cell set, band,

46 mask schema and provenance shape) is preserved; the pixel

47 values and WCS are *not* physically meaningful. This is

48 the "morph cells in place" fallback: it sidesteps the

49 outer-ring problem (inputs/PSFs that overlap kept cells)

50 by rebuilding a self-consistent miniature coadd rather

51 than trying to carve an accurate subset out of the real

52 one. An accurate per-cell subset would inline several

53 150x150 planes per cell and produce multi-megabyte JSON,

54 which defeats the purpose of a fixture.

56Run interactively (CellCoadd works with just this package installed;

57VisitImage needs a full Rubin environment so the real PSF can be read)::

59 python -c "

60 from lsst.images.tests._minify_for_fixtures import minify

61 minify('cell_example.fits', 'tests/data/schema_v1/legacy/cell_coadd.json')

62 minify('dp1.fits', 'tests/data/schema_v1/legacy/visit_image_dp1.json')

63 minify('dp2.fits', 'tests/data/schema_v1/legacy/visit_image_dp2.json')

64 "

66The helper is invoked manually by developers when they have a real

67on-disk file to derive from; it is not exercised by CI.

68"""

70from __future__ import annotations

72__all__ = ("minify",)

74import os

75from collections.abc import Callable

76from typing import Any

78import numpy as np

80from .. import DifferenceImage, VisitImage

81from .._cell_grid import CellGrid, CellGridBounds, CellIJ, PatchDefinition

82from .._geom import YX, Box

83from .._image import Image

84from .._mask import Mask

85from .._transforms import SkyProjection, TractFrame, Transform

86from .._transforms._ast import PolyMap

87from ..cells import CellCoadd, CellField, CellPointSpreadFunction, CoaddProvenance

88from ..psfs import PiffWrapper

89from ..serialization import read

90from ._creation import make_random_sky_projection

92# Default morph parameters for CellCoadd. ``CELL_SIZE`` should divide the

93# native cell size evenly; ``KERNEL_SIZE`` must be odd. ``MAX_INPUTS`` caps

94# the provenance ``inputs`` table (a real coadd has hundreds of visits); the

95# full provenance schema is already exercised by the ``coadd_provenance``

96# fixture, so here we keep just enough rows to be representative.

97_CELL_SIZE = 6

98_KERNEL_SIZE = 5

99_MAX_INPUTS = 6

100

101# Default trim parameters for VisitImage. Amplifiers and aperture-correction

102# entries are homogeneous collections, so a couple of each cover the schema

103# just as well as the full set (a real detector has 16 amplifiers and dozens

104# of aperture corrections).

105_MAX_AMPLIFIERS = 2

106_MAX_APERTURE_CORRECTIONS = 2

107

108# Field-interpolation order to truncate a Piff PSF to (the solution table of a

109# real order-4 PixelGrid PSF dominates the fixture at ~225 KB). Order 0 is the

110# field-averaged PSF; set to `None` to leave the PSF untouched.

111_PSF_INTERP_ORDER = 0

112

113# Maximum permitted deviation (radians) when approximating a projection's

114# pixel->sky mapping with an affine one. Over a fixture's tiny box the real

115# mapping is linear well below this, so the fit always succeeds.

116_PROJECTION_LINEAR_APPROX_TOL = 1e-8

117

118

119def minify(in_path: str, out_path: str) -> None:

120 """Read a real archive at ``in_path``, take a small subset, and write JSON.

121

122 Parameters

123 ----------

124 in_path

125 Path to a FITS (``.fits`` / ``.fits.gz``) or NDF (``.sdf`` / ``.ndf``)

126 file to read.

127 out_path

128 Path to the output file to write. The parent directory is

129 created if it does not exist. Can be any supported file format.

130 File extension controls the output format.

131

132 Raises

133 ------

134 ValueError

135 If the file extension is not recognised.

136 NotImplementedError

137 If the top-level type is not one this helper knows how to subset.

138 """

139 obj = read(in_path)

140 subsetter = _dispatch(obj)

141 subset = subsetter(obj)

142

143 os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)

144 subset.write(out_path)

145

146

147def _dispatch(image: Any) -> Callable[[Any], Any]:

148 """Return the relevant subsetter for this image object."""

149 match image:

150 case VisitImage() | DifferenceImage(): 150 ↛ 151line 150 didn't jump to line 151 because the pattern on line 150 never matched

151 return _subset_visit_image

152 case CellCoadd(): 152 ↛ 153line 152 didn't jump to line 153 because the pattern on line 152 never matched

153 return _subset_cell_coadd

154 case _:

155 raise NotImplementedError(f"No minify rule for image of type {type(image)}.")

156

157

158# -- VisitImage ------------------------------------------------------------

159

160

161def _subset_visit_image(

162 visit_like_image: VisitImage | DifferenceImage,

163 *,

164 size: int = 16,

165 max_amplifiers: int = _MAX_AMPLIFIERS,

166 max_aperture_corrections: int = _MAX_APERTURE_CORRECTIONS,

167 linearize_projection: bool = True,

168 projection_tol: float = _PROJECTION_LINEAR_APPROX_TOL,

169 psf_interp_order: int | None = _PSF_INTERP_ORDER,

170) -> VisitImage | DifferenceImage:

171 """Crop a VisitImage's pixel planes to a small corner and trim its

172 homogeneous collections.

173

174 The detector frames are a single structure carried through unchanged by

175 ``__getitem__``. The detector's amplifiers and the aperture-correction map

176 are repeated, schema-identical entries, so they are trimmed to a

177 representative few. The projection's pixel->sky mapping is replaced by its

178 affine approximation over the kept box (see ``_linear_approx_projection``)

179 unless ``linearize_projection`` is false. A Piff PSF's field interpolation

180 is truncated to ``psf_interp_order`` (see ``_simplify_piff_psf``) unless

181 that is `None`.

182 """

183 bbox = visit_like_image.bbox

184 y0 = bbox.y.start

185 x0 = bbox.x.start

186 y1 = min(y0 + size, bbox.y.stop)

187 x1 = min(x0 + size, bbox.x.stop)

188 subset = visit_like_image[Box.factory[y0:y1, x0:x1]]

189

190 # ``subset`` is a fresh throwaway object whose detector amplifier list,

191 # aperture-correction map and PSF are live, mutable components. Trim them

192 # in place through the public accessors rather than reaching for private

193 # attributes.

194 del subset.detector.amplifiers[max_amplifiers:]

195 aperture_corrections = subset.aperture_corrections

196 for key in list(aperture_corrections)[max_aperture_corrections:]:

197 del aperture_corrections[key]

198 if psf_interp_order is not None and isinstance(subset.psf, PiffWrapper):

199 _simplify_piff_psf(subset.psf, order=psf_interp_order)

200

201 if not linearize_projection or subset.sky_projection is None:

202 return subset

203

204 # The pixel planes carry the projection immutably (there is no public

205 # setter for it), so install the affine approximation by rebuilding the

206 # VisitImage from its public components with re-viewed planes. Only the

207 # image plane's projection is actually serialized, but keeping all three

208 # consistent avoids surprises.

209 linear = _linear_approx_projection(subset.sky_projection, subset.image.bbox, tol=projection_tol)

210 return type(visit_like_image)(

211 subset.image.view(sky_projection=linear),

212 mask=subset.mask.view(sky_projection=linear),

213 variance=subset.variance.view(sky_projection=linear),

214 sky_projection=linear,

215 psf=subset.psf,

216 obs_info=subset.obs_info,

217 bounds=subset.bounds,

218 summary_stats=subset.summary_stats,

219 detector=subset.detector,

220 photometric_scaling=subset.photometric_scaling,

221 aperture_corrections=subset.aperture_corrections,

222 backgrounds=subset.backgrounds,

223 band=subset.band,

224 metadata=subset.metadata,

225 )

226

227

228def _linear_approx_projection(sky_projection: SkyProjection, bbox: Box, *, tol: float) -> SkyProjection:

229 """Return a copy of ``sky_projection`` whose pixel->sky mapping is replaced

230 by its best linear (affine) approximation over ``bbox``.

231

232 Real WCS mappings (e.g. TAN-SIP) serialize as large AST polynomial dumps.

233 Over the small box of a fixture they are linear to far below a pixel, so

234 an affine approximation is schema-identical but orders of magnitude

235 smaller. The result carries no FITS approximation (the affine is itself

236 trivially FITS-representable).

237

238 This is written as a self-contained ``sky_projection -> sky_projection``

239 transform so it can be promoted to a public

240 ``SkyProjection.linear_approx(bbox, tol)`` method later with essentially

241 no change. It assumes a 2-D pixel->sky

242 mapping.

243

244 Parameters

245 ----------

246 sky_projection

247 The projection to approximate.

248 bbox

249 Box (in pixel coordinates) over which the approximation must hold.

250 tol

251 Maximum permitted deviation from linearity, as a Cartesian

252 displacement in the output (sky, radians) coordinates. AST raises

253 ``RuntimeError`` if no fit within ``tol`` exists.

254 """

255 transform = sky_projection.pixel_to_sky_transform

256 mapping = transform._ast_mapping

257 lbnd = [bbox.x.start, bbox.y.start]

258 ubnd = [bbox.x.stop, bbox.y.stop]

259 # linearApprox yields [offsets; Jacobian] as a (1 + n_out, n_in) array on

260 # both AST backends (astshim returns the flat buffer in the same order, so

261 # the reshape recovers the same layout the starlink-pyast bridge returns).

262 fit = np.asarray(mapping.linearApprox(lbnd, ubnd, tol), dtype=float).reshape(3, 2)

263 offset = fit[0] # (lon0, lat0), radians

264 jacobian = fit[1:] # jacobian[i, j] = d(out_i) / d(in_j), in = (x, y)

265 jacobian_inv = np.linalg.inv(jacobian)

266 forward = _affine_polymap_coeffs(jacobian, offset)

267 inverse = _affine_polymap_coeffs(jacobian_inv, -jacobian_inv @ offset)

268 affine = Transform(

269 transform.in_frame,

270 transform.out_frame,

271 PolyMap(forward, inverse),

272 in_bounds=sky_projection.pixel_bounds,

273 )

274 return SkyProjection(affine)

275

276

277def _affine_polymap_coeffs(matrix: np.ndarray, offset: np.ndarray) -> np.ndarray:

278 """Build AST ``PolyMap`` coefficients for ``out = matrix @ in + offset``.

279

280 Each row is ``[coefficient, output_axis (1-based), power_of_in_1, ...]``;

281 one constant row plus one row per input per output axis. Returned as a

282 float array, which is the form both AST backends require.

283 """

284 n = len(offset)

285 coeffs: list[list[float]] = []

286 for i in range(n):

287 coeffs.append([float(offset[i]), i + 1, *([0] * n)])

288 for j in range(n):

289 powers = [1 if k == j else 0 for k in range(n)]

290 coeffs.append([float(matrix[i][j]), i + 1, *powers])

291 return np.array(coeffs, dtype=float)

292

293

294def _simplify_piff_psf(psf: PiffWrapper, *, order: int) -> None:

295 """Truncate a Piff PSF's field interpolation to ``order``, in place.

296

297 A real Piff PSF interpolates a per-pixel model across the focal plane with

298 a high-order 2-D polynomial; that solution table dominates the serialized

299 size (a 25x25 PixelGrid x order-4 polynomial is ~225 KB). Truncating to

300 ``order`` keeps only the lowest-order field terms -- order 0 is the

301 field-averaged PSF -- which is schema-identical but far smaller, and needs

302 no stars or refit (the fitted ``stars`` are already dropped on serialize).

303

304 Only ``BasisPolynomial``-interpolated PSFs are handled; anything else (a

305 higher-order model already at/under ``order``, a non-polynomial interp) is

306 left untouched.

307

308 ``piff`` is imported lazily because it is an optional dependency; this is

309 only ever reached when the PSF being simplified is itself a Piff PSF.

310 """

311 interp = getattr(psf.piff_psf, "interp", None)

312 if interp is None or type(interp).__name__ != "BasisPolynomial" or interp.q is None:

313 return

314 if order >= max(interp._orders):

315 return

316

317 from piff import BasisPolynomial

318

319 # ``q`` has one column per active basis term; the terms are the True cells

320 # of ``_mask`` in row-major (i, j) order (see BasisPolynomial.basis). Make

321 # the same ordering for a lower-order interp and copy the shared columns.

322 def _terms(orders: tuple[int, ...], mask: np.ndarray) -> list[tuple[int, ...]]:

323 grids = np.meshgrid(*[np.arange(o + 1) for o in orders], indexing="ij")

324 return list(zip(*(grid[mask].tolist() for grid in grids)))

325

326 old_terms = _terms(interp._orders, interp._mask)

327 truncated = BasisPolynomial(order, keys=list(interp._keys))

328 new_terms = _terms(truncated._orders, truncated._mask)

329 column_of = {term: index for index, term in enumerate(old_terms)}

330 truncated.q = np.ascontiguousarray(interp.q[:, [column_of[term] for term in new_terms]])

331 psf.piff_psf.interp = truncated

332

333

334# -- CellCoadd -------------------------------------------------------------

335

336

337def _subset_cell_coadd(

338 cell_coadd: CellCoadd,

339 *,

340 cell_size: int = _CELL_SIZE,

341 kernel_size: int = _KERNEL_SIZE,

342 max_inputs: int = _MAX_INPUTS,

343) -> CellCoadd:

344 """Crop a CellCoadd to a small block of cells and morph it onto a tiny

345 grid (see the module docstring for the rationale).

346 """

347 if kernel_size % 2 == 0:

348 raise ValueError(f"kernel_size must be odd, got {kernel_size}.")

349

350 # 1. Pick a block of (up to) 2x2 cells, preferring one that contains a

351 # missing cell so the sparse-grid path is exercised. Falls back to the

352 # first available block when the coadd is fully dense.

353 block = cell_coadd[_choose_block_bbox(cell_coadd)]

354

355 grid = block.grid

356 cs = grid.cell_shape

357 start = block.bounds.subgrid_start

358 stop = block.bounds.subgrid_stop

359 n_i = stop.i - start.i

360 n_j = stop.j - start.j

361

362 # 2. Build a tiny full-patch grid with the same cell *count* as the

363 # original patch but ``cell_size`` pixels per cell, anchored at (0, 0).

364 full_shape = grid.grid_size

365 new_grid = CellGrid(

366 bbox=Box.factory[0 : full_shape.i * cell_size, 0 : full_shape.j * cell_size],

367 cell_shape=YX(y=cell_size, x=cell_size),

368 )

369 new_block_bbox = _scale_box_to_grid(block.bbox, grid, cell_size)

370 new_bounds = CellGridBounds(grid=new_grid, bbox=new_block_bbox, missing=block.bounds.missing)

371

372 # 3. Decimate each plane. Because the block's planes tile the kept cells

373 # contiguously, a uniform stride that maps one native cell onto

374 # ``cell_size`` samples is equivalent to per-cell decimation.

375 step_y = max(1, cs.y // cell_size)

376 step_x = max(1, cs.x // cell_size)

377 ny = n_i * cell_size

378 nx = n_j * cell_size

379

380 def shrink2d(array: np.ndarray) -> np.ndarray:

381 return np.ascontiguousarray(array[::step_y, ::step_x][:ny, :nx])

382

383 def shrink3d(array: np.ndarray) -> np.ndarray:

384 return np.ascontiguousarray(array[::step_y, ::step_x, :][:ny, :nx, :])

385

386 # 4. Synthetic-but-valid sky_projection over the tiny tract frame.

387 rng = np.random.default_rng(0)

388 tract_frame = TractFrame(skymap=cell_coadd.skymap, tract=cell_coadd.tract, bbox=new_grid.bbox)

389 sky_projection = make_random_sky_projection(rng, tract_frame, new_block_bbox)

390

391 unit = cell_coadd.unit

392 image = Image(shrink2d(block.image.array), bbox=new_block_bbox, unit=unit, sky_projection=sky_projection)

393 mask = Mask(shrink3d(block.mask.array), schema=block.mask.schema, bbox=new_block_bbox)

394 variance = Image(shrink2d(block.variance.array), bbox=new_block_bbox, unit=unit**2)

395 mask_fractions = {

396 name: Image(shrink2d(plane.array), bbox=new_block_bbox)

397 for name, plane in block.mask_fractions.items()

398 }

399 noise_realizations = [

400 Image(shrink2d(plane.array), bbox=new_block_bbox) for plane in block.noise_realizations

401 ]

402

403 # 5. Crop the PSF kernels to a small odd window about their centre,

404 # keeping the (n_i, n_j) per-cell structure and NaN-for-missing cells.

405 psf_array = block.psf._array

406 ky, kx = psf_array.shape[2:]

407 half = kernel_size // 2

408 cy, cx = ky // 2, kx // 2

409 psf_array = np.ascontiguousarray(psf_array[:, :, cy - half : cy + half + 1, cx - half : cx + half + 1])

410 psf = CellPointSpreadFunction(psf_array, bounds=new_bounds)

411

412 # 6. Patch geometry scaled onto the tiny grid; provenance and backgrounds

413 # are reused as-is (provenance is cell-indexed and already subset).

414 patch = PatchDefinition(

415 id=block.patch.id,

416 index=block.patch.index,

417 inner_bbox=_scale_box_to_grid(block.patch.inner_bbox, grid, cell_size),

418 cells=new_grid,

419 )

420

421 provenance = block._provenance

422 if provenance is not None:

423 provenance = _trim_provenance(provenance, max_inputs=max_inputs)

424

425 # Aperture corrections are not subset when CellCoadd is subset with a

426 # bounding box, because they're always tiny. But that makes setting

427 # up a consistent new grid for them tricky.

428 aperture_corrections = {}

429 new_apcorr_bounds = None

430 for i, (name, field) in enumerate(block.aperture_corrections.items()):

431 if new_apcorr_bounds is None:

432 new_apcorr_bounds = CellGridBounds(

433 grid=new_grid,

434 bbox=_scale_box_to_grid(field.bounds.bbox, grid, cell_size),

435 missing=cell_coadd.bounds.missing,

436 )

437 aperture_corrections[name] = CellField(new_apcorr_bounds, field._array)

438 if i >= 2:

439 break

440

441 return CellCoadd(

442 image,

443 mask=mask,

444 variance=variance,

445 mask_fractions=mask_fractions,

446 noise_realizations=noise_realizations,

447 sky_projection=sky_projection,

448 band=block.band,

449 psf=psf,

450 patch=patch,

451 provenance=provenance,

452 backgrounds=block._backgrounds,

453 aperture_corrections=aperture_corrections,

454 )

455

456

457def _trim_provenance(provenance: CoaddProvenance, *, max_inputs: int) -> CoaddProvenance:

458 """Cap the provenance ``inputs`` table to ``max_inputs`` rows and drop any

459 contributions that reference the removed inputs.

460

461 The two-table structure, polygon arrays and string dictionary-compression

462 paths are all preserved; only the number of contributing visits shrinks.

463 """

464 inputs = provenance.inputs

465 if len(inputs) <= max_inputs:

466 return provenance

467 kept_inputs = inputs[:max_inputs]

468 keys = {(str(row["instrument"]), int(row["visit"]), int(row["detector"])) for row in kept_inputs}

469 contributions = provenance.contributions

470 mask = np.array(

471 [

472 (str(instrument), int(visit), int(detector)) in keys

473 for instrument, visit, detector in zip(

474 contributions["instrument"], contributions["visit"], contributions["detector"]

475 )

476 ],

477 dtype=bool,

478 )

479 return CoaddProvenance(inputs=kept_inputs, contributions=contributions[mask])

480

481

482def _choose_block_bbox(cell_coadd: CellCoadd) -> Box:

483 """Return the pixel bbox of a (up to) 2x2 block of cells to keep.

484

485 Prefers a block containing a missing cell; otherwise the block anchored at

486 the start of the populated region. Never raises if there is no missing

487 cell.

488 """

489 bounds = cell_coadd.bounds

490 grid = bounds.grid

491 start = bounds.subgrid_start

492 stop = bounds.subgrid_stop

493 span_i = min(2, stop.i - start.i)

494 span_j = min(2, stop.j - start.j)

495

496 target = next(iter(sorted(bounds.missing)), None)

497 if target is not None:

498 # Anchor the block so it includes the missing cell, clamped to the

499 # populated index range.

500 i0 = min(max(target.i, start.i), stop.i - span_i)

501 j0 = min(max(target.j, start.j), stop.j - span_j)

502 else:

503 i0 = start.i

504 j0 = start.j

505

506 lo = grid.bbox_of(CellIJ(i=i0, j=j0))

507 hi = grid.bbox_of(CellIJ(i=i0 + span_i - 1, j=j0 + span_j - 1))

508 return Box.factory[lo.y.start : hi.y.stop, lo.x.start : hi.x.stop]

509

510

511def _scale_box_to_grid(box: Box, grid: CellGrid, cell_size: int) -> Box:

512 """Map a grid-aligned box onto a grid with ``cell_size`` pixels per cell,

513 anchored at the origin.

514 """

515 cs = grid.cell_shape

516 s = grid.bbox.start

517 iy0 = (box.y.start - s.y) // cs.y

518 iy1 = (box.y.stop - s.y) // cs.y

519 ix0 = (box.x.start - s.x) // cs.x

520 ix1 = (box.x.stop - s.x) // cs.x

521 return Box.factory[iy0 * cell_size : iy1 * cell_size, ix0 * cell_size : ix1 * cell_size]

Coverage for python/lsst/images/tests/_minify_for_fixtures.py: 20%

175 statements