Coverage for python/lsst/images/tests/_minify_for_fixtures.py: 20%
175 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-24 08:46 +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.
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
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
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
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
119def minify(in_path: str, out_path: str) -> None:
120 """Read a real archive at ``in_path``, take a small subset, and write JSON.
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.
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)
143 os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
144 subset.write(out_path)
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)}.")
158# -- VisitImage ------------------------------------------------------------
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.
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]]
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)
201 if not linearize_projection or subset.sky_projection is None:
202 return subset
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 )
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``.
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).
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.
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)
277def _affine_polymap_coeffs(matrix: np.ndarray, offset: np.ndarray) -> np.ndarray:
278 """Build AST ``PolyMap`` coefficients for ``out = matrix @ in + offset``.
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)
294def _simplify_piff_psf(psf: PiffWrapper, *, order: int) -> None:
295 """Truncate a Piff PSF's field interpolation to ``order``, in place.
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).
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.
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
317 from piff import BasisPolynomial
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)))
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
334# -- CellCoadd -------------------------------------------------------------
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}.")
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)]
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
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)
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
380 def shrink2d(array: np.ndarray) -> np.ndarray:
381 return np.ascontiguousarray(array[::step_y, ::step_x][:ny, :nx])
383 def shrink3d(array: np.ndarray) -> np.ndarray:
384 return np.ascontiguousarray(array[::step_y, ::step_x, :][:ny, :nx, :])
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)
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 ]
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)
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 )
421 provenance = block._provenance
422 if provenance is not None:
423 provenance = _trim_provenance(provenance, max_inputs=max_inputs)
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
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 )
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.
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])
482def _choose_block_bbox(cell_coadd: CellCoadd) -> Box:
483 """Return the pixel bbox of a (up to) 2x2 block of cells to keep.
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)
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
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]
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]